Study of Advective Energy Transport in the Inflow and the Outflow of Super-Eddington Accretion Flows

We consider a steady($\partial/\partial t=0$) and axisymmetric($\partial/\partial\phi=0$) accretion flow onto a non-spinning BH at the origin of spherical coordinates ($r\theta\phi$), neglecting self-gravity of accreted matter. The equations of continuity and motion can be respectively written as

where $\Psi$ is the gravitational potential, T the tensor of viscous stress, $\boldsymbol{F}$ the radiative flux, and $\kappa$ the total opacity. We assume that for the accretion flow, the $r\phi$-component of the viscous stress tensor, $t_{r\phi}$, is dominant, and adopt the $\alpha$ prescription of viscosity (Shakura & Sunyaev, 1973), $t_{r\phi}=-\alpha p$, where $p=p_{\mathrm{gas}}+p_{\mathrm{rad}}$ is the total pressure. We adopt the Newtonian gravitational potential, $\Psi=-GM/r$, in this paper. Note that close to the central BH, relativistic effects become strong and self-similarity may no longer hold, so we keep our study beyond that region and focus on the region in the range of intermediate radii away from the inner and the outer boundaries, where self-similarity is a good approximation (Narayan & Yi, 1994; Narayan et al., 1997; Jiao et al., 2015). Thus Newtonian potential is good enough in our calculation.

For a typical super-Eddington accretion flow, both steady models and numerical simulations have shown that the principal part of the flow should be both optically thick and extremely radiation-pressure dominated (the ratio of gas pressure to total pressure $\beta\to 0$), ranging from the inner edge of the flow to several hundred or thousand times of Schwarzschild radius $r_{\mathrm{S}}$ $(\equiv 2GM/c^{2})$, depending on the accretion rate. The accretion flow beyond this range generally resembles the structure of a SSD, though it should still be optically thick and radiation-pressure dominated, up to a much larger radius (Kato et al., 2008, p.112). Here we focus the study on the typical region of a super-Eddington accretion flow. As the flow is optically thick, we can adopt the Eddington approximation, $p_{\mathrm{rad}}=E_{\mathrm{rad}}/3$, where $E_{\mathrm{rad}}$ is the radiation energy density (per unit volume). According to the first moment of the radiative transfer equation, for a steady radiation field, we have (Kato et al., 2008, p.490)

so in Equation (2) the radiation force term can actually be incorporated into the gas-pressure gradient term to form the gradient of the total pressure, $p=p_{\mathrm{gas}}+p_{\mathrm{rad}}$, which is the form of the equation of motion adopted in Xue & Wang (2005) and Jiao & Wu (2011).

The energy equation can be written as (Ohsuga et al., 2005; also see Jiao et al., 2015),

where $E=E_{\mathrm{gas}}+E_{\mathrm{rad}}$ is the total energy density of the coupled matter and radiation. As mentioned above, gas pressure can be neglected in a typical super-Eddington accretion flow, so here we actually use $p=p_{\mathrm{rad}}$ and $E=E_{\mathrm{rad}}=3p$, and the left hand side (LHS) of Equation (4) essentially represents the advection of radiation field. On the right hand side (RHS), $q_{\mathrm{vis}}$ is the viscous heating rate, which can be expressed as (Kato et al., 2008, p.239)

and $\nabla\boldsymbol{\cdot}\boldsymbol{F}$ represents the change of the energy density in a fixed region due to the transportation via the radiative flux $\boldsymbol{F}$, which is usually referred to as "radiative cooling" term in height-integrated accretion disk models, though it could actually be either positive or negative here.

To close this set of equations, we need to know the opacity $\kappa$. We consider the low metallicity case in this paper. The scattering opacity comes from electron scattering, so we have $\kappa_{\mathrm{sca}}=\kappa_{\mathrm{es}}\approx 0.4$ cm² g^-1. The absorption opacity is dominated by the free-free absorption, $\kappa_{\mathrm{ff}}\propto\rho T^{-7/2}$. In super-Eddington accretion flows, temperature is relatively high, and $\kappa_{\mathrm{es}}$ generally dominates over $\kappa_{\mathrm{ff}}$, so here we can safely neglect the absorption opacity and set $\kappa=0.4$ cm² g^-1. Note that for the high metallicity case, the bound-free absorption opacity can play an important role (Ohsuga et al., 2005), especially in the outflow region, which could be the cause of blue-shifted absorption lines in some broad absorption line quasars. This will be studied in our future work.

With some mathematical deduction, Equations (1), (2) and (4) can be transformed to 5 PDEs of $r$ and $\theta$ with 5 unknowns, $\rho$, $v_{r}$, $v_{\theta}$, $v_{\phi}$, and $p$. We adopt self-similar assumptions in the radial direction, so that the solution takes the form

where $v_{\mathrm{K}}=\sqrt{GM/r}$ is the Keplerian velocity. With these assumptions, Equations (1) and (2) are reduced to four ODEs independent of $r$ [see equations (16)–(19) in Jiao & Wu, 2011 for the detailed forms]. However, the $r$-dependency is not immediately eliminated in Equation (4). Here we write the equation as

where $q_{\mathrm{adv}}$ is the advective energy transport, defined as

and $q_{\mathrm{rad}}$ is the radiative energy transport, defined as

Under self-similar assumptions, we find that $q_{\mathrm{adv}}\propto r^{-5/2-n}$, $q_{\mathrm{vis}}\propto r^{-5/2-n}$, and $q_{\mathrm{rad}}\propto r^{-3}$. Thus in the energy equation, we cannot eliminate the $r$-dependency as in the continuity and momentum equations, unless we set $n=1/2$. Note that the mass inflow rate obeys $\dot{M}_{\mathrm{in}}\propto r^{3/2-n}$ according to Equation (6) & (7). Theoretical studies of Begelman (2012) proposed that hypercritical accretion should have $\dot{M}\propto r$, so that $n=0.5$ (but see Jiao, 2022 for a caveat in their work). As for numerical simulations, most works have shown that outflow becomes negligibly small close to the BH (generally for $r<10r_{\mathrm{S}}$ while the exact location is dependent on the specific simulation), which is not the region we study here. Outside this region, some simulations have presented the power-law fits of the density or mass inflow rate profiles. For example, Ohsuga et al. (2005) performed 2D radiation-hydrodynamic (RHD) simulations of super-Eddington accretion flows and found that the mass accretion rate is roughly in proportion to radius in the region which achieves a quasi-steady state, so that $n\sim 0.5$. Yang et al. (2014) expanded their work with an improved treatment of the viscous stress and found $n\approx 0.55$ (denoted $p$ in their paper). Kitaki et al. (2018) performed a series of 2D RHD simulations with various BH masses and accretion rates for super-Eddington accretion covering the space of 2–3000$r_{\mathrm{S}}$, and achieved a relatively large $r_{\mathrm{qss}}\sim 200r_{\mathrm{S}}$ inside which the simulation establishes a quasi-steady state, and they found $\rho\propto r^{-0.73}$ so that $n=0.73$. On the other hand, the three-dimensional general relativistic radiation-magnetohydrodynamic (3D GRRMHD) simulation performed by McKinney et al. (2014) found that $\dot{M}\propto r^{1.7}$ so that $n=-0.2$, though their simulation is for super-Eddington accretion flow onto a fast-spinning BH ($a/M=0.9375$) and the power-law fit includes the region that has not achieved inflow equilibrium, so their result may not be applicable to our study here. Overall, we can see that the value of $n$ is still an open question. Note that while simulations with more physically-realistic assumptions keep being developed, they are also more expensive in computation and generally achieve a quasi-steady state with inflow equilibrium at a smaller radius (e.g., the simulation in McKinney et al., 2014 has $r_{\mathrm{qss}}\sim 14r_{\mathrm{S}}$), so their power-law fits of radial density profiles are not necessarily better. In this paper, we will just set $n=1/2$, and leave the discussion of $n$ values to future works.

For $n=1/2$, the full set of equations becomes a set of ODEs independent of $r$, which can then be solved with proper boundary conditions and three input parameters [$\alpha$, $M$, $\rho(\pi/2)$]. Here $\alpha$ is the viscous parameter, $M$ is the mass of the central BH, and $\rho(\pi/2)$ is the value of $\rho(\theta)$ in cgs units on the equatorial plane in Equation (6), which is determined eventually by the accretion rate. So the input parameters are essentially identical to conventional height-integrated disk models, which are ($\alpha$, $M$, $\dot{M}$). Compared with the input parameters in our previous self-similar solution (Jiao & Wu, 2011), which are ($\alpha$, $f$, $\gamma_{\mathrm{equ}}$, $n$), the advective factor, $f=q_{\mathrm{adv}}/q_{\mathrm{vis}}$, is no longer a constant input parameter, but instead can be calculated from the equations; the equivalent heat capacity ratio, $\gamma_{\mathrm{equ}}$, is omitted, as gas pressure is negligible in super-Eddington accretion flows; the radial density index, $n$, is now set to be 1/2. On the other hand, the calculation of radiative flux involves the actual values of density and radius, unlike Jiao & Wu (2011) where only the relative value of density matters, so $M$ and $\rho(\pi/2)$ are now added as two new input parameters.

We assume that the accretion flow has reflection symmetry about the equatorial plane, so that

As mentioned in Section 2.1, we also set $\rho(\pi/2)$ as a boundary condition (given as an input parameter), which is eventually determined by the mass accretion rate of the flow. Note that in the energy Equation (4), the term $\nabla\boldsymbol{\cdot}\boldsymbol{F}$ contains 2nd-order derivatives of $p$, which must be dealt with. Under self-similar assumptions, the radiative flux $\boldsymbol{F}$ is expressed as

neglecting gas pressure. We can see that $r^{2}F_{r}$ actually becomes independent of $r$, so when calculating $\boldsymbol{\nabla\boldsymbol{\cdot}F}$, the contribution of $\boldsymbol{F}_{r}$ is 0. Physically this represents that the radiative flux flowing into and out of a fixed region in $r$ direction just cancels out each other. Thus we only need to consider $p^{\prime\prime}(\theta)$ as the 2nd-order term in the energy equation. Here we define $z_{p}(\theta)\equiv p^{\prime}(\theta)$, so that $p^{\prime\prime}(\theta)=z_{p}^{\prime}(\theta)$, transferring the 2nd-order differential equation into two 1st-order ones, so that the problem becomes a set of six ODEs with six unknowns.

To start the ODE solver from the equatorial boundary, we need to calculate the boundary values of the six unknowns, namely $\rho(\pi/2)$, $v_{r}(\pi/2)$, $v_{\theta}(\pi/2)$, $v_{\phi}(\pi/2)$, $p(\pi/2)$, and $z_{p}(\pi/2)$. The value of $\rho(\pi/2)$ is given as an input parameter, and Equation (14) provides the other five boundary conditions, although some of them are derivatives of the unknowns on the equatorial boundary. Ideally, we should be able to calculate the equatorial values of the six unknowns with the six ODEs at $\theta=\pi/2$. However, the momentum equation in $\theta$ direction, which is

under self-similar assumptions, becomes the form "0=0" on the equatorial plane. To solve this problem, we differentiate this equation with respect to $\theta$, which yields

As we have $v_{\theta}^{\prime\prime}(\pi/2)=0$ according to the reflection symmetry, Equation (2.2) can be used in the place of Equation (16) to solve for values of the six unknowns on the equatorial boundary. Note that Equation (2.2) is only used for calculating the boundary values of the unknowns, and we still use Equation (16) for the ODE solver, so the problem is still a set of six first-order ODEs and the equatorial boundary conditions are sufficient for solving the problem.

A typical solution of our model, with input parameters of $\alpha=0.1$, $M=10M_{\odot}$, and $\rho(\pi/2)=0.5$ g cm^-2.5 [for $n=1/2$, see Equation (6)], is presented in Figure 1. The corresponding density on the equatorial plane (denoted $\rho_{\mathrm{eq}}$), mass inflow rate (denoted $\dot{M}_{\mathrm{in}}$), and mass outflow rate (denoted $\dot{M}_{\mathrm{out}}$) are $\rho_{\mathrm{eq}}=9.2\times 10^{-6}$ g cm^-3, $\dot{M}_{\mathrm{in}}=6800L_{\mathrm{Edd}}/c^{2}$, and $\dot{M}_{\mathrm{out}}=853L_{\mathrm{Edd}}/c^{2}$ respectively at $r=1000r_{\mathrm{S}}$, which scale with radius in the form $\rho\propto r^{-1/2}$ and $\dot{M}\propto r$ (in a more general case, $\rho\propto r^{-n}$ and $\dot{M}\propto r^{3/2-n}$; see the discussion on mass conservation of steady solutions in Section 3.2). Here we set $1000r_{\mathrm{S}}$ as the outer boundary of the self-similar region, beyond which the mass inflow and outflow rates no longer change with radius, so $6800L_{\mathrm{Edd}}/c^{2}$ can also be regarded as the mass supply rate from the environment. The top three panels show the $\theta$ profiles of the three components of velocity, normalized by $v_{\mathrm{K}}=\sqrt{GM/r}$. The bottom three panels show the $\theta$ profiles of density, pressure, and radiative fluxes, from left to right, normalized by $\rho_{\mathrm{eq}}$, $\rho_{\mathrm{eq}}v_{\mathrm{K}}^{2}$, and $GMcr^{-2}\kappa^{-1}$, respectively. In the bottom right panel, the solid line represents $-F_{\theta}$ while the dashed line represents $F_{r}$. The values of $F_{\theta}$ are negative as they are in the $\theta$-decreasing direction towards the polar axis. Note that all these profiles are independent of $r$.

We can see that the accretion flow consists of an inflow region around the equatorial plane and an outflow region beyond the inflow, with $v_{r}=0$ at a certain inclination angle (for the solution presented in Figure 1, $\theta=1.03$ rad) between the inflow and the outflow. This can be more clearly seen in the velocity field diagram for the solution, Figure 2, where the inflow region is below the solid line (due to reflection symmetry, it is actally within 0.54 rad around the equatorial plane) and the outflow region is between the solid and the dashed lines. The dashed line represents the calculation boundary of our self-similar solution, which is further explained as follows.

The calculation starts from the equatorial plane, moves towards the polar axis, and stops at the inclination where density and pressure drops very low and numerical error begins to rise so steeply that even reducing the integration step to machine accuracy fails to find a next-step solution within the error control setting. We call this inclination angle the calculation boundary of our solution. Mathematically, the calculation boundary represents the place where the ODE solver cannot proceed due to a steep increase in numerical errors, and should be distinguished from the equatorial bounday which provides boundary conditions for our calculation. We have tried different ODE solvers (see, e.g., Hosea & Shampine 1996; Shampine & Reichelt 1997; Shampine et al. 1999; Shampine 2002) with different error control settings, and find that the calculation boundaries and the solutions are basically the same. For example, when we use the same ODE solver, changing the relative error tolerance from $10^{-12}$ to $10^{-3}$ only changes the calculation boundary for less than $10^{-6}$ rad. Thus we expect that this calculation boundary is intrinsic rather than dependent on the numerical method. The physical interpretation is that self-similarity cannot well describe the whole space (which it should not due to mass conservation, see Section 3.2 for more details), so that above the self-similar accretion flow that our solution describes, there exists non-self-similar flow in the region around the polar axis. When our self-similar solution approaches the boundary between the self-similar and the non-self-similar regions, numerical errors will inevitably undergo a steep rise, regardless of the ODE solver and the error control setting that we choose. So the calculation boundary of our solution is a good indication of the physical boundary between the self-similar and the non-self-similar regions in the accretion flow, which is also the surface of the self-similar accretion flow that we solves. Nevertheless, our solution solves the structure of the accretion flow in the self-similar region with a good numerical accuracy, and our analyses below based on this solution should be trustworthy within the framework of our assumptions. We have also calculated the optical depth $\tau$ from the surface, and find that the position of $\tau=1$ is very close to the surface ($\Delta\theta=0.01$ rad at $r=10r_{\mathrm{S}}$, and even smaller for larger radius). So the assumption of optical thickness is well held in most of the solution, even when the optical depth in the non-self-similar region is not counted.

It is also worth noticing in Figure 1 that the rotational speed $v_{\phi}$ is only $0.8964v_{\mathrm{K}}$ on the equatorial plane and becomes super-Keplerian at higher latitudes. The reason of $v_{\phi}<v_{\mathrm{K}}$ near the equatorial plane is that there is significant pressure gradient pointing outward in the $r$ direction, so that the required centrifugal force to resist the gravity is reduced. We also calculated the slim disk model solution for the mass supply rate $6800L_{\mathrm{Edd}}/c^{2}$, and found that the value of $v_{\phi}/v_{\mathrm{K}}$ on the equatorial plane inside $1000r_{\mathrm{S}}$ is in the range 0.6480–0.8995, so this is to be expected in conventional models for super-Eddington accretion. The equatorial value of $v_{\phi}$ in our solution also agrees with numerical simulations which have discussed this value (e.g., Ohsuga et al., 2005; Ohsuga & Mineshige, 2011). Note that while in figure 5 of Jiang et al. (2014) $v_{\phi}$ is very close to $v_{\mathrm{K}}$ outside of 10$r_{\mathrm{S}}$, their $v_{\phi}$ is a vertically averaged value, and cannot be compared directly with our equatorial value of $v_{\phi}$. The reason that $v_{\phi}$ could become super-Keplerian at higher latitudes is related to the conservation of angular momentum and the geometry of spherical coordinates (i.e., the distance to the axis of rotation decreases as $\theta$ decreases), which is further discussed in Section 3.3.

Compared with our previous self-similar solution (Jiao & Wu, 2011), a large improvement is that now the model includes the calculation of the radiative energy transport $q_{\mathrm{rad}}$, so that the advective factor $f$ is calculated self-consistently instead of being set as an input parameter. The energy mechanisms are presented in the left panel of Figure 3, where the solid line represents the heating rate of viscous dissipation $q_{\mathrm{vis}}$, the dashed line represents the radiative energy transport $q_{\mathrm{rad}}$, the dotted line represents the advective energy transport in $r$ direction $q_{\mathrm{adv}}^{r}$, and the dot-dashed line represents the advective energy transport in $\theta$ direction $q_{\mathrm{adv}}^{\theta}$, all of which are normalized by $(GM)^{3/2}r^{-n-5/2}$. The values of $q_{\mathrm{rad}}$, $q_{\mathrm{adv}}^{r}$, and $q_{\mathrm{adv}}^{\theta}$ can be either positive or negative, while positive values represent energy loss (cooling) and negative values represent energy gain (heating), and the sum of these values should equal $q_{\mathrm{vis}}$ in order to maintain energy conservation. We can see that $q_{\mathrm{rad}}$ is always a cooling mechanism in our solution, so that the term "radiative cooling" is actually accurate. It also gets stronger closer to the surface, where photons can escape more easily. Advection, on the other hand, is somewhat complex. While the radial advection in the slim disk model is an important or even dominant cooling mechanism for super-Eddington accretion, it is now a heating mechanism in the inflow of our solution. As shown in Equation (12), the advective energy transport can be divided to contributions of advected internal energy (here mostly radiative energy), $\boldsymbol{\nabla\boldsymbol{\cdot}}(E\boldsymbol{v})$, and compression heating/expansion cooling, $p\boldsymbol{\nabla\boldsymbol{\cdot}v}$. Under self-similar assumptions, their contributions in $r$ direction can be further expressed as (see Jiao, 2022, for more details)

where

and $q_{\mathrm{adv}}^{r1}$ and $q_{\mathrm{adv}}^{r2}$ denote the advected internal energy and compression heating/expansion cooling in $r$ direction, respectively. As we set $n=1/2$ in the solution, we actually have $q_{\mathrm{adv}}^{r1}=0$, which means that the incoming and outgoing internal energy in $r$ direction for a fixed region in the rest frame just cancel out each other and do not influence the energy balance. The radial advection is thus dominated by compression heating in the inflow and expansion cooling in the outflow, as shown by the dotted line in the left panel of Figure 3. The advection in $\theta$ direction $q_{\mathrm{adv}}^{\theta}$, on the other hand, acts as a cooling mechanism in the inflow and carries the left-over entropy that cannot be cooled by the local radiation upwards, which is eventually cooled by radiation or carried away by the outflow (i.e. cooled by expansion) at higher latitude. In the outflow, $q_{\mathrm{adv}}^{r}$ becomes a cooling mechanism due to expansion. The advection in $\theta$ direction $q_{\mathrm{adv}}^{\theta}$ switches to heating when the sum of radiative and expansion cooling becomes large enough to offset the viscous heating, which happens somewhere in the outflow. The radiative cooling $q_{\mathrm{rad}}$ gets much stronger closer to the surface, where $q_{\mathrm{adv}}^{\theta}$ becomes the main heating mechanism, which is essentially the accumulated entropy that cannot be sufficiently cooled in the inflow and the lower part of the outflow.

The $\theta$ profiles of the local advective factor $f$ (solid line) and the cumulative advective factor $f_{\mathrm{cum}}$ (dashed line) are presented in the right panel of Figure 3, where $f_{\mathrm{cum}}$ represents the ratio of advective energy transport to viscous heating integrated from the equatorial plane to the current $\theta$ and is calculated as

We present this value because as density and pressure drops close to 0 near the surface, the local viscous heating $q_{\mathrm{vis}}$, which is the denominator in the calculation of $f$, also approaches 0 and causes $f$ to become unusually large in absolute value, so that $f_{\mathrm{cum}}$ actually describes the energy mechanism better. We can see that as $\theta$ decreases, radiation gets stronger, and consequently $f_{\mathrm{cum}}$ keeps decreasing, indicating that radiation becomes more important in the total cooling. The mean advective factor of the accretion flow (denoted $\bar{f}$), which is $f_{\mathrm{cum}}$ calculated at the surface, is only 0.29, much lower than the local advective factor on the equatorial plane, 0.81, and also smaller than that expected by the slim disk model, which is generally advection-dominated ($>0.5$).

We have also made multiple calculations with different input parameters, and find that outflow always exists as long as the input parameters are reasonable values under our assumptions. Note that our model itself does not consider the effects of gas pressure and absorption opacity, so that we have to rely on other theoretical studies or numerical simulations to ensure that neglecting these effects are feasible. We have discussed in Section 2.1 that these assumptions are justified in the inner region of a super-Eddington accretion flow, ranging from the inner edge of the flow to several hundred or thousand times of Schwarzschild radius $r_{\mathrm{S}}$, depending on the accretion rate. So a simple way to ensure that our model is applicable is to ensure that the flow is super-Eddington. For observations, if the observed source is known to have an Eddington ratio larger than 1, then the application of our model is feasible for its inner region. For theoretical studies with a given set of input parameters, we can calculate the effective accretion rate of our model as

where $\theta_{\mathrm{s}}$ is the value of $\theta$ at the surface of our solution. As the black hole mass is known as an input parameter, we can calculate $\dot{M}_{\mathrm{Edd}}$ (here we adopt the definition of $\dot{M}_{\mathrm{Edd}}=16L_{\mathrm{Edd}}/c^{2}$, see, e.g., Abramowicz et al. 2010; we try to avoid a direct calculation of the radiative efficiency in our model, which would require additional assumptions, see Section 3.5) and compare these two values. The application of our model would be well justified in the inner region of the accretion flow when $\dot{M}_{\mathrm{eff}}>\dot{M}_{\mathrm{Edd}}$. Note that $\dot{M}_{\mathrm{eff}}\propto r$, and the eventual mass accretion rate that falls into the accretor can be estimated as $\dot{M}_{\mathrm{eff}}$ at the inner edge of the self-similar region, inside which outflow becomes negligible. While we do not know the exact place of the inner edge of the self-similar region, it is generally considered to be away from the inner boundary of the accretion flow (Narayan & Yi, 1994; Narayan et al., 1997; Jiao et al., 2015), so we can be certain that this inner edge is larger than the innermost circular orbit (ISCO) at $r_{\mathrm{ISCO}}=3r_{\mathrm{S}}$. Therefore, $\dot{M}_{\mathrm{eff}}$ calculated at $r_{\mathrm{ISCO}}$ is an lower limit (hereafter denoted $\dot{M}_{\mathrm{ll}}$), and if it is larger than $\dot{M}_{\mathrm{Edd}}$, the accretion flow can surely be considered super-Eddington (the real accretion rate is likely much larger than this lower-limit value).

The above criterion also incorporates the effects of changing input parameters. For example, with fixed values of $\alpha=0.1$ and $M=10M_{\odot}$, the solution for $\rho(\pi/2)=0.46$ g cm^-2.5 has a lower-limit of the effective accretion rate $\dot{M}_{\mathrm{ll}}=\dot{M}_{\mathrm{Edd}}$, and the corresponding boundary between inflow and outflow (hereafter denoted $\theta_{0}$) is $\theta_{0}=1.036$ rad, the surface of the solution is $\theta_{\mathrm{s}}=0.764$ rad, and the mean advective factor is $\bar{f}=0.266$. As a comparison, the solution presented in Figure 1 for $\rho(\pi/2)=0.5$ g cm^-2.5 has $\dot{M}_{\mathrm{ll}}=1.12\dot{M}_{\mathrm{Edd}}$, $\theta_{0}=1.03$ rad, $\theta_{\mathrm{s}}=0.75$ rad, and $\bar{f}=0.29$. We can see that for $\rho(\pi/2)=0.46$ g cm^-2.5, the sizes of the inflow and the outflow regions and the mean advective factor have all decreased, corresponding to a reduced accretion rate. The effective accretion rate is not linear to $\rho(\pi/2)$, because the sizes of the inflow and the outflow regions also decrease as $\rho(\pi/2)$ decreases.

The influence of changing black hole mass $M$ on the solution is not so intuitive as changing $\rho(\pi/2)$. When we compare the structure of different solutions, we usually study the region with the same values of $r$ in unit of $r_{\mathrm{S}}$ ($\equiv 2GM/c^{2}$). Note that the value of $r$ in Equation (24) is in cgs units, so that $r\propto r_{\mathrm{S}}\propto M$. Then, according to Equation (24), we can see that $\dot{M}_{\mathrm{eff}}\propto M^{3/2}$ when the differences in the profiles of $v_{r}(\theta)$ and $\rho(\theta)$ are ignored. On the other hand, the Eddington accretion rate $\dot{M}_{\mathrm{Edd}}\propto M$, so eventually the accretion rate in unit of $\dot{M}_{\mathrm{Edd}}$ is $\dot{m}=\dot{M}_{\mathrm{eff}}/\dot{M}_{\mathrm{Edd}}\propto M^{1/2}$ (when the change in vertical structure is considered, $\dot{m}$ should change faster than $M^{1/2}$, while a positive correlation still holds). Thus when $\alpha$ and $\rho(\pi/2)$ are fixed, decreasing the black hole mass $M$ also decreases the accretion rate (in unit of $\dot{M}_{\mathrm{Edd}}$). For example, with $\alpha=0.1$ and $\rho(\pi/2)=0.5$ g cm^-2.5, the solution for $M=5M_{\odot}$ has $\dot{M}_{\mathrm{ll}}=0.69\dot{M}_{\mathrm{Edd}}$, $\theta_{0}=1.05$ rad, $\theta_{\mathrm{s}}=0.84$ rad, and $\bar{f}=0.18$. We can see that the sizes of the inflow and the outflow regions and the mean advective factor have all decreased again, corresponding to a reduced accretion rate (in unit of $\dot{M}_{\mathrm{Edd}}$).

Zahra Zeraatgari et al. (2020) also established a self-similar solution for super-Eddington accretion. Compared with our solution, they considered all components of the viscous stress tensor, so that more boundary conditions are required and they applied the symmetric boundary conditions on both the equatorial plane and the polar axis. Consequently, the self-similar assumptions are adopted in the full space in their solution. Note that for a steady solution, the total accretion rate should be constant, which with reflection symmetry can be expressed as

If self-similarity holds for the whole range of $\theta$, we can further get,

As Zahra Zeraatgari et al. (2020) also set $n=1/2$, their solution must have $\dot{M}_{\mathrm{total}}=0$. However, the direct applications of super-Eddington accretion models in astronomy are to observations with super-Eddington luminosities, such as ULXs, whose energy eventually comes from the gravitational energy of accreted matter. The solution of Zahra Zeraatgari et al. (2020) requires that the mass outflow exactly cancels the mass inflow, so that no (or negligible) mass are actually accreted, which cannot explain the high luminosities in these observations.

On the other hand, our self-similar solution does not cover the whole range of $\theta$, and we expect that there exists a non-self-similar region around the polar axis, so that the total accretion rate in our solution is

where $\dot{M}_{\mathrm{ss}}$ is the net accretion rate in the self-similar region and $\dot{M}_{\mathrm{polar}}$ is the accretion rate of the non-self-similar region around the polar axis. As $\dot{M}_{\mathrm{polar}}$ does not follow the self-similar form, our solution no longer requires a total accretion rate of 0 as in Zahra Zeraatgari et al. (2020), so that a large total accretion rate can be allowed in our solution. The trade-off is that our solution can only describe a limited range of $\theta$. The direction of velocity on the surface of our solution in Figure 1 is still pointing upward and outward, so that there is outflow blowing from the surface into the non-self-similar region, and we speculate that the accretion flow in the non-self-similar region is also in the form of outflow, which helps to maintain a constant total accretion rate. This configuration is consistent with the structure obtained in numerical simulations of super-Eddington accretion (see references in Section 1), though the exact structure in the non-self-similar region is beyond our calculation.

Another difference is that we only consider the $r\phi$-component of the viscous stress tensor in the calculation, while Zahra Zeraatgari et al. (2020) considered all the components. However, the inclusion of other components of the viscous stress does not change the formula of advective energy transport, and only changes the formula of the viscous heating, which would still be a heating term anyway, so the general results in this paper should still hold.

In the classical picture of viscous accretion disk theory, the viscosity caused by differential rotation transfers angular momentum outward and disk material inward (Lynden-Bell & Pringle, 1974). When we focus on a fixed region in the rest frame, the viscous torque reduces the angular momentum in the region, but the material flowing into the region carries more angular momentum than that flowing out of the region, so that the angular momentum in the fixed region can be conserved for a steady solution. We can regard this gain (or loss) of angular momentum due to flow motion as the advection of angular momentum. In conventional height-integrated models, vertical motion is ignored, so that the radial advection of angular momentum exactly balances the viscous torque, and the angular momentum per unit mass relative to the axis of rotation would remain constant in the vertical direction. For a solution in cylindrical coordinates, this would result in a constant $v_{\phi}$ in $z$ direction. However, in spherical coordinates, the angular momentum per unit mass is $l=r\sin\theta v_{\phi}$, and a constant $l$ means that $v_{\phi}$ would increase as $\theta$ decreases in the form of $v_{\phi}\propto\csc\theta$, becasue the distance to the axis of rotation is proportional to $\sin\theta$. The real case is more complicated when vertical motion and advection are considered, especially in the outflow region where the radial advection also tends to decrease the angular momentum in the fixed region, and the conservation of angular momentum is upheld by the vertical advection.

To study the conservation of angular momentum relative to the polar axis (which is the axis of rotation here), we calculate $\boldsymbol{r}\times$ Equation (2), project the result onto the polar axis, and divide the equation by $\rho$ (to study the angular momentum per unit mass), which eventually yields

where $N_{\phi}$ is the $\phi$-component of the viscous force

The RHS of Equation (28) is the total torque per unit mass, contributed completely by the viscous force as gravity and pressure gradient exert no torque about the polar axis due to axisymmetry, while the two terms on the LHS are the radial and vertical advection of angular momentum, respectively. The two components of the advection of angular momentum and the viscous torque per unit mass, calculated from our solution, are displayed in the left panel of Figure 4. The solid line represents the total torque, the dotted line represents the total advection of angular momentum, and the dashed and dot-dashed lines represent the advection of angular momentum in $r$ and $\theta$ directions, respectively. The advection of angular momentum is calculated as the difference between the outgoing and the incoming angular momentum, so that positive values here represent loss, and negative values represent gain, of angular momentum. The total torque is always negative, indicating that viscosity tends to reduce the angular momentum in the fixed region. The total advection of angular momentum coincides with the total torque, indicating that the angular momentum conservation is maintained, as expected for a steady solution. The radial advection starts out as angular momentum gain on the equatorial plane, decreases as $\theta$ decreases, and becomes angular momentum loss in the outflow region, where the angular momentum conservation is upheld solely by the vertical advection. The vertical advection always increases the angular momentum in the fixed region, which indicates that the angular momentum per unit mass tends to decrease in the direction of $v_{\theta}$. This can be clearly seen in the right panel of Figure 4, where the solid line represents angular momentum per unit mass normalized by $r$, $l/r=\sin\theta v_{\phi}$, and the dotted line represents $v_{\phi}$, shown here as a comparison. We can see that although the angular momentum per unit mass decreases as $\theta$ decreases for a fixed radius, $v_{\phi}$ keeps increasing as the distance to the axis of rotation decreases faster than the angular momentum.

We expand the equation of motion, Equation (2), in the $r$ and $\theta$ directions, and get

where the last terms on the LHS are the $r$ and $\theta$ components of the centrifugal force per unit mass, $f_{\mathrm{cen}}=v_{\phi}^{2}/(r\sin\theta)$, whose direction is perpendicular to the axis of rotation. Note that in spherical coordinates, gravity has no component in the vertical ($\theta$) direction while centrifugal force has a component, $\cos\theta\cdot f_{\mathrm{cen}}=v_{\phi}^{2}\cot\theta/r$, and the motion in the vertical ($\theta$) direction is thus determined by the pressure gradient and the centrifugal force. This is quite different from the analyses of accretion flow solutions in cylindrical coordinates, where the motion in the vertical ($z$) direction is determined by the $z$-components of gravity and pressure gradient while centrifugal force has no influence. As we ignored gas pressure due to the domination of radiation pressure in our solution, the pressure gradient is actually the radiative force, whose components in $r$ and $\theta$ directions agree with the profiles of $F_{r}$ and $F_{\theta}$ (the bottom-right panel of Figure 1).

Components of forces (per unit mass) in $r$ and $\theta$ directions are shown in Figure 5. In each panel, the solid line represents the total force, the dotted line represents the radiative force (i.e., pressure gradient), the dashed line represents the centrifugal force, and the dot-dashed line represents the gravitational force. The directions are in typical definition of spherical coordinates (except for the gravitational force which we present the absolute value here), so that positive values in $r$ direction points outward away from the accretor, and positive values in $\theta$ direction points downward away from the polar axis.

In $r$ direction, the motion is determined by the combined effects of radiative, centrifugal, and gravitational forces. The total force in $r$ direction is marginally negative (i.e., pointing inward) near the equatorial plane, but it needs to provide the inward acceleration of $v_{r}$ in $r$ direction, which means that the first term in Equation (30) must be negative. The total effect for the acceleration of $v_{r}$ in $\theta$ direction is thus always equivalent to an outward-pointing force, i.e., $\partial v_{r}/\partial\theta$ is always negative (except for the equatorial plane where $\partial v_{r}/\partial\theta=0$ due to reflection symmetry). Therefore, as $\theta$ decreases, $v_{r}$ tends to increase outward, and the inward radial motion on the equatorial plane will eventually turn outward and outflow is developed, as the combination of radiative and centrifugal forces overcomes the gravity. We can also see that the $r$-component of the radiative force maintains a relatively high level in the majority of the flow except for the region very close to the surface, and plays an important role in resisting the gravity and driving the outflow, although the dominant outward-pointing force is still the centrifugal force. This agrees with numerical simulations of super-Eddington accretion that have presented the analyses of forces in $r$ direction (e.g., fig. 5 of Ohsuga & Mineshige 2007; figure 5 of Yang et al. 2014). Some simulations have stated that close to the polar axis, the radiative force becomes dominant in $r$ direction (e.g., Ohsuga et al., 2005, 2009; Ohsuga & Mineshige, 2007, 2011; Yang et al., 2014; Sądowski et al., 2014), but our solution can not reach there as the region around the polar axis is non-self-similar.

In $\theta$ direction, the radiative force always points away from the equatorial plane, as the radiation tries to escape the accretion flow and the $\theta$-component of radiative flux, $F_{\theta}$, always points away from the equatorial plane and keeps increasing as $\theta$ decreases. The $\theta$-component of centrifugal force, on the other hand, pushes the material in the opposite direction. The total force is always negative in $\theta$ direction (except for on the equatorial plane where the $\theta$-components of both the total force and the velocity are 0), so that $v_{\theta}$ starts as 0 on the equatorial plane and becomes negative (pointing away from the equatorial plane) while its absolute value keeps increasing as $\theta$ decreases, as shown in Figure 1. The force analysis in $\theta$ direction agrees with the numerical simulation of Yang et al. (2014) which has presented this analysis in their figure 6.

Note that the definition of outflow and consequently the discussion of outflow-driving mechanism may differ in different papers. For example, in Jiang et al. (2014) whose calculation is performed in cylindrical coordinates, the material that escapes vertically is also considered outflow, and the outflow-driving mechanism is discussed based on vertical ($z$) components of forces (their figure 8). If we make a similar discussion in $\theta$ direction, we can also say that our flow is radiation-driven because the radiative force in $\theta$ direction is the dominant force that drives the flow away from the equatorial plane.

The relatively small value of mean advective factor indicates that compared with the conventional slim disk model, a larger portion of viscous heating is converted into radiation. As mentioned in Section 3.1, the advection in $\theta$ direction carries the left-over entropy at lower latitudes that cannot be radiated away upward, part of which is converted into radiation at higher latitudes. This effect is not considered in the slim disk model which neglects vertical motion. However, our solution cannot resolve the structure in the non-self-similar region, and the net accretion rate in the self-similar region of our solution decreases as radius decreases, which reduces the total viscous heating. So the consideration of vertical motion and outflow in our model actually influences the radiative efficiency $\eta$ in two opposite ways, when it is calculated according to the mass supply rate. Thus, it is necessary to compare the emergent radiative fluxes and luminosities of our model and the slim disk model.

The radiative flux in our model has two components, $F_{r}$ and $F_{\theta}$, under axisymmetry. As shown in the bottom right panel of Figure 1, $F_{\theta}$ starts out as 0 on the equatorial plane and increases towards the polar axis as the local heating is converted into radiation, reaching its maximum value at the surface of the accretion flow (negative values represent the direction). On the other hand, $F_{r}$ does not contribute to the "radiative cooling" as mentioned in Section 2.2, and can be perceived as being generated in the inner part of the accretion flow close to the central BH, which then propagates outward in $r$ direction following the inverse-square law. We can actually calculate the contribution of $F_{r}$ to the total luminosity as

where $\theta_{\mathrm{s}}$ is the value of $\theta$ at the surface. For the solution presented in Figure 1, we get $L_{r}=0.124L_{\mathrm{Edd}}$. In height-integrated models in cylindrical coordinates, such as the SSD and the slim disk, people usually calculate the radiative cooling per unit area on the equatorial plane (denoted $Q_{\mathrm{rad,cyl}}$ with a subscript "cyl" to distinguish from our model in spherical coordinates) as

where $H$ is the half thickness of the disk in cylindrical coordinates, and $q_{\mathrm{rad,cyl}}$ represents the radiative cooling per unit volume in these models. In this way, the emergent radiative flux $F=Q_{\mathrm{rad,cyl}}/2$. The total luminosity is calculated as

Apparently the emergent radiative flux as well as the luminosity in height-integrated models does not consider the contribution of the radiative flux in $r$ direction, so it corresponds to the emergent $F_{\theta}$ at the surface (denoted $F_{\theta\mathrm{s}}$) in our model. We can also calculate the radiative cooling per unit area of our model in a similar way,

However, the relation between $F_{\theta\mathrm{s}}$ and $Q_{\mathrm{rad}}$ is actually $F_{\theta\mathrm{s}}=Q_{\mathrm{rad}}/(2\sin{\theta_{\mathrm{s}}})$. The reason is that the emitting surface is smaller by a factor of $\sin{\theta_{\mathrm{s}}}$ compared with the height-integrated models in cylindrical coordinates. The luminosity emergent in $\theta$ direction is thus

which is basically the same formula for $Q_{\mathrm{rad}}$ as Equation (34). Note that self-similar assumptions are only applicable to a certain range of radii in the accretion flow, and does not hold in, e.g., the region very close to the BH where general relativistic effects are strong, while the contribution to luminosity there cannot be neglected. So the luminosity emergent in $\theta$ direction cannot be obtained in this paper without further assumptions. However, it can still be compared indirectly by comparing the radiative fluxes in the self-similar region.

The comparison of emergent radiative fluxes are shown in Figure 6. The solid line represents $F_{\theta\mathrm{s}}$, the dotted line represents $F_{\theta\mathrm{s}}\sin{\theta_{\mathrm{s}}}$, with a factor $\sin{\theta_{\mathrm{s}}}$ to keep the luminosity formula the same as models in cylindrical coordinates [Equations (34) and (36)], and the dashed line represents the emergent radiative flux of the slim disk for the mass supply rate from the environment, $\dot{M}_{\mathrm{supply}}=6800L_{\mathrm{Edd}}/c^{2}$, calculated with the code for traditional slim disk model in Jiao et al. (2009) (not the revised slim disk model whose accretion rate has an upper limit beyond which outflow will be inevitably developed), where we set $\kappa_{\mathrm{es}}=0.4$ cm² g^-1 in accordance with the low metallicity case in this paper, instead of $\kappa_{\mathrm{es}}=0.34$ cm² g^-1 in Jiao et al. (2009) corresponding to solar elemental abundances. Note that although we plot the lines all the way to the ISCO at $r=3r_{\mathrm{S}}$, the self-similar assumptions are generally no longer reliable close to the BH. We can see that both $F_{\theta\mathrm{s}}$ and $F_{\theta\mathrm{s}}\sin{\theta_{\mathrm{s}}}$ are significantly larger than $F_{\mathrm{slim}}$, indicating that the radiative efficiency of our model is larger than that of the slim disk model, even when calculated based on the mass supply rate. This result agrees with the numerical simulations of Jiang et al. (2014, 2019).

As mentioned above, $F_{\theta}$ keeps increasing as $\theta$ decreases because the local radiative cooling adds into the flux, which can be clearly seen in Equation (35). As the non-self-similar region beyond our calculation boundary is an outflow region, the photons will not be carried into the central BH and we expect that the total radiation will continue to increase in this region, so that the radiative flux obtained in our solution is only a lower limit. This further enhances our conclusion that the radiative efficiency of super-Eddington accretion should be larger than that obtained from the slim disk model.