Component of energy flow from supercritical accretion disks around rotating stellar mass black holes

By performing two-dimensional axisymmetric general relativistic radiation magnetohydrodynamics simulations with spin parameter $a^*$ varying from -0.9 to 0.9, we investigate the dependence on the black hole spin of the energy flow from supercritical accretion disk around stellar mass black hole. It is found that optically and geometrically thick disks form near the equatorial plane, and a part of the disk matter is launched from the disk surface in all models. The gas ejection is mainly driven by the radiative force, but magnetic force cannot be neglected, when $|a^*|$ is large. The energy outflow efficiency (total luminosity normalized by $\dot{M}_{\rm in} c^2 $; $\dot{M}_{\rm in}$ and $c$ are the mass accretion rate at the event horizon and the light speed) is larger for rotating black holes than for non-rotating black holes. This is $0.7\%$ for $a^*=-0.7$, $0.3\%$ for $a^*=0$, and $5\%$ for $a^*=0.7$ for $\dot{M}_{\rm in} \sim 100L_{\rm Edd}/c^2$ ($L_{\rm Edd}$ is Eddington luminosity). Also, although the energy is mainly released by radiation when $a^* \sim 0$, the Poynting power increases with $|a^*|$ and exceeds the radiative luminosity for models with $a^* \geq 0.5$ and $a^* \leq -0.7$. The more the black hole rotates, the larger the power ratio of the kinetic luminosity to the isotropic luminosity tends to be. This implies that objects with large (small) power ratio may have rapidly (slowly) rotating black holes. Among ultraluminous X-ray sources, IC342 X-1, is a candidate with a rapidly rotating black hole.


INTRODUCTION
Luminous compact objects such as X-ray binaries, microquasars, active galactic nuclei, and gamma-ray bursts are considered to be powered by the gas accretion on black holes. Black hole accretion is known to exhibit several different modes. If the mass accretion rate comparable to or less that the Eddington limit (L Edd /c 2 where L Edd and c being the Eddington luminosity and speed of light), geometrically thin and optically thick disk (socalled standard accretion disk) is formed (Shakura & Sunyaev 1973). When the mass accretion rate is much smaller than the Eddington limit, radiatively inefficient accretion flows (RIAF), which consists of the hot rarefied plasma appear (Narayan & Yi 1994, 1995aAbramowicz et al. 1995). Supercritical accretion flows, on the other hand, are geometrically and optically thick disks, whose luminosity exceeds the Eddington luminosity and form when the accretion rate exceeds the Eddington limit.
The supercritical disks are thought to be the energy source for extremely bright compact objects. Ultraluminous X-ray sources (ULXs) are off-nuclear sources with X-ray luminosities of 10 39 erg s −1 and may be undergoing supercritical accretion into stellar mass black holes. Other possibilities have been suggested as subcritical accretion onto intermediate-mass black holes and supercritical accretion onto neutron stars. For the ULX, for which X-ray pulses have been detected, a neutron star is believed to exist in the center, but the central objects of the other ULXs are unknown. In this study, we investigate the case of supercritical accretion onto the stellar mass black hole.
Radiation hydrodynamics (RHD) is necessary in the study of supercritical accretion because the interaction between the radiation field and the material is important. Ohsuga et al. (2005) was the first to succeed in reproducing the steady-state of supercritical accretion flow using RHD simulations (see also Eggum et al. 1985Eggum et al. , 1988. Subsequently, the radiation magnetohydrodynamics (RMHD) simulations were developed in order to account for the effects of magnetic fields, which are the origin of disk viscosity (e.g., Ohsuga et al. 2009;Ohsuga & Mineshige 2011;Jiang et al. 2014a). It has been reported by these simulations that a part of the disk matter is blown away by the strong radiation pressure force. The bubble structure detected in some ULXs (Pakull & Mirioni 2002Grisé et al. 2006;Abolmasov et al. 2007;Cseh et al. 2012) may have been formed by these radiatively driven outflows (Hashizume et al. 2015).
However, a non-rotating black hole is supposed in these simulations and supercritical accretion on a rotating black hole is not yet well understood even though the black hole may be rotating. ULXs would be kind of X-ray binaries (Komossa & Schulz 1998), but there is a possibility of an isolated black hole that plunges into an interstellar matter (Mii & Totani 2005;Kobayashi et al. 2019). In the latter case, the rotation axis of the black hole can be misaligned with that of the disk. It is also possible that the rotation direction of the black hole and the rotation direction of the disk are opposite. It is known that the SS433 has a precessed jet (Crampton et al. 1980;Cherepashchuk 1981) and one possible interpretation of this is a misaligned system.
In order to investigate the effect of black hole spin, it is necessary to perform RMHD simulations that incorporate general relativistic effects, i.e., general relativistic (GR) RMHD simulations. GR-RMHD simulations of supercritical accretion flows around the non-rotating black hole have been performed by Sądowski et al. (2014) (see also Sądowski & Narayan 2015a). Sądowski et al. (2014) has been carried out the simulation for the case of a rotating black hole, but only the spin parameter of 0.9 is employed (see also Takahashi et al. 2016).
In the present study, we investigate the effect of the black hole spin on the supercritical accretion flows around the black holes by performing two-dimensional GR-RMHD simulations with spin parameter a * varying from −0.9 (retrograde) to 0.9 (prograde). The disk luminosity, the power of the ejected gas, their dependence on the spin parameter in the quasi-steady state are investigated and compared with the observations of some ULXs. The paper is organized as follows: we describe the methods of calculations and initial conditions in Section 2. In Section 3, we show our simulation results. Section 4 is devoted to discussion. This section contains comparisons with observations of ULXs. We also discuss about the magnetically arrested disks (MADs) (e.g., Igumenshchev et al. 2003;Narayan et al. 2003;Tchekhovskoy et al. 2011McKinney et al. 2012), since our present study focuses on the standard and normal evolution (SANE) disks and does not deal with the MADs. Finally, we summarize our main findings in Section 5.

BASIC EQUATIONS AND NUMERICAL METHOD
We numerically solve GR-RMHD equations using the code developed by Takahashi et al. (2016), in which the zeroth and first moment equations of the radiation are solved with applying a M1 closure (Levermore 1984;Kanno et al. 2013;Sądowski et al. 2013b). In the following, the Greek suffixes indicate spacetime components, and the Latin suffixes indicate space components. Length and time are normalized by gravitational radius (r g = GM/c 2 ) and its light crossing time (t g = r g /c), where G is the gravitational constant and M is the black hole mass. Hereafter c and G are taken as unity.
The mass conservation equation is given by where ρ is the proper mass density and u µ is the fluid four velocity. The energy momentum conservation equations for the magnetofluid and the radiation are given by and Here, the energy momentum tensor for fluid T µν is given by where e is the gas internal energy, p g is the gas pressure, δ ν µ is the Kronecker delta. The gas internal energy is related to the gas pressure by e = (Γ − 1) −1 p g with Γ = 5/3 being the specific heat ratio. The energy momentum tensor for electromagnetic field M µν is given by where b µ is the magnetic four vector and p m = b µ b µ /2 is the magnetic pressure. The magnetic four vector b µ is related to its three vectors B i through where g µν is the metric. The radiation energy momentum tensor is given by where p rad is the radiation pressure and u µ rad is the radiation frame's four velocity. The radiation four force G µ is given by where κ abs and κ sca are the opacities for free-free absorption and Thomson-scattering, κ abs = 6.4 × 10 22 ρT −3.5 where T gas is the gas temperature, which is related to the gas pressure as Here k B is the Boltzmann constant, µ = 0.5 is the mean molecular weight, and m p is the proton mass. The blackbody intensity is given by B = a rad T 4 gas /4π with a rad being the radiation constant. We included the thermal Comptonization as follows: where T e is the electron temperature,Ê rad is the comoving frame radiation energy density, T rad = (Ê rad /a rad ) 1/4 is the radiation temperature, and m e is the electron rest mass (Sądowski & Narayan 2015b). We take T e = T g for simplicity. The induction equation is given by where g = det(g µν ). We set Γ = 5/3 for a gas, and 4/3 for a radiation (see Equation (8)). The specific heat ratio of the mixed fluid is determined by their temperature. The assumption of Γ = 5/3 could be violated for the relativistically hot gas. We find that the T g is less than m p c 2 /k B in the whole simulation box after simulation. We believe the Γ = 5/3 for the gas is a good approximation, even though more realistic equation of state is needed to implement for a detailed analysis (Mignone & McKinney 2007;Chattopadhyay & Ryu 2009;Sądowski et al. 2017;Dihingia et al. 2019). We solve these equations in polar coordinate (t, r, θ, φ) with the Kerr-Schild metric by assuming axisymmetry with respect to the rotation axis (θ = 0 and π). The computational domain consists of r = [R in , R out ] = [r H , 250r g ], θ = [0, π], where r H = 1 + [1 + (a * ) 2 ] 1/2 is a horizon radius with a * being the normalized spin parameter of black hole. Number of numerical grid points is (N r , N θ , N φ ) = (264, 264, 1).
A grid spacing exponentially increases with radius in the radial direction. The radius of the j-th radial grid is given as r(j) = exp [ln r H + j ln(R out /R in )/N r ]. A polar angle of k-th polar grid is given by θ(k) = kπ/N θ + 0.25 sin [2πk/N θ ]. We adopted the outgoing boundary at inner and outer radius, and the reflective boundary is adopted at θ = 0 and π.
We start simulations from an equilibrium torus (Fishbone & Moncrief 1976) in which weak poloidal magnetic fields are initially embedded. In the torus, the magnetic flux vector A φ is proportional to ρ and the plasma beta β = (p g + p rad )/p m is set to be 100 at most. The inner edge of the initial torus is situated at r = 20r g , while its pressure maximum is situated at r = 33r g . We take ρ 0 = 1.4 × 10 −2 g cm −3 as maximum mass density.
We employ M = 10M throughout the present study, where M is the solar mass. In total, we calculate nine models with a * = 0, ±0.3, ±0.5, ±0.7, ±0.9 until t = 5000t g . We note that the difference in the structure of the initial torus due to differences in a * is negligibly small.

Overview of Simulations
At the beginning, we introduce in Figure 1 the global structure of supercritical accretion flows, jets, and winds (we will define later) for the case with a * = 0.7 (top), 0 (middle), and −0.7 (bottom). The color contour on the left side represents the radiative energy density, and the solid lines indicate the magnetic field lines. The radiation energy density is expressed as where n α = (−α, 0) is the normal observer's four velosity with α = (−g tt ) −1/2 being the lapse function. The mass density is shown by the color contour on the right side. The streamlines, which are calculated from the velocity field in the zero angular momentum observer (ZAMO) frame, are overlaid with the arrows. In this figure, all physical quantities are time-averaged over t = 3000 − 5000t g . Here we note that the flow is basically quasisteady during this time-averaged period (see below), but there are some fluctuations. In particular, somewhat larger fluctuations appear in the region where the mass density is much less than 10 −2 ρ 0 (where the gas flows outward). Hereafter, the time averaged quantities over t = 3000 − 5000t g is denoted with . The white dashed lines show the photosphere (τ tot = 1), where τ tot is the total optical depth measured from the polar axis, Figure 1. The structure of time-averaged (t = 3000 − 5000tg) accretion disks, jets, and winds for a * = 0.7 (top), a * = 0 (middle), a * = −0.7 (bottom). In each panel, the left half shows the radiative energy density (colors) and the magnetic field lines (contours), and the right half shows the mass density (colors), Be = 0(magenta solid contours), Be = 0.05 (cyan solid contours), the photosphere (white dashed contours) and the mass flux direction (stream lines). Here, γ = αu t is the Lorentz factor. The magenta and cyan solid lines indicate the surface where Bernoulli parameter, B e , is 0 and 0.05, respectively. Here B e is given by (Sądowski et al. 2013a;Sądowski & Narayan 2015a. Note that B e ≥ 0 is the condition for the steady flow to reach infinity, and the flows with B e ≥ 0.05 can have ∼ 30% of the speed of light at infinity. In all models, we distinguish the whole structure into three regions. One of them is the disk region with B e < 0, and another is the jet region with B e ≥ 0.05 around the rotation axis. The remaining one is the wind region, which appears between these two regions (0 ≤ B e < 0.05). We must note that the boundary between the jet and the wind is not well defined, but they are smoothly connected. Thus, the jet boundary defined as B e ≥ 0.05 might not agree with the observed one. Despite this, we think the definition of the jet boundary by B e ≥ 0.05 is useful to distinguish the fast jet from the slower wind. In the disk region, we can see that the gas flows toward the black hole with the turbulent motion. The large number of photons is accumulated since the disk is very optically thick. In the jet region, the rarefied matter is blown away in all models except very vicinity of the rotation axis when a * = 0. The radial velocity, which is measured in the ZAMO frame, in that region reaches a maximum of ∼ 60% of the speed of light for the case with a * = 0.7 and a * = −0.7, although it is less than ∼ 30% of the speed of light when a * = 0. The magnetic field lines in the jet region have a radial structure along the gas flow although the complex structure arises in the disk region. In the wind region, where the Bernoulli parameter is roughly 0 ≤ B e < 0.05, relatively dense matter (ρ/ρ 0 ∼ 10 −5 − 10 −2 ) moves outward at ∼ 10 − 20% of the speed of light, regardless of a * . Although not shown in Figure  1, accretion disk, winds, and jets also appear in all other models.
In the present study, only the model with a * = 0.9 eventually reaches MAD limit  and shows violent behavior such as a rapid increase or decrease in accretion rate and outflow rate. Therefore, we mainly compare the models with a * = ±0.7 with the model with a * = 0.

Mass Flow and Luminosity
Next, we measure the mass accretion rate at the event horizon, asṀ The mass outflow rate is calculated at r = 200r g aṡ where θ Be (θ Be ) is the angle to identify the jet and wind region for θ < π/2 (θ > π/2). We set the angle of B e = 0.05, which indicates the boundary of jet and wind regions. We also set the angle of B e = 0, which separates the wind and disk regions. The angle of B e = 0.05, which corresponds to the opening angle of the jets, is slightly larger for the models with a * ± 0.7 than the model with a * = 0. The jet opening angle is around 35 • (a * = −0.7), 20 • (a * = 0), and 30 • (a * = 0.7) at r = 200r g . Furthermore, the opening angle of the jet roughly corresponds to the opening angle of the photosphere. (see Figure 1) The total luminosity is evaluated at r = 200r g by where L rad is the radiative luminosity, L mag is the Poynting power (luminosity, which is calculated by Poynting flux), L kin is the kinetic power (luminosity, which is calculated by kinetic flux), at r = 200r g , Here −K t r is the energy momentum tensor of the fluid minus the thermal energy, the rest mass energy, and the potential energy components, , and becomes the kinetic flux in the non-relativistic limit. Figure 2 shows the time evolution of the mass accretion rate (top), the mass outflow rate (middle), and the total luminosity (bottom) for a * = 0, ±0.7 normalized by Eddington luminosity given by As shown in the top panel, the mass accretion rate rapidly increases at around t = 1000t g since the gas from the initial torus reaches the event horizon in all models. At t ≥ 1000t g , a quasi-steady supercritical disk is formed, and the accretion rate always exceeds the Eddington limit (L Edd ). The average mass accretion rate is ∼ 470, 240, and 150L Edd for a * = −0.7, 0, and 0.7, respectively. Although there are some fluctuations, the mass outflow rate remains above L Edd in all models at t ≥ 3000t g , meaning that the powerful winds are launched. It is also found that the outflow rate in the region of B e ≥ 0 (solid line, jet+wind region) is much higher than that in the region of B e ≥ 0.05 (dashed line, jet region). Thus, we find that the matter is mainly blown passing through the wind region. We can see in the bottom panel that the total luminosity in the region of B e ≥ 0 (solid line) is kept above L Edd in all models at t ≥ 3000t g . This panel indicates that the total luminosity for B e ≥ 0.05 (dashed line) is slightly smaller than that in the region of B e ≥ 0 (solid line), but not by much. Thus, we can conclude that the energy is mainly released through the jet region unlike the ejecting matter. Here, the luminosity, accretion rate, and mass outflow rate are measured at 200r g . We note that they are almost constant with radius far from the black hole for r > 150r g , and the mass accretion rate, along with outflow rate, reach steady states at t ≥ 3000r g . Thus, these results do not depend on the choice of the radius for r > 150r g . Hereafter, we take time average between t = 3000t g and 5000t g since we can recognize that the disk, wind, and jet reach a steady state.
In Figure 3, we plot Ṁ in /L Edd (panel a), L tot /L Edd (panel b1), L tot / Ṁ in (panel b2) as a function of a * from top to bottom. The mass accretion rate tends to decrease with increasing a * (see panel a). For instance, the mass accretion rate is ∼ 470L Edd for a * = −0.7, ∼ 240L Edd for a * = 0, and ∼ 150L Edd for a * = 0.7.
It is found in the panel b1 that the total luminosity tends to increase with |a * |. We attribute this result to the effect of the Blandford-Znajek (BZ) mechanism (Blandford & Znajek 1977) (we will discuss later). Since the total luminosity is more sensitive to the spin parameter in the range of a * > 0, L tot is higher for the prograde case than the retrograde case with same |a * |. The total luminosity in the region of B e ≥ 0 for a * = 0.7 is ∼ 9L Edd , that is higher than the luminosity for a * = 0 (∼ 3L Edd ) and for a * = −0.7 (∼ 5L Edd ). As noted before, the opening angle of the jet is larger for a rapidly rotating case, but the size of the solid angle of the jet is only 1.5 times different between a * = 0 and a * = 0.7. Thus, the higher luminosity for a rapidly ro- tating black hole originates from the higher flux rather than the larger solid angle. Note that the difference between the luminosity in the region of B e ≥ 0 (red) and that in the region of B e ≥ 0.05 (blue) is very small, which means that the energy is mainly released through the jet region as we have already mentioned in Figure 2.
Since the mass accretion rate tends to be higher for a * < 0 than for a * > 0, the enhancement of the energy outflow efficiency ( L tot / Ṁ in ) is more pronounced when a * is large (see panel b2). The panel b2 shows that the energy outflow efficiency becomes 0.05 only when a * ≥ 0.7. Figure 4 shows the mass outflow rate Ṁ out (panel c1 and c2), and the mass outflow efficiency Ṁ out / Ṁ in (panel c3 and c4). The outflow rate is measured where B e ≥ 0.05 (panel c1 and c3), and where B e ≥ 0 (c2 and c4). As we have mentioned in Figure 2, the mass is mainly ejected through the wind region, so that the outflow rate for B e ≥ 0.05 (blue line) is much lower than that for B e ≥ 0 (red line) (see panels c1, c2, c3, and c4). It is found that Ṁ out and Ṁ out / Ṁ in for B e ≥ 0.05 tend to increase with increasing |a * |. Such a trend may be related to the BZ effect. On the other hand, although the mass outflow rate in the region of B e ≥ 0 is very high when a * = 0.9, there is no clear dependence of Ṁ out on a * in panel c2. The ratio of Ṁ out to Ṁ in slightly increases with a * (see panel c4). The mechanism of gas ejection from the surface of the accretion disk is mainly radiative force, but magnetic force also comes into play when |a * | is large (we will discuss in Section 3.3).
The top panel in Figure 5 represents the radiative flux. Color contour shows the magnitude of the radial component of the radiative flux |R t r | and arrows denotes the direction of the radiative flux on the poloidal plane. We find that the magnitude and direction of the radiative flux do not depend much on the black hole spin. In all models, the radiative flux near the equatorial plane is the highest, and the radiative flux is also strong in the jet region near the rotation axis. The direction of the radiative flux is basically outward, but in the disk, there are places where the flux turns inward or toward the equatorial plane. This is thought to be due to the turbulent motion of the optically thick matter within the disk.
In the bottom three panels of Figure 5, we show the magnitude of the radial component of the Poynting flux  t r | by color contour and the direction of the Poynting flux by arrows. As with |R t r |, |M t r | is higher in the disk region in all models. This is thought to be due to the amplification of the magnetic field inside the disk. In addition, the enhanced outward Poynting flux appears in the regions of jet and wind extending diagonally upward from the black hole in models with a * = ±0.7, except very vicinity of the rotation axis. Although the outward Poynting flux in the wind region is slightly enhanced in the model with a * = 0, the magnitude of it is much smaller than the other models. Therefore, the enhanced Poynting fluxes in the jet and wind regions that appear in models of rotating black holes may be due to the BZ mechanism. Indeed, we can confirm that the resulting Poynting flux at the event horizon, −M r t (r = r H ), is consistent with the BZ flux F BZ | r=rH derived in McKinney & Gammie (2004), Here Ω H ≡ a * /2r H and ω are the rotation frequency of the black hole and the magnetic field at the event horizon. For example, the magnitudes calculated at θ = π/4 are ∼ 6 × 10 −4 ρ 0 (−M r t ) and ∼ 4 × 10 −4 ρ 0 (F BZ ) for a * ± 0.7. We note that direction of black hole spin, accretion flow, and magnetic field rotation are aligned even for a * < 0 at the event horizon.
It is clearly understood in Figure 6 that the BZ effect enhances the Poynting flux. The top panel in the figure shows the energy outflow efficiency via the radiative flux, L rad / Ṁ in , the Poynting flux, L mag / Ṁ in , and the kinetic flux, L kin / Ṁ in , as a function of a * . In the bottom panel, we also plot L rad , L mag , and L kin normalized by L tot . Here, L rad , L mag , and L kin in the figure indicates the luminosities in the jet regions (see Equations (21)-(23)).
The luminosity due to the Poynting flux tends to increase with an increase of |a * |. Although L mag / Ṁ in is only ∼ 10 −4 in the model with a * = 0, it exceeds ∼ 0.01 in the case of a * = −0.9 and a * ≥ 0.5. The luminosity calculated from the Poynting flux is responsible for more than 50% of the total luminosity for a * ≥ 0.5 and a * ≤ −0.7 (see bottom panel).
The energy outflow efficiency via the radiative flux exhibits a slightly increasing trend with increasing a * , but is not as sensitive as that via the Poynting flux. Although L rad is higher than L mag for a * ∼ 0, we find L rad < L mag for the case of a * ≥ 0.5 and a * ≤ −0.7. As a result, the radiative luminosity becomes dominant in the range of a * = −0.5 and a * = 0.3 as shown in the bottom panel. Also, the kinetic power is enhanced in the models with large a * . However, L kin is never the most dominant at r = 200r g Our results suggest that the black hole spin enhances the energy flux, especially the Poynting flux. The cause is probably the BZ effect. Supercritical flows release the energy mainly by radiation for slowly rotating black hole, and by Poynting flux for rapidly rotating black hole.

Mass Outflow Rate
In the present simulations, the matter ejected from the disk surface is mainly blown away through the wind region as we have mentioned in Section 3.2. Here, we discuss about forces, which launch the matter from the disk.
The general relativistic force per unit mass, which is modified to more precisely describe the forces related to the magnetic field in the previous study (Moller & Sadowski 2015), is shown below. Assuming a steady state, the equation of motion is obtained from the energy momentum conservation for magnetofluids as follows: Figure 7. The time averaged force per unit mass in z(= r cos θ) direction as a function of z for a * = −0.7 (left), a * = 0 (middle), a * = 0.7 (right). Red, blue, and black lines show the radiative force, the magnetic force, and the total force of radiation plus magnetic force. Black dashed lines show gravitational force, and dashed lines indicate force in the negative direction. The radius at which these are measured is R(= r sin θ) = 50rg. Forces are measured in Boyer-Lindquist coordinates (BL. co.) and in the ZAMO frame.
where w = ρ + e + p g + 2p m denotes the relativistic enthalpy, and Γ λ µν is the christoffel symbol. Based on this equation, the gas pressure gradient force, the radiative force, the magnetic force, and the metric force are defined by All terms related to the magnetic field are included in the magnetic force, f mag ν . We use the gravity in the non-relativistic limit, since we don't investigate the very vicinity of the black hole in this section. Figure 7 shows the component of the force perpendicular to the equatorial plane as a function of z(= r cos θ) at R(= r sin θ) = 50r g . Here, these forces are transformed to the ZAMO frame in Boyer-Lindquist coordinates. At the deep inside the disk (z 20r g ), the radiative force is more dominant than the magnetic force and is balanced with gravity in all models. Around the disk surface (z ∼ 20 − 30r g and B e ∼ 0), the radiative force is stronger than the magnetic force, and the sum of the radiative and magnetic forces (f tot ) slightly exceeds gravity. Thus, the gas is launched from the disk surface . Here we note that the gas pressure gradient force and metric force are negligibly small compared to radiative force or magnetic force.
Since f tot is stronger than the gravity force even at the region above the disk (z 30r g ), the gas is accelerated. However, whether the radiative force or the magnetic force is stronger depends on the black hole spin. The magnetic force steeply increases from inside the disk z < 20r g to wind region, and exceeds the radiative force at around z 35r g for the models with a * = ±0.7.

Astrophysical Implications
Now, we compare our results with ULX observations. Figure 8 shows the spin dependence of isotropic X-ray luminosity L iso rad , where As Kitaki et al. (2021) noted, the ratio of kinetic luminosity L kin and isotropic X-ray luminosity is ∼ 0.04 − 0.14 for ULX Holmberg II X-1, and ∼ 0.15−0.3 for ULX IC342 X-1 (Kaaret et al. 2004;Abolmasov et al. 2007;Shidatsu et al. 2017). These values are specified with red and blue hatches in Figure 8. Even though model a * = −0.3 is an exception, the value of ULX Holmberg II X-1 is consistent with our models of between a * = −0.7 and 0.3. Also, there is not an inconsistency between ULX IC342 X-1 and our models of a * ≥ 0.5 and a * = −0.9. These results indicate that the Holmberg X II X-1 originates from a slowly, or mildly rotating black hole, and IC342 X-1 contains relatively a rapidly rotating black hole. Here, we note that the ratio of L kin /L iso rad for Holmberg X II X-1 is three times larger based on the jet power reported by Cseh et al. (2014). In this case, a rotating black hole is preferred as the central object.
Regarding V404 Cygni, the ratio of kinetic to isotropic X-ray luminosities is 10 −4 −10 −3 , which is much smaller than our results. The X-ray luminosity at the outburst is about 10 39 erg s −1 (Motta et al. 2017) and jet power is 10 35 − 10 36 erg s −1 (Tetarenko et al. 2019) in V404 Cygni. This X-ray luminosity is close to the Eddington luminosity of the black hole with 10M (Shahbaz et al. 1994;Khargharia et al. 2010) so that the radiation force is not strong enough to efficiently accelerate the gas. In the case of ULX Holmberg II X-1 and ULX IC342 X-1, the X-ray luminosity is about 10 times higher than the Eddington luminosity, so the strong radiation can help the mass ejection.
The L kin /L iso rad of GRS 1915+105 at high state is, on the other hand, about 20 (Done et al. 2004 ;Fender & Belloni 2004), which is much higher than our results'. This fact indicates that the GRS 1915+105 might originate from the rapidly rotating black hole (McClintock et al. 2006), or from an efficient energy conversion mechanism. If the radiation energy and/or the magnetic field energy are converted to the kinetic energy at a large distance, the large value of L kin /L iso rad might be explained since we measure L kin and L iso rad at r = 200r g in the present study. The energy conversion from magnetic to kinetic energy could occur at a large distance (Vlahakis & Königl 2003;Sapountzis & Janiuk 2019). Also, it has been reported that the radiation energy is converted to the kinetic energy in the distant region (Sądowski & Narayan 2015a). If this is the case, even the supercritical disk is identified as an object exhibiting strong jets/winds. If the energy conversion efficiencies from the radiation energy to the kinetic energy, and from the magnetic energy to the kinetic energy are understood, it may be possible to restrict the spin parameter from the ratio of L kin /L iso rad . For example, if the conversion from radiation energy to kinetic energy is inefficient, sources with large L kin /L iso rad could not be explained by a radiatively driven model, but a magnetically accelerated model would be preferred. Such a source would originate from the rapidly rotating black hole. For a more detailed study, the simulations with a larger computational domain are required.
The large L kin /L iso rad may be due to observer's viewing angle (angle between the rotation axis and the line of sight). The X-ray luminosity estimated from the observation is approximately L iso rad for a face-on observer, but is expected to decrease as the observer's viewing angle increases. Indeed, a relatively large angle, ∼ 60 − 70 • , is suggested by observations of GRS 1915+105 (Mirabel & Rodríguez 1994;Fender et al. 1999;Blum et al. 2009). To elucidate this, it is necessary to know the large-scale structure. The simulations with a larger computational domain are left as important future work.

Limitations and Future Work
We showed that the energy outflow efficiency L tot /Ṁ in of supercritical accretion disk increases with |a * |. A similar result is obtained by Sądowski et al. (2014). The L tot /Ṁ in is 0.003 for a * = 0 in our model is comparable to that of Sądowski et al. (2014). For a * = 0.9, our result is ∼ 10 times larger than that of Sądowski et al. (2014). This big difference would be due to an occurrence of transition from SANE to MAD in our simulation. We computed the normalized magnetic flux threading the horizon φ BH , and confirmed that the flux exceeds the critical value φ BH = 50 at t ∼ 4500r g , only for a * = 0.9, before the end of the simulations. Thus, the disk is in MAD state and the strong jet emerges at t 4500r g Narayan et al. 2022), leading to an increase of L tot /Ṁ in . When the time domain of t > 4500r g is excluded, L tot /Ṁ in is 0.09 for a * = 0.9. This value is consistent with that by Sądowski et al. (2014), L tot /Ṁ in ∼ 0.02. We also note that the transition occurs at about t = 4500t g in our model, which is earlier than Sądowski et al. (2014) (t = 12000t g ). We think the reason might be the plasma beta at the initial state, which is smaller in our model than in theirs. The stronger magnetic field at the initial state leads to a faster transition to MAD. We aim to study the MAD state in the supercritical accretion in the forthcoming paper, but our discussion is concentrated on mainly the SANE state in this manuscript.
Enhancement of the energy outflow efficiency via the black holes spin is also shown by GR-MHD calculations, which is used in the study of disks with very small accretion rates Narayan et al. 2022). However, the absolute values of the energy outflow efficiencies are very different, while our results are less than 30% at maximum, and their results can exceed 100%. This may be due to the difference in the mass accretion rate and/or magnetic field strength. They investigate the MAD, which has strong magnetic fields, while our simulations investigate SANE. The problem can be solved by performing simulations of the MAD with the supercritical accretion rate. Such simulations are left as an important future work.
Regarding the radiation luminosity, L rad /Ṁ in is almost consistent, but is slightly larger (by a factor of 3) than that of Sądowski et al. (2014). The reason why their results are slightly smaller, may be due to a higher mass accretion rate (by a factor of 10) than ours. Most of the photons generated inside the disk are swallowed by the black hole due to the photon trapping effect. This mechanism is prominent for a higher mass accretion rate. Even though the mass accretion rate increases, the radiation luminosity does not rise so much in supercritical accretion disks. The resulting L rad /Ṁ in would be smaller for a higher mass accretion rate. Further investigation about the dependency of the efficiency on the accretion rate is needed for more details.
We start the simulations from the rotating equilibrium torus with the outer edge r < 100r g and study the accretion disk near the black hole and the gas ejection from it. However, Kitaki et al. (2018) showed by the two-dimensional RHD calculations that the gas is ejected from just inside the trapping radius, which is the boundary between the supercritical and standard disks, ∼ (Ṁ in /L Edd )r g . Thus, it is necessary to perform long term GR-RMHD simulations with a larger compu-tational domain in order to reveal overall structure of the supercritical flows.
The investigation of dependence of the initial magnetic field is also a future work. Although the single loop of poloidal magnetic field is used as the initial magnetic field in the present study, it has been reported that the magnetic flux at the event horizon differs depending on the initial magnetic field Sądowski et al. 2013aSądowski et al. , 2015. Moreover, we need to perform three-dimensional calculations in order to understand more realistic inflow-outflow structure, although the two-dimensional simulations are performed in this study. For the cases of the supercritical disks, it has been reported that higher variability in the radiative flux appears in the two-dimensional simulations . Also the disk dynamo does not work in two-dimensions, so that it affects on the accretion rate and radiative efficiency. Global three-dimensional simulations can solve these problem.
We use M1 closure to approximately calculate the radiation field in this study. It has been pointed out that the M1 closure does not give correct results where the radiation fields is highly anisotropic region like the vicinity of the rotation axis (Asahina et al. 2020). To solve this problem, we need to solve the radiation transfer equation without using the M1 closure approximation (Jiang et al. 2014b(Jiang et al. , 2019aAsahina et al. 2020).
In this work, we evaluated the electron temperature T e appeared in the opacity by assuming T e = T g . It is well known that this assumption is violated for lowdensity flows, e.g, low luminosity accretion disks (RIAF) (Narayan et al. 1995;Manmoto et al. 1997;Mościbrodzka et al. 2016Mościbrodzka et al. , 2021. The supercritical accretion disk targeted in this study is much denser than RIAF, so that the Coulomb interaction would work efficiently. The electron temperature is close to the ion temperature. In fact, we confirmed from our numerical results that the time scale of the Coulomb interaction is shorter than the dynamical time scale of the fluid in the entire computational domain, excluding the very vicinity of the rotation axis. The electron temperature is then determined by the balance between radiative cooling, viscous heating, adiabatic heating/cooling, and advective heating/cooling. Although ions and electrons are considered to decouple near the axis, the internal energy of the gas is sufficiently smaller than the magnetic and radiation energy. Thus, the effect of this decoupling on the dynamics would be negligible, despite the determination of the electron temperature being important to evaluate the emergent spectra. Two-temperature GR-RMHD simulation has already modelled the low luminosity disk (Sądowski et al. 2017;Ryan et al. 2018). We leave the more detailed analysis as an important future work.
Radiation transfer calculations using the results of two-temperature simulations are also an important future work to restrict the black hole spin. This is because our simulations show that the strength of the Poynting flux around the rotation axis depends on the black hole spin. Since synchrotrons and synchrotron self-comptons are very sensitive to magnetic fields, differences in black hole spin are expected to induce differences in the spectral energy distribution. Polarization degree and polarization distributions, which are strongly affected by Faraday rotation and conversion, may also depend on black hole spin (Tsunetoe et al. 2020(Tsunetoe et al. , 2021Mościbrodzka et al. 2017;Event Horizon Telescope Collaboration et al. 2021). Polarized radiation transfer calculations are also important future work.

SUMMARY
We have performed two-dimensional axisymmetric GR-RMHD simulations of the supercritical flows in the SANE state around stellar mass black holes with a * varying from −0.9 to 0.9. Our simulations show optically and geometrically thick disks near the equatorial plane and the powerful gas ejection from the disk surface in all models. The gas ejection is mainly induced by the radiative force, but magnetic force also works to accelerate when |a * | is large.
The energy outflow efficiency, L tot / Ṁ in , is larger for rotating black holes than for the non-rotating black hole. In particular, this is found to be larger for the prograde models. We find L tot / Ṁ in = 0.7% for a * = −0.7, 0.3% for a * = 0, and 5% for a * = 0.7 foṙ M in ∼ 100L Edd . Although the disks around the nonrotating black holes release the energy mainly by radiation, the Poynting power is enhanced more than the other luminosities in the case of a rotating black hole.
The power ratio of the kinetic and isotropic radiative luminosity L kin /L iso rad tends to be larger the case for rotating the black hole than the case for not rotating black hole. This result suggests that objects with large L kin /L iso rad ( 0.15), such as ULX IC342 X-1, may have a rapidly rotating black hole. On the other hand, a slowly rotating or non-rotating black hole would be suitable for an central object for ULXs with small L kin /L iso rad ( 0.15), such as ULX Holmberg II X-1. For more detailed comparison with observations, GR-RMHD simulations with larger computational domain and radiation transfer simulations are needed.
Numerical computations were performed with computational resources provided by the Multidisciplinary Cooperative Research Program in the Center for Computational Sciences, University of Tsukuba, Oakforest-PACS operated by the Joint Center for Advanced High-Performance Computing (JCAHPC) and, Cray XC50 at the Center for Computational Astrophysics (CfCA) of the National Astronomical Observatory of Japan (NAOJ). This work was supported by JST, the establishment of university fellowships towards the creation of science technology innovation,