Accretion flows with comparable radiation and gas pressures

TL;DR

This paper models the vertical structure of accretion flows where radiation and gas pressures are comparable, revealing high energy advection and increased disc thickness with higher gas pressure, using self-similar solutions.

Contribution

It introduces a detailed model considering photon absorption and separate energy equations for matter and radiation in accretion flows with comparable pressures.

Findings

High energy advection occurs even below Eddington accretion rates.

Greater accretion rates and lower gas pressure portions lead to more energy being advected inward.

The accretion disc remains convectively stable across its structure.

Abstract

By taking into account photon absorption, we investigate the vertical structure of accretion flows with comparable radiation and gas pressures. We consider two separate energy equations for matter and radiation in the diffusion limit. In order to solve the set of radiation hydrodynamic equations in steady state and axisymmetric configuration, we employ self-similar technique in the radial direction. We need the reflection symmetry about the mid-plane to find gas density at the equator. For a typical solution, we assume that the gas pressure has 10-50\% portion of the total pressure. In this paper, since the radiation energy is involved directly, we are able to estimate how much energy of viscous heating is transported in the radial direction and advected towards the central object. Our results show that although the mass accretion rate does not approach the Eddington limit, the energy advection is rather high. Moreover, in a disc with greater accretion rate and less portion of gas pressure at the total pressure, more energy is advected to its center. In addition, as we expect the accretion flow becomes thicker with greater values of gas pressure. Based on Solberg-H\~A\v{z}iland conditions, we notice that the flow is convectively stable in all parts of such a disc.

Figures

Tables

Table 1. Table 1: Equatorial values of quantities

${\beta }_0$	0.1	0.3	0.5
$T_0/{10}^9\left(K\right)$	0.5	2	5
${\rho }_0/{10}^{-6}(g.c{m}^{-3})$	2.64	1.65	1.37
$v_{\varphi 0}/v_K$	0.93	0.90	0.85
$-v_{r0}/v_K$	0.02	0.03	0.05
$p_{g0}/{10}^{11}(dyn.c{m}^{-2})$	2.18	5.46	11.29
$p_{r0}/{10}^{12}(dyn.c{m}^{-2})$	1.96	1.27	1.13

Table 2. Table 2: Results

${\beta }_0$	0.1	0.3	0.5
$\dot{m}$	0.20	0.28	0.44
$\Delta {\theta }^{\circ }$	13.17	25.46	39.72
H/r	0.23	0.43	0.64
$f_{adv}$	0.15	0.14	0.13
$Q_{vis}/{10}^{21}\left(ergc{m}^{-3}{s}^{-1}\right)$	0.86	1.17	1.80
$Q_{adv}/{10}^{20}\left(ergc{m}^{-3}{s}^{-1}\right)$	1.33	1.72	2.37
$Q_c/{10}^{20}\left(ergc{m}^{-3}{s}^{-1}\right)$	3.14	5.25	10.5

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \text{v})=0,$$

$$\rho \frac{D\text{v}}{Dt}=-\nabla p_g-\rho \nabla \Psi -\nabla \cdot {\text{T}}^{\nu }+\frac{\chi }{c}{\text{F}}_0$$

$$\chi =\frac{{\sigma }_T}{m_p}\rho +\kappa$$

$${\kappa }_{ff}=1.7\times 10^{-25}{T}^{-7/2}(\frac{\rho }{m_p}{)}^{2}c{m}^{-1},$$

$${\kappa }_{bf}=4.8\times 10^{-24}{T}^{-7/2}\frac{Z}{Z_{\odot }}(\frac{\rho }{m_p}{)}^{2}c{m}^{-1},$$

$$\frac{\partial e}{\partial t}+\nabla \cdot (e\text{v})=-p_g\nabla \cdot \text{v}+{\Phi }_{vis}-4\pi \kappa B+c\kappa E_0,$$

$$\frac{\partial E_0}{\partial t}+\nabla \cdot (E_0\text{v})=-\nabla \text{v:}{\text{P}}_0-\nabla \cdot {\text{F}}_0+4\pi \kappa B-c\kappa E_0,$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\rho v_r)=0,$$

$${\text{P}}_{0ij}=\{\begin{array}{ll}{p}_r=E_0/3& \text{if}i=j\\ 0& \text{if}i\ne j\end{array}$$

$${\text{F}}_0=-\frac{c}{\chi }\nabla p_r,$$

$$F_{0r}=-\frac{c}{\chi }\frac{\partial p_r}{\partial r},F_{0\theta }=-\frac{c}{\chi }\frac{\partial p_r}{r\partial \theta },F_{0\varphi }=0,$$

$$v_r\frac{\partial v_r}{\partial r}-\frac{v_{\varphi }^2}{r}=-\frac{GM}{r^2}-\frac{1}{\rho }\frac{\partial }{\partial r}(p_g+p_r),$$

$$v_{\varphi }^2\cot \theta =\frac{1}{\rho }\frac{\partial }{\partial \theta }(p_g+p_r)$$

$$v_r\frac{\partial }{\partial r}(r{v}_{\varphi })=\frac{1}{\rho r^2}\frac{\partial }{\partial r}(r^3{T}_{r\varphi })$$

$$T_{r\phi}=\eta r\frac{\partial}{\partial r}\big{(}\frac{v_{\phi}}{r}\big{)}$$

$$\eta =\alpha \frac{p_g+p_r}{{\Omega }_K}$$

$$e=\frac{p_g}{\gamma -1}$$

$$\Phi_{vis}=\eta\big{[}r\frac{\partial}{\partial r}\big{(}\frac{v_{\phi}}{r}\big{)}\big{]}^{2}.$$

$$B=\frac{\sigma T^4}{\pi },$$

$$T=\frac{\mu m_p}{k_B}\frac{p_g}{\rho },$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}(\frac{p_g}{\gamma -1}{r}^2{v}_r)=-\frac{p_g}{r^2}\frac{\partial }{\partial r}(r^2{v}_r)-4\kappa \sigma T^4,$$

$$\hskip 8.5359pt+3c\kappa p_{r}+\alpha\frac{p_{g}+p_{r}}{\Omega_{K}}\big{[}r\frac{\partial}{\partial r}\big{(}\frac{v_{\phi}}{r}\big{)}\big{]}^{2},$$

$$\frac{1}{r^{2}}\frac{\partial}{\partial r}(3r^{2}p_{r}v_{r})=\frac{c}{r^{2}}\frac{\partial}{\partial r}\bigg{(}\frac{r^{2}}{\chi}\frac{\partial p_{r}}{\partial r}\bigg{)}+\frac{c}{r\sin\theta}\frac{\partial}{\partial\theta}\bigg{(}\frac{\sin\theta}{\chi}\frac{\partial p_{r}}{r\partial\theta}\bigg{)}$$

$$-\bigg{(}\frac{\partial v_{r}}{\partial r}+2\frac{v_{r}}{r}\bigg{)}p_{r}+4\kappa\sigma T^{4}-3c\kappa p_{r},$$

$$v_i\propto r^{-1/2},\rho \propto r^{-3/2},p_i\propto r^{-5/2},T\propto r^{-1},{\text{F}}_{0i}\propto r^{-2}$$

$$v_{\varphi }^2=v_K^2-\frac{v_r^2}{2}-\frac{5}{2\rho }(p_g+p_r),$$

$$\rho v_{\varphi }^2\cot \theta =\frac{d}{d\theta }(p_g+p_r),$$

$$v_r=-\frac{3\alpha }{2\rho v_K}(p_g+p_r),$$

$$\frac{1}{2}{v}_r(3\rho v_{\varphi }^2-\frac{5-3\gamma }{\gamma -1}{p}_g)=3rc\kappa [p_r-\frac{4\sigma }{3c}{T}^4]$$

$$\frac{3}{2}{v}_r{p}_r+\frac{c}{r}(\cot \theta +\frac{d}{d\theta })[\frac{1}{\chi }\frac{d{p}_r}{d\theta }]=3rc\kappa [p_r-\frac{4\sigma }{3c}{T}^4]$$

Full text

††pagerange: Accretion flows with comparable radiation and gas pressures–LABEL:lastpage††pubyear: 2018

1 INTRODUCTION

In active galactic nuclei (AGN), the major part of energy originates from the vicinity of the central black hole which is occupied by an optically thick accretion flow. Any realistic accretion flow has sufficient angular momentum to form a rotationally supported disc. Due to the liberation of gravitational energy and the existance of turbulent viscosity, the disc can become visible for us. Modeling of such a disc which surrounds a black hole started by pioneers like Zeldovich (1964), Salpeter (1964) and Zeldovich & Novikov (1967). Lynden-Bell (1969) explained that the energy source of quasars is a viscous accretion disc. The standard accretion disc model proposed by Shakura & Sunyaev (1973) became very successful to describe optically thick and geometrically thin accretion discs observed in X-ray binaries. However, this model turned out unable to clarify some observational feauters arising from very hot plasma around black holes. After that popular model, slim disc model were introduced by Abramowicz et al. (1988), presented the properties of those optically thick discs with high mass accretion rates close to its critical values. A few years before it, another popular model had been gradually forming (Ichimaru 1977, Rees et al. 1982 ), and finally became comprehensive and identified as radiatevely inefficient accretion flows (RIAFs) (Narayan & Yi 1994, 1995a,b; Abramowicz et al. 1995; Chen et al. 1997). This model was applied to explain very hard radiation, such as hard X-rays and gamma-rays (up to $100keV$) which were not perceptible by the standard disc model. In the slim discs and RIAFs, the heat energy released by viscous dissipation is stored in the gas and it is advected onto the disc’s center. Therefore, they both belong to advection-dominated accretion flows and can be appropriately described by the self-similar solutions in the radial direction (Narayan & Yi 1994). On the other hand, the presence of radiation can influence significantly the dynamics of the accretion system. Moreover, we expect the radiation pressure exceeds the gas pressure in the central regions, since the accretion rate is greater than about 1 $\%$ of the Eddington value (Turner 2004; Shakura & Sunyaev 1973; Pringle 1981). It seems that without including radiation pressure in the definition of $\alpha-$prescription, we cannot explain observed X-ray luminosities of some compact objects. Near Eddington limit, photons carry a significant fraction of the total momentum, hence radiation also associates with matter for transportation of the angular momentum out of the disc (Jiang et al. 2013). So far, some simulations have been done by several athours (e.g. Eggum et al. 1987, 1988; Okuda et al. 1997; Fujita & Okuda 1998; Kley & Lin 1999; Takeuchi et al. 2009; Sadowski et al. 2014) to study slim discs especially in supercritical accretion regime. Fujita & Okuda (1998) performed 2D calculations to examine subcritical accretion discs with a luminosity comparable with Eddington one. Okuda (2002) noticed that it is possible for a flow which has supercritical accretion rate in the outer boundary layers to become subcritical near the black hole. An analytical investigation for a time-independent case has been done by Gu (2012) (Gu12) which was drawn attention to the vertical structure of such a flow with dominant radiation pressure. In the pioneer theoretical works of Shakura-Sunyaev and slim disks, it is common to consider asymptotic states of the accreting systems which one of the two pressures of gas or radiation is dominant and also they assume just electron scattering has the main role in the opacity of disc. In the standard model, in the region close to the central mass, the radiation pressure reaches to its ultimate limit, and the high amount of luminosity comes out from this part of disc. Further that luminous part, the main assumption changes to the other ultimate limit of pressure where the gas pressure becomes more significant than radiation one. Therefore, we do not expect these ultimate limits occur between two completely discrete neighborhood regions, it does not sound unreasonable to have an expectation to find a transition area between these two territories with both considerable pressures of photons and real particles. On the other hand, we think of another probable situation which is related to being an independent area extended from $r=10r_{g}$ to several tens of Schwarzschild radius with $p_{r}\approx p_{g}$. Another point related to luminosity of the mentioned region might come to our mind: whether the high luminosity can also be produced by a part of the flow where both gas and radiation have comparable portions and a non-negligible number of photons are trapped and therefore this might affect the flow’s dynamics. Moreover, one of basic assumptions of the standard disc is the local energy balance between viscous heating and cooling via radiation but a bit deviation from it leads to extra energy for advection. If we consider conservation energy of matter and radiation separately, we will be able to investigate the variation of energy advection at the intermediate area with mass accretion rate.

In this paper, we follow the similar procedure of Gu12 and employ self-similar solutions to remove $r$- dependency of the system quantities but we aim to study the mentioned different situation where both radiation and gas have comparable portions to form the total pressure of system. The accretion disc’s thickness has been examined by Gu et al. (2009) with assumption of semi-polytropic relationship, $p(\theta)=K\rho(\theta)^{\gamma}$, where $p,\rho$ and $\gamma$ are gas pressure, density and the ratio of specific heats. The important result of that work was to show that the accretion flows in advection-dominated regime tend to be quite thick. In our earlier studies (Samadi et al. 2014 and Samadi & Abbassi 2017), we focused on the effect of magnetic field with toroidal and poloidal configurations, separately. For both cases, the magnetic force acted in the opposite direction of the gas pressure force and made the disc thinner. In the current work, we might expect the radiation pressure behaves like gas pressure and decreases towards the polar axis. The main purpose of this paper is to revisit vertical structure of optically thick accretion flows with rather high accretion rates (saying, several tens of percent of Eddington mass accretion rate) by taking into account photon absorption. The outline of this paper is as follows. In Section 2 we present the basic radiation hydrodynamics equations, which include the two separate equations for energies of gas and radiation. Self-similar equations are employed for a typical radius of $r=30r_{g}$ in Section 3. Section 4 is devoted to boundary conditions and the permissible values of density at the equatorial plane regarding the gas temperature. The results of the numerical solution and a derivation of the disc thickness are presented in Section 5 and, finally, discussions and conclusions are given in Section 6.

2 Equations and assumptions

We study an accreting flow around a black hole with mass $M$, in steady-state regime which is assumed to be axisymmetric and non-relativistic. We define the system in spherical coordinates,$(r,\theta,\phi)$ and we use Newtonian potential, that is, $\Psi=GM/r$. The basic equations to describe such a system, are consist of

1. continuity equation,

[TABLE]

where $\rho$ is density and v is the time-averaged flow velocity.

2. the momentum conservation,

[TABLE]

where $D/Dt=\partial/\partial t+\textbf{v}.\nabla$, and, $p_{g}$ is gas pressure, $\textbf{T}^{\nu}$ is stress tensor, $c$ is light speed,$\textbf{F}_{0}$ is the radiation flux. The total opacity $\chi$ has two parts including Thomson scattering and absorption

[TABLE]

where $\sigma_{T}=6.652\times 10^{-25}\hskip 5.69046ptcm^{2}$ is the Thomson cross-section, $m_{p}$ is the proton mass and $\kappa$ is the absorption opacity. In this problem, we consider two kinds of absorption, saying free-free and bound-free, as

[TABLE]

[TABLE]

So we can find the total absorption opacity as $\kappa=\kappa_{ff}+\kappa_{bf}$.

3. the energy equation of gas,

[TABLE]

where $e$ is the internal energy density, $\Phi_{vis}$ is the viscous dissipative function, $B$ is the blackbody intensity and $E_{0}$ is the radiation energy density.

4. and finally based on conservation of radiation energy density (Mihalas & Mihalas

1984. we have,

[TABLE]

where $\textbf{P}_{0}$ is the radiation pressure tensor. It might be helpful to point out the physical interpretation of new terms which have appeared in Eq.(5) due to considering radiation and also the interaction between photons and gas particles. For instance, $\nabla\textbf{v:}\textbf{P}_{0}$ shows the rate at which radiation pressure is doing work, $\nabla\cdot\textbf{F}_{0}$ specifies the rate of transport of radiation energy, $4\pi\kappa B$ is the rate at which energy of matter is being added to the radiation field (the same term in Eq.(4) appears with negative sign that means decreasing energy of matter due to its radiation) and finally $-c\kappa E_{0}$ denotes the amount of radiation energy which is removed from the total radiation energy and transported to the energy of matter in the unit of time (all terms are per unit volume). In the spherical coordinates, the velocity field has these three components, $(v_{r},v_{\theta},v_{\phi})$ but we follow Narayan & Yi (1995), and assume $v_{\theta}=0$111This assumption is often used for all kinds of accretion discs.. With steady-state, axisymmetric flow and neglecting $v_{\theta}$, Eq. (1) yields,

[TABLE]

For simplicity, we use Eddington approximation for the radiation tensor, that is,

[TABLE]

this means that $\textbf{P}_{0}$ behaves like a scalar and simply equals to radiation pressure, $p_{r}$. The radiation flux vector, $\textbf{F}_{0}=(F_{0r},F_{0\theta},F_{0\phi})$, in diffusion approximation has a relationship with gradient $E_{0}$ as $\textbf{F}_{0}=-c\nabla E_{0}/(3\chi)$, using Eq.(7) we will find:

[TABLE]

Therefore, the components of radiation flux become,

[TABLE]

Now we can find three scalar equations based on Eq.(2) as,

[TABLE]

[TABLE]

[TABLE]

where we have assumed that the $r\phi-$ component of the viscous stress tensor is dominant and we neglect the other components of $T^{\nu}$,

[TABLE]

where $\eta$ is the dynamical viscosity coefficient. We apply $\alpha$-prescription (Shakura & Sunyaev 1973) to determine the dynamical viscosity, means

[TABLE]

where $\Omega_{K}$ is the Keplerian angular velocity. To determine the internal energy density, we can use,

[TABLE]

The last term in Eq. (4) that we need to know is the heating rate of viscosity, $\Phi_{vis}$, given by

[TABLE]

Moreover, $B$ is the blackbody intensity and defined as,

[TABLE]

where $T$ is the gas temperature and $\sigma$ is Stefan-Boltzmann constant. To calculate the gas temperature, we can use the following relationship,

[TABLE]

where $k_{B}$ is the Boltzmann constant and is the mean molecular weight; in the present study we assume $\mu=0.5$. Now substituting Eq.(13)-(17) in Eq.(4) besides our basic assumptions, we obtain,

[TABLE]

[TABLE]

Using Eq.(7), (8), (15) and (16) in the last basic equation of (5) leads us to find,

[TABLE]

[TABLE]

(refer to appendix to see the details of calculation of $\nabla\textbf{v}\textbf{:}\textbf{P}_{0}$).

3 Self-similar solutions

In this work, we are interested in studying vertical structure of the disc with comparable pressures of gas and radiation, so we need to use self-similar technique in the radial direction to remove $r$- dependency of the quantities and find equations with respect to just one variable, means $\theta$. Based on self-similarity we know the dependency of quantities are power law as,

[TABLE]

Employing these solutions in Eq. (9)-(11), (18) and (19) gives,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now, we have six equations (17) and (21-25) including six unknowns ($T,\rho,v_{r},v_{\phi},p_{r}$ and $p_{g}$) that we should find their dependencies to $\theta$. In order to solve the system of ODEs, we should first find boundary conditions that comes in the following.

4 Solutions at the equator

So far, we have obtained a set of six equations that only two of them are differential equations, means Eq.(22) and (25). From these three algebraic Eq. (17), (21) and (23-24) which consist of $T,\rho,v_{r},v_{\phi},p_{g}$ and $p_{r}$, so we can determine only four unknowns. Since we are somewhat familiar with the temperature of gas in optically thick discs, we treat this quantity as a free parameter. We also apply $\beta_{0}$ parameter (zero index denotes the value of each quantity at $\theta=\pi/2$) as the ratio of gas pressure to the total pressure at the equatorial plane,

[TABLE]

It is clear that for radiation pressure supported discs we expect $\beta_{0}\sim 0$ and for gas pressure dominated discs we have $\beta_{0}\sim 1$. Here, it is convenient to fix the equatorial value of gas temperature and calculate the corresponding $p_{r}$ and then find the proper density for several $\beta_{0}$. The solutions of the three-term algebraic equations system are displayed in figures 1 and 2 for $M=10M_{\odot}$, $\mu=0.5$ at the radius $r=30r_{g}$, where $r_{g}(=2GM/c^{2})$ is the Schwarzschild radius. According to Fig.(1), the maximum permissible density at the equator roughly reaches $3\times 10^{-6}\ g\ cm^{-3}$ for $T_{0}\approx 10^{9}\ K$ and $\beta_{0}=0.1$. Moreover, the maximum available density at the mid-plane becomes smaller with greater $\beta_{0}$. The variation of rotational velocity at $\theta=90^{\circ}$ in the second panel of Fig.1 is presented and shows that $v_{\phi 0}$ in fairly cold discs (with temperature less than $10^{8}K$) is approximately Keplerian and in discs with greater radiation pressure it drops more rapidly. Here, we need to take care of the sign of $\nabla\cdot\textbf{F}_{0}$, because it is expected to be positive and show the transportation of energy outwards. The first panel of Fig.(2) reveals that the mentioned term becomes zero at $T_{0}=3.5\times 10^{9}K$ for $\beta_{0}=0.1$. Comparing the second and third panels with the first one, we notice that with greater $\beta_{0}$, the radiation cooling transforms to heating term at larger gas temperatures. At the limit of low temperatures of all panels, the viscous heating is equal to the radiation cooling ($\Phi_{vis0}=(\nabla\cdot\textbf{F}_{0})_{0}$) and energy advection is zero.

5 Vertical Structure

In this section, we refer to the two differential equations of our system, i.e.. equations (22) and (25). For solving Eq.(22), we need to specify one boundary condition. From the previous section, density and gas and radiation pressures corresponding to temperature and a certain value of $\beta_{0}$ are known at the equatorial plane. Similarly for Eq.(25) which is second order ODE, two boundary conditions are required. Here beside $\rho_{0},T_{0},p_{g0}$ and $p_{r0}$, we employ the reflection symmetry about the disc’s equator to assume the derivatives of density and pressures are equal to zero, means

[TABLE]

We start the integrations of ODEs from the equatorial plane and stop our calculations when gas pressure becomes zero. For these sets of input parameters: $M=10M_{\odot},\alpha=0.3,\gamma=1.5$ and for $\beta_{0}=0.1,0.3,0.5$ with these temperatures: $T_{0}/10^{8}K=0.5,2,5$ at $r=30r_{g}$, the integrations might be continued and end at the poles for some other cases. With smaller gas temperature like $5\times 10^{8}K$ that we have used for $\beta_{0}=0.1$, we stop integrations at $\theta_{s}\approx 77^{\circ}$ (we call $\theta_{s}$, the surface angle) because the gas pressure, $p_{g}$ becomes zero at this angle.

Figure (3) shows the angular variation of dimensionless density, radial and azimuthal components of the flow’s velocity and also the two pressures of gas and radiation. All quantities in this figure have been scaled by their boundary values at the equatorial plane which are listed in table (1). According to Fig.(3a), the density of gas starts from a maximum value at the mid-plane and tends to its minimum at the surface angle. Comparing the minimum value of our density here with most of those in ADAFs we have faced with in our previous works, we notice this minimum cannot reach to zero and just becomes about 20% of its equatorial value. However, this result is not strange for theoretical solutions, since similar little ascending trend can be seen in Fig.3 Jiao & Wu 2011 (the minimum value of $\rho$ is about $0.4\rho_{0}$ and also in Fig.1 of Narayan & Yi 1995a, $\rho$ varies by only 10% from the polar axis to the mid-plane with $\epsilon^{\prime}=0.1$ and for $\epsilon^{\prime}=1$ we see $\rho(\theta=90^{\circ})\sim 2\rho(\theta=0)$. In profile (b), it can be seen that the rotational velocity is decreasing outwards, whereas the accretion velocity is minimum at the mid-plane and ends up maximum at the surface. The gas temperature, $T$ (equivalently the sound speed squared $c_{s}^{2}=k_{B}T/\mu m_{p}$ ) peaks at $\theta=90^{\circ}$, that is the same as ADAF’s solutions in other works (such as Fig.1,4 of Samadi et al. 2013, Fig.1 of Gu et al. 2009). The two pressures appear somehow different, $p_{r}$ does not show noticeable change inside the disc especially for lower value of $\beta_{0}$, i.e.. 0.1. As we expect for gas pressure, it faces with a rather sharp decrease towards the surface of disc.

To estimate boundary conditions we have set $\beta_{0}=0.1$ and according to Fig.4, $\beta$ is not constant and goes down and will tend to be zero at the surface, thus radiation is dominant at the outer layers. Here we can substitute the density and radial velocity in the following equation and evaluate the mass accretion rate,

[TABLE]

at $r=30r_{g}$ with $\beta_{0}=0.3$, we find $\dot{M}=0.20\dot{M}_{Edd}$ (see table 2 for two other values of $\beta_{0}$) that $\dot{M}_{Edd}$ is Eddington accretion rate and has this definition: $\dot{M}_{Edd}=4\pi GM/(\eta c\kappa_{esc}$, where $\eta$ is the radiative efficiency of the flow and we chose $\eta=1/16$ since it is comparable to the Schwarzschild black hole efficiency of 0.057 (Gu12). In the presence of radiation flux, the flow’s dynamics is affected by the two effective components of radiation force. According to Fig.(5), the both components, radial $F_{r}=+5cE_{0}/(6\chi r)$ and meridional $F_{\theta}=-c/(3\chi r)dE_{0}/d\theta$ are directed out of the flow. Moreover, the absolute value of the radiation flux force is increasing outwards in the current case. An important aspect of hot accretion flows is their ability to store some of dissipated energy (via turbulent viscosity) as entropy and it is advected with the flow. This advection energy can be written in the following form,

[TABLE]

where $S_{t}$ and $p_{t}$ are the total entropy and total pressure, respectively. Using the definition of Lagrangian time derivative, Eq.(15) and also assumptions of steady-state, $v_{\theta}=0$, self-similar radial solutions in Eq.(27), we obtain

[TABLE]

The advection energy besides two other energies of $\Phi_{vis}$ and $\nabla\cdot\textbf{F}_{0}$ are plotted in Fig.6 for three values of $\beta_{0}$. For all these $\beta_{0}$’s, we see that the amount of $\Phi_{vis}$ is larger than $q_{adv}$ and divergence of radiation flux. As mentioned before, in common solutions the sign of $\nabla\cdot\textbf{F}_{0}$ is expected to be positive and show the cooling term. This expectation is satisfied in the typical solution presented here but not all over the regions of the accretion flow as seen in Fig.6 and it becomes zero at $81.5^{\circ}$ for the disc with$\beta_{0}=0.1$. We speculate some reasons for this issue. First of all, the radiation energy should be transported outside, along the vertical direction over the photosphere, here with just $\theta-$dependency of quantities, we can be sure about agreement between our results in spherical coordinates and those in cylindrical one from the equatorial plane up to just the small distances from the mid-plane. On the other hand, we have ignored $v_{\theta}$ and this causes $v_{r}$ remains negative and even more negative towards the rotation axes. The negative value of $v_{r}$ besides self-similarity solutions leads us to evaluate the left-hand side of both Eq.(4) and (5) with positive sign ($ev_{r}\propto r^{-3}\rightarrow\nabla\cdot(e\textbf{v})=1/r^{2}d/dr(r^{2}ev_{r})=-ev_{r}/r>0$) which means that both the internal gas and radiation density energy are increasing with time in the fixed frame. The main difference between these two equations is that the viscous heating assure us about the presence of a source to increase $e$ but in Eq.(5) instead of that, there is $-\nabla\cdot\textbf{F}_{0}$ which is expected to be negative and act as a cooling term. At regions close enough to the equatorial plane, the last two terms in Eq.(5) by exchanging energy between particles and photons compensate for lack of immediate or direct source of heating to enhance the radiation energy density. The other reason for negative sign of $\nabla\cdot\textbf{F}_{0}$ might be unsatisfactory of self-similar technique near the surface of the disc. Despite being negative of divergence $\textbf{F}_{0}$, the meridional component of radiation flux remains in the same direction as before according to the right panel of figure 5. Equation (28) says that the advection energy depends on both pressures of gas and radiation and radial velocity. Although both gas and radiation pressures are minimum at the surface, the profile of advection energy $q_{adv}$ in Fig.6 does not diminish because of considerable growth in the absolute value of radial velocity. The total energy of advection ($Q_{adv}$), viscosity ($Q_{vis}$) and cooling ($Q_{c}$) at a certain radius can be calculated as,

[TABLE]

[TABLE]

[TABLE]

For the typical solution at $r=30r_{g}$ with $\beta_{0}=0.3$, we find $Q_{adv}=1.33\times 10^{20}\ erg\ cm^{-3}s^{-1}$, $Q_{vis}=8.65\times 10^{20}\ erg\ cm^{-2}s^{-1}$ and $Q_{c}=3.14\times 10^{20}\ erg\ cm^{-3}s^{-1}$ (see table 2 for two other values of $\beta_{0}$). Therefore, the vertically averaged advection factor becomes,

[TABLE]

which is approximately equal to $0.15$.

In figure (7), we have plotted (a) the half-thickness of disc (b) the total mass accretion rate and (c) advection parameter versus radius for three values of $\beta_{0}$. As seen in panels of this figure, all these quantities increase at larger radii. Furthermore, in the last figure of this section, fig.8, we see a direct effect from mass accretion rate on the advection factor. The similar result have been presented by Gu (2012) in figure 4 of his paper.

5.1 Convective Stability

The last thing which would be our interest to examine is the convection instability in the optically thick accretion discs with comparable gas and radiation pressures. To do this, we should refer to Solberg-Høiland conditions for dynamical stability which are defined in cylindrical coordinates $(R,\phi,z)$ as following,

[TABLE]

[TABLE]

where $j=Rv_{\phi}=r\sin\theta v_{\phi}$ is the specific angular momentum per unit mass and the differential of entropy is defined by $dS=d\ln(p_{t}/\rho^{\gamma})$. On the other hand, the oscillation frequencies related to the internal hydrodynamic modes are driven by equilibrium gradients. In this way, the square of the epicyclic frequency is obtained as,

[TABLE]

The two components of the Brunt–Väisälä frequency are given by

[TABLE]

[TABLE]

Substituting the last three equations in Eq.(31) gives the first condition of a convectively stable flow as,

[TABLE]

where we have named the total frequency as $N_{t}^{2}$. Furthermore, the second condition is simplified as

[TABLE]

as seen, we have omitted the factor $\partial p_{t}/\partial z$ because this term is usually negative in accretion discs, in addition our solutions of $p_{r}$ and $p_{g}$ in figure 3 are in agreement with it. Since we have found the angular dependency of the quantities, it is convenient to transform derivatives with respect to $R$ and $z$ into the spherical coordinates, by means of the chain rule,

[TABLE]

[TABLE]

By using the self-similar solutions and our results from the previous section we have calculated these frequencies and shown in figure 9. The radial component of Brunt–Väisälä frequency, $N_{R}^{2}$ is negative and the other component of it is positive. As it is seen in fig.9, the total frequency of $N_{t}^{2}$ and $\Delta_{ls}$ for all three are positive, thus the current disc is convectively stable.

6 Summary and conclusion

We studied the vertical structure of optically thick accretion flows in the presence of both radiation and gas pressures. We did not employ radiative transfer, instead of that we adopted the separated energy equation similar to the internal energy of gas based on diffusion approximation in radiation field. The self-similar solutions in the radial direction helped us to convert the partial differential equations into a set of algebraic equations and ODE’s. Solving equations on the equator and considering some basic assumptions like high optical depth, we could obtain proper values of density and gas temperature. At the typical radius of $r=30r_{g}$, we assumed the gas pressure had 10-50% portion of the total pressure and calculated the $\theta-$dependency of physical quantities. With a permissible gas temperature at mid-plane ($T_{0}/10^{9}K=0.5,2,5$), we calculated corresponding gas pressure which yields $\beta_{0}=0.1,0.3,0.5$. The velocity’s components at the equatorial plane were found from two algebraic equations. With these known quantities, we could start integrating the differential equations and then stopped calculation because of gas pressure which tended to negative values. In fact, the gas pressure showed a descending behaviour and it became gradually zero at a certain angle ($0<\theta_{0}<90$) before approaching the vertical axis, whereas the radiation pressure faced with a little change at this angular distance. Therefore, we expect to find accretion discs with two vertical parts, one is located inside from the mid-plane to $\theta_{0}$, and affected by both $p_{g}$ and $p_{r}$, and the other one as a surrounding layer of the equatorial part whose radiation pressure is dominant and hence it has different features. Although the temperature of gas at the inside region was rather high at the equator, it became negligible at the surface of the disc, meanwhile the temperature at the outer layer is determined based on radiation and it must be very high. Furthermore, we found out with greater accretion rate and less portion of gas pressure in the total pressure, more energy was advected to the central object. In our solutions, we realized that despite of the outward direction of radiation flux, it is possible that divergence of the radiation flux to become negative. However, our results for other input parameters are not practical at outer regions (especially for hotter discs) because we neglected the meridional component of the velocity and it is just ignorable near the mid-plane but $v_{\theta}$ is very effective to make $v_{r}$ positive and change significantly the angular profile of energy transportation outward the flow. Moreover, the typical accretion flows presented here with relatively we saw them convectively stable. Our results here are useful to apply for intermediate cases between slim and thin discs which include two parts, one part near the mid-plane which has nonignorable gas pressure and the other part with dominant radiation pressure and negligible gas pressure (about 3 order of magnitude smaller). The vertical structure of inner part can be found by the presented method but it is incomplete. The outer part can be described by omitting the gas pressure and use one energy equation for the total energy of the flow like which has been done by Gu12. We can use the surface value of quantities as boundary conditions for solving differential equations in the outer part.

Acknowledgements

We are grateful to the anonymous referee for his/her thoughtful and constructive comments, which helped us to improve the first version of this paper. Maryam Samadi also wishes to thank Faride Danesh-Manesh for her kind helps in some parts of this project. This work has been supported financially by Research Institute for Astronomy & Astrophysics of Maragha (RIAAM) under research project No. 1/6025-24.

Appendix A Formulas of Tensor production

To simplify the tensor production of $\nabla\textbf{v:}\textbf{P}_{0}$ in the radiation energy equation (Eq.5) in the spherical coordinate, we present the basic formulas here to calculate it conveniently. The double dot (or scalar, or inner) product produces a scalar, (Phan-Thien 2013)

[TABLE]

The gradient of a vector is given by

[TABLE]

[TABLE]

[TABLE]

So the production of $\nabla\textbf{v}$ and $\textbf{P}_{0}$ becomes,

[TABLE]

[TABLE]

with the assumptions of $v_{\theta}=0$, $\partial/\partial\phi=0$ and $\textbf{P}_{0ii}=E_{0}/3$ we have,

[TABLE]