Long-term evolution of supercritical black hole accretion with outflows: a subgrid feedback model for cosmological simulations

We study the long-term evolution of the global structure of axisymmetric accretion flows onto a black hole (BH) at rates substantially higher than the Eddington value ($\dot{M}_{\rm Edd}$), performing two-dimensional hydrodynamical simulations with and without radiative diffusion. In the high-accretion optically-thick limit, where the radiation energy is efficiently trapped within the inflow, the accretion flow becomes adiabatic and comprises of turbulent gas in the equatorial region and strong bipolar outflows. As a result, the mass inflow rate decreases toward the center as $\dot{M}_{\rm in}\propto r^{p}$ with $p\sim 0.5-0.7$ and a small fraction of the inflowing gas feeds the nuclear BH. Thus, super-Eddington accretion is sustained only when a larger amount of gas is supplied from larger radii at $>100-1000~\dot{M}_{\rm Edd}$. The global structure of the flow settles down to a quasi-steady state in millions of the orbital timescale at the BH event horizon, which is $>10-100$ times longer than that addressed in previous (magneto-)RHD simulation studies. Energy transport via radiative diffusion accelerates the outflow near the poles in the inner region but does not change the overall properties of the accretion flow compared to the cases without diffusion. Based on our simulation results, we provide a mechanical feedback model for super-Eddington accreting BHs. This can be applied as a sub-grid model in large-scale cosmological simulations that do not sufficiently resolve galactic nuclei, and to the formation of the heaviest gravitational-wave sources via accretion in dense environments.

Several formation channels of high-redshift SMBHs have been proposed (see Inayoshi et al. 2020;Volonteri et al. 2021, for recent comprehensive reviews). One pos-sible approach is to initiate with seed BHs with mass of ∼ 10 − 100 M , formed after the collapse of the first generation stars (Abel et al. 2002;Bromm et al. 2002;Yoshida et al. 2008;Hirano et al. 2014). It is possible for these seeds to grow into SMBHs within ∼ 1 Gyr if the BHs continuously accrete at or slightly above (Tanaka & Haiman 2009;Madau et al. 2014) the Eddington rate ofṀ Edd [≡ L Edd /(0.1c 2 )] with a 10% radiative efficiency, where L Edd is the Eddington luminosity, c is the speed of light. An alternative approach is to begin with heavier BH seeds with ∼ 10 4 −10 6 M via collapse of supermassive stars (SMSs) or runaway stellar mergers in a dense star cluster (Bromm & Loeb 2003;Omukai et al. 2008;Devecchi & Volonteri 2009;Inayoshi et al. 2014;Hirano et al. 2017;Chon & Omukai 2020;Li et al. 2021). Under special circumstances with strong H 2 photo-dissociating radiation produced from nearby galaxies, violent successive halo mergers, or high baryon-dark matter streaming motion, massive seed BHs are expected to form efficiently and get a head start toward the SMBH regime (Omukai 2001;Dijkstra et al. 2008;Tanaka & Li 2014;Inayoshi et al. 2018a;Wise et al. 2019).
In either case, efficient subsequent growth of the seed BHs is required to explain the 10 9 M SMBHs at early times. The nature of rapid mass accretion onto BHs has been extensively investigated. Observations of ultra-luminous X-ray sources (King, A. R. et al. 2001;Watarai et al. 2001) and narrow-line Seyfert-1 galaxies (Wang et al. 1999;Mineshige et al. 2000) have provided several lines of evidence that both stellar mass BHs and SMBHs are able to accrete gas at rates exceeding the Eddington rates.
Super-Eddington accretion onto stellar-mass BHs in dense gaseous environments (e.g., accretion disks of active galactic nuclei; hereafter AGNs) has been further explored (McKernan et al. 2012;Stone et al. 2017;Yang et al. 2019;Secunda et al. 2019;Tagawa et al. 2020Tagawa et al. , 2021 since the discoveries of intermediate-mass BHs detected as gravitational-wave (GW) sources (Abbott et al. 2021a). The mass spectrum of those merging BHs shows a continuous distribution, namely a single power-law + a Gaussian peak or a double power-law (Abbott et al. 2021b), without a discontinuity that is expected to appear owing to (pulsational) pair-instability supernova (PISN) explosions of massive stellar progenitors (Barkat et al. 1967;Fowler & Hoyle 1964;Heger & Woosley 2002;Heger et al. 2003;Vink et al. 2012;Belczynski et al. 2016;Abbott et al. 2021b). Indeed, several GW events such as GW190521, GW190602 175927, and GW190519 153544 are inferred to have their primary BH mass in the expected PISN mass gap. The existence of such massive BHs would not be preferred in the isolated field-binary scenario due to PISNe caused by single massive stars with M ∼ 50−120 M . Therefore, other formation channels through hierarchical BH mergers (Tagawa et al. 2021) and/or the subsequent growth of stellar remnant BHs via super-Eddington accretion have been considered as a solution for the existence of massive GW sources (Safarzadeh & Haiman 2020).
Rapid mass accretion onto BHs has been investigated by analytical (e.g. Begelman 1979) and numerical work. In recent decades, radiation magneto-hydrodynamical (RMHD) simulations of accreting BHs have shown that the BH feeding rate can exceed the critical value as long as a dense gaseous torus already exists or a sufficient amount of gas supply is maintained from larger radii to the vicinity of the BH (see also Ohsuga & Mineshige 2011;Jiang et al. 2014;Sadowski et al. 2015).
Indeed, they found that BHs can be fed at mildly super-Eddington accretion rates of ∼ 2 − 40Ṁ Edd at least for short durations of ≈ 10 4−5 t 0 , where t 0 (≡ GM BH /c 3 ) is the light-crossing timescle at half of the Schwarzschild radius (r Sch ) around a BH with a mass of M BH . However, the gravitational energy released as radiation (and/or in outflows) heats the gas surrounding the BH and thus can quench the gas supply from the BH's gravitational influence radius. As a result of radiative feedback, when the density of the ambient gas is low, the BH feeding rate is strongly suppressed and the time-averaged rate is limited below the Eddington accretion rate (e.g. Ciotti & Ostriker 2001;Alvarez et al. 2009;Milosavljević et al. 2009b,a;Park & Ricotti 2011Jeon et al. 2012). This discrepancy between the results regarding the BH feeding rates and impacts of radiative feedback is owing to inconsistent treatments of numerical simulations, most of which have focused on the dynamics of gas and radiation at either small or large scales. This fact emphasizes the importance of understanding the global structure of accretion flows (see a recent work by Lalakos et al. 2022, where both the BH event horizon scale and Bondi scale for hot accreting gas are simultaneously resolved.).
Recently, the global properties of super-Eddington accretion flows have been extensively studied with radiation hydrodynamical simulations with a sufficiently large computation domain that enables us to capture the multi-scale physics properly (e.g., Inayoshi et al. 2016;Takeo et al. 2018;Toyouchi et al. 2018;Takeo et al. 2020). Unlike the cases mentioned above, when the BH is embedded in sufficiently dense gas, the emergent radiation flux is substantially reduced by photon trapping and dust absorption in the flow. Radiative feedback does not prevent the gas inflow and instead the ionized region surrounding the BH collapses. As a result, rapid quasi-steady mass accretion exceeding the Eddington rate ( 500Ṁ Edd ) through a dense, optically thick disk is achieved. The extremely high accretion rate is expected to take place in dense regions such as the centers of high-redshift protogalaxies and/or the interiors of AGN disks even at lower redshifts. If such intense accretion flows are maintained down to the BH event horizon and can directly feed the BH without significant mass loss at larger radii, this process potentially solves the problems regarding (1) SMBH formation and (2) the lack of the PISN mass gap on the mass spectrum of merging binary BHs.
In this paper, we investigate the properties of accretion flows onto a BH at very high rates ofṀ ∼ 100 − 2, 000Ṁ Edd with two-dimensional, axisymmetric radiation hydrodynamical simulations extending to ∼ 10 3 r Sch . In particular, the following two issues are addressed. First, we study the long-term evolution of highrate accretion flows over t > 3×10 6 t 0 , which is required to reach a true steady state of the flow; namely, one viscous timescale at ∼ 10 3 r Sch . We note that the previous simulations that studied mildly super-Eddington accreting flows were able to reach 10% of this duration due to computational limitations such as short characteristic timescales of radiation transport and requirement of a more accurate numerical algorithm to solve the radiation-transfer equations for those systems. On the other hand, in our simulations where the inflow rate is much higher and the radiation energy is efficiently trapped within the inflow, radiation transport is well approximated by diffusion and the timescale of radiation transport becomes longer in the dense accreting matter, enabling longer-term simulations. Secondly, we consider various types of initial conditions and boundary conditions to check their effects on the resulting accretion flow.
The rest of this paper is organized as follows. In §2, we describe the methodology of our simulations. In §3, we discuss the simulation results and present the selfsimilar behavior of supercritical accretion flows. In §4 we present a subgrid model of AGN feedback for large scale simulations in the presence of outflows. In §5, we summarize the main conclusions of this paper.

Basic equations
We study the gas dynamics at the vicinity of a BH, performing axisymmetric two-dimensional radiation hydrodynamical simulations with the open source code PLUTO (Mignone et al. 2007). The mass flux for the ideal fluid is computed using the Harten-Lax-vanLeer Riemann solver (Harten et al. 1983), and the second-order accuracy in space and time is ensured. The diffusion terms (viscosity and radiation transport) are treated with the explicit time integration (see below). The basic equations in spherical coordinates (r, θ, φ) are the equation of continuity, and the equation of motion, where ρ is the density, v is the velocity, p is the total pressure, and Φ = −GM BH /(r − r Sch ) is the gravitational potential of the BH (Paczyńsky & Wiita 1980), and σ is the viscous stress tensor. We solve the energy where e = 1 2 ρ|v| 2 + ε is the total energy density, ε is the internal energy density, F 0 is the radiative flux in the comoving frame, and the right-hand-side represents viscous heating.
In our axisymmetric simulations without magnetohydrodynamical (MHD) effects, angular momentum transport in the accretion flow is calculated by imposing viscosity. The viscous stress tensor is given by where ν is the shear viscosity and the bulk viscosity is neglected. To mimic angular momentum transport associated with MHD turbulence driven by the magnetorotational instability (MRI) in a disk (Stone et al. 1996;Balbus & Hawley 1991, 1998, we assume the azimuthal components of the shear tensor are non-zero. In spherical polar coordinates, the tensor components are given by (e.g., Stone et al. 1999;D'Orazio et al. 2013). The strength of anomalous shear viscosity is calculated with the α-prescription (Shakura & Sunyaev 1973), where α is the viscous parameter, Ω K ≡ (GM BH /r 3 ) 1/2 is the Keplerian angular frequency, and c s = γp/ρ is the sound speed. In our simulations, we adopt a single value of α = 0.01 and activate the viscous process within the flow using the function f (θ) defined by The choice of this functions is to mimic viscosity restricted inside the dense region (or disk). Since the gas density near the poles is several orders of magnitude lower than that near the equator, weak viscosity in the polar region hardly affects our conclusions.
In our simulations, we assume that the radiation pressure dominates the total pressure (Begelman 1978). Note-This table shows the initial setups and some results for 2D simulations. Density is in g cm −3 .Ṁ0 is the inflow rate at the outer boundary.ṀBH is the net accretion onto the BH at the inner boundary. r tr,eff is the effective trapping radius inferred from the rest-frame luminosity. Here the 'H', 'L', 'N', 'I', 'D' represent the setup of boundary conditions and 'LR' is short for low resolution. For Case-I, a constant radial inflow velocity is imposed at the outer boundary in the equatorial region; otherwise the outer boundary conditions depend on the direction of the radial velocity (see the text). As discussed in §3.3, the mass inflow and BH accretion rate show the self-similar nature and thus those absolute rates in the cases without radiative diffusion can be renormalized to other values as long asṀin Ṁ Edd and τ c/v hold. Note that the photon trapping radius is well-defined only with radiative diffusion, but we show the values for the cases without radiative diffusion for reference values.
This approximation is valid when the mass accretion rate significantly exceeds the Eddington rate and the flow is highly optically thick. The internal energy density is then given by where E 0 is radiation energy density, P 0 ( p) is the radiation pressure, and the adiabatic index corresponds to γ 4/3. We also adopt the flux-limited diffusion (FLD) approximation, where the radiation flux is calculated with where κ es is the opacity due to electron scattering and the flux limiter λ is given by where R = |∇E 0 | / (ρκ es E 0 ). In the optically thick limit (R → 0), the diffusion coefficient approaches λ = 1/3. The Lorentz transformation of the radiation flux between the rest frame F and the comoving frame F 0 is given by Note that using this transformation, the radiation moments in the divergence term of Eq.
(3) are replaced by the radiation flux F in the observer's rest-frame. The first and second term corresponds to the diffusion and advection term, respectively. Within dense and fast accretion flows, the two terms are comparable and the restframe flux decreases because |F − F 0 | τ |(v/c) · F 0 |, where τ is the optical depth within a scale where the radiation energy density varies. In fact, the ratio of these terms in the spherically symmetric flow (along the radial direction) is whereṀ = 4πρvr 2 is the mass flow rate. This means that the radiative flux in the rest frame has a positive (outward) value at r > r tr and a negative (inward) value at r < r tr , respectively. We note that the expression of the so-called photon-trapping radius (r tr ≡ κ esṀ /4πc) is valid for spherically symmetric flows (Begelman 1978). Nevertheless, we refer to this as a reference in the following discussion.
In what follows, we find that when the mass inflow rate in the computational domain is substantially higher than the Eddington value (Ṁ in Ṁ Edd ), the simulation result with radiative diffusion is qualitatively similar to those without radiative diffusion because most of the photons are transported via advection with the flow rather than via diffusion. Moreover, since the radiative diffusion time becomes substantially shorter than the dynamical time, we need a long computational time to reach a steady state of the flow. Therefore, to avoid this issue, we first focus on simulations without radiative diffusion and discuss how the outer boundary conditions affect the flow properties (see §3.1 and §3.2) and next show the case with radiative diffusion for comparison (see §3.3).

Initial & Boundary Conditions
In this section, we present the initial and boundary conditions for our simulations. We set a computational region of r min ≤ r ≤ r max and 0 + ≤ θ ≤ π − with r min = 3 r Sch , r max 1500 r Sch and = 0.01 to avoid a numerical singularity at the poles. We assume a BH with mass of M BH to be located at the coordinate center (r = 0). We set logarithmically spaced grids in the rdirection to capture the gas behavior at high density and set uniformly spaced grids in the polar direction. The number of grid points is (N r , N θ ) = (256, 256), except for low resolution cases where (N r , N θ ) = (128, 128). We run six simulations with the same initial conditions but with different outer boundary conditions. In Table 1, the simulation setups are summarized.
As initial conditions, we adopt power-law distributions of the gas density, pressure, and radial velocity that follow a Bondi-like inflow; namely, ρ(r) = ρ out (r/r max ) −3/2 , p(r) = p out (ρ/ρ out ) 4/3 , where the adiabatic index is assumed to be γ = 4/3, v r (r) = v out (r/r max ) −1/2 , and v θ = 0. The normalization factors of those profiles are set to ρ out = 7.5×10 −9 g cm −3 , p out = 2 × 10 8 erg cm −3 , v out = −5 × 10 8 cm s −1 . This choice is adopted so that the mass inflow rate from r max becomes as high as ∼ O(100 − 1000)Ṁ Edd and the flow is highly optically thick. We add gas rotation to this Bondi-like profile, assuming the initial specific angular momentum to be the Keplerian value of j z = √ GM BH r, at R(≡ r sin θ) ≥ r max ; otherwise no initial rotation. Note that the initial distribution of the rotational velocity does not affect the flow properties in the late time owing to angular momentum injection through the outer boundary (see below).
We impose injection boundary conditions at r = r max , unlike previous studies where the BH is fed by a dense torus that is initially set at the vicinity of the BH (e.g. Jiang et al. 2014;Sadowski et al. 2015). In our simulations, we consider four different outer-boundary conditions as summarized in Table 1. First, we divide the outer boundary layer into an equatorial (|θ − π/2| < θ 0 ) and a polar region (|θ − π/2| > θ 0 ), where θ 0 is a parameter to be determined in each case. In the polar regions, zero gradients across the boundary are imposed on physical quantities of the inflowing gas using ghost cells (numerically, boundary conditions are set two cells outside the grid referred to as ghost cells) and inflows are prohibited (i.e., v r ≥ 0 is imposed). In the equatorial region, we set two different boundary conditions depending on the direction of the radial velocity. Specifically, the zero-gradient boundary conditions are imposed as in the polar region if v r ≥ 0 (i.e., outflows), and constant values of the density (ρ 0 ) and rotational velocity (v rot ) are set at the outer boundary if v r < 0 (i.e., inflows). This treatment allows us to maintain a high accretion rate through the outer boundary when the gas is inflowing. For Case-L (Low), we adopt θ 0 = π/4, ρ 0 = 6 × 10 −9 g cm −3 , and v rot = 5 × 10 8 cm s −1 . For Case-N (Narrow), we set a smaller injection angle of θ 0 = π/12. For Case-H (High), we increase the inflow density from the outer boundary by a factor of ten.
As for "Injection" boundary conditions, the polar region conditions are the same as for the "Low" boundary condition, while in the equatorial region, the gas is constantly injected into the domain with fixed density (ρ 0 = 6×10 −9 g cm −3 ), radial velocity (v r = −3×10 8 cm s −1 ) and rotation velocity (v rot = 5 × 10 8 cm s −1 ). As for inner boundary conditions, sink boundary conditions are adopted since the inner radius is roughly the "inner-most stable circular orbit" (ISCO) radius. For the boundaries near the polar axes, periodic boundary conditions are applied to the rotational velocity and reflection boundary conditions are imposed to all the other physical quantities. Since we do not include cooling/heating or chemical reactions, the Courant time (the Courant number is set to 0.3) is sufficiently accurate to solve the hydrodynamical equations.

2D HYDRODYNAMICAL SIMULATIONS
In this section, we present the simulation results with different boundary conditions. For all cases, we set M BH = 10 3 M . In §3.1, we show the results for Case-H (regarded as the"Fiducial" case) and discuss the general behaviour of gas dynamics. In §3.2, we discuss how the different outer boundary conditions which provide different gas supply rates affect the properties of the accreting gas and the BH feeding rate, In §3.3, we show the results with radiative diffusion in the dense, optically thick inflows and discuss the effect on supercritical accretion. We also discuss the luminosity as well as photon-trapping effects. The curves indicate the accretion onto the BH (red solid), and the inflow and outflow rates at the outer boundary (blue solid and dotted, respectively). While in the early stage, both the inflow and outflow rate at the outer boundary oscillate, they approach a constant 1 × 10 6 2 × 10 6 3 × 10 6

The Fiducial Case
Net accretion Inflow at Rmax Outflow at Rmax Figure 1. Time evolution of the accretion rate onto the BH (red curve), the mass inflow (blue solid curve), and the outflow rate (blue dotted curve) from the outer boundary for the fiducial case (the Case-H). The accretion rates and timescales are normalized by the Eddington accretion rate and t0, respectively. The simulation settles down to a quasisteady state atṀBH 100Ṁ Edd after t 2.3 × 10 6 t0 1.2 × 10 4 s (MBH/10 3 M ).
value of ∼ 2000Ṁ Edd at t 2.3 × 10 6 t 0 . On the other hand, the accretion rate onto the BH becomes as small as ∼ 100Ṁ Edd in the quasi-steady state. Fig. 2 shows the time-averaged radial structure of the angle-integrated mass inflow (red) and outflow (blue) rates in the fiducial case when the system is in a quasisteady state (3 t/(10 6 t 0 ) 4). These rates are defined byṀ andṀ where the net accretion rate is calculated byṀ net = M in −Ṁ out with mass loading factor β = |Ṁ out (r = r max )/Ṁ BH | andṀ BH =Ṁ net (r = r min ). Note thaṫ M out in Eq. (15) is not necessarily the same as the outflow rate of gas escaping to infinity, since it includes all gas with v r > 0. At larger radii (r > 30 r Sch ), the inflow and outflow rates are almost balanced and thus the inflow rate decreases toward the center follow-ingṀ in (r) ∝ r 0.6 . Within r 20 − 30 r Sch , the outflow ceases and the inflowing gas dominates. As a result, the accretion rate is reduced significantly from the inflow rate at r r max to 100Ṁ Edd . We note that the . Time-averaged radial profiles of the net accretion rate (black), the mass inflow rate (red), and the mass outflow rate (blue) in the fiducial case (the Case-H). These profiles are time-averaged over 3 t/(10 6 t0) 4. While the inflow rate balances the outflow rate at larger radii, the inflow rate dominates at the vicinity of the BH. As a result of mass removal by outflows, the inflow rate decrease toward the center followingṀin(r) ∝ r 0.6 , and 5% of the inflowing gas from the outer boundary feeds the BH.
net accretion profile is constant within ∼ 300 r Sch 1 , indicating that the gas over a wide range of the spatial scales sufficiently approaches a quasi-steady state. The reduction of the inflow rate (i.e.,Ṁ in ∝ r p ; p > 0) is generally seen in most numerical simulations of radiatively inefficient accretion flows (RIAF, p ∼ 0.75 − 1; Stone et al. 1999;Igumenshchev & Abramowicz 1999). The power-law index p 0.6 is smaller than that predicted in the convection-dominated accretion flow (CDAF) model Quataert & Gruzinov 2000;Abramowicz et al. 2002;Igumenshchev 2002;Yuan & Narayan 2014). The value of p can be used to characterize the outflow strength. In different accretion models, the outflow strength has been found in the range 0.5 p 1. Fig. 3 shows the two-dimensional distribution of the sound speed and gas density at two different scales for the fiducial case when the simulation terminates. On the larger scale (left panel), the overall flow pattern is governed by convective motion and shows that the gas is inflowing through the equatorial region and a sim-1 Outside this radius, a fountain-like structure forms in the accretion flow and thus the flow does not settle down completely yet.
The time-dependent fountain-like motion carries a small fraction (∼ 3%) of the total gas mass outward from the domain in a 10 6 t 0 duration. Nevertheless, the flow properties in the quasisteady state would not be much different from those we discuss here. Small scale (~120 r Sch ) Large scale (~1000 r Sch ) Density Density Sound speed Sound speed radius (10 9 cm) radius (10 9 cm) 2.0e+8 1.4e+9 1.0e+10 2.0e+8 6.3e+8 2.0e+9 1.0e-9 1.0e-7 1.0e-5 1.0e-8 3.1e-7 1.0e-5 cm/s cm/s g/cm 3 g/cm 3   Figure 4. Time-averaged profiles of the density (left), sound speed (middle), and rotation velocity (right) along the equator (red) and the directions at θ = 30 • (blue) and 150 • (orange), respectively. The density profile along the equator is characterized by a broken power-law profile that follows ρ ∝ r −1/2 in the inner region of r < 30 r Sch and ρ ∝ r −1 in the outer region. The former slope agrees with that expected in a CDAF solution (Quataert & Gruzinov 2000). The accretion flow on the equatorial plane is rotationally (pressure) supported inside (outside) r ∼ 100 r sch , but the gravitational force exerted by the central BH marginally dominates at all radii over the sum of the pressure-gradient force and centrifugal force.
ilar amount of gas is outflowing to the polar regions. On the smaller scale (right panel), the inflows dominate over the outflows owing to strong gravity. The sound speed at the inner and polar region is higher than elsewhere indicating that radiation is trapped in the inner region but can escape through polar region and heat the ambient gas. The accretion flow becomes highly turbulent due to energy generated via viscosity. This behavior is consistent with those found in simulations of RIAFs (Stone et al. 1999;Igumenshchev & Abramowicz 1999). However, in reality, turbulent motion is essentially a three-dimensional (3D) effect that would break the axisymmetric flow structure imposed in our simulations (but see also Igumenshchev et al. 2000). Exploring the nature of multi-dimensional effects is left for future investigation. Fig. 4 presents the time-averaged radial profiles of the density (left), sound speed (middle), and rotational velocity (right). The sound speed and rotational velocity is normalized by the speed of light and the Keplerian velocity, respectively. The three different curves show the profiles along three different directions: θ = 30 • (blue), θ = 90 • (equator; red), and θ = 150 • (orange). The density profile along the equator follows a power-law of ρ ∝ r −1/2 at r 30 r Sch and agrees with those expected in a CDAF solution (Quataert & Gruzinov 2000), while the slope becomes as steep as ρ ∝ r −1 in the outer region. Strong and fast outflows driven by the radiation pressure gradient force create low-density cavities in the bipolar directions, where the outflowing matter obeys ρ ∝ r −1/2 (for > 30 r Sch ). The sound speed profile is approximated by a single power-law of c s ∝ r −1/2 for all the directions, indicating that the pressure gradient force partially balances with the BH gravity; namely the ratio is < 0.3 at r 100 r Sch . The rotational velocity is close to the Keplerian value and the ratio of v rot /v Kep ∼ 0.87 is consistent with the analytical expression for γ = 4/3 (Quataert & Gruzinov 2000) with a gradual increase toward the center. Indeed, the accretion flow on the equatorial plane is rotationally (pressure) supported inside (outside) r ∼ 100 r Sch , but the gravitational force exerted by the central BH marginally dominates at all radii over the sum of the pressuregradient force and centrifugal force.
In summary, due to the presence of strong outflows driven by the pressure gradient force, the inflow rate decreases toward the center asṀ in ∝ r 0.6 and thus the net accretion rate onto the BH is significantly reduced from the gas supply rate injected at the outer boundary. Since the efficiency of mass accretion defined by η [≡ M in (r min )/Ṁ in (r max )] is kept as high as η 0.05 (or mass loading factor β ∼ 19, in our cases) even under strong outflows, the BH feeding rate reaches 100Ṁ Edd in the fiducial case. We note that the global structure of the flow is found to settle down to this quasi-steady state after t 2.3 × 10 6 t 0 , which is ∼ 10 − 100 times longer than the computational timescales addressed in previous (magneto-)radiation hydrodynamical simulation studies (e.g. Ohsuga et al. 2005;Jiang et al. 2014;Sadowski et al. 2015;Jiang et al. 2019).

Simulations with Different Boundary Conditions
Next, we discuss the effect of the outer boundary conditions on the properties of the accretion flows. In the left panel of Fig. 5, we show the time evolution of the net accretion rate (solid) and the mass inflow rate (dotted) for the four cases; the fiducial Case-H, Case-L, Case-N, and Case-I. Depending on the outer boundary conditions, the BH accretion rates range oveṙ M BH 3 − 100Ṁ Edd and each of them approaches a constant value at t 2 × 10 6 t 0 . The reduction factor of the net accretion rate compared to the fiducial case (Case-H) is consistent with the reduction in the inflow rate due to the differences at the outer boundary. Regardless of the different gas supply rates, the accretion efficiency for all cases is η ∼ 0.05. The overall behavior of the mass inflow rate from the outer boundary is also qualitatively similar except in Case-N, which shows periodic oscillations with larger amplitudes. This is because a large-scale circulation crossing the equator disturbs the gas inflowing from ∼ r max through a narrow injection angle.
In the right panel of Fig. 5, we present the timeaveraged radial profiles of the net accretion and mass inflow rates in all four cases. These rates are normalized by the BH accretion rate in the fiducial case and a scale factor F ≡ ρ 0 v r Ω/(ρ 0,Case−H v r,Case−H Ω Case−H ), which is the ratio of the injection density, solid angle, and velocity in each case relative to the fiducial case. Namely, F = 10.0 (Case-L), 27.7 (Case-N), and 9.2 (Case-I), respectively. The radial profiles of these normalized rates clearly show the self-similar nature of the radiation-dominated accretion flows, i.e.,Ṁ in (r) ∝ r p with p ∼ 0.5 − 0.7. Similarly, Fig. 6 shows the selfsimilarity of the gas density (left), sound speed (middle), and rotational velocity (right) both along the equator (solid) and the polar angle θ = 30 • (dotted). Note that the density is normalized by that injected through the outer boundary. Reflecting the scale-free nature of the flows, the profiles of these quantities are similar to those in the fiducial case, without introducing any new characteristic physical scale.
Recently, Kitaki et al. (2021) (hereafter, K21) conducted simulations similar to our Case-I run, where gas  is injected at a constant inflow rate of ∼ 300Ṁ Edd through a narrow funnel around the equator (∆θ = 0.02 π). In their simulations, K21 found that almost all the gas accretes onto the BH with an accretion efficiency of η 0.88 and produces weak outflows with a mass loading factor β 0.14. While the outflows are launched from the interior of the photon-trapping radius (r tr 450 r Sch in their simulations), most of the outflowing gas fails to escape from the outer boundary and falls back onto the disk. In contrast, we find a large mass loading factor for the outflows with β 20 on large scales and thus only a small fraction of the inflowing gas feeds the BH. The large discrepancy in the bulk properties of the flow and mass loading factor is attributed to the following. First, K21 adopted a computational domain that cover a hemisphere and impose equatorial mirror-symmetry of the accretion flow across the equatorial plane. The imposed symmetry tends to suppress global convective motion that crosses the equator and is known to produce outflows coherently through the equator rather than the polar regions (see Li et al. 2013). Under K21's setup, the equatorial outflow dissipates its kinetic energy by forming shocks at the interface with injected gas, which pushes the gas inward, and thus the production of outflows is suppressed. Secondly, they assume a larger viscous parameter of α = 0.1, which is known to cease convective/turbulent motion of RIAFs and thus to reduce the outflow rate (e.g., Inayoshi et al. 2018b). This would yield a flatter radial profile of the mass inflow rate, consistent with the results in K21. Finally, it is worth noting that in K21 the dense disk region is surrounded by a hot optically thin atmosphere set in the initial configuration, while dense optically thick outflows fill the polar regions in our cases. We speculate that the difference might be caused by the initial and outer-boundary conditions. We leave further investigation about this to future work.

Simulations with Radiative Diffusion
Finally, we discuss the effect of radiative diffusion. Note that the simulation setup with radiative diffusion (Case-D) is the same as the fiducial case (Case-H) except that a lower resolution is adopted (N r × N θ = 128 × 128; see Table 1) to avoid a significantly shorter diffusion timescale, i.e., t diff t dyn . For comparison, we also show simulation results of the fiducial model with the same lower grid resolution (Case-H-LR). Fig. 7 presents the time evolution of the mass flow rates (left panel) and the time-averaged radial profiles (right panel) of these rates for the diffusion case (Case-D; dashed curves) and Case-H-LR (solid curves). The red and blue curves correspond to the net accretion rate and the mass inflow rate, respectively. Overall, the mass flow rates in the two cases are similar, although the net accretion rate increases slightly due to energy transport by radiative diffusion. This suggests that energy transport via radiative diffusion plays a minor role in determining the bulk properties of highly optically thick accretion flows at the vicinity of the BH (r < 10 3 r Sch ), where energy advection owing to photon trapping (the angle-averaged trapping radius for these two cases are r tr,eff 50 r Sch , inferred from rest-frame luminosity) Adv. Lum. rest Kinetic(out) Figure 10. Time-averaged radial profiles of three luminosities for the case with radiative diffusion (Case-D): the radiation luminosity in the rest (observed) frame (black), the advection luminosity (blue), the radiation luminosity in the fluid (comoving) frame (red), and the kinetic luminosity of outflows with vr > 0 (purple). The rest-frame luminosity overall follows the advection component and thus it becomes negative due to inward advection (i.e., photon trapping) within r 40 r Sch .
dominates as expected in previous numerical and analytical studies. The left panel of Fig. 8 presents the angular dependence of the density at two radii (r = 10 and 1000 r Sch ) for Case-H-LR (solid) and Case-D (dashed). The profiles in the two cases are similar near the equator, but the density near polar regions at r 10 r Sch in Case-D is an order of magnitude lower than that in Case-H-LR. The low-density regions form due to mass loading into bipolar outflows accelerated by the strong pressure gradient force due to radiation energy leakage from the photospheric surface. As a reference, we overlay the analytical result for the density distribution of a CDAF solution with γ = 4/3 (solid black curve; Quataert & Gruzinov 2000). In both cases, the level of density concentration near the equator is different from the selfsimilar solution predicted by the CDAF solution. That is attributable to the fact that the radial density profile is as steep as ρ ∝ r −1 at large radii and energy leakage due to radiative diffusion at small radii.
The right panel of Fig. 8 shows the radial profiles of the density along the equator (red) and the polar angle θ = 30 • (blue) for the two cases. Along the equator, the profiles are consistent with each other because radiative diffusion does not play an important role in the highly log 10 ( | F adv /F rad | ) outflow inflow Figure 11. Two-dimensional distribution of the ratio of the advection flux (blue for outflow, orange for inflows) to the radiative flux (log 10 (|F adv /F rad |)) in the comoving frame at t = 3.0×10 6 t0 for Case-D. The two panels show the distribution in the whole region (left) and in the inner region within the effective trapping radius (right). The effective trapping radius inferred from the rest frame luminosity is overlaid (∼ 40 r Sch , the grey half circle in the right panel). For reference, the spherical photon-trapping radius is rtr ∼ 6000 r Sch . The radiative flux dominates in the polar regions and the advective flux dominate in the equatorial regions. The 2D distribution does not show a distinct spherical surface characterizing the photon trapping radius.
optically thick dense region. On the other hand, the density profile near the polar region at r 10 r Sch is affected by radiative diffusion because a larger fraction of radiation energy is transported via diffusion into the polar region. The radiative diffusion effect in mildly optically thick layers creates a strong pressure gradient, which pushes the gas outward, and thus the outflowing gas piles up near ∼ 30 r Sch . Nevertheless, the overall impact by radiative diffusion is limited near the inner and polar regions.
In Fig. 9, we present the 2D distribution of the sound speed for Case-H-LR (left) and Case-D (right) when the systems reach a quasi-steady state (t = 2.4 × 10 6 t 0 ). This clearly shows that the radiation energy is mainly transported from the disk and increases the sound speed in the polar regions significantly. We also find that the bipolar outflows are accelerated from the vicinity of the BH horizon owing to strong pressure gradient and lead to a steeper density gradient (see also Fig. 8). Note that although the outflow velocity in the diffusion case reaches 0.1 c 2 , this effect increases the outflow rate by a factor of ∼ 1.3 within the trapping radius compared to the case without radiative diffusion (see the right panel of Fig. 7). The faster outflow at the vicinity of the launching scale does not increase the outflow momentum output at larger radii because the acceleration mechanism owing to radiative diffusion operates only at the small radii.
For the highly accreting gas at rates ofṀ in Ṁ Edd , the radiative and mechanical outputs onto larger scales are crucially important to decide the efficiency of feedback, which in turn governs the gas supply rate onto the nuclear region of interest. In a super-Eddington accreting system, photon trapping is considered to be an important effect to reduce the radiative output and thus moderate the feedback strength. With a steady, spherically symmetric flow, the size of the trapping radius is estimated as r tr /r Sch 5000 (Ṁ in /10 3Ṁ Edd ). However, this simplest assumption is broken in our axisymmetric 2D simulations owing to strong bipolar (non-spherical) outflows, which reduce the mass inflow rate toward the center. Therefore, we need to quantify the nature of photon trapping via relative contributions of different energy transport mechanisms.
In Fig. 10, we show the radial profiles of four kinds of luminosities for Case-D: the radiation luminosity in the rest (observed) frame (black), the advection luminosity (blue), the radiation luminosity in the fluid (comoving) frame (red), and the kinetic luminosity of outflows with v r > 0 (purple) Note that if the gas were static, i.e. v r = 0, the comoving luminosity would be identical to the rest-frame luminosity. In the optically thick flow, the comoving luminosity calculated by the FLD method is almost constant at ∼ 0.6 L Edd over the simulation domain. This value is 2 The maximum outflow speed reaches at most 0.3 c, which corresponds to a Lorentz factor of 1.05. Thus, the relativistic effects do not change our simulation results significantly. Figure 12. The relations between the inflow at the outer boundary (Ṁ0 ≡Ṁin(r = rmax)) and the net BH accretion rate (left panel) and the momentum ejection rate of the outflow (right panel) in our six simulations with different inflow rates. The solid lines in both panels present the fitted linear relations given in Eqs. (20) and (21). Each quantity is time-averaged and normalized by its Eddington value. Note that the results with a lower resolution for Case-D and Case-H-LR yield BH feeding rates nearly half of the high resolution case (Case-H). This small but systematic offset is comparable to those for other cases around the linear relation and thus does not affect our conclusion. The blue-shaded regions mark the 1 σ deviation relative to the scaling relations. We also note that despite radiative diffusion effects in accelerating the outflow in the inner region near the poles (see Fig. 9), the outflow momentum flux at outer boundary for Case-D (compared with Case-H-LR) is not significantly affected due to the fact that radiative diffusion has no effects on large scale. consistent with the radiation luminosity that diffuses out in a spherically symmetric flow (Begelman 1979). On the other hand, the advection luminosity becomes negative (dotted) in the inner region of r 40 r Sch , while it becomes positive (solid) at larger radii owing to strong outflows. The rest-frame luminosity basically follows the advection component and thus it becomes negative within r 40 r Sch , which we can regard as the effective photon-trapping radius. The radiation-dominated outflow also carries kinetic energy and the kinetic luminosity reaches ∼ 3 L Edd , which is consistent with that obtained by an analytical model for an optically thick outflow (see Piran et al. 2015).
The effective photon-trapping radius defined in Fig. 10 is an angle-averaged value. As an alternative way to demonstrate the photon trapping effect, Fig. 11 presents the 2D distribution of the ratio of the advective flux for the inflow component to the comoving radiative flux at t = 3.0×10 6 t 0 on the largest scale (left) and in the inner region of r 40 r Sch (the outflow regions are shown by blue). Overall, the radiative flux dominates the advective flux for outflows near the polar region, as indicated by light blue regions in the right panel. The advection flux dominates in dense and optically thick regions, as indicated by regions with dark blue and orange colors (right panel in Fig. 8 and equatorial region in Fig. 11). Importantly, the 2D distribution does not show a dis-tinct characteristic spherical surface, within which the advection flux overcomes the diffusive flux. It is reasonable that the angle-averaged (effective) trapping radius (30−50 r Sch ) is much smaller than the spherical trapping radius (220 − 7000 r Sch ). As shown in Table 1, the effective trapping radius is always near r tr,eff ∼ 30 − 50 r Sch with a weak dependence on the gas supply rate.

SUBGRID MODELS FOR COSMOLOGICAL SIMULATIONS
In general, the nature of BH feeding and feedback is determined by the global response of the accreting and outflowing gas over a wide range of scales, from the vicinity of BHs (r max ∼ 10 −7 pc for M BH = 10 3 M ) to galactic scales (∼ 1 kpc). The required dynamical range of the computational domain is therefore about 10 orders of magnitude in radius, which makes such simulations computationally challenging. Recently, Lalakos et al. (2022) have conducted MHD simulations of radiatively inefficient accretion flows (i.e., substantially low accretion rates) that cover 5 orders of magnitude in space to study the jet launching mechanism associated with MHD effects and its impact on the gas at larger scales. Similar types of RHD simulations in a suf- = 20 = 4 outflow Figure 13. The angular distribution of the time-averaged radial velocity at the outer-most radius (orange curve) with black dots denoting the outflow components (vr > 0). The outflow velocity in the analytical expression with Eqs. (19) and (20) is shown; N = 0 (red), N = 4 (black), and N = 20 (blue), where N characterizes the degree of anisotropy of the outflow. The case with N = 20 describes the simulation results near the poles well, with the velocity reaching the escape velocity at r = rmax (horizontal gray line). ficiently long term 3 are required to better understand the global interaction between super-Eddington accretion flows and outflows. An alternative approach for this problem is to separate the computation domain into several regions and to treat the dynamics of accretion and feedback crossing those regions in a self-consistent way. For example, Botella et al. (2022) showcase the impact of outflows originating from supercritical accretion, using nested large-scale simulations. However, they limit their study to a specific configuration, and do not provide prescriptions that can be used in other large scale simulations. For this purpose, we provide a subgrid feedback model for super-Eddington accretion flows accompanied by strong outflows that carry mass and momentum.
In Fig. 12, we present the BH accretion rate (left) and the outflow momentum flux (right) as a function of the inflow rate (each point corresponds to the case denoted in the figure). These three quantities show nearly linear correlations as 4 P out = 1.7 +0.9 −0.6 × 10 −3 cṀ 0 whereṀ 0 ≡Ṁ in (r = r max ) and p = 0.5 is adopted. For reference, the values ofṀ 0 for the six cases with different outer boundary conditions are summarized in Table 1. From Eqs. (20) and (21), one can obtain the mass loading factor and the angular-averaged radial velocity of the outflow v wind (≡Ṗ out /Ṁ out ), respectively, as v wind = 1.7 +0.9 −0.6 × 10 −3 c r max 1500 r Sch −0.5 . (23) It is worth noting that the outflow is collimated to the polar regions and the outflow velocity is higher than the escape velocity within the bi-conical regions. Here, we model the angular profile of the outflow velocity as v wind = v wind (N + 1) cos N θ.
where N is a free parameter that characterizes the anisotropy of the outflow, and is determined so that the outflow velocity near the poles exceeds the escape velocity. Fig. 13 presents the radial velocity in Case-H (orange curve; the dots represent the outflowing component), the angle-averaged wind velocity (red curve), and the models from Eq. (24). This figure suggests that the case with N = 20 agrees well with the simulation results near the poles and the outflow velocity near the pole exceeds the escape velocity (grey line). We note that the relations shown above are valid when the mass inflow rate through the outer-most radius is substantially higher than the Eddington value so that the spherical photon-trapping radius for the givenṀ 0 is larger than the simulation domain. Therefore, if the spherical photon-trapping radius is resolved (though this is not the case in our simulations), the outer-most radius would be replaced with the trapping 4 Without imposing the linear relation, the best-fit slopes are slightly different from unity; namelyṀ BH ∝Ṁ 0.95 0 andṖout ∝ M 0.86 0 . Note that Case-D and Case-H-LR are not taken into account for the fitting data, because the two cases are lowerresolution runs and yield nearly twice lower BH accretion rates. We here assume that the power-law distribution of the mass inflow rate (Ṁ ∝ r p ) is maintained until the exterior of the computational domain (presumably, until the spherical photon-trapping radius). radius (r tr = 5ṁ 0 r Sch , whereṁ 0 ≡Ṁ 0 /Ṁ Edd ). Substituting this r tr for r max , we can rewrite the above four equations asṀ BH = 17.1Ṁ Edd ṁ 0 300 P out = 0.51 cṀ Edd ṁ 0 300 β = 17.5 ṁ 0 300 and v wind = 1.7 × 10 −3 c ṁ 0 300 −0.5 . (28) Note that those quantities are supposed to be the values at larger distances outside the trapping radius, where no further production and acceleration of the outflow is considered. Therefore, assuming that the mass accretion rate through sink cells (or onto sink particles) is given byṀ 0 , we propose to apply Eqs. (24)-(28) to large-scale cosmological simulations that hardly resolve the nuclear region but take into account mechanical feedback accompanied by super-Eddington accreting BHs (ṁ 0 10). As an important caveat, our simulations focus on axisymmetric accretion flows without taking into account MHD effects explicitly, but adopting the α-viscosity prescription (α = 0.01) to treat angular momentum transport. With a toroidal magnetic field or multiple poloidal loops configuration in the initial conditions or in injected gas, the accretion flow is dominated by turbulent motion driven by MRI and the saturated magnetic energy density is limited to less than 0.1 − 1% of thermal energy density. As a result, this type of accretion flow is qualitatively similar to that obtained from non-MHD simulations. In contrast, when a poloidal magnetic field configuration is assumed in an accretion disk around a rapidly spinning BH, strong outflows and/or jets tend to be launched owing to amplified magnetic fields near the BH (Tchekhovskoy et al. 2011;McKinney et al. 2014;Sadowski et al. 2015;Tchekhovskoy & Bromberg 2016), which may further suppress the BH feeding rate. We leave a study of strongly magnetized outflows from super-Eddington accretion to a future investigation.

SUMMARY
In this paper, we study the long-term evolution of the global structure of axisymmetric accretion flows onto a black hole fed on large scales at rates substantially higher than the Eddington value, performing 2D RHD simulations. The global structure of the flow settles down to a quasi-steady state in millions of the orbital timescale at the BH event horizon, which is 10 − 100 times longer than that addressed in previous (magneto-) RHD simulation studies.
In the high-accretion optically thick limit, where the radiation energy is efficiently trapped within the flow, the accretion flow becomes adiabatic and thus comprises of turbulent gas in the equatorial region and strong bipolar outflows. As a result, the mass inflow rate decreases toward the center asṀ in ∝ r p with p ∼ 0.5 − 0.7, and only a small fraction of the inflowing gas ends up accreted by BH (see Figs. 2, 5 and 7). In our simulations, super-Eddington accretion is sustained only when a larger amount of gas is supplied from larger radii at 100 − 1000Ṁ Edd . Despite different gas supply rates owing to different outer-boundary conditions, the bulk dynamical properties of the gas inflow in all cases we studied show a self-similarity (see Figs. 5 and 6). Even in the super-Eddington accretion regime, a fraction of the radiation energy in the flow can escape by diffusing out. Energy transport via radiative diffusion accelerates the outflow near the poles in the inner region but does not affect the overall properties of the accretion flow.
As for the strength of the outflows, a small value of p (≤ 0.5 − 0.7) value is expected for supercritical accretion flows since most of the radiation is thought to be trapped near the BH. Simulations with different scales and initial setups show that the value of p is near 0.5−1, with the uncertainties coming from the treatment of radiation feedback and the strength of the gas supply. Further studies of the global properties of outflows and inflows are necessary to determine the outflow strength in a self-consistent manner.
Cosmological and galactic scale simulations cannot resolve the inner regions of the BH accretion flows and need to implement BH growth and feedback with subgrid descriptions (e.g. Massonneau et al. 2022). In this paper, we provide simple formulae (see Eqs. 24-28), which can be used in simulations of super-Eddington accretion at their resolution limit, to predict the BH accretion rate, as well as the mass loading factor, velocity, and momentum injection rate of the outflow, as functions of the mass inflow rate at the resolution limit. This physically motivated implementation of outflow models can help us understand the co-evolution between SMBH and host galaxies in a self-consistent way, connecting large and small scale simulations.