Vertical Structure of Radiation-pressure-dominated Thin Disks: Link between Vertical Advection and Convective Stability

The standard thin accretion disk model under the alpha description was constructed by Shakura and Sunyaev (Shakura & Sunyaev 1973), and has been successfully applied to X-ray binaries and active galactic nuclei (for a review, see Frank et al. 2002; Kato et al. 2008). However, when we compare the theoretical results with observations and simulations, some basic problems become apparent, such as thermal stability and convective stability. Recent simulations (Hirose et al. 2009; Jiang et al. 2013) on thin disks based on a shearing box showed that, for high-mass accretion rates around $0.1{\dot{M}}_{\mathrm{Edd}}$, where ${\dot{M}}_{\mathrm{Edd}}$ is the Eddington accretion rate, the inner disk is radiation-pressure-dominated, and the vertical advection can be an efficient process for energy transport (this is probably related to the magnetic buoyancy; Hirose et al. 2009; Jiang et al. 2013). In addition, the local convective stability criterion, ${dS}/{dz}\gt 0$, where S is the entropy, is well satisfied according to the simulation results. On the other hand, according to the classic theory of thin disks, the vertical energy transport is completely dominated by the diffusion process. Moreover, when the radiation pressure dominates the gas pressure, the local convective stability criterion ${dS}/{dz}\gt 0$ is not satisfied, and therefore the disk may suffer convective instability (Sa̧dowski et al. 2011). The question arises: is there any link between the vertical advection and the convective stability?

In the present work, we will revisit the vertical structure and energy transport of radiation-pressure-dominated thin disks by including the possible vertical advection process. In addition, based on the simulation results, we will take into account the local convective stability criterion in order to modify the set of equations. The remainder of the paper is organized as follows. Equations and boundary conditions are described in Section 2. Numerical results are shown in Section 3, and conclusions and discussion are presented in Section 4.

The set of equations for the vertical structure of thin disks is based on the alpha-stress assumption. The gas pressure ${P}_{\mathrm{gas}}$ and the radiation pressure ${P}_{\mathrm{rad}}$, and their derivatives, take the following expressions:

$$\begin{eqnarray}&&{P}_{\mathrm{gas}}=\displaystyle \frac{\rho {k}_{{\rm{B}}}T}{\mu {m}_{{\rm{p}}}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{rad}}=\displaystyle \frac{1}{3}{{aT}}^{4},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dP}}_{\mathrm{rad}}}{{dz}}=-\displaystyle \frac{{\kappa }_{\mathrm{es}}\rho {F}_{z}}{c},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dP}}_{\mathrm{gas}}}{{dz}}=-\rho {{\rm{\Omega }}}_{{\rm{K}}}^{2}z+\displaystyle \frac{{\kappa }_{\mathrm{es}}\rho {F}_{z}}{c},\end{eqnarray} \tag{ 4 }$$

where ρ is the mass density, T is the temperature, ${k}_{{\rm{B}}}$ is the Boltzmann constant, and ${m}_{{\rm{p}}}$ is the mass of a proton. The Keplerian angular velocity ${{\rm{\Omega }}}_{{\rm{K}}}$ is written as ${{\rm{\Omega }}}_{{\rm{K}}}\,=\sqrt{{{GM}}_{\mathrm{BH}}/R{(R-{R}_{{\rm{g}}})}^{2}}$ under the well-known Paczyński–Wiita potential ${{\rm{\Psi }}}_{\mathrm{PW}}=-{{GM}}_{\mathrm{BH}}/(R-{R}_{{\rm{g}}})$ (Paczyński & Wiita 1980), where the gravitational radius is defined as ${R}_{{\rm{g}}}\equiv 2{{GM}}_{\mathrm{BH}}/{c}^{2}$.

In the standard disk model, the angular velocity is assumed to be Keplerian, i.e., ${\rm{\Omega }}={{\rm{\Omega }}}_{{\rm{K}}}$. Moreover, the local energy balance at each radius is such that the viscous heating rate equals the radiative cooling rate, and the vertical energy transport is dominated by the photon diffusion process. Following on from the simulation results, the vertical advection may have an additional contribution to the energy transport. Thus, the vertical flux F_z due to the diffusion is expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{{dF}}_{z}}{{dz}}=-{\tau }_{r\phi }{{\rm{\Omega }}}_{{\rm{K}}}{g}_{* }\,-\,\displaystyle \frac{{{dF}}_{\mathrm{adv}}}{{dz}},\end{eqnarray} \tag{ 5 }$$

where ${F}_{\mathrm{adv}}$ is the vertical flux owing to the vertical advection process, and the parameter g_* takes the form ${g}_{* }\,\equiv -d\mathrm{ln}{{\rm{\Omega }}}_{{\rm{K}}}/d\mathrm{ln}R=3/2+1/(R/{R}_{{\rm{g}}}-1)$. The $r\phi$-component of the shear stress ${\tau }_{r\phi }$ is assumed to be proportional to the total pressure, i.e., ${\tau }_{r\phi }=-\alpha ({P}_{\mathrm{gas}}+{P}_{\mathrm{rad}})$, where α is a constant parameter.

The total entropy S is the sum of gas and radiation, which is written as

$$\begin{eqnarray*}&&S={S}_{\mathrm{gas}}+{S}_{\mathrm{rad}}=\displaystyle \frac{3}{2}\displaystyle \frac{{k}_{{\rm{B}}}}{\mu {m}_{{\rm{p}}}}\mathrm{ln}\left(\displaystyle \frac{{P}_{\mathrm{gas}}}{{\rho }^{5/3}}\right)+\displaystyle \frac{4}{3}\displaystyle \frac{{{aT}}^{3}}{\rho }+\mathrm{const}.\end{eqnarray*}$$

It is well known that the local convective stability criterion can be expressed as

$$\begin{eqnarray*}&&\displaystyle \frac{{dS}}{{dz}}\gt 0\qquad (z\gt 0),\end{eqnarray*}$$

which is equivalent to the relationship between the radiative gradient ${{\rm{\nabla }}}_{\mathrm{rad}}$ and the adiabatic gradient ${{\rm{\nabla }}}_{\mathrm{ad}}$:

$$\begin{eqnarray*}&&{{\rm{\nabla }}}_{\mathrm{rad}}\lt {{\rm{\nabla }}}_{\mathrm{ad}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{{\rm{\nabla }}}_{\mathrm{rad}}\equiv {\left(\displaystyle \frac{\partial \mathrm{ln}T}{\partial \mathrm{ln}P}\right)}_{\mathrm{rad}},\qquad {{\rm{\nabla }}}_{\mathrm{ad}}\equiv {\left(\displaystyle \frac{\partial \mathrm{ln}T}{\partial \mathrm{ln}P}\right)}_{S}.\end{eqnarray*}$$

With the above expression of the total entropy S, we can derive the explicit form of ${{\rm{\nabla }}}_{\mathrm{ad}}$ as

$$\begin{eqnarray*}&&{{\rm{\nabla }}}_{\mathrm{ad}}=\displaystyle \frac{4-3\beta }{1.5{\beta }^{2}+12\beta {(1-\beta )+(4-3\beta )}^{2}},\end{eqnarray*}$$

where β is defined as $\beta \equiv {P}_{\mathrm{gas}}/({P}_{\mathrm{gas}}+{P}_{\mathrm{rad}})$.

As mentioned in the first section, simulations have not shown convective instability, which indicates that the relationship ${{\rm{\nabla }}}_{\mathrm{rad}}\leqslant {{\rm{\nabla }}}_{\mathrm{ad}}$ may be well satisfied in the disk. Following on from this, we assume the thermodynamical gradient to be

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}T}{d\mathrm{ln}P}=\min ({{\rm{\nabla }}}_{\mathrm{rad}},\lambda {{\rm{\nabla }}}_{\mathrm{ad}}),\end{eqnarray} \tag{ 6 }$$

where $\lambda \leqslant 1$ can guarantee the disk to be convectively stable.

The system consists of six equations, Equations (1)–(6), for the six unknown variables ρ, T, ${P}_{\mathrm{gas}}$, ${P}_{\mathrm{rad}}$, F_z, and ${F}_{\mathrm{adv}}$. There are four first-order differential equations in this system. In addition, the position of photosphere ${H}_{\mathrm{phot}}$ is unknown. Thus, five boundary conditions are required to solve the system between the equatorial plane (z = 0) and the photosphere ($z={H}_{\mathrm{phot}}$).

On the equatorial plane there exist two natural boundary conditions:

$$\begin{eqnarray}&&{F}_{z}=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{adv}}=0.\end{eqnarray} \tag{ 8 }$$

At the photosphere, the other three boundary conditions can be derived:

$$\begin{eqnarray}&&-{\kappa }_{\mathrm{es}}{\rho }^{2}{\left(\displaystyle \frac{d\rho }{{dz}}\right)}^{-1}=1,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{F}_{z}=2\sigma {T}^{4},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{F}_{z}+{F}_{\mathrm{adv}}=\displaystyle \frac{1}{4\pi }\dot{M}{{\rm{\Omega }}}_{{\rm{K}}}^{2}{f}_{* }{g}_{* },\end{eqnarray} \tag{ 11 }$$

where Equation (9) can be regarded as the definition of the photosphere position. The five boundary conditions ((7)–(11)) together with Equations (1)–(6) enable us to solve the system and derive the vertical structure.

Following Hirose et al. (2009), we define a vertical thickness H as

$$\begin{eqnarray*}&&H=\displaystyle \frac{3}{8\pi }\displaystyle \frac{{\kappa }_{\mathrm{es}}\dot{M}}{c},\end{eqnarray*}$$

which is used as a scale of the height. For simplicity, we choose $\lambda =1$ in Equation (6) for our numerical calculations. The other parameters are ${M}_{\mathrm{BH}}=10{M}_{\odot }$, $\dot{M}=0.1{\dot{M}}_{\mathrm{Edd}}$, and $\alpha =0.01$.

In order to directly compare the structure, including the vertical advection process, with that in the classic theory, we made numerical calculations for both of these two models. Figure 1 shows the vertical density profile, where the solid line corresponds to the results including vertical advection, and the dashed line corresponds to the results of the classic model. The solid line has a peak at z = 0, whereas the dashed line has two peaks at $z/H\approx \pm 0.3$. The peculiar shape of the dashed line, i.e., increasing ρ with z in the range $0\lt z/H\lt 0.3$, indicates that the disk suffers convective instability. The well-known local convective stability criterion, ${dS}/{dz}\gt 0$ for $z\gt 0$ (or ${dS}/{dz}\lt 0$ for $z\lt 0$), may work well in geometrically thin disks, where the radial velocity is low and therefore the advection effects may be negligible. Even though the entropy profile is not plotted in Figure 1, it is obvious that the entropy S will decrease with increasing z in the range $0\lt z/H\lt 0.3$ due to the peculiar density profile. Thus, the local convective stability criterion is not satisfied and therefore the disk is likely to be convectively unstable, as previously investigated by Sa̧dowski et al. (2011). On the contrary, the profile of the density in the case including vertical advection (solid line) is similar to that of the simulations (Hirose et al. 2009), which shows a continuously decreasing density with increasing z for $z\gt 0$. Thus, a convectively stable disk is quite possible.

Zoom In Zoom Out Reset image size

Figure 2 shows the vertical profiles of flux; the black dashed line corresponds to the diffusive flux F_z in the classic thin disk model, and the two solid lines correspond to the case including the vertical advection, where the black and red lines show the variations of the diffusive flux F_z and the advection flux ${F}_{\mathrm{adv}}$, respectively. The vertical advection makes a significant contribution to the total radiation flux, in particular for the region near the equatorial plane ($0\lt z/H\lt 0.1$). In some regions on the right part ($z\gt 0$), the red line is even higher than the black solid line, which means that ${F}_{\mathrm{adv}}$ can dominate over F_z. In addition, the profile of ${F}_{\mathrm{adv}}$ in Figure 2 is quite similar to that in the simulations (e.g., Hirose et al. 2009). Since our numerical calculation is based on a convectively stable disk, the results may indicate a link between the energy transport due to vertical advection and the convectively stable property. This link is probably due to the additional energy transport that decreases the entropy in the region near the equatorial plane, therefore allowing the entropy to increase with increasing vertical height, thus satisfying the local convective stability criterion.

Zoom In Zoom Out Reset image size

Taking $\dot{M}=0.1{\dot{M}}_{\mathrm{Edd}}$ as a typical accretion rate, we also investigate the strength of vertical advection for different radii. Obviously, the ratio of the gas pressure to the total pressure β will increase with an increasing radius for a fixed $\dot{M}$. In the outer part, where gas pressure is dominant, the vertical advection may be negligible. In other words, the diffusion process is the dominant energy transport and the disk is well convectively stable. However, in the inner part, where the radiation pressure dominates the gas pressure, the diffusion and the vertical advection may both be of importance for energy transport.

Figure 3 shows the profiles of diffusive flux (black) and advective flux (red) at the three radii $R=10{R}_{{\rm{g}}}$ (solid lines), $100{R}_{{\rm{g}}}$ (dashed lines), and $1000{R}_{{\rm{g}}}$ (dotted line). Here, we take the height of photosphere ${H}_{\mathrm{phot}}$ as the length unit. The red solid line shows that at $R=10{R}_{{\rm{g}}}$ the advective flux ${F}_{\mathrm{adv}}$ covers the range from the equatorial plane (z = 0) to the photosphere ($z={H}_{\mathrm{phot}}$). At the larger radius $R=100{R}_{{\rm{g}}}$, the red dashed line shows a smaller vertical range for ${F}_{\mathrm{adv}}$ ($\lesssim 2{H}_{\mathrm{phot}}/3$). Moreover, for the sufficiently large radius $R=1000{R}_{{\rm{g}}}$, Figure 3 shows that the advective flux disappears and diffusion is the only mechanism for the vertical energy transport, as shown by the dotted line. Thus, for thin disks, the two conditions for the occurrence of vertical advection are high accretion rates $\dot{M}\gtrsim 0.1{M}_{\odot }\,{{\rm{s}}}^{-1}$ and small radii $R\lesssim 100{R}_{{\rm{g}}}$.

Zoom In Zoom Out Reset image size

In this work, we have revisited the vertical structure of radiation-pressure-dominated thin disks by taking into account the role of the vertical advection process. Our study has shown that vertical advection not only significantly contributes to the energy transport, but it also plays an important role in making the disk convectively stable. Thus, a link may exist between the vertical advective energy transport and the convectively stable property. This link is probably due to the additional energy transport that decreases the entropy in the region near the equatorial plane, therefore allowing the entropy to increase with increasing vertical height, thus satisfying the local convective stability criterion. Our work may be helpful for understanding the discrepancy between classic theory and simulation results.

It is important to note that a detailed study of convective stability may require global, rather than local, stability analyses. For example, Abramowicz et al. (1993) demonstrated that when the viscosity is fully taken into account, stability analyses cannot be discussed within the framework of a local analysis, and a fully global treatment is required. Moreover, global stability analyses of the vertical convection of a thin gaseous disk were performed by several groups (e.g., Ruden et al. 1988). On the other hand, it is known that the advection may play an important role in stabilizing the disk against dynamic, thermal, and viscous perturbations. For instance, a geometrically thick disk without radial motion (Abramowicz et al. 1980; Paczyński & Wiita 1980) may suffer the Papaloizou-Pringle instability (Papaloizou & Pringle 1984), which is a dynamic instability based on the acoustic perturbations propagating between two boundaries in a differential rotating system. By including the advection terms, Blaes (1987) found that all of the unstable modes for the purely rotating flow are quickly stabilized by the advection process. Furthermore, Abramowicz et al. (1988) proposed the well-known slim disk model (or the optically thick, advection-dominated accretion disk) for super-Eddington accretion systems. The slim disk was found to be dynamically stable, thermally stable, and viscously stable, which may be related to the radial advection process. In the present work, we focus on the stability of geometrically thin disks, where the radial velocity (${v}_{R}\approx -\alpha {(H/R)}^{2}{v}_{{\rm{K}}}$, where ${v}_{{\rm{K}}}$ is the Keplerian velocity) is quite low, and therefore the effects of advection may also be quite weak.

In a previous work, Gu (2012) showed that, for high-mass accretion rates $\dot{M}\gtrsim {\dot{M}}_{\mathrm{Edd}}$ where the radiation pressure completely dominates the gas pressure, the disk is likely to be convectively stable without including the energy transport through the vertical advection. The physics of the stable property is probably related to the radial advection effects and the geometrically thick structure. On the other hand, in recent years many global simulation studies have been conducted on the super-Eddington accretion flows (e.g., Ohsuga et al. 2005; Jiang et al. 2014; Sa̧dowski & Narayan 2015). The simulations of Jiang et al. (2014) revealed the importance of the vertical advection, which can essentially enhance the radiative efficiency. In addition, outflows may play another important role in such flows (Jiang et al. 2014; Sa̧dowski & Narayan 2015). In our opinion, the theory of super-Eddington accretion flows is worth further investigation.

The authors thank Yan-Fei Jiang for providing the vertical profile of entropy in the simulations, and thank the referee for helpful comments that improved the paper. This work was supported by the National Basic Research Program of China (973 Program) under grants 2014CB845800, the National Natural Science Foundation of China under grants 11573023, 11333004, and 11222328, and the CAS Open Research Program of Key Laboratory for the Structure and Evolution of Celestial Objects under grant OP201503.