© ESO, 2011

1. Introduction

Modeling accretion flows onto black holes is crucial for understanding the energetic emissions observed from many sources, both Galactic and extragalactic (e.g., X-ray binaries, ULXs, AGNs). The accreting matter commonly settles into a disk-like configuration in which angular momentum is transported outwards as a result of shear in the differentially rotating fluid, allowing mass to be accreted onto the central compact object. Of particular interest are sources such as the low-mass X-ray binaries in which matter is transferred onto the X-ray source from a binary companion, usually a late-type star. When containing a black hole, these sources undergo irregular outbursts (lasting several months or years) that are analogous to the outbursts observed in dwarf novae, and are thought to be related to the viscous-thermal instability, which leads to an enhanced transfer of angular momentum in the accreting fluid (Lasota 2001).

In outburst, the source is a luminous X-ray emitter and as the luminosity decays from the outburst maximum to minimum, the disk is thought to follow a sequence of quasi-stationary states, whose spectra can be derived from steady-state models of accretion disks. At high luminosities describing the inner regions of such accretion disks by the most-often used “thin disk” models is no longer valid. A description of such disks in terms of “slim disk” models is more appropriate.

Following the approach of Shakura & Sunyaev (1973), in the “standard" discussion of stationary disk structure one assumes that there is no radial advection of heat (the flow is “radiatively efficient”), the radial pressure gradients are negligible, and the distribution of angular momentum is Keplerian. All of these assumptions hold for very thin disks, including the original Shakura & Sunyaev (1973) model (henceforth SS) and its general-relativistic version (NT), which was developed for the Kerr geometry by Novikov & Thorne (1973). However, it has long been realized that not all disk-like accretion flows satisfy these assumptions In particular, slim disks (e.g., Abramowicz et al. 1988) form a stable branch of accretion, in which advection of entropy is important and the disk is not necessarily geometrically thin. These solutions are particularly relevant to sources with high accretion rates, and they converge to the standard case in the low accretion rate limit.

In the case of thin disks, methods for treating the radiative transfer and computing the vertical structure of the disk have been developed by Shaviv & Wehrse (1986) and Hubeny (1998) and applied by, e.g., Davis & Hubeny (2006), Idan et al. (2008), and Różańska & Madej (2008). These have not yet been applied to slim disks. Until recently, the properties of slim disk models have been represented with quantities averaged over the vertical structure of the disk. This paper presents improved advective steady-state models of accretion, in which the vertical structure of the “slim disk” is explicitly taken into account.

Recently, Sądowski (2009) revisited slim disk models with an improved numerical code, while Sądowski et al. (2009) performed a preliminary study of the slim disk vertical structure. While our present work is based on these two papers, it breaks with the long tradition of computing the radial dependence of disk quantities, including the emergent flux F(r), before the vertical structure is analyzed. In this paper, by closely coupling the radial properties of the disk to its vertical structure, we offer a more consistent treatment of the slim disk. A similar approach has been taken by Dotan & Shaviv (2010) in their discussion of winds from Super-Eddington slim disks in a pseudo-Newtonian potential.

The high luminosities observed in the bright X-ray sources imply there is an effective mechanism of angular momentum transport in the accreting fluid, and much work has been devoted to MHD simulations of one such mechanism, the magneto-rotational instability (MRI). For thin disks, the effectiveness of angular momentum transport is conveniently parametrized by the alpha parameter that was introduced by Shakura & Sunyaev (1973). Observations of systems in which variability is driven by accretion disk instabilities suggest rather high values of α for fully ionized disks, much higher than the value α = 0.001 that was adopted in the original slim-disk article (Abramowicz et al. 1988). Smak (1999) has already showed that one has to take α ≈ 0.2 to describe hot disks in dwarf novae. As reviewed by King et al. (2007), the value of α in hot, fully ionized disks has to be  ~0.1−0.4, whereas in cold protostellar and FU Urionis disks, much lower values of this parameter are required  ~0.01 and  ~0.001−0.003, respectively. Numerical simulations give a value of α ~ 0.01 for ionized disks, ten times lower than suggested by observations, and it probably reflects the limitations of the present MRI calculations (King et al. 2007). However, the apparently observational determinations of α mentioned above are in fact strongly model-dependent. In the case of dwarf-novae, for example, these determinations assume that outbursts are described well by the thermal-viscous instability model. However, the case of the epitome dwarf-nova SS Cyg shows that such a description might be inadequate (Schreiber & Lasota 2007; Smak 2010), putting presumably observational determinations into doubt. Therefore, one cannot exclude that the numerical simulations are correct after all.

In this paper we adopt the value α = 0.01. With such a low value of α, the resulting disk models are optically thick (τ_eff > 1), thus allowing a simplified treatment of radiative transfer. We treat the vertical energy transport in the diffusive approximation for the radiative flux (and account for vertical convection in the mixing-length approximation).

We begin the body of the paper with a brief discussion of the model, then both the radial and the vertical equations of disk structure are presented. In Sect. 3 we describe the numerical method used to solve the problem. The radial structure of the solutions is described in Sect. 4.1, while the vertical profiles of various disk quantities are discussed in Sect. 4.2. By comparing our results with previous work, we find that, with a judicious choice of parameters, a standard polytropic slim disk model can be made to approximate our full solution. Analytical formulae for the appropriate values of the parameters can be found in Sect. 5. Finally, in the last section we discuss the limitations of the model and its possible applications.

2. Slim-disk structure

Slim disks are alpha disks. Like Shakura-Sunyaev disks they are based on the assumption that the dissipation mechanisms operating in accretion flows may be described by a viscous stress tensor whose leading component is proportional to the pressure.

Slim disk solutions, again just like the SS and NT solutions, are obtained by solving vertically averaged (or height-integrated), radial equations of motion. Thus, steady slim disks, like the SS and NT ones, neglect the vertical structure of flow, and describe essentially flat fluid configurations. Although an expression for the radial dependence of disk thickness is obtained, the slope of the disk surface is usually neglected in the discussion of emergent spectra, where a plane-parallel atmosphere is typically considered for the purposes of computing radiative transfer in the vertical direction.

In the case of slim disks, the radial equations of structure are a set of ordinary differential equations (ODEs). For a detailed account of these Kerr-metric equations (and state-of-the-art traditional slim-disk models), see Sądowski (2009), who followed previous work beginning with Lasota (1994), and considered subsequent improvements (Abramowicz et al. 1996, 1997; Gammie & Popham 1998). The set of equations is closed by including relations describing vertical hydrostatic equilibrium and vertical transport of energy. By solving these equations with appropriate boundary (or regularity) conditions, one obtains the radial profiles of the central temperature, T_c(r), surface density, Σ(r), the height-integrated pressure, P(r), the half-thickness of the disk, h(r), a radial velocity V(r), the emerging flux of radiation, F(r), and certain other physical quantities.

In this paper we are taking a first step towards constructing truly three-dimensional slim disk models, by solving a set of differential equations describing the vertical structure of the disk. The resulting z-dependence of physical quantities is used to compute certain coefficients that enter the radial equations and that up till now have been estimated with algebraic expressions. We refer to the resulting models as “two-dimensional (2-D) slim disk” solutions, although it has to be understood that in contrast to the thin-disk solutions of Urpin (1984) and Kluźniak & Kita (2000), the actual meridional flow has not been computed in this paper, but the average radial velocity alone has been considered in the structure equations. However, the models are now self-consistent in the sense that the vertical averages of physical quantities that form the coefficients of the radial ODEs do correspond to the vertical structure considered in the radiative transfer calculation – in previous work the vertical structure of the radiative atmosphere was considered a posteriori, and it had no influence on the radial structure of the disk.

2.1. Basic assumptions, parameters, and coefficients

We assume an axially symmetric, stationary fluid configuration in the Kerr metric, with fixed values of the black hole (BH) mass, M, and spin, a, parameters and the fluid disk is symmetric under reflection in the equatorial plane of the metric. Matter is supplied at a steady rate, Ṁ, through a boundary “at infinity” and angular momentum is removed through the same boundary (in practice, we use the NT solution for the outer boundary condition), whereas zero torque is assumed at the BH horizon. We assume that no mass or angular momentum crosses the disk surface. We neglect the loss of angular momentum to both wind and radiation. We do not consider self-irradiation of the disk and assume that the magnetic pressure may be neglected. Neglecting the incoming radiation may not be justified for super-Eddington accretion rates for which the disk is geometrically thick.

A fraction, (1 + f^adv(r))^-1, of the entropy generated locally by dissipative processes is released into the radiation field, while the remainder is advected by the gas.

A unique solution to the slim-disk model can only be found if certain additional assumptions are made. We make the following arbitrary choice. We neglect the vertical variation (z dependence) of the velocity field, considering only its height-averaged value; thus, the velocity is always directed radially inwards and is a function of the radial coordinate alone. Similarly, we assume there is no z variation of the advection factor f^adv(r). Dissipation and angular momentum transport are given by the alpha prescription (Shakura & Sunyaev 1973), with a constant value of α. We assume that the dissipation rate is proportional to the total pressure, p. For a more detailed statement, see Eq. (13) and the comment following it. Calculations are carried out for the value α = 0.01.

We are looking for 2-D slim-disk solutions at a definite value of mass accretion rate for a given Kerr metric. Thus, for a fixed value of α, there are three fundamental parameters describing a given slim-disk solution: M, a, and Ṁ.

In the structure equations, we take G = c = 1 and make use of the following expressions involving the BH spin: $\begin{array}{ccc}\mathrm{\Delta }& \mathrm{=}& {\mathit{r}}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{Mr}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathit{,}\\ \mathit{A}& \mathrm{=}& {\mathit{r}}^{\mathrm{4}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{Mr}{\mathit{a}}^{\mathrm{2}}\mathit{,}\\ \mathrm{𝒞}& \mathrm{=}& \mathrm{1}\mathrm{-}\mathrm{3}{\mathit{r}}_{\mathrm{\ast }}^{-1}\mathrm{+}\mathrm{2}{\mathit{a}}_{\mathrm{\ast }}{\mathit{r}}_{\mathrm{\ast }}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\\ \mathrm{𝒟}& \mathrm{=}& \mathrm{1}\mathrm{-}\mathrm{2}{\mathit{r}}_{\mathrm{\ast }}^{-1}\mathrm{+}\mathrm{2}{\mathit{a}}_{\mathrm{\ast }}^{\mathrm{2}}{\mathit{r}}_{\mathrm{\ast }}^{-2}\\ \mathrm{\mathscr{H}}& \mathrm{=}& \mathrm{1}\mathrm{-}\mathrm{4}{\mathit{a}}_{\mathrm{\ast }}{\mathit{r}}_{\mathrm{\ast }}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{3}{\mathit{a}}_{\mathrm{\ast }}^{\mathrm{2}}{\mathit{r}}_{\mathrm{\ast }}^{-2}\end{array}$with a_∗ = a/M and r_∗ = r/M.

We use the Boyer-Lindquist system of coordinates and introduce the vertical coordinate z = rcosθ. The metric near the equatorial plane takes the form (Abramowicz et al. 1996), $\mathrm{d}{\mathit{s}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{{\mathit{r}}^{\mathrm{2}}\mathrm{\Delta }}{\mathit{A}}\mathrm{d}{\mathit{t}}^{\mathrm{2}}\mathrm{+}\frac{\mathit{A}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{(}\mathrm{d}\mathit{\phi }\mathrm{-}\mathit{\omega }\mathrm{d}\mathit{t}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{r}}^{\mathrm{2}}}{\mathrm{\Delta }}\mathrm{d}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{d}{\mathit{z}}^{\mathrm{2}}\mathit{,}$(3)where ω = 2Mar/A is the angular velocity of the frame dragging.

The radial gas (three-)velocity, V, as measured by an observer co-rotating with the fluid at a fixed value of r, is given by the relation (Abramowicz et al. 1996) $\mathit{V}\mathit{/}\sqrt{\mathrm{1}\mathrm{-}{\mathit{V}}^{\mathrm{2}}}\mathrm{=}{\mathit{u}}^{\mathit{r}}{\mathit{g}}_{\mathit{rr}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\frac{\mathit{r}{\mathit{u}}^{\mathit{r}}}{{\mathrm{\Delta }}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,}$(4)where u^r is the contravariant radial component of the fluid 4-velocity and has the dimension of physical velocity. The disk surface density is , while the height-integrated pressure is . The total pressure is the sum of gas and radiation pressures, p = p_gas + p_rad. We adopt an equation of state corresponding to the choice p_gas = kρT/(μm_p), and p_rad = aT⁴/3, with k the Boltzmann constant, m_p the proton mass, and a the radiation constant (no confusion with the spin parameter may arise). The mean molecular weight is taken to be μ = 0.62, but see the comment following Eq. (21). The height-integrated energy density is $\mathit{E}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}\left(\frac{{\mathit{p}}_{\mathrm{gas}}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{+}\mathrm{3}{\mathit{p}}_{\mathrm{rad}}\right)\mathrm{d}\mathit{z,}$(5)where γ = 5/3.

The following averages (moments) enter the radial equations as coefficients of certain terms: $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{1}}& \mathrm{\equiv }& \frac{\mathrm{1}}{{\mathit{T}}_{\mathrm{c}}^{\mathrm{4}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathit{h}}{\mathit{T}}^{\mathrm{4}}\mathrm{d}\mathit{z,}\\ {\mathit{\eta }}_{\mathrm{2}}& \mathrm{\equiv }& \frac{\mathrm{2}}{\mathrm{\Sigma }{\mathit{T}}_{\mathrm{c}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathit{h}}\mathit{\rho T}\mathrm{d}\mathit{z,}\\ {\mathit{\eta }}_{\mathrm{3}}& \mathrm{\equiv }& \mathit{E}\mathit{/}\mathit{P,}\\ {\mathit{\eta }}_{\mathrm{4}}& \mathrm{\equiv }& \frac{\mathrm{1}}{\mathrm{\Sigma }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{h}}\mathit{\rho }{\mathit{z}}^{\mathrm{2}}\mathrm{d}\mathit{z}\mathit{.}\end{array}$All integrals in this section are taken at a fixed value of r. Here, T(z) is the gas temperature, and T_c is its value at the equatorial plane, T(0) = T_c.

The vertical epicyclic frequency squared, which can be thought of as the vertical component of gravity, is (Kato 1993), ${\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{\equiv }\frac{\mathit{M}}{{\mathit{r}}^{\mathrm{3}}}\frac{\mathrm{\mathscr{H}}}{\mathrm{𝒞}}\mathrm{·}$(10)We also define $\begin{array}{}\mathrm{˜}\\ {\mathrm{\Gamma }}_{\mathrm{1}}\end{array}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{3}}}\mathrm{·}$(11)

2.2. Vertical structure equations

We describe the vertical structure of an accretion disk in the optically thick regime by the following equations:

• (i)

  Hydrostatic equilibrium (Katoet al. 2008),$\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathit{z,}$(12)where Ω _⊥  is defined in Eq. (10). Although other expressions for the right-hand side of Eq. (12) can be found in the literature (e.g., Abramowicz et al. 1997, includes v^z ≠ 0), the form above is appropriate for our scheme, in which the vertical structure is precalculated before any information about the radial variables becomes available. Thus, in Eqs. (10) and (12), we assume Keplerian angular velocity , as well as v_z = 0 in hydrostatic equilibrium.

• (ii)

  The energy generation equation. We assume that the vertical flux of energy inside the disk ℱ is generated according to $\frac{\mathrm{dℱ}}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathrm{3}\mathrm{𝒟}}{\mathrm{2}\mathrm{𝒞}}\left(\frac{\mathit{\alpha p}}{\mathrm{1}\mathrm{+}{\mathit{f}}^{\mathrm{adv}}}\right){\left(\frac{\mathit{M}}{{\mathit{r}}^{\mathrm{3}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(13)Strictly speaking, this does not correspond to a constant α prescription, as the term (3/2)(M/r³)^1/2, which is derived from Keplerian strain, departs somewhat from the value that would follow from the actually computed distribution of angular momentum (cf., Fig. 1). However, as the departure for sub-Eddingtonian accretion rates is small, we expect Eq. (13) to afford a good approximation. Note that f^adv = 0 corresponds to the NT disk (going over into the Shakura-Sunyaev disk in the nonrelativistic limit of thin disks), f^adv > 1 characterizes advection-dominated disks, while f^adv < 0 describes those disk regions where the advected heat is being released. The amount of heat advected Q^adv = f^advℱ(h).

• (iii)

  Energy transport. The structure of the disk has to be such that the actual value of the divergence of the flux corresponds to Eq. (13). Radiative transport is computed in the diffusive approximation $\mathrm{ℱ}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{16}\mathit{\sigma }{\mathit{T}}^{\mathrm{3}}}{\mathrm{3}\mathit{\kappa \rho }}\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{z}}\mathit{,}$(14)while convective transport is computed in a mixing-length approximation. Energy is transported in the vertical direction through diffusion of radiation or convection according to the value of the thermodynamical gradient, which can be either radiative or convective. Accordingly, we take $\frac{\mathrm{dln}\mathit{T}}{\mathrm{dln}\mathit{p}}\mathrm{=}\{\begin{array}{c}\\ {\mathrm{\nabla }}_{\mathrm{rad}}\mathit{,}& \mathrm{for}& {\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{\le }{\mathrm{\nabla }}_{\mathrm{ad}}\\ {\mathrm{\nabla }}_{\mathrm{conv}}\mathit{,}& \mathrm{for}& {\mathrm{\nabla }}_{\mathrm{rad}}\mathit{>}{\mathrm{\nabla }}_{\mathrm{ad}}\end{array}$(15)with the adiabatic gradient given by a derivative at constant entropy: ∇_ad ≡ (∂lnT/∂lnp)_S. The radiative gradient ∇_rad is calculated in the diffusive approximation, ${\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{=}\frac{\mathrm{3}\mathit{p}{\mathit{\kappa }}_{\mathrm{R}}\mathrm{ℱ}}{\mathrm{16}\mathit{\sigma }{\mathit{T}}^{\mathrm{4}}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathit{z}}\mathit{,}$(16)where κ_R is the Rosseland mean opacity, and σ is the Stefan-Boltzmann constant. At the equatorial plane we apply the boundary conditions described in point (iv) below. When the temperature gradient exceeds the value of the adiabatic gradient, we have to consider the convective energy flux. The convective gradient ∇_conv is calculated using the mixing length theory introduced by Paczyński (1969). We take the following mixing length, ${\mathit{H}}_{\mathrm{ml}}\mathrm{=}\mathrm{1.0}{\mathit{H}}_{\mathrm{p}}\mathit{,}$(17)with pressure scale height H_p defined as (Hameury et al. 1998) ${\mathit{H}}_{\mathrm{p}}\mathrm{=}\frac{\mathit{p}}{\mathit{\rho }{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathit{z}\mathrm{+}\sqrt{\mathit{p\rho }}{\mathrm{\Omega }}_{\mathrm{\perp }}}\mathrm{·}$(18)The convective gradient is defined by the formula ${\mathrm{\nabla }}_{\mathrm{conv}}\mathrm{=}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{+}\mathrm{(}{\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}\mathit{y}\mathrm{(}\mathit{y}\mathrm{+}\mathit{w}\mathrm{)}$(19)where y is the solution of the equation $\frac{\mathrm{9}}{\mathrm{4}}\frac{{\mathit{\tau }}_{\mathrm{ml}}^{\mathrm{2}}}{\mathrm{3}\mathrm{+}{\mathit{\tau }}_{\mathrm{ml}}^{\mathrm{2}}}{\mathit{y}}^{\mathrm{3}}\mathrm{+}\mathit{w}{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{w}}^{\mathrm{2}}\mathit{y}\mathrm{-}\mathit{w}\mathrm{=}\mathrm{0}\mathit{,}$(20)with the typical optical depth for convection τ_ml = ρκ_RH_ml, and w given by $\frac{\mathrm{1}}{{\mathit{w}}^{\mathrm{2}}}\mathrm{=}{\left(\frac{\mathrm{3}\mathrm{+}{\mathit{\tau }}_{\mathrm{ml}}^{\mathrm{2}}}{\mathrm{3}{\mathit{\tau }}_{\mathrm{ml}}}\right)}^{\mathrm{2}}\frac{{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathit{z}{\mathit{H}}_{\mathrm{ml}}^{\mathrm{2}}{\mathit{\rho }}^{\mathrm{2}}{\mathit{C}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{512}{\mathit{\sigma }}^{\mathrm{2}}{\mathit{T}}^{\mathrm{6}}{\mathit{H}}_{\mathrm{P}}}{\left(\frac{\mathit{\partial }\ln \mathit{\rho }}{\mathit{\partial }\mathrm{lnT}}\right)}_{\mathrm{p}}\left({\mathrm{\nabla }}_{\mathrm{rad}}\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\right)\mathit{.}$(21)The thermodynamical quantities C_p, ∇_ad and (∂lnρ/∂lnT)_p are calculated using standard formulae (e.g., Chandrasekhar 1967) assuming solar abundances (X = 0.70, Y = 0.28) and, when necessary, taking the effect of partial ionization of gas on the gas mean molecular weight into account . We use Rosseland mean opacities κ_R (including the processes of absorption and scattering) taken from Alexander et al. (1983) and Seaton et al. (1994). Following other authors (e.g., Idan et al. 2008), we neglect expansion opacities, in agreement with our neglect of vertical velocity gradients.

• (iv)

  We set the following boundary conditions. At the equatorial plane (z = 0) we set ℱ(0) = 0, in accordance with the assumption of reflection symmetry, while at the disk surface, τ(h) = 0, we follow the Eddington approximation (Mihalas 1982) and require . In practice, for a fixed r, and prescribed values of T_c and f^adv, a trial value of the central density, ρ_c, is assumed and the equations are integrated in z until ρ(z_∗) = 10^-16gcm^-3 (as a stand-in for the disk surface, z_∗ = h). If ℱ(h) and T(h) fail to satisfy the surface boundary condition, the assumed value of ρ_c is adjusted, and the integration is repeated, until the condition ℱ(h) = 2σT⁴(h) is met. Convergence is usually attained in a few iterations. The emergent flux of radiation at any given r is then F = ℱ(h).

2.3. Radial structure equations

The radial sector of the model is described by four laws of conservation and a regularity condition:

• (i)

  Mass conservation, $\mathit{Ṁ}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi }\mathrm{\Sigma }\mathit{r}{\mathit{u}}^{\mathit{r}}\mathit{.}$(22)

• (ii)

  Conservation of angular momentum, $\frac{\mathit{Ṁ}}{\mathrm{2}\mathit{\pi }}\mathrm{(}\mathrm{ℒ}\mathrm{-}{\mathrm{ℒ}}_{\mathrm{in}}\mathrm{)}\mathrm{=}\frac{{\mathit{A}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{\Delta }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\Gamma }}{\mathit{r}}\mathit{\alpha P,}$(23)where ℒ = u_φ is the specific angular momentum, ℒ_in is a constant, whose value is to be specified later, and Γ is the Lorentz factor (Gammie & Popham 1998): ${\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}{\mathit{V}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathrm{ℒ}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}}{\mathit{A}}\mathrm{·}$

• (iii)

  Conservation of radial momentum, $\frac{\mathit{V}}{\mathrm{1}\mathrm{-}{\mathit{V}}^{\mathrm{2}}}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathrm{𝒜}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{\Sigma }}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathit{,}$(24)where $\mathrm{𝒜}\mathrm{=}\mathrm{-}\frac{\mathit{MA}}{{\mathit{r}}^{\mathrm{3}}\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{k}}^{\mathrm{+}}{\mathrm{\Omega }}_{\mathrm{k}}^{\mathrm{-}}}\frac{\mathrm{(}\mathrm{\Omega }\mathrm{-}{\mathrm{\Omega }}_{\mathrm{k}}^{\mathrm{+}}\mathrm{)}\mathrm{(}\mathrm{\Omega }\mathrm{-}{\mathrm{\Omega }}_{\mathrm{k}}^{\mathrm{-}}\mathrm{)}}{\mathrm{1}\mathrm{-}\begin{array}{}\mathrm{˜}\\ {\mathrm{\Omega }}^{\mathrm{2}}\end{array}\begin{array}{}\mathrm{˜}\\ {\mathit{R}}^{\mathrm{2}}\end{array}}\mathit{,}$(25)and Ω = u^φ/u^t is the angular velocity with respect to a stationary observer, is the angular velocity with respect to an inertial observer, are the angular frequencies of the corotating and counterrotating Keplerian orbits and is the radius of gyration. The value of P is taken from vertical structure solutions for given values of Σ and T_c, $\mathit{P}\mathrm{=}{\mathit{\eta }}_{\mathrm{2}}\frac{\mathit{k}}{\mathit{\mu }{\mathit{m}}_{\mathrm{p}}}\mathrm{\Sigma }{\mathit{T}}_{\mathrm{c}}\mathrm{+}{\mathit{\eta }}_{\mathrm{1}}\frac{\mathrm{2}}{\mathrm{3}}\mathit{a}{\mathit{T}}_{\mathrm{c}}^{\mathrm{4}}\mathrm{;}$(26)hence, the radial derivative of the height-integrated pressure takes the form $\begin{array}{ccc}\frac{\mathrm{d}\ln \mathit{P}}{\mathrm{d}\mathit{r}}& \mathrm{=}& \mathrm{(}\mathrm{4}\mathrm{-}\mathrm{3}\mathit{\beta }\mathrm{)}\frac{\mathrm{d}\ln {\mathit{T}}_{\mathrm{c}}}{\mathrm{d}\mathit{r}}\mathrm{+}\mathit{\beta }\frac{\mathrm{d}\ln \mathrm{\Sigma }}{\mathrm{d}\mathit{r}}\\ & & \mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\frac{\mathrm{d}\ln {\mathit{\eta }}_{\mathrm{1}}}{\mathrm{d}\mathit{r}}\mathrm{+}\mathit{\beta }\frac{\mathrm{d}\ln {\mathit{\eta }}_{\mathrm{2}}}{\mathrm{d}\mathit{r}}\mathit{,}\end{array}$(27)with β = η₂(k/μm_p)ΣT_c/P.

• (iv)

  Energy conservation. The advective cooling is defined following Kato et al. (2008) in terms of the vertically integrated quantities, as ${\mathit{Q}}^{\mathrm{adv}}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\mathrm{\left(}\mathit{r}{\mathit{u}}^{\mathit{r}}\mathrm{\right(}\mathit{E}\mathrm{+}\mathit{P}\mathrm{\left)}\mathrm{\right)}\mathrm{-}{\mathit{u}}^{\mathit{r}}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}\mathrm{-}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}{\mathit{u}}^{\mathit{z}}\frac{\mathit{\partial p}}{\mathit{\partial z}}\mathrm{d}\mathit{z}\mathit{.}$(28)Using mass conservation (Eq. (22)) and hydrostatic equilibrium (Eq. (12)), the expression can be rewritten as $\begin{array}{ccc}{\mathit{Q}}^{\mathrm{adv}}& \mathrm{=}& \mathrm{-}\frac{\mathit{Ṁ}}{\mathrm{2}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}}({\mathit{\eta }}_{\mathrm{3}}\frac{\mathit{P}}{\mathrm{\Sigma }}\frac{\mathrm{dln}\mathit{P}}{\mathrm{dln}\mathit{r}}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\eta }}_{\mathrm{3}}\mathrm{)}\frac{\mathit{P}}{\mathrm{\Sigma }}\frac{\mathrm{dln}\mathrm{\Sigma }}{\mathrm{dln}\mathit{r}}\mathrm{+}\\ & & \mathrm{+}{\mathit{\eta }}_{\mathrm{3}}\frac{\mathit{P}}{\mathrm{\Sigma }}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{3}}}{\mathrm{dln}\mathit{r}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{4}}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{4}}}{\mathrm{dln}\mathit{r}})\mathrm{·}\end{array}$(29)Just like P, the advective cooling term, Q^adv, and the coefficients η₃, η₄ are all determined from the vertical structure solutions. The manipulation of the last term in Eq. (28) was as follows: $\begin{array}{ccc}& & \mathrm{-}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}{\mathit{u}}^{\mathit{z}}\frac{\mathit{\partial p}}{\mathit{\partial z}}\mathrm{d}\mathit{z}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}{\mathit{u}}^{\mathit{z}}\mathit{\rho z}\mathrm{d}\mathit{z}\\ & & \mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}\frac{\mathit{\partial }}{\mathit{\partial z}}\mathrm{\left(}{\mathit{u}}^{\mathit{z}}\mathit{\rho }\mathrm{\right)}{\mathit{z}}^{\mathrm{2}}\mathrm{d}\mathit{z}\\ & & \mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathit{r\rho }{\mathit{u}}^{\mathit{r}}\mathrm{\right)}{\mathit{z}}^{\mathrm{2}}\mathrm{d}\mathit{z}\\ & & \mathrm{=}\mathrm{-}\frac{\mathit{Ṁ}}{\mathrm{2}\mathit{\pi r}}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\frac{\mathrm{d}}{\mathrm{d}\mathit{r}}\left(\frac{\mathrm{1}}{\mathrm{\Sigma }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{h}}\mathit{\rho }{\mathit{z}}^{\mathrm{2}}\mathrm{d}\mathit{z}\right)\\ & & \mathrm{=}\mathrm{-}\frac{\mathit{Ṁ}}{\mathrm{2}\mathit{\pi r}}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{4}}}{\mathrm{d}\mathit{r}}\mathrm{·}\end{array}$

3. Numerical method

3.1. Vertical structure

The set of ordinary differential equations describing the vertical structure, i.e., Eqs. (12), (13) and (15) together with appropriate boundary conditions are solved for a given BH spin on a three-dimensional grid spanned by the radius r, the central temperature T_c, and the advection factor f^adv. For a given set of these parameters, we start the integration from the equatorial plane (z = 0), and the solution satisfying the outer boundary condition is found as described at the end of Sect. 2.2. The resulting quantities describing the vertical structure (T_c, Σ, Q^adv, P, η₁, η₂, η₃, and η₄), together with r, are printed out to tables for subsequent use in interpolation routines. As it turns out, the first two of these parameters, T_c and Σ, can be used to uniquely determine all the other quantities characterizing the vertical structure, including f^adv (see Fig. 2).

Calculating the full grid of vertical solutions for a single value of BH spin takes about 5 h on a 4-CPU workstation.

3.2. Radial structure

By a series of algebraic manipulations of Eqs. (22)–(24), and (29), we obtain the following set of two ordinary differential equations for V(r) and T_c(r) $\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}{\mathit{V}}^{\mathrm{2}}}\frac{\mathrm{d}\ln \mathit{V}}{\mathrm{d}\ln \mathit{r}}\mathrm{=}\frac{{\mathrm{𝒩}}_{\mathrm{1}}}{{\mathrm{𝒟}}_{\mathrm{0}}}\mathit{,}$(30)$\begin{array}{ccc}\mathrm{(}\mathrm{4}\mathrm{-}\mathrm{3}\mathit{\beta }\mathrm{)}\frac{\mathrm{d}\ln {\mathit{T}}_{\mathrm{c}}}{\mathrm{d}\ln \mathit{r}}\mathrm{=}\left(\frac{{\mathrm{𝒩}}_{\mathrm{1}}}{{\mathrm{𝒟}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{r}\mathrm{(}\mathit{r}\mathrm{-}\mathit{M}\mathrm{)}}{\mathrm{\Delta }}\right)\mathrm{(}\mathit{\beta }\mathrm{-}\begin{array}{}\mathrm{˜}\\ {\mathrm{\Gamma }}_{\mathrm{1}}\end{array}\mathrm{)}& & \\ \mathrm{-}\frac{\mathrm{2}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}}{\mathit{Ṁ}{\mathit{\eta }}_{\mathrm{3}}}\frac{\mathrm{\Sigma }}{\mathit{P}}{\mathit{Q}}^{\mathrm{adv}}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{1}}}{\mathrm{dln}\mathit{r}}\mathrm{-}\mathit{\beta }\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{2}}}{\mathrm{dln}\mathit{r}}& & \\ \mathrm{-}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{3}}}{\mathrm{dln}\mathit{r}}\mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{\Sigma }}{\mathit{P}}\frac{{\mathit{\eta }}_{\mathrm{4}}}{{\mathit{\eta }}_{\mathrm{3}}}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{4}}}{\mathrm{dln}\mathit{r}}& & \end{array}$(31)with and given by $\begin{array}{ccc}{\mathrm{𝒩}}_{\mathrm{1}}& \mathrm{=}& \mathrm{𝒜}\mathrm{+}\frac{\mathrm{2}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}}{\mathit{Ṁ}{\mathit{\eta }}_{\mathrm{3}}}{\mathit{Q}}^{\mathrm{adv}}\mathrm{+}\frac{\mathit{P}}{\mathrm{\Sigma }}(\frac{\mathit{r}\mathrm{(}\mathit{r}\mathrm{-}\mathit{M}\mathrm{)}}{\mathrm{\Delta }}\begin{array}{}\mathrm{˜}\\ {\mathrm{\Gamma }}_{\mathrm{1}}\end{array}\\ & & \mathrm{+}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{3}}}{\mathrm{dln}\mathit{r}})\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\frac{{\mathit{\eta }}_{\mathrm{4}}}{{\mathit{\eta }}_{\mathrm{3}}}\frac{\mathrm{dln}{\mathit{\eta }}_{\mathrm{4}}}{\mathrm{dln}\mathit{r}}\end{array}$(32)$\begin{array}{ccc}{\mathrm{𝒟}}_{\mathrm{0}}& \mathrm{=}& {\mathit{V}}^{\mathrm{2}}\mathrm{-}\begin{array}{}\mathrm{˜}\\ {\mathrm{\Gamma }}_{\mathrm{1}}\end{array}\frac{\mathit{P}}{\mathrm{\Sigma }}\mathrm{·}\end{array}$(33)Typically, vanishes close to the black hole, as ℒ approaches ℒ_in.

To obtain a solution one has to solve this system of two ordinary differential equations, together with the following regularity conditions at the sonic radius r_S, defined by the same conditions: ${\mathrm{𝒩}|}_{{\mathit{r}}_{\mathrm{S}}}\mathrm{=}{{\mathrm{𝒟}}_{\mathrm{0}}|}_{{\mathit{r}}_{\mathrm{S}}}\mathrm{=}\mathrm{0}\mathit{,}$(34)as well as outer boundary conditions given at some large radius r_out. The solution between the outer boundary and the sonic point is found using the relaxation technique (Press 2002), with ℒ_in treated as the eigenvalue of the problem. The sonic point for a_∗ = 0 is located at 5.9M for accretion rate 0.01Ṁ_Edd and at 5.0M for 2.0Ṁ_Edd. To start the relaxation process we have to provide a trial solution that is obtained by a method similar to the one described in Sądowski (2009) assuming the Novikov & Thorne (1973) outer boundary conditions. Once the trial solution is found, one can start the relaxation process with a free inner boundary corresponding to the location of the sonic point. To find the solution inside the sonic point we make a small step inward to cross the critical point and then integrate down to BH horizon using a Runge-Kutta method of the fourth order. In this work we use 25 mesh points spaced logarithmically in the radius on the section between the sonic point and r_out = 1000M. This particular number of grid points is enough to resolve all disk features. We have verified that the results are accurately reproduced with a denser grid.

The parameters linked to the vertical structure (P, Q^adv, η₁, η₂, η₃, and η₄) for given Σ and T_c are linearly interpolated from pre-calculated tables of the vertical structure solutions (for any value of V, Σ is determined directly from mass conservation, Eqs. (4), (22)). The radial derivatives dlnη₁/dlnr, dlnη₂/dlnr, dlnη₃/dlnr, and dlnη₄/dlnr are evaluated numerically from the η₁, η₂, η₃, and η₄ profiles in the previous iteration step. A relaxed solution is obtained in a few iteration steps. Once a solution outside the sonic point is found we numerically estimate the radial derivatives of V and T_c at the sonic point using values given at r > r_S and use these derivatives to start direct integration inside the sonic point.

The solution thus obtained may then be used as a trial solution when looking for the relaxation solution of another slim disk, i.e., when one of the three fundamental parameters (Ṁ, a, M) has a slightly different value. Each relaxation step takes approximately 5 seconds on a single-CPU workstation.

4. Results

In the following two sections we present and discuss both the radial and vertical structure of slim accretion disks. All the solutions, if not stated otherwise, were computed assuming α = 0.01 and M = 10   M_⊙.

4.1. Disk radial structure

Angular momentum

The angular momentum profiles for our accretion disk solutions near a nonrotating BH are presented in Fig. 1. Results for two values of mass accretion rate are shown, Ṁ = 0.1M_Edd and 2.0M_Edd, where Ṁ_Edd = 16L_Edd/c² is the critical accretion rate that for a disk around a nonrotating BH approximately corresponds to the Eddington luminosity, L_Edd. For the lowest accretion rates the profiles follow the Keplerian profile and reach its minimal value (ℒ_in in Eq. (23)) at the marginally stable orbit (ISCO, Bardeen et al. 1972). The higher the accretion rate, the stronger the deviation from the Keplerian profile. The disk is sub-Keplerian at large distances and super-Keplerian at moderate radii. The Keplerian profile is crossed again at a point located inside the marginally stable orbit, and corresponding to what is usually called “the cusp” or “the potential spout”. For a detailed study of the physics of the inner edge of a see Abramowicz et al. (2010).

S-curves

Figure 2 presents slim disk solutions at r = 20M on the T_c-Σ plane, for a nonrotating BH. Solutions of the polytropic, height-averaged models are presented for comparison; for detailed discussion see Sect. 5. The locus of solutions for various values of the mass accretion rate has the shape of the so-called “S-curve” (Abramowicz et al. 1988). The lower, gas-pressure dominated branch accurately follows the track of radiatively efficient solutions (f^adv = 0). The middle, radiation-pressure dominated branch is reached at Ṁ ≈ 0.1Ṁ_Edd. As advection becomes significant, the slim-disk solution leaves the f^adv = 0 track and moves to higher advection rates (S-curve). Around Ṁ = 5Ṁ_Edd, the solutions enter the upper advection-dominated branch corresponding to f^adv > 1.0 (more than 50% of heat stored in the accreted gas). At Ṁ = 20Ṁ_Edd this rate almost increases up to 80% (f^adv ≈ 4.0).

Surface density

Profiles of the surface density for 2-D slim-disk solutions for a non-rotating BH are presented in the upper panel of Fig. 3. Different regimes, corresponding to different branches of the “S-curve” on the (Σ, T_c) plane are visible. For large radii the surface density increases with increasing accretion rate (the lower gas-pressure dominated branch), while this relation is opposite for moderate radii (the middle radiation-pressure dominated branch). For accretion rates Ṁ > 5.0Ṁ_Edd the upper advection-dominated branch would be reached. The local maxima in the surface density profiles (discussed in detail in, e.g., Sądowski 2009) are visible for moderate accretion rates (~0.5 Ṁ_Edd). The bottom panel of Fig. 3 presents corresponding profiles of the radial velocity V as measured by an observer corotating with the fluid.

The surface density dependence on BH rotation is presented in Fig. 4. The profiles are shifted to lower radii as the inner edge of the disk moves inward for higher BH spins. The outer parts of the accretion disk are insensitive to the metric.

Optical depth

In the top panel of Fig. 5 we plot the total optical depth of the vertical slim disk solutions for different accretion rates (a_∗ = 0). The total optical depth

${\mathit{\tau }}_{\mathrm{tot}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{h}}{\mathit{\kappa }}_{\mathrm{R}}\mathit{\rho } \mathrm{dz}\mathit{,}$where the total opacity coefficient κ_R, which includes the processes of absorption and scattering, is closely related to the surface density. Indeed, the radial profiles of the optical depth shown in Fig. 5 follow the corresponding profiles of surface density. Any differences in the profiles come from the dependence of the opacity coefficient on local density and temperature. Outside the ISCO the total optical depth is always large (τ_tot > 10³).

The diffusive approximation for radiative transport may only be used if photons are absorbed, otherwise LTE cannot be established. In a scattering-dominated atmosphere, the effective optical depth is then the relevant quantity to be used in checking for the self-consistency of the diffusive approximation. The bottom panel of Fig. 5 presents corresponding profiles of the effective optical depth, which is estimated in the following way:

${\mathit{\tau }}_{\mathrm{eff}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{h}}\sqrt{\mathrm{(}{\mathit{\kappa }}_{\mathrm{R}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{es}}\mathrm{)}{\mathit{\kappa }}_{\mathrm{R}}} \mathit{\rho }\mathrm{dz}$where κ_es = 0.34   cm²g^-1 is the mean opacity for Thomson electron scattering. For Ṁ > 0.3Ṁ_Edd the inner region of the disk becomes radiation-pressure dominated (it enters the middle branch on the corresponding S-curve), the surface density decreases with increasing accretion rate, and electron scattering begins to dominate absorption. Therefore, the effective optical depth decreases with increasing accretion rate and reaches values τ_eff < 1 in the inner parts of the disk for accretion rates above 1.0Ṁ_Edd. As a result, for Ṁ > 1.0Ṁ_Edd, the diffusive approximation can no longer be applied, and our present approach to solving for the disk structure breaks down. In Fig. 6 we exhibit the dependence of the effective optical depth on the value of α. In general, τ_eff is inversely proportional to α (as is the surface density). At lower accretion rates (Ṁ < 0.1Ṁ_Edd) the effective optical depth remains large even for high values of α.

Flux profiles

Profiles of the flux emitted from the disk surface () in the case of a nonrotating BH are presented in Fig. 7. Results corresponding to accretion rates from 0.01 up to 5.0Ṁ_Edd are shown. For the lowest rates the emission from inside the marginally stable orbit is negligible as expected in the standard accretion disk models. This is no longer true for higher accretion rates, and the advection of energy causes significant emission from smaller radii (Abramowicz et al. 2010). For super-Eddington accretion rates the emitted flux continues to grow with a decreasing radius even inside the marginally stable orbit. Radiation coming from the direct vicinity of the black hole is suppressed by the gravitational redshift (the g-factor). Therefore, an observer at infinity will observe a maximum in the profile of the effective temperature even for the highest accretion rates.

The increase in the advective flux with increasing accretion rate is clearly visible in Fig. 8. The ratio of the heat advected to the amount of energy emitted is presented for different accretion rates. For very low accretion rates these profiles approach the limit of a radiatively efficient disk (f^adv ≡ 0), and the advection component becomes significant for higher accretion rates. Some part (up to 30% at r = 20M for 2.0Ṁ_Edd) of the energy generated at moderate radii is advected with matter and radiated away at r < 10M. This causes the significant change in the emitted flux profile at the higher accretion rates visible in Fig. 7.

In Fig. 9 we present the emitted flux profiles for different BH angular momenta at a constant accretion rate Ṁ = 0.1Ṁ_Edd. These profiles coincide at large radii where the influence of the BH rotation is negligible; however, the higher the BH spin, the closer to the horizon the marginally stable orbit. Therefore, in the case of rotating BHs, the accreting matter can move much deeper into the gravity well, compared to nonrotating BHs. This effect leads to an increase in the disk luminosity and hardening of its spectrum, which can be inferred from Fig. 9 – the higher the spin, the higher the disk luminosity, and the higher the flux (which corresponds to the effective temperature).

Photosphere location

The flux observed at infinity may be obtained by performing ray tracing of photons emitted from the accretion disk (see Sect. 5). Scattering in the layer above the photosphere must also be taken into account. An accurate calculation requires detailed knowledge of atmospheric properties, including the location of the photosphere. This is particularly important when accretion rates are high and the disk is no longer geometrically thin (Sądowski et al. 2009).

In Fig. 10 we plot the profiles for the z-location of the photosphere, H_phot, obtained in our model at different accretion rates for a nonrotating black hole. Clearly, for high accretion rates, Ṁ > 0.1Ṁ_Edd, the inner regions become thicker (effects of radiation pressure). For Ṁ = 1.0Ṁ_Edd, the ratio H_phot/r reaches a value as high as 0.25. Close to the ISCO the height of the photosphere rapidly decreases because of vigorous cooling (compare Fig. 8). Although the rapid change in disk thickness violates the assumption of the hydrostatic equilibrium that we make when solving for the disk vertical structure, the accelerations connected with the vertical motions involved are much lower than the vertical component of gravity. Indeed, the vertical accelerations are close to . In deriving this estimate we liberally assumed that dH/dr ~ 1 and used the fact that the rapid decrease in disk height occurs near the sonic point, while the speed of sound v_s is approximately rΩ _⊥ (H/r). Thus the acceleration terms modifying Eq. (12) would be smaller than the gravitational acceleration by a factor of a few percent: (H/r) ~ 10^-1. In Fig. 12 we plot the dynamical and gravitational components of the vertical equilibrium equation at the photosphere. It is clear that the former is at least 10 times smaller than the latter at the sonic radii for Ṁ ≤ Ṁ_Edd. Therefore, in all likelihood, our results correctly describe the disk structure that would be obtained without assuming strict hydrostatic equilibrium.

In Fig. 11 we present radial photosphere profiles, at a fixed accretion rate and different values of the BH spin. The photospheric heights coincide for large radii. In the inner regions of the disk, the height of the photosphere increases with BH spin, reflecting the increased luminosity and radiation pressure.

4.2. Vertical structure

In Fig. 13 we present a vertical cross-section of a Schwarzschild slim disk for Ṁ = 0.01Ṁ_Edd. At this accretion rate the disk is radiatively efficient and no advection of entropy is expected. The top panel presents the radial profiles of the photosphere and the disk surface (defined as a layer with ρ = 10^-16g/cm³).

The total optical depth reaches values as high as  ~20000 on the equatorial plane and decreases monotonically towards the photosphere at τ = 2/3 (Fig. 13, second panel from the top). Within the marginally stable orbit, the total optical depth significantly decreases, as shown in Fig. 5.

The third panel of Fig. 13 presents the temperature, which decreases with height from T_c at the equatorial plane to T_eff at the photosphere. The maximum value of the temperature is attained on the equatorial plane at r ~ 10M.

Density is presented in the fourth panel, and the next panel presents the vertical flux that is generated inside the disk according to Eq. (13). It is set to zero on the equatorial plane by reflection symmetry, and then rapidly increases, because in an alpha disk the dissipation is proportional to the pressure, which reaches its maximum on the z = 0 plane. Close to the disk surface, where pressure is almost negligible, the flux slowly settles down to the emitted value. At the accretion rate chosen for the figure the flux rapidly decreases inside the marginally stable orbit.

The bottom panel of Fig. 13 presents the nonmonotonic distribution of the termodynamical gradient (Eq. (15)), which ranges between 0.2 and 0.4. Therefore, the disk’s vertical structure cannot be described by a simple polytropic relation. Moreover, in the region of the highest temperature (close to the equatorial plane at moderate radii, r ≈ 20M), the heat is transported upward through convection.

The vertical structure of an accretion disk with ten times higher accretion rate, Ṁ = 0.1Ṁ_Edd, is presented in Fig. 14. The general picture remains the same, since the advection of heat is still insignificant. However, as the temperatures increase, the inner disk regions become dominated by radiation pressure. For 0.1Ṁ_Edd the convective region extends from the marginally stable orbit up to 300M and covers more than half of the disk thickness.

The disk structure is significantly different in the case of a high accretion rate (e.g., 1.0Ṁ_Edd), with a significant amount of advection. The vertical cross-sections of the slim disk are presented in Fig. 15. The inner regions are dominated by radiation pressure, so the disk geometrically thickens and the photosphere is now higher. The total optical depth mostly follows the surface density (compare Fig. 5), and therefore decreases considerably towards the black hole in the inner parts of the disk.

The temperature maximum is again located at the equatorial plane close to r ≈ 10M. The maximum of the effective temperature (corresponding to the vertical flux shown in the fifth panel (compare also Fig. 7) is shifted inwards, down to r ≈ 8M.

The fourth panel of Fig. 15 presents the density distribution. Despite the fact that the surface density monotonically increases outwards (see Fig. 3), ρ has two maxima in the equatorial plane: at r ≈ 200M and r ≈ 6M. Finally, the bottom panel presents the termodynamical gradient distribution. For such a high accretion rate, the convective zone extends nearly to the photosphere for r < 200M and is present up to r = 2000M.

In Fig. 16 we present the density in the meridional plane for three different accretion rates. The equatorial plane lies in the middle of each plot. The violet boundaries denote the photosphere. The two maxima of density at the equatorial plane for Ṁ = 1.0Ṁ_Edd) are clearly visible also in this representation, as are some other features discussed above.

5. Comparison with height-averaged slim disk solutions

In this section we compare our 2-D slim disk solutions, where the radial structure equations are coupled to those for the vertical structure, with the standard polytropic slim disk model, in which the slim-disk equations and properties are averaged over the thickness of the disk (e.g., Kato et al. 2008). We note here that the slim disk solutions presented in one of our previous papers (Sądowski 2009) did not follow the polytropic formalism, assuming different relations between the vertically integrated quantities and their values on the equatorial plane (e.g., Σ = 2ρ₀H instead of Σ = 2I_Nρ₀H, where I_N is defined in Eq. (37)) In this work, we also use a more general form of the energy equation (compare our Eq. (29) with Eq. (6) in Sądowski 2009).

Let us assume the polytropic equation of state with the polytropic index N: p = Kρ^1 + 1/N. The vertical integration of the hydrostatic equilibrium formula, Eq. (12), gives $\mathit{\rho }\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\mathrm{1}\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{H}}^{\mathrm{2}}}\right)}^{\mathit{N}}\mathrm{·}$(35)We also assume $\mathit{T}\mathrm{=}{\mathit{T}}_{\mathrm{0}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{H}}^{\mathrm{2}}}\right)\mathrm{·}$(36)One can now calculate analytical formulae for η₁ to η₄ (Sect. 2.3): $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{T}}_{\mathrm{0}}^{\mathrm{4}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{H}}{\mathit{T}}^{\mathrm{4}}\mathrm{d}\mathit{z}\mathrm{=}{\mathit{I}}_{\mathrm{4}}\mathit{H}\\ {\mathit{\eta }}_{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{2}}{\mathrm{\Sigma }{\mathit{T}}_{\mathrm{0}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{H}}\mathit{\rho T}\mathrm{d}\mathit{z}\mathrm{=}{\mathit{I}}_{\mathit{N}\mathrm{+}\mathrm{1}}\mathit{/}{\mathit{I}}_{\mathit{N}}\\ {\mathit{\eta }}_{\mathrm{3}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{P}}\left(\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\frac{\mathit{k}}{\mathit{\mu }}\frac{{\mathit{I}}_{\mathit{N}\mathrm{+}\mathrm{1}}}{{\mathit{I}}_{\mathit{N}}}\mathrm{\Sigma }{\mathit{T}}_{\mathrm{C}}\mathrm{+}\mathrm{2}{\mathit{I}}_{\mathrm{4}}\mathit{a}{\mathit{T}}_{\mathrm{C}}^{\mathrm{4}}\mathit{H}\right)\\ {\mathit{\eta }}_{\mathrm{4}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{\Sigma }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{H}}\mathit{\rho }{\mathit{z}}^{\mathrm{2}}\mathrm{d}\mathit{z}\mathrm{=}{\mathit{J}}_{\mathit{N}}{\mathit{H}}^{\mathrm{2}}\mathit{,}\end{array}$where ${\mathit{I}}_{\mathit{N}}\mathrm{=}\frac{\sqrt{\mathit{\pi }}}{\mathrm{2}}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{N}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{+}\mathit{N}\mathrm{)}}\stackrel{\mathit{N}\mathrm{\in }N}{\mathrm{=}}\frac{\mathrm{(}{\mathrm{2}}^{\mathit{N}}\mathit{N}\mathrm{!}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{N}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}$(37)${\mathit{J}}_{\mathit{N}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{+}\mathit{N}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{2}\mathrm{+}\mathit{N}\mathrm{)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{6}\mathrm{+}\mathrm{4}\mathit{N}}\mathrm{·}$(38)In this approach we do not solve the vertical structure consistently, so we need to make some additional assumptions about the vertical equilibrium of forces and radiation transfer. Following other authors, we simplify the hydrostatic equilibrium (Eq. (12)) by applying a finite difference approximation and write ${\mathit{H}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{\perp }}^{\mathrm{2}}\mathrm{=}\mathrm{(}\mathrm{2}\mathit{N}\mathrm{+}\mathrm{3}\mathrm{)}\frac{\mathit{P}}{\mathrm{\Sigma }}\mathrm{·}$(39)One has to remember that disk thickness, H, defined in this way is not the exact location of the photosphere. Therefore we introduce a factor f_H relating these quantities, ${\mathit{H}}_{\mathrm{phot}}\mathrm{=}{\mathit{f}}_{\mathit{H}}\mathit{H}\mathit{.}$(40)Assuming that radiation is transported in the vertical direction through diffusion, the radiative flux is given by $\mathrm{ℱ}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{16}\mathit{\sigma }{\mathit{T}}^{\mathrm{3}}}{\mathrm{3}\mathit{\kappa \rho }}\frac{\mathrm{d}\mathit{T}}{\mathrm{d}\mathit{z}}\mathrm{·}$(41)Under the one-zone approximation one obtains the following formula for the total flux emitted from disk surface,

$\mathit{F}\mathrm{=}{\mathit{f}}_{\mathit{F}}\frac{\mathrm{64}\mathit{\sigma }{\mathit{T}}_{\mathrm{C}}^{\mathrm{4}}}{\mathrm{3}\mathrm{\Sigma }\mathit{\kappa }}\mathrm{·}$(42)where factor f_F has been introduced to account for inaccuracies arising from this approximation, as well as from the dominance of the disk convection in certain regions (as discussed in Sect. 4.2). Now, advective cooling takes the form ${\mathit{Q}}^{\mathrm{adv}}\mathrm{=}{\mathit{f}}^{\mathrm{adv}}\mathit{F}\mathrm{=}\mathrm{-}\mathit{\alpha P}\frac{\mathit{A}{\mathit{\gamma }}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{3}}}\frac{\mathrm{d\Omega }}{\mathrm{d}\mathit{r}}\mathrm{-}{\mathit{f}}_{\mathit{F}}\frac{\mathrm{64}\mathit{\sigma }{\mathit{T}}_{\mathrm{C}}^{\mathrm{4}}}{\mathrm{3}\mathrm{\Sigma }\mathit{\kappa }}\mathrm{·}$(43)Usually, the following values of N, f_H and f_F are assumed: $\begin{array}{ccc}\mathit{N}& \mathrm{=}& \mathrm{3.0}\\ {\mathit{f}}_{\mathit{H}}& \mathrm{=}& \mathrm{1.0}\\ {\mathit{f}}_{\mathit{F}}& \mathrm{=}& \mathrm{1.0}\end{array}$(44)In Fig. 17 we compare radial profiles of the flux, photospheric height and surface density of the 2-D slim disk model described in this paper (consistently taking its vertical structure into account) with profiles obtained for the polytropic, height-averaged, slim disk, as presented in this section with parameters defined in Eq. (44). The comparison was carried out for two accretion rates (0.1 and 1.0Ṁ_Edd), with α = 0.01. For the lower accretion rate (0.1Ṁ_Edd), the disk is radiatively efficient (advective cooling is negligible), and therefore both flux profiles almost coincide and correspond to the Novikov & Thorne solutions. However, the disk photosphere location in the polytropic model turns out to be overestimated by more than 20% in the region corresponding to the maximal emission. The profiles of the surface density do not coincide either, the 2-D solutions giving values  ~25% lower than the corresponding polytropic solutions (compare Fig. 2).

For the higher accretion rate (1.0Ṁ_Edd), the flux profiles remain similar (up to 1%). However, as advection becomes important, the emission is shifted inwards with respect to the Novikov & Thorne profile. The photosphere location in the polytropic model is overestimated by  ~30% and the surface density by  ~20%.

A question arises as to whether such differences in the flux, photosphere, and surface density profiles affect the resulting disk spectrum. In Fig. 18 we present spectral profiles and their ratios (2-D to polytropic) for two accretion rates. The spectra were calculated with ray-tracing routines (Bursa 2006) using the BHSPEC package (Davis et al. 2005), assuming the inclination angle i = 70° and distance to the observer d = 10kpc. As BHSPEC gives tabulated solutions of the full, frequency-dependent, radiative transfer equations for the disk vertical structure taking the Compton scatterings into account in the disk atmosphere, the spectra presented in Fig. 18 are not those of a simple multi-color blackbody. However, one should be aware that using BHSPEC for calculating spectral color correction is not consistent with the assumed vertical structure (calculated or height-averaged) because it is based on a stand-alone disk model.

The general shape of the spectra is similar for both types of slim disk models, because the emission profiles nearly coincide. However, the spectra are not identical. The bottom panel of Fig. 18 presents ratios of the spectral profiles of the corresponding solutions for each accretion rate. They all coincide att low energies (<0.1keV), while for higher energies the discrepancies are as large as 3% for 1.0Ṁ_Edd at 5   keV. These differences are atributed to (slightly) different profiles of the flux, the photosphere, and the surface density. We conclude that the proper treatment of the vertical structure hardly affects the spectra of slim disks for this range of accretion rates (Ṁ < Ṁ_Edd) and the viscosity parameter (α ≤ 0.01).

We end this section by giving fitting formulae for N, f_H, and f_F, which approximate the full numerical 2-D solutions described in this paper with a polytropic slim disk model. An advantage of using these formulae lies in avoiding the need to perform time-consuming calculations of the vertical structure and avoiding numerical problems connected to interpolation in the vertical solutions grid. The formulae for the polytropic model parameters for α = 0.01 are $\begin{array}{ccc}\mathit{N}& \mathrm{=}& \mathrm{3.25}\mathrm{\times }{\mathrm{𝒮}}_{\mathit{N}}\\ {\mathit{f}}_{\mathit{H}}& \mathrm{=}& \mathrm{0.63}\mathrm{\times }{\mathrm{𝒮}}_{\mathit{H}}\\ {\mathit{f}}_{\mathit{F}}& \mathrm{=}& \mathrm{0.94}\mathrm{\times }{\mathrm{𝒮}}_{\mathit{F}}\mathit{,}\end{array}$(45)where the spin correction coefficients , , and are given by $\begin{array}{ccc}{\mathrm{𝒮}}_{\mathit{N}}& \mathrm{=}& \mathrm{1}\mathrm{+}\mathrm{0.002} \mathrm{(}\mathrm{6}\mathrm{-}{\mathit{r}}_{\mathrm{ms}}\mathit{/}\mathit{M}\mathrm{)}\mathit{,}\\ {\mathrm{𝒮}}_{\mathit{H}}& \mathrm{=}& \mathrm{1}\mathrm{+}\mathrm{0.003} \mathrm{(}\mathrm{6}\mathrm{-}{\mathit{r}}_{\mathrm{ms}}\mathit{/}\mathit{M}\mathrm{)}\mathit{,}\\ {\mathrm{𝒮}}_{\mathit{F}}& \mathrm{=}& \mathrm{1}\mathrm{+}\mathrm{0.064} \mathrm{(}\mathrm{6}\mathrm{-}{\mathit{r}}_{\mathrm{ms}}\mathit{/}\mathit{M}\mathrm{)}\mathit{.}\end{array}$Here, r_ms is the radius of the marginally stable orbit. For a nonrotating BH, at the radius r = 7M (corresponding to the highest disk effective temperature), the fitting formulae are accurate to 1% for the emitted flux, the photospheric height and the surface density.

6. Discussion

Motivated by a desire to explain and fit the observed spectra of accreting black holes in binary systems to theoretical models, we have developed a 2-D model of optically thick slim disks. These should be particularly relevant to transient binaries. In quiescence, their inner disk regions are described well by optically thin advection-dominated accretion flows (ADAFs; see e.g., Lasota et al. 1996; Dubus et al. 2001). However, a few black hole systems (e.g., GRS 1915+105 and LMC X-3) have been observed in thermal states corresponding to disk luminosities higher than 0.3L_Edd (McClintock & Remillard 2003; Steiner et al. 2010), and modeling these require optically thick models going beyond the standard thin disks.

In this work we present a 2-D slim disk model, in which the radial and vertical structures are coupled. Such an approach eliminates arbitrary factors that influence solutions of the usual polytropic slim disk model. The results were obtained under two key assumptions: an alpha disk was assumed (dissipation proportional to pressure), with a uniform value of α, and the fraction of the generated entropy that is advected was computed at every radius under the assumption that this fraction does not vary with the height above the disk plane (Eq. (13)). Both of these assumptions seem arbitrary, and we can offer no physical motivation for the (conventional) choice we made.

Under these assumptions and for the value α = 0.01 of the viscosity parameter, we computed and presented the detailed structure of 2-D slim disks, parametrized by the mass accretion rate, and the two Kerr metric parameters, M and a. Somewhat surprisingly, the spectra observed at infinity from such disks differ by only a few percent from those obtained from previously considered slim disk models (in which the equations and structure correspond to a height average over a polytropic atmosphere). Such differences are unlikely to introduce any large corrections to spin measurements based on X-ray continuum fits made with corresponding height-averaged polytropic models of slim disks. However, already the latter produce significantly softer spectra in the sub-Eddington regime than the Novikov & Thorne model. For high luminosities, fits based on slim disk models may therefore provide higher values of the black hole spin parameter than corresponding fits based on the kerrbb model (e.g., Shafee et al. 2006). This issue will be discussed in detail in a forthcoming paper (Bursa & Sądowski, in prep.).

One has to be aware that the model of vertical structure presented here is only an approximation of the real physical processes taking place in disk interiors. The diffusion approximation and the convection treatment in the mixing length approach are known to successfully describe media with large effective optical depths but break down when the disk becomes optically thin. We have shown that the effective optical depth of slim accretion disks may drop below unity for super-Eddington luminosities and sufficiently high values of α. For such conditions, a more sophisticated model of radiation transfer should be implemented. However, for α ≤ 0.01 and L ≤ L_Edd the assumptions of this work are self-consistent. For higher values of α their range of applicability is limited to lower luminosities (e.g., to 0.5L_Edd for α = 0.1). Kerr-metric slim disks with low effective optical depth were discussed by Beloborodov (1998), who finds them to be significantly hotter than the optically thick ones. We have already started implementing a radiative transfer scheme valid for disks with small effective optical depths into the scheme introduced in this work. It will be presented and discussed in a future paper.

Another remark is connected to the fact that one can expect winds to be blown out of the disk surface at super-Eddington luminosities. Such a phenomenon may significantly change the disk structure, e.g., its thickness. This feature of slim disks, not described in our calculations, has recently been cleverly modeled by Dotan & Shaviv (2010).

Acknowledgments

This work was supported in part by Polish Ministry of Science grants N203 0093/1466, N203 304035, N203 380336, N N203 381436. J.P.L. acknowledges support from the French Space Agency CNES, MB from ESA PECS project No. 98040. We thank the anonymous referee for valuable comments.

References

All Figures