© ESO, 2012

1. Introduction

There is a consensus that a standard model does exist for geometrically thin accretion disks (see e.g. Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Pringle 1981). Their structure is obtained assuming steady accretion and hydrodynamic equilibrium. While the viscosity is needed to ensure accretion, it can be completely eliminated in favour of the mass accretion rate – once the steady disk solution is known to exist – from the explicit formulation of structure equations. The luminosity can be inferred from the disk geometry and the mass accretion rate, but it does not influence the structure of steadily accreting disks.

In contrast to that, there is no generally accepted model for geometrically thick disks. We review just a sample of the related literature. Paczyński (1978) did not solve the full system of thick-disk equations, but assumed elements of thin-disk approximation in determining the vertical structure. Paczyński & Wiita (1980) imposed some ingredients of thin-disk models onto thick disks. The luminosity was obtained analogously to thin disks – it was deduced from the already known disk structure. In a similar analysis Abramowicz et al. (1980) assumed test gas approximation and adapted the mass accretion rate in a way suitable to ensure the energy conservation. In a scenario discussed in Paczyński & Abramowicz (1982) the mass accretion rate has been given as part of primary data and the bulk of the disk was convective. This description has been applied to supercritical accretion, but the selfgravity of the disk was not included. Hashimoto et al. (1993) and Hashimoto et al. (1995) considered optically and geometrically thick toroidal stars and Keplerian disks, with the selfgravity taken into account. They observed that the luminosity can be significantly higher than the Eddington limit. In these constructions the mass accretion rate has been estimated a posteriori.

The main goal of this investigation is to find a simple selfconsistent model of steady accretion with radiation and to study its structure. We considered a scenario with a geometrically thick self-gravitating Newtonian disk (as in Hashimoto et al. 1995), in which the mass accretion rate is prescribed a priori (as in Paczyński & Abramowicz 1982). We assumed, following Paczyński & Abramowicz (1982), that matter spirals inwards in the immediate vicinity of the central disk plane z = 0, where radiation is generated, by a phenomenological mechanism that is inessential for the purpose of this analysis. This radiation interacts with the gas in the bulk of the disk through Thompson scattering (i.e., we employed the thin-gas approximation). The structure of the bulk is then obtained from the hydrodynamic equations (that include a radiation force) supplemented by the energy conservation equation. We demonstrate that a self-contained description of accretion disks emerges – a simple variant of radiation hydrodynamics. The whole problem is then analysed numerically, without any other simplifications.

Our approach is inspired by the classical accretion model of Bondi (1952), in which one prescribes the mass accretion rate Ṁ. The mass accretion rate is not arbitrary, it has to agree with other accretion data: the asymptotic mass density ρ_∞, the speed of sound a_∞ and the equation of state. In the original model of Bondi the selfgravity of gas has been ignored, but an analogous construction gives a model that also includes selfgravity (Karkowski et al. 2006). As a consequence, the structure of the accreting spherical ball of fluid depends on Ṁ. It has recently been discovered that these models are stable, even if selfgravity is taken into account (Mach & Malec 2008).

The content of the remainder of this paper is following. The notation is explained and the model is formulated in Sect. 2. Section 3 deals with radiating test fluids. We explicitly show that the equations can be solved by a Born-type approximation scheme, which appears to be convergent when the mass transfer rate is sufficiently low. Section 4.1 describes an iterative numerical scheme that has been specifically invented to solve this problem. Section 4.2 presents the numerical results. The imposed conditions – thin-gas approximation, steady accretion and approximately stationary disks – appear quite stringent. We find disk configurations with the central mass of the order of 10⁷ M_⊙. Obtained disks are geometrically thick, and become even thicker for higher accretion rates, but the effect of increasing accretion rate and radiation is not very strong. This is different from features known in “Polish donuts” (Abramowicz et al. 1978, 1980; Paczyński & Wiita 1980; Paczyński & Abramowicz 1982). There is an upper limit for the luminosity, and it is significantly lower than the Eddington luminosity. These numerical solutions pass the recently formulated virial test described in Mach (2012).

2. Description of the model

2.1. Notation and equations

We assume that a steady accretion disk rotates around a central mass M_c. The fluid is barotropic, i.e., p = p(ρ), where p is the pressure and ρ is the mass density. Symbols U and Φ will be used to denote the velocity of the fluid and the gravitational potential, respectively.

We assume that there exists a phenomenological mechanism that produces radiation in the zone occupied by a radial infall flow (this coincides with the support of the mass accretion function Ṁ – see a definition in Sect. 2.2), near the z = 0 plane, but we do not specify this mechanism. Pringle remarked that a steady disk can be constructed for almost any combination of viscosity and radiation process (Pringle 1981). We believe that one can find – by trial and error – a suitable (point-dependent) viscosity near the surface z = 0 that can produce the radiation that is obtained in the way described below.

The radiation is produced with the emissivity per unit mass j_ν and undergoes the Thompson scattering. Our description of the radiative transfer follows Padmanabhan (2000). Neglecting the change of the frequency of the scattered photon, we can write the scattering cross section as

$\frac{\mathrm{d}\mathit{\sigma }\left(\hat{{k}}\mathrm{\to }{k̂}\mathrm{\prime }\right)}{\mathrm{d\Omega }}\mathrm{=}\mathit{\sigma \varphi }\left(\hat{{k}}\mathit{,}{k̂}\mathrm{\prime }\right)\mathit{,}$where σ denotes the total Thompson scattering cross section, ϕ describes the angular dependence, and

$\mathrm{\int }\mathrm{d\Omega }\mathit{\varphi }\left(\hat{{k}}\mathit{,}\hat{{k}}\mathrm{\prime }\right)\mathrm{=}\mathrm{1.}$Here $\hat{{k}}$ and $\hat{{k}}{}^{\mathrm{\prime }}$ denote appropriate directions in which the radiation propagates.

The radiation transport can be described in terms of its intensity I_ν. For a stationary process we have $\hat{{k}}\mathrm{·}\mathrm{\nabla }{\mathit{I}}_{\mathit{\nu }}\mathrm{=}\frac{\mathit{\rho }{\mathit{j}}_{\mathit{\nu }}}{\mathrm{4}\mathit{\pi }}\mathrm{-}\mathit{c\rho \kappa }{\mathit{I}}_{\mathit{\nu }}\mathrm{+}\mathit{c\rho \kappa }\mathrm{\int }\mathrm{d}{\mathrm{\Omega }}^{\mathrm{\prime }}\mathit{\varphi }\left(\hat{{k}}\mathit{,}\hat{{k}}\mathrm{\prime }\right){\mathit{I}}_{\mathit{\nu }}\left(\hat{{k}}\mathrm{\prime }\right)\mathit{.}$(1)Here κ is the scattering opacity divided by the speed of light c. For the fully ionised hydrogen κ = σ/(m_pc).

Introducing the radiative flux

${F}\mathit{\nu }\mathrm{=}\mathrm{\int }\mathrm{d\Omega }\hat{{k}}{\mathit{I}}_{\mathit{\nu }}\left(\hat{{k}}\right)$and integrating Eq. (1) over the solid angle, we obtain ∇F_ν = ρj_ν, where we assumed that j_ν is independent of the direction $\hat{{k}}$. The scattering terms that are present in Eq. (1) cancel after integration.

In the following we will use the frequency integrated radiative flux or the radiation momentum density

${j}\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{\nu }{F}\mathit{\nu }\mathit{,}$so that $\mathrm{\nabla }{j}\mathrm{=}\mathit{\rho }\mathrm{\int }\mathrm{d}\mathit{\nu }{\mathit{j}}_{\mathit{\nu }}\mathit{.}$(2)The density and the velocity of the fluid must obey the continuity equation $\mathrm{\nabla }\left(\mathit{\rho }{U}\right)\mathrm{=}\mathrm{0.}$(3)The momentum conservation – stationary Euler equations with the radiation term – is given by $\mathrm{(}{U}\mathrm{·}\mathrm{\nabla }\mathrm{)}{U}\mathrm{=}\mathrm{-}\mathrm{\nabla \Phi }\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{\nabla }\mathit{p}\mathrm{+}\mathit{\kappa }{j}\mathit{,}$(4)where Φ is the gravitational potential. The last term describes the interaction of gas and radiation in the thin-gas approximation; again, only Thomson scattering is taken into account.

The energy conservation equation states that the radiation energy flux j is correlated with the energy flux of the infalling gas: $\mathrm{\nabla }\mathrm{·}\left({U}\mathit{\rho }\left(\mathit{h}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{U}\mathrm{2}\mathrm{+}\mathrm{\Phi }\right)\right)\mathrm{+}\mathrm{\nabla }\mathrm{·}{j}\mathrm{=}\mathrm{0}\mathit{,}$(5)where h is the specific gas enthalpy (dh = dp/ρ).

Finally, the gravitational potential is given by the Poisson equation $\mathrm{\Delta \Phi }\mathrm{=}\mathrm{4}\mathit{\pi G\rho ,}$(6)where G is the gravitational constant.

It is now clear that Eqs. (3)–(6) decouple from Eq. (2) and any other possible equations describing the radiation transport.

Let L denote the total luminosity: L = ∫_SdS·j, where S is the surface of a disk and dS denotes an oriented normal surface element. Integration of Eq. (5) over the disk volume leads to the approximate expression (notice that the neglected term with the enthalpy is comparatively small on the boundary) $\mathit{L}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathit{S}}\mathrm{d}{S}\mathrm{·}{U}\mathit{\rho }\left(\frac{\mathrm{1}}{\mathrm{2}}{U}\mathrm{2}\mathrm{+}\mathrm{\Phi }\right)\mathit{.}$(7)The luminosity depends on the shape of the disk, boundary values of the gravitational potential, rotation velocity, and the mass accretion rate. All these quantities are intricately related according to the differential equations; the shape itself can be determined a posteriori, after solving all of them.

2.2. Axially symmetric equations

We will seek axially symmetric solutions of Eqs. (3)–(6). The adopted coordinates are cylindrical variables (r,φ,z). Here r is the cylindrical radius, while R denotes the polar radius, so that $\mathit{R}\mathrm{=}\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$.

The condition of stationarity essentially implies U^z = 0; we assume that from now on. The velocity vector can be written as U = U∂_r + ω∂_φ, where we shall demand |U| ≪ ωr.

Let us define the mass accretion function $\mathit{Ṁ}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi Ur\rho }\mathit{.}$(8)Equation (3) implies that Ṁ depends only on the variable z, i.e., Ṁ = Ṁ(z).

The Euler equations read $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{r}}\mathrm{\Phi }\mathrm{+}\frac{\mathrm{1}}{\mathit{\rho }}{\mathit{\partial }}_{\mathit{r}}\mathit{p}\mathrm{+}\mathit{U}{\mathit{\partial }}_{\mathit{r}}\mathit{U}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\mathit{r}& \mathrm{=}& \mathit{\kappa }{\mathit{j}}^{\mathit{r}}\mathit{,}\\ {\mathit{\partial }}_{\mathit{z}}\mathrm{\Phi }\mathrm{+}\frac{\mathrm{1}}{\mathit{\rho }}{\mathit{\partial }}_{\mathit{z}}\mathit{p}& \mathrm{=}& \mathit{\kappa }{\mathit{j}}^{\mathit{z}}\mathit{,}\\ \mathit{U}\left({\mathit{\partial }}_{\mathit{r}}\mathrm{\left(}\mathit{r\omega }\mathrm{\right)}\mathrm{+}\mathit{\omega }\right)& \mathrm{=}& \end{array}$(9)For barotropes we have ∇p/ρ = ∇h, so that ∇(h + Φ) =  −(U·∇)U + κj. The expression on the right-hand side has a vanishing curl. Writing this explicitly, one finds that the consistency condition is ∂_z(ω²r) − ∂_z(U∂_rU) =  −κ(∂_zj^r − ∂_rj^z). We search for solutions with small U, and therefore we neglect the term ∂_z(U∂_rU). If in addition the rotation velocity ω depends only on the radius, i.e., ω = ω(r), then ∂_zj^r − ∂_rj^z = 0, and we can assume that j^r = ∂_rΨ, j^z = ∂_zΨ, where Ψ(r,z) is a scalar function. It is also convenient to introduce . We shall refer to Ψ (or ) as to the radiation potential.

Let us take the full divergence of both sides of Eq. (9). After trivial rearrangements and neglecting the small term with U one obtains $\mathrm{\Delta }\mathrm{(}\mathit{h}\mathrm{+}\mathrm{\Phi }\mathrm{-}\mathrm{\Psi ̂}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}}{\mathit{\partial }}_{\mathit{r}}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(10)Introducing the centrifugal potential

${\mathrm{\Phi }}_{\mathrm{c}}\mathrm{=}\mathrm{-}{\mathrm{\int }}^{\mathit{r}}\mathrm{d}{\mathit{r}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{\prime }}{\mathit{\omega }}^{\mathrm{2}}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$and integrating Eq. (10), we obtain $\mathit{h}\mathrm{+}\mathrm{\Phi }\mathrm{+}{\mathrm{\Phi }}_{\mathrm{c}}\mathrm{-}\mathrm{\Psi ̂}\mathrm{=}\mathit{C,}$(11)where C is a constant. Equation (11) can be also obtained directly from the Euler Eq. (4). From Eq. (5) we have, again neglecting terms with U, $\mathrm{\Delta }\mathrm{\Psi ̂}\mathrm{=}\frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{2}\mathit{\pi r}}{\mathit{\partial }}_{\mathit{r}}\left(\mathit{h}\mathrm{+}\mathrm{\Phi }\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\right)\mathit{.}$(12)Combining Eqs. (11) and (12), one arrives at $\begin{array}{ccc}\mathrm{\Delta }\mathrm{\Psi ̂}& \mathrm{=}& \frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{2}\mathit{\pi r}}{\mathit{\partial }}_{\mathit{r}}\left(\mathrm{-}{\mathrm{\Phi }}_{\mathrm{c}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{\Psi ̂}\right)\\ & \mathrm{=}& \end{array}$(13)We have to keep in mind that the right-hand side should be evaluated only in the region occupied by the disk. This is a linear equation that can be solved iteratively. There exists a unique solution that is regular in the open space R³.

We observe that the whole problem reduces to Eqs. (6), (11) and (13). Equations (6) and (11) can be transformed into a nonlinear integral equation. The gravitational potential Φ can be expressed using the Green function formula $\mathrm{\Phi }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{\mathit{R}}\mathrm{+}\mathrm{4}\mathit{\pi G}{\mathrm{\int }}_{\mathit{V}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}^{\mathrm{\prime }}\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}{x}\mathrm{-}{x}\mathrm{\prime }\mathrm{)}\mathit{\rho }\mathrm{\left(}{x}\mathrm{\prime }\mathrm{\right)}\mathit{,}$(14)where  − GM_c/R is the contribution due to the central mass, V ⊂ R³ is the region occupied by the disk and $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{G}\end{array}$ denotes the Green function of the laplacian in the open space R³, i.e., $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{G}\end{array}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{\left(}\mathrm{4}\mathit{\pi }\mathrm{\right|}{x}\mathrm{\left|}\mathrm{\right)}$. Equation (14) can be combined with Eq. (11), and it is convenient to exploit the fact that the density ρ is connected with the specific enthalpy h by the assumed equation of state. For the polytropic equation of state p = Kρ^Γ the specific enthalpy is h = (KΓ/(Γ − 1))ρ^Γ − 1, and we obtain $\begin{array}{ccc}\mathit{h}& \mathrm{-}& \frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{\mathit{R}}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{c}}\mathrm{-}\mathrm{\Psi ̂}\\ & & \end{array}$(15)In summary, we need the following elements to describe the hydrostationary equilibrium of a Newtonian radiating disk: (i) the accretion mass function Ṁ(z), the rotation curve ω(r), the central gravitational potential  − GM_c/R and the equation of state p = p(ρ); (ii) the radiation potential , which can be determined from the linear Eq. (13). It is convenient to write it down using the integral equation

$\mathrm{\Psi ̂}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathit{V}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}^{\mathrm{\prime }}\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}{x}\mathrm{-}{x}\mathrm{\prime }\mathrm{)}\frac{\mathit{\kappa }\mathit{Ṁ}}{{\mathit{r}}^{\mathrm{\prime }}}{\mathit{\partial }}_{{\mathit{r}}^{\mathrm{\prime }}}\left(\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}{{\mathit{r}}^{\mathrm{\prime }}}^{\mathrm{2}}\mathrm{-}{\mathrm{\Phi }}_{\mathrm{c}}\mathrm{+}\mathrm{\Psi ̂}\right)\mathrm{;}$(iii) the distribution of the density or the specific enthalpy. This can be obtained from Eqs. (6) and (11), or from an equation similar to Eq. (15).

In all above formulae the rotation law ω = ω(r) is treated as known a priori. A popular choice for test fluid solutions is to assume a (modified) Keplerian rotation ${\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(16)where z₀ is a constant. In this case the centrifugal potential is Φ_c = GM_c/$\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}$. Other simple possibilities are the rigid rotation ω = const., the rotation law ω = v₀/r, where v₀ = const. is the linear angular velocity, and ω = j₀/r², where j₀ = const. is the specific angular momentum. We will refer to the last relation as to the j-const. rotation law.

The j-const. rotation is exceptional. In this case we have

$\mathrm{2}\mathit{r}{\mathit{\omega }}^{\mathrm{2}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{\omega }{\mathit{\partial }}_{\mathit{r}}\mathit{\omega }\mathrm{=}\mathrm{0}\mathit{,}$and the radiation equation takes the form

$\mathrm{\Delta }\mathrm{\Psi ̂}\mathrm{=}\frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{2}\mathit{\pi r}}{\mathit{\partial }}_{\mathit{r}}\mathrm{\Psi ̂}\mathit{.}$Assuming that is constant at the spatial infinity, we can conclude that everywhere. For the proof, notice that it suffices to deal with the homogeneous case where as R → ∞. Multiply the last equation by and integrate over R³. Integrating by parts and employing the fact that ∂_rṀ = 0 one arrives at

$\mathrm{-}{\mathrm{\int }}_{R^{\mathrm{3}}}{\mathrm{d}}^{\mathrm{3}}\mathit{x}{\left(\mathrm{\nabla }\mathrm{\Psi ̂}\right)}^{\mathrm{2}}\mathrm{=}\mathit{\kappa }{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{d}\mathit{z}\mathit{Ṁ}{\mathrm{\left(}\mathrm{\Psi ̂}\mathrm{\right)}}^{\mathrm{2}}\mathit{.}$The left-hand side is nonpositive, while the right-hand side is nonnegative. Therefore . That is a torus with a j-const. rotation law is not emitting any radiation in r- and z-directions. Since the j^φ component of the radiation flux vector vanishes identically for the j-const. rotation, we conclude that this rotation law is not compatible with the radiation. This is consistent with the picture emerging from the analysis of Lynden-Bell & Pringle (1974); the radiated luminosity balances the energy budget of the accreting matter and it is accompanied by the shedding of the angular momentum.

The quantity C in Eq. (11) is a free parameter. The boundary of the disk is defined as a closed two-surface on which the specific enthalpy vanishes; thus it cannot be defined a priori; it is free. The exception is the test gas approximation without radiation, where the shape of a boundary is completely dictated by the central potential and the rotation curve ω(r). For radiating disks, the shape of a boundary also depends on the luminosity, which in turn is related to the mass accretion rate. To uniquely define a disk, one needs additional information. These questions will be discussed in forthcoming sections.

2.3. Luminosity and fluxes

Two of the flux densities are given by j^r = ∂_rΨ, j^z = ∂_zΨ. The third component j^φ can be obtained from the φth Euler equation.

The formula (7) yields for a disk with the minimal and maximal radial extensions r_in and r_out, respectively, $\mathit{L}\mathrm{\equiv }{\mathrm{\int }}_{\mathit{S}}\mathrm{d}{S}\mathrm{·}{j}\mathrm{\approx }\mathrm{\int }\mathrm{d}\mathit{z}\mathit{Ṁ}\left(\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{(}\mathit{\omega }{\mathit{r}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{+}\mathrm{\Phi }\right)\begin{array}{}{\mathit{r}}_{\mathrm{out}}\\ {\mathit{r}}_{\mathrm{in}}\end{array}\mathit{,}$(17)where we employed the condition |U| ≪ ωr. This agrees with the formula derived by Lynden-Bell & Pringle (1974) for the central star that is co-rotating with the disk, if r_in ≪ r_out. The local formulae for the flux densities are different. In our case the flux density j is defined uniquely, whereas in the standard approach it is given up to the total divergence (Lynden-Bell & Pringle 1974). The accretion flow originates outside of the disk and falls onto the central body. It is concentrated close to the plane z = 0. The quantity Ṁ does not depend on r and thus the mass density cannot vanish at the boundary; that means that the actual shape of the disk is not well defined near z = 0. Nevertheless, in formula (17) we assume that the disk extends from a definite exterior cylinder (r_out) to a definite inner cylinder (r_in).

3. Test fluid solutions

Assume that the gas density is low and the gravitational potential is dominated by the central term  − GM_c/R³. It is reasonable to expect (and in fact this expectation can be proved) that there are solutions of Eq. (10) that are well approximated by solutions of the linear inhomogeneous equation $\mathrm{\Delta }\mathit{h}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}}{\mathit{\partial }}_{\mathit{r}}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(18)Consider the rotation law given by Eq. (16). One can check that $\mathit{h}\mathrm{=}\mathit{G}{\mathit{M}}_{\mathrm{c}}\left(\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}\right)$(19)solves Eq. (18), with the boundary given by two planes |z| = z₀. For small z₀ one recovers the solution of Paczyński (1978):

$\mathit{h}\mathrm{=}\mathit{G}{\mathit{M}}_{\mathrm{c}}\left(\frac{\mathrm{1}}{\mathrm{2}{\mathit{r}}^{\mathrm{3}}}\mathrm{(}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{z}}^{\mathrm{2}}\mathrm{)}\right)\mathit{.}$Notice that the boundary is not closed. We show below that radiation causes the closure of the external end of the disk, but it is the influence of selfgravity that can make a finite disk.

3.1. Radiation

Let the mass of the disk be negligible and the potential be dominated by the central term Φ =  −GM_c/R. Assume a modified Keplerian rotation curve (16). The solution of Eqs. (10) and (13) can be found perturbatively, treating the mass accretion term as a perturbation.

The zeroth-order term h₀ is given by Eq. (19), that is ${\mathit{h}}_{\mathrm{0}}\mathrm{=}\mathit{G}{\mathit{M}}_{\mathrm{c}}\left(\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}\right)\mathit{,}$(20)and the h₁ term is $\begin{array}{ccc}{\mathit{h}}_{\mathrm{1}}\mathrm{(}{x}\mathrm{)}& \mathrm{=}& \int \mathit{V}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}_{\mathrm{1}}\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}{x}\mathrm{1}\mathrm{-}{x}\mathrm{)}\frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{1}}}{\mathit{\partial }}_{{\mathit{r}}_{\mathrm{1}}}\left(\mathrm{\Phi }\mathrm{+}{\mathit{h}}_{\mathrm{0}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{r}}_{\mathrm{1}}^{\mathrm{2}}\right)\\ & \mathrm{=}& \mathrm{-}{\mathrm{\int }}_{\mathit{V}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}_{\mathrm{1}}\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}{x}\mathrm{1}\mathrm{-}{x}\mathrm{)}\frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{1}}}\mathit{G}{\mathit{M}}_{\mathrm{c}}\\ & & \end{array}$(21)Higher order terms h_k (k = 2,3,...   ) are defined by ${\mathit{h}}_{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{V}}{\mathrm{d}}^{\mathrm{3}}{\mathit{x}}_{\mathit{k}}\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}{x}\mathit{k}\mathrm{-}{x}\mathrm{)}\frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathit{k}}}{\mathit{\partial }}_{{\mathit{r}}_{\mathit{k}}}{\mathit{h}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{(}{x}{\mathit{k}}^{\mathrm{)}}\mathit{.}$(22)Here, for k = 1,2,...   , x_k = (r_k,z_k,φ_k); r_k, z_k and φ_k are the cylindrical components of x_k, and the integration volume element d³x_k equals d³x_k = dr_kdz_kdφ_kr_k. The unlabelled quantities r, z and φ are reserved herein and in what follows for the components of the vector x = (r,z,φ).

The specific enthalpy corresponding to the full solution of Eqs. (10) and (13) is h = h₀ + h₁ + δh, where δh = h₂ + h₃ + ...    It is easy to check that all individual functions h₁, h₂, h₃, ...are non-positive. Assume furthermore that the disk extends from the inner cylinder of a radius r_in to an outer cylinder of a radius r_out. Then it is a simple exercise to show that $\mathrm{\left|}{\mathit{\partial }}_{{\mathit{r}}_{\mathit{k}}}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathrm{\le }\frac{\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}}{{\mathit{r}}_{\mathrm{in}}}\mathit{,}$(23)for any k = 1,2,...   

Stationary disks can exist when the above expansion scheme is convergent. The application of standard criteria for convergence of series leads to a bound onto the mass accretion rate Ṁ and on the total luminosity. These statements can be proved quite generally. We would like to point out that the rotation law (16) is not compatible with the geometry of a closed disk. While it allows for a closure at a finite distance from the centre, it should be modified in the interior region to give a disk configuration. Let us point out, however, that the inclusion of self-gravity should allow for a closed disk configuration.

3.2. Radiation and the boundary of the disk

The boundary of a disk is a smallest surface where the specific enthalpy vanishes: h = h₀ + h₁ + δh ≡ 0. In the absence of accretion, i.e., when Ṁ = 0, we have h_k = 0 for k = 1,2,...    Thus at z = z₀ there exists a horizontal part of the boundary.

If the mass accretion term Ṁ is positive, then it is clear from the inspection of Eqs. (21) and (22) that h₁ and δh are non-negative. This means that the boundary is pushed inward to a location with a value |z| < z₀.

The directional derivative of the specific enthalpy h vanishes along the boundary: dh = 0, dh = ∂_rhdr + ∂_zhdz. Notice that $\mathit{\kappa }{\mathit{j}}_{\mathit{r}}\mathrm{=}{\mathit{\partial }}_{\mathit{r}}{}^{\mathrm{(}}\mathit{h}_{\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{h}}^{\mathrm{)}}$ and $\mathit{\kappa }{\mathit{j}}_{\mathit{z}}\mathrm{=}{\mathit{\partial }}_{\mathit{z}}{}^{\mathrm{(}}\mathit{h}_{\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{h}}^{\mathrm{)}}$. This leads to the equation $\frac{\mathrm{d}\mathit{z}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathit{\kappa }{\mathit{j}}_{\mathit{r}}\mathrm{+}{\mathit{\partial }}_{\mathit{r}}{\mathit{h}}_{\mathrm{0}}}{\mathit{\kappa }{\mathit{j}}_{\mathit{z}}\mathrm{+}{\mathit{\partial }}_{\mathit{z}}{\mathit{h}}_{\mathrm{0}}}\mathrm{=}\mathrm{-}\frac{\mathit{\kappa }{\mathit{j}}_{\mathit{r}}\mathrm{+}{\mathit{\partial }}_{\mathit{r}}{\mathit{h}}_{\mathrm{0}}}{\mathit{\kappa }{\mathit{j}}_{\mathit{z}}\mathrm{-}{\mathit{\partial }}_{\mathit{z}}\mathrm{\Phi }}\mathrm{·}$(24)Values of j_r and j_z are positive on the upper face of the disk (z > 0) and negative on the lower face (z < 0). For vanishing radiation, Eq. (24) reduces to the well-known formula $\mathrm{d}\mathit{z}\mathit{/}\mathrm{d}\mathit{r}\mathrm{=}{}^{\mathrm{(}}\mathit{\omega }^{\mathrm{2}}\mathit{r}\mathrm{-}{\mathit{\partial }}_{\mathit{r}}{\mathrm{\Phi }}^{\mathrm{)}}\mathit{/}{\mathit{\partial }}_{\mathit{z}}\mathrm{\Phi }$.

3.3. A family of analytic solutions

Assume Ṁ ≡ F(r)δ(z), where F(r) ≡ F = const. > 0 for r_in ≤ r ≤ r_out and F(r) = 0 for r < r_in and r > r_out. The delta-like distribution can be easily handled analytically, and a close inspection shows that some aspects of the disk solutions weakly depend on the specific form of the mass accretion function Ṁ. We assume that F ≪ ∫dzωr²ρ, because the radial drift should be negligible in comparison to the rotational motion. For r_in < b < r < r_out the rotation curve is assumed to be given by Eq. (16), which yields the h₀ enthalpy term according to Eq. (20). With this assumption one can obtain the cusp-like end of the external part of the disk. This rotation law does not allow to close the internal part of the disk. To achieve this, we have to modify the rotation curve in a transient region just above r_in. This transient zone is not particularly important, since it does not seriously impact the structure of the disk nor its luminosity. For instance, one can take in the interval r_in < r < b, ${\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{c}}}{{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\left(\mathrm{2}\mathrm{+}\left(\mathit{\alpha }\mathrm{-}\mathrm{1}\mathrm{+}\mathit{\alpha }{\left(\frac{{\mathit{z}}_{\mathrm{0}}}{\mathit{r}}\right)}^{\mathrm{2}}\right){\left(\frac{\mathit{r}}{\mathit{b}}\right)}^{\mathit{\alpha }}\right)\mathrm{·}$(25)The constant F describes the intensity of the radial baryonic drift, and it can be found from the condition that h(r_out) = 0. The parameter α can be specified by the demand that h(r_in) = 0. Notice that the rotation curve (25) yields the specific enthalpy

${\mathit{h}}_{\mathrm{0}}\mathrm{=}\mathit{G}{\mathit{M}}_{\mathrm{c}}\left(\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}\frac{\mathrm{2}\mathrm{-}\mathrm{(}\mathit{r}\mathit{/}\mathit{b}{\mathrm{)}}^{\mathit{\alpha }}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}\right)$in the region (r_in,b). It is continuous everywhere, albeit it is not differentiable at r = b.

Let us define the elliptic function $\begin{array}{ccc}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathit{z}}^{\mathrm{)}}& \mathrm{=}& \mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}{\mathit{\phi }}_{\mathit{k}}{\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}{x}\mathrm{-}{x}{\mathit{k}}^{\mathrm{)}}|}_{{\mathit{z}}_{\mathit{k}}\mathrm{=}\mathrm{0}}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}{\mathit{\phi }}_{\mathit{k}}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{r}}_{\mathit{k}}\mathit{r}\cos {\mathit{\phi }}_{\mathit{k}}}}\end{array}$(here k = 1,2,...   ). The term h₁ is now given by $\begin{array}{ccc}{\mathit{h}}_{\mathrm{1}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathit{\kappa FG}{\mathit{M}}_{\mathrm{c}}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathrm{1}}{\mathit{r}}_{\mathrm{1}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathrm{1}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\\ & & \mathrm{\times }\left(\frac{\mathrm{1}}{{\mathrm{(}{\mathit{r}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathrm{(}{\mathit{r}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{5}\mathit{/}\mathrm{2}}}\right)\\ & & \mathrm{+}\frac{\mathit{\kappa FG}{\mathit{M}}_{\mathrm{c}}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{\mathit{b}}\mathrm{d}{\mathit{r}}_{\mathrm{1}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathrm{1}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\\ & & \mathrm{\times }{\mathit{\partial }}_{{\mathit{r}}_{\mathrm{1}}}\left(\frac{\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\alpha }}^{\mathrm{)}}{\mathrm{(}{\frac{{\mathit{r}}_{\mathrm{1}}}{\mathit{b}}}^{\mathrm{)}}}^{\mathit{\alpha }}}{\mathrm{2}{\mathrm{(}{\mathit{r}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{+}\frac{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}\mathrm{-}{\mathrm{(}{\frac{{\mathit{r}}_{\mathrm{1}}}{\mathit{b}}}^{\mathrm{)}}}^{\mathit{\alpha }}\right)}{\mathrm{2}{\mathrm{(}{\mathit{r}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\right)\mathrm{·}\end{array}$(26)One can specify the parameters b and α of the transition so that the second integral in Eq. (26) can be neglected, at least for thin disks.

The next expansion terms h_k’s (k = 2,3,...   ) are given by

${\mathit{h}}_{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{\kappa F}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathit{k}}{\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathit{z}}^{\mathrm{)}}{\mathit{\partial }}_{{\mathit{r}}_{\mathit{k}}}{\mathit{h}}_{\mathit{k}\mathrm{-}\mathrm{1}}|}_{{\mathit{z}}_{\mathit{k}}\mathrm{=}\mathrm{0}}\mathit{.}$We can use estimate (23) to obtain the inequality $\begin{array}{}\sup \end{array}\mathrm{\left|}{\mathit{h}}_{\mathit{k}}\mathrm{\right|}\mathrm{\le }\frac{\mathit{\kappa F}}{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{in}}}\begin{array}{}\sup \end{array}\left({\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathit{k}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\right)\begin{array}{}\sup \end{array}\mathrm{\left|}{\mathit{h}}_{\mathit{k}\mathrm{-}\mathrm{1}}\mathrm{\right|}\mathit{.}$(27)The elliptic function E⁽r,r_k,z⁾ is integrable on any finite disk. The inequality (27) means that the series expansion is convergent for sufficiently low values of the product

$\frac{\mathit{\kappa F}}{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{in}}}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathit{k}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\mathit{.}$It suffices that the mass current constant F is small enough, which means – notice the conservation law (17) – that stationary disks cannot have arbitrarily high luminosity.

The mass accretion rate influences the disk geometry, as seen from the following discussion. At the boundary we have h₀ + h₁ ≈ 0; thus (keeping only the leading term of h₁) we have $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}& \mathrm{=}& \frac{\mathit{\kappa F}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathrm{1}}{\mathit{r}}_{\mathrm{1}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathrm{1}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\\ & & \end{array}$(28)for r > b and $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}\frac{\mathrm{2}\mathrm{-}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{b}}^{\mathrm{)}}}^{\mathit{\alpha }}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}& \mathrm{=}& \frac{\mathit{\kappa F}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathrm{1}}{\mathit{r}}_{\mathrm{1}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathrm{1}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\\ & & \end{array}$(29)for r ≤ b. Assuming that the outer part of the disk extends up to r_out, one obtains the value of F$\frac{\mathit{\kappa F}}{\mathrm{4}\mathit{\pi }}\mathrm{\approx }\frac{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{\mathit{r}}_{\mathrm{out}}\mathit{I}}\mathit{,}$(30)where

$\mathit{I}\mathrm{=}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}\mathit{/}{\mathit{r}}_{\mathrm{out}}}^{\mathrm{1}}\mathrm{d}\mathit{xxE}\mathrm{(}\mathrm{1}\mathit{,x,}{\mathrm{0}}^{\mathrm{)}}\left(\frac{\mathrm{1}}{{\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{{\mathit{r}}_{\mathrm{out}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{+}\frac{\mathrm{3}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{r}}_{\mathrm{out}}^{\mathrm{2}}}{{\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{{\mathit{r}}_{\mathrm{out}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{5}\mathit{/}\mathrm{2}}}\right)\mathrm{·}$One can check that $\mathit{I}\mathrm{\ge }{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathit{k}}\mathit{E}\mathrm{(}{\mathit{r}}_{\mathrm{out}}\mathit{,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathrm{0}}^{\mathrm{)}}\mathrm{\ge }\mathit{\gamma }\mathrm{(}{\mathit{r}}_{\mathrm{in}}\mathit{/}{\mathit{r}}_{\mathrm{out}}\mathrm{)}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathit{k}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathit{z}}^{\mathrm{)}}$(31)for small z₀. Here γ(r_in/r_out) is a coefficient, with values depending on the ratio r_in/r_out, ranging between 0.01 and 1 for r_in/r_out = 10^-9,...,1. It is found numerically as the largest coefficient ensuring that the inequality

${\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}\mathit{/}{\mathit{r}}_{\mathrm{out}}}^{\mathrm{1}}\mathrm{d}\mathit{xE}\mathrm{\left(}\mathrm{1}\mathit{,x,}\mathrm{0}\mathrm{\right)}\mathrm{\ge }\mathit{\gamma }\mathrm{(}{\mathit{r}}_{\mathrm{in}}\mathit{/}{\mathit{r}}_{\mathrm{out}}\mathrm{)}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}\mathit{/}{\mathit{r}}_{\mathrm{out}}}^{\mathrm{1}}\mathrm{d}\mathit{xE}\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{out}}\mathit{,x,}\mathrm{0}\mathrm{)}$is satisfied. Inequality (31) combined with Eq. (30) yields

$\frac{\mathit{\kappa F}}{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{in}}}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{in}}}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{d}{\mathit{r}}_{\mathit{k}}\mathit{E}\mathrm{(}\mathit{r,}{\mathit{r}}_{\mathit{k}}\mathit{,}{\mathit{z}}^{\mathrm{)}}\mathrm{\le }\frac{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\gamma }\mathrm{(}{\mathit{r}}_{\mathrm{in}}\mathit{/}{\mathit{r}}_{\mathrm{out}}\mathrm{)}{\mathit{r}}_{\mathrm{in}}{\mathit{r}}_{\mathrm{out}}}\mathrm{·}$Obviously, for thin disks the right-hand side is much smaller than one; this implies the convergence of our approximation scheme. The consistency condition F ≪ ∫dzωr²ρ, discussed earlier, can be checked only a posteriori, after solving all equations.

Equation (29) allows, putting r = r_in, to specify the parameter α that appears in the rotation curve in the transient zone (r_in,b). By a proper choice of b, r_in, r_out and z₀ one can always achieve α ≪ 1, which gives the total luminosity close to the value corresponding to the Keplerian rotation curve. Formula (17) now yields

$\mathit{L}\mathrm{\approx }\frac{\mathit{FG}{\mathit{M}}_{\mathrm{c}}}{\mathrm{2}}\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{in}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{out}}}\right)\mathrm{·}$Other quantities, including the gas density, can be found from the Euler equations and the radiation transport equations.

4. Heavy self-gravitating disks

4.1. Numerical methods

Solutions corresponding to heavy self-gravitating disks can be obtained numerically by solving Eqs. (6), (11) and (13). In the following we assume the polytropic equation of state p = Kρ^Γ, so that h = (KΓ/(Γ − 1))ρ^Γ − 1. The rotation law is fixed up to a multiplicative constant. Numerical examples will be given for Keplerian-type rotation ω ~ r ^− 3/2.

Assume that the disk spreads up to R = R_out at the equatorial plane, and that the maximal density of the gas is ρ_max. The quantity $\mathit{u}\mathrm{=}\mathit{G}{\mathit{R}}_{\mathrm{out}}^{\mathrm{2}}{\mathit{\rho }}_{\max }$ has the dimension of the potentials Φ, Φ_c and . It can be used to define following dimensionless quantities $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Phi }\end{array}\mathrm{=}\mathrm{\Phi }\mathit{/}\mathit{u}$, $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathrm{\Phi }}_{\mathrm{c}}\end{array}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{c}}\mathit{/}\mathit{u}$, $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Psi }\end{array}\mathrm{=}\mathrm{\Psi ̂}\mathit{/}\mathit{u}$, $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{K}\end{array}\mathrm{=}\mathit{K}{\mathit{\rho }}_{\max }^{\mathrm{\Gamma }\mathrm{-}\mathrm{1}}\mathit{/}\mathit{u}$, $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\omega }\end{array}\mathrm{=}\mathit{\omega }{\mathit{R}}_{\mathrm{out}}\mathit{/}$$\sqrt{\mathit{u}}$. In similar fashion we introduce $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathit{\rho }\mathit{/}{\mathit{\rho }}_{\max }$ and $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}\mathrm{=}{\mathit{M}}_{\mathrm{c}}\mathit{/}\mathrm{\left(}{\mathit{\rho }}_{\max }{\mathit{R}}_{\mathrm{out}}^{\mathrm{3}}\mathrm{\right)}$. Dimensionless spatial coordinates are defined as $\begin{array}{}{}_{\mathrm{˜}}\\ {x}\end{array}\mathrm{=}{x}\mathit{/}{\mathit{R}}_{\mathrm{out}}$. The relevant equations of the model can be rewritten in the form $\begin{array}{ccc}& & \begin{array}{}\mathrm{˜}\\ \mathit{K}\end{array}\mathrm{\Gamma }\mathit{/}\mathrm{(}\mathrm{\Gamma }\mathrm{-}\mathrm{1}\mathrm{)}\stackrel{}{{{\mathit{\rho }}^{\mathrm{\Gamma }\mathrm{-}\mathrm{1}}}_{\mathrm{˜}}}\mathrm{+}\begin{array}{}\mathrm{˜}\\ \mathrm{\Phi }\end{array}\mathrm{+}\begin{array}{}\mathrm{˜}\\ {\mathrm{\Phi }}_{\mathrm{c}}\end{array}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathrm{\Psi }\end{array}\mathrm{=}\begin{array}{}\mathrm{˜}\\ \mathit{C}\end{array}\mathit{,}\\ & & \begin{array}{}\mathrm{˜}\\ \mathrm{\Delta }\end{array}\begin{array}{}\mathrm{˜}\\ \mathrm{\Phi }\end{array}\mathrm{=}\mathrm{4}\mathit{\pi }\begin{array}{}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathit{,}\\ & & \begin{array}{}\mathrm{˜}\\ \mathrm{\Delta }\end{array}\begin{array}{}\mathrm{˜}\\ \mathrm{\Psi }\end{array}\mathrm{=}\begin{array}{}\mathrm{˜}\\ \mathit{S}\end{array}\mathit{,}\end{array}$\arraycolsep1.75ptwhere $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{C}\end{array}\mathrm{=}\mathit{C}\mathit{/}\mathit{u}$, $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Delta }\end{array}$ is the laplacian with respect to the rescaled coordinates $\begin{array}{}{\mathrm{˜}}_{}\\ {x}\end{array}$, and

$\begin{array}{}\mathrm{˜}\\ \mathit{S}\end{array}\mathrm{=}\frac{\mathit{\kappa }\mathit{Ṁ}}{\mathrm{2}\mathit{\pi }\begin{array}{}\mathrm{˜}\\ \mathit{r}\end{array}}\mathrm{(}{\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{r}\end{array}}\begin{array}{}\mathrm{˜}\\ \mathrm{\Psi }\end{array}\mathrm{+}\mathrm{2}\begin{array}{}\mathrm{˜}\\ \mathit{r}\end{array}\stackrel{}{{{\mathit{\omega }}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{+}\stackrel{}{{{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\begin{array}{}\mathrm{˜}\\ \mathit{\omega }\end{array}{\mathit{\partial }}_{\begin{array}{}\mathrm{˜}\\ \mathit{r}\end{array}}{\begin{array}{}\mathrm{˜}\\ \mathit{\omega }\end{array}}^{\mathrm{)}}\mathit{.}$The system of Eqs. (32) and (33) with $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Psi }\end{array}\mathrm{=}\mathrm{0}$, i.e., without radiation, corresponds to a simple model of a rotating star or a toroid. There is a long tradition in the numerical analysis of these equations (see e.g. Stoeckly 1965; Ostriker & Mark 1968; Eriguchi 1978; Eriguchi & Müller 1985). In particular, it is known that they can be solved iteratively starting with some initial guess of the density distribution. The very fruitful iterative scheme of Ostriker & Mark (1968) based on the Green function expression for the gravitational potential is known as the self-consistent field method. It has been used successfully by many authors and in many variants, including computation of general-relativistic solutions (cf. Clement 1974; Blinnikov 1975; Komatsu et al. 1989; Nishida et al. 1992). An analytic expansion scheme has been also proposed by Odrzywołek (2003). Since in our case Eq. (34) for the radiation potential is also Poisson-like, the numerical method of this paper follows the self-consitent field pattern.

Knowing the density distribution, one can compute the gravitational potential $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Phi }\end{array}$ from a rescaled version of Eq. (14), i.e., $\begin{array}{}\mathrm{˜}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {x}\end{array}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ {\mathit{M}}_{\mathrm{c}}\end{array}}{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}\mathrm{+}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathit{V}}{\mathrm{d}}^{\mathrm{3}}\mathit{x}\begin{array}{}\mathrm{˜}\\ \mathit{G}\end{array}\mathrm{(}\begin{array}{}\mathrm{˜}\\ {x}\end{array}\mathrm{-}{x}\mathrm{)}\begin{array}{}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}$(35)Equation (35) is clearly not suited for a direct numerical implementation due to the singularity at $\begin{array}{}{}_{\mathrm{˜}}\\ {x}\end{array}\mathrm{=}{x}$ present in the Green function $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{G}\end{array}\mathrm{(}\begin{array}{}{}_{\mathrm{˜}}\\ {x}\end{array}\mathrm{-}{x}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{\left(}\mathrm{4}\mathit{\pi }\mathrm{\right|}\begin{array}{}{}_{\mathrm{˜}}\\ {x}\end{array}\mathrm{-}{x}\mathrm{\left|}\mathrm{\right)}$. This difficulty can be avoided by a regularisation procedure. The standard trick is to expand the Green function in terms of spherical harmonics. Assuming axial symmetry, and spherical coordinates $\mathrm{(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}\mathrm{=}\cos \begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\theta }\end{array}\mathit{,}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\phi }\end{array}\mathrm{)}$ we can write $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ \mathrm{\Phi }\end{array}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ {\mathit{M}}_{\mathit{c}}\end{array}}{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}\mathrm{-}\mathrm{4}\mathit{\pi }\sum _{\mathit{j}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathit{P}}_{\mathrm{2}\mathit{j}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{R}{\mathit{R}}^{\mathrm{2}}{\mathit{F}}_{\mathrm{2}\mathit{j}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{,R}\mathrm{\right)}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{\mu }{\mathit{P}}_{\mathrm{2}\mathit{j}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathrm{\left(}\mathit{R,\mu }\mathrm{\right)}\mathit{,}\end{array}$(36)where

${\mathit{F}}_{\mathit{j}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{,R}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{1}}{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}{\mathrm{(}{\frac{\mathit{R}}{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}}^{\mathrm{)}}}^{\mathit{j}}\mathit{,}& \begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathrm{\ge }\mathit{R,}\\ \frac{\mathrm{1}}{\mathit{R}}{\mathrm{(}{\frac{\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}}{\mathit{R}}}^{\mathrm{)}}}^{\mathit{j}}\mathit{,}& \begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{<}\mathit{R,}\end{array}$and P_j denotes the jth Legendre polynomial. From the numerical point of view the regularisation of the discretised integral in expression (35) consists of taking a finite number of terms in the series in Eq. (36).

The equation for the radiation potential $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$ can be solved in a similar fashion. Assuming that $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{S}\end{array}$ is known, we can obtain $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$ as $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ \mathrm{\Psi }\end{array}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{,}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\sum _{\mathit{j}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathit{P}}_{\mathrm{2}\mathit{j}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{R}{\mathit{R}}^{\mathrm{2}}{\mathit{F}}_{\mathrm{2}\mathit{j}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{R}\end{array}\mathit{,R}\mathrm{\right)}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{d}\mathit{\mu }{\mathit{P}}_{\mathrm{2}\mathit{j}}\mathrm{\left(}\mathit{\mu }\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{S}\end{array}\mathrm{\left(}\mathit{R,\mu }\mathrm{\right)}\mathit{.}\end{array}$(37)Knowing $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$ is necessary to compute $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{S}\end{array}$, which in turn is needed to find $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{S}\end{array}$. We demonstrate that in the case investigated here an iterative process can be applied that eventually yields an appropriate solution both for $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{S}\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$. The expression for $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{S}\end{array}$ is valid in the interior of the disk, i.e., where $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\ne \mathrm{0}$. The integral in Eq. (37) has to be evaluated only over this region, hence $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ is implicitly assumed to be known.

Conversely, once $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Phi }\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$ are known, we can compute the density $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ directly from Eq. (32). This equation, as it is written, allows negative values of enthalpy. In this case one has to modify the physical boundary of the disk. We will only search for $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ in domains where the enthalpy given by Eq. (32) is positive; otherwise we will assume $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{0}$.

The strategy of solving Eqs. (32), (33) and (34) is as follows. We fix the value of the central mass $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{M}}_{\mathrm{c}}\end{array}$, and the ratio of the inner and outer radii of the boundary of the disk R_in/R_out. In the next step we assume temporarily that $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Psi }\end{array}\mathrm{=}\mathrm{0}$, and construct an initial density distribution in the form of a toroid with assumed ratio R_in/R_out and the maximum density $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{1}$. For this toroid, the gravitational potential can be obtained form Eq. (36), and the new density distribution can be computed from Eq. (32). This procedure is iterated until a desired convergence is reached. After establishing a hydrostatic structure of the self-gravitating disk, we “switch on” the radiation. We fix the function $\mathit{\kappa }\mathit{Ṁ}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{z}\end{array}\mathrm{\right)}$ and obtain the radiation potential coming from the already found hydrostatic structure. This can be performed by solving iteratively Eq. (37). Finally, all three Eqs. (36), (37), and (32) are solved iteratively in the same fashion. This eventually produces a configuration with the hydrostationary structure that is influenced by radiation and mass accretion.

A subtle point in this procedure is that the appropriate values of constants $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{C}\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$, as well as a multiplicative constant appearing in the rotation law, are not known a priori. One has to perform a kind of renormalisation of these constants during the iterative process described above. From the condition that the inner and outer edges of the disk are located at R_in and R_out, respectively (i.e., the density must vanish there), we can compute the values of $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{C}\end{array}$ and the multiplicative constant in the rotation law, provided that the potentials $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Phi }\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$ are known. Expressing the rotation law as $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\omega }\end{array}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{r}\end{array}\mathrm{\right)}\mathrm{=}\mathit{\omega ̅}\mathit{f}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{r}\end{array}\mathrm{\right)}$ with denoting the appropriate multiplicative constant, we can write $\begin{array}{ccc}\mathit{\omega ̅}\mathrm{2}& \mathrm{=}& \frac{\mathrm{1}}{{\mathrm{\int }}_{{\mathit{R}}_{\mathrm{in}}\mathit{/}{\mathit{R}}_{\mathrm{out}}}^{\mathrm{1}}\mathrm{d}\mathit{rr}{\mathit{f}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{(}\begin{array}{}\mathrm{˜}\\ \mathrm{\Phi }\end{array}\mathrm{(}\begin{array}{}\mathrm{˜}\\ \mathit{r}\end{array}\mathrm{=}\mathrm{1}\mathit{,}\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathrm{\Psi }\end{array}\mathrm{(}\begin{array}{}\mathrm{˜}\\ \mathit{r}\end{array}\mathrm{=}\mathrm{1}\mathit{,}\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}\mathrm{=}\mathrm{0}{\mathrm{)}}^{}\\ & & \end{array}$(38)The value of $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{C}\end{array}$ can be now obtained from Eq. (32) taken at R_in or R_out with $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{=}\mathrm{0}$. The constant $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$ is given a value that ensures that the density corresponding to the found maximal value of enthalpy h is precisely ρ_max. These computations have to be repeated each time a new density distribution is computed.

Convergence properties of our numerical scheme are dependent on the spatial resolution of the numerical grid. In this paper we follow the optimisation techniques for the computation of integrals appearing in Eqs. (36) and (37) that are described by Müller & Steinmetz (1995). They allow one to achieve spatial resolutions of the order of 5000 × 5000 on a parallel computer consisting of 64 processors, with a total computing time of several minutes only. We truncate the expansion in Legendre polynomials P_l in Eqs. (36) and (37) around l_max = 20. The convergence is mainly controlled by computing the maximum norm of the two subsequent density distributions in the iteration scheme, i.e., we compute $\mathit{ϵ}\mathrm{=}{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\stackrel{}{{{\mathit{\rho }}^{\mathit{i}\mathrm{+}\mathrm{1}}}_{\mathrm{˜}}}\mathrm{-}\stackrel{}{{{\mathit{\rho }}^{\mathit{i}}}_{\mathrm{˜}}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\max }$, where the index i numbers subsequent iterations. We continue to iterate until ϵ reaches the level of 10^-6, or even 10^-8. Convergence in $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Phi }\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$ is controlled in a similar fashion, the difference being that the maximum absolute values of these quantities are not known a priori.

4.2. Numerical examples

Fig. 1 Density $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane. The density is greyscale-coded.

Fig. 2 Part of the gravitational potential $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Phi }\end{array}$ due to the accretion disk only, i.e., $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Phi }\end{array}\mathrm{+}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}\mathit{/}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}$, obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane.

Fig. 3 Radiation potential obtained for solution (a). The quantity ${\log }_{\mathrm{10}}\mathrm{\left|}\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Psi }\end{array}\mathrm{\right|}$ is plotted in greyscale. The graph shows a section of the upper hemisphere in a meridian plane.

We discuss two sample solutions (a) and (b), obtained with the method described in the preceding section. They are computed for the polytropic exponent Γ = 5/3, and the central mass $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}\mathrm{=}\mathrm{2}$. We assume the quasi-Keplerian rotation law $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\omega }\end{array}\mathrm{=}\mathit{\omega ̅}\mathit{/}\stackrel{}{{{\mathit{r}}^{\mathrm{3}\mathit{/}\mathrm{2}}}_{\mathrm{˜}}}$. We refer to this rotation curve as “quasi-Keplerian” because $\mathit{\omega ̅}\ne \sqrt{\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}}$, although the actual difference between and $\sqrt{\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{M}}_{\mathrm{c}}\end{array}}$ appears small. Notice, for the subsequent discussion, that setting $\mathit{\omega ̅}\mathrm{=}\sqrt{\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}}$ is equivalent to $\mathit{\omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{c}}}\mathit{/}{\mathit{r}}^{\mathrm{3}\mathit{/}\mathrm{2}}$. We take κṀ distribution in the form $\mathit{\kappa }\mathit{Ṁ}\mathrm{=}\mathit{A}\exp \mathrm{(}\mathrm{-}{\left(\mathrm{50}\begin{array}{}\mathrm{˜}\\ \mathit{z}\end{array}\right)}^{\mathrm{2}}\mathrm{)}\mathit{.}$(39)Solution (a) is obtained for R_in/R_out = 10^-4 and A = 10^-4. Solution (b) is computed assuming R_in/R_out = 10^-3 and A = 10^-2. The density $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$, gravitational potential of the disk $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Phi }\end{array}\mathrm{+}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}\mathit{/}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}$, and the radiation potential $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Psi }\end{array}$, obtained for solution (a) are shown in Figs. 1–3, respectively. We do not display solution (b), because its shape is not much different from that of solution (a).

The obtained value of is almost that of a strictly Keplerian rotation law. Expressing as $\mathit{\omega ̅}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\alpha }\mathrm{)}\sqrt{\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}}$ we obtain α ≈ 4.7 × 10^-6 for solution (a) and α ≈ 4.0 × 10^-5 for solution (b). This can be intuitively understood by the careful inspection of Eq. (38); for elongated disks (i.e., for low values of R_in/R_out), the total gravitational potential is dominated by the divergent term  − GM_c/R near the centre.

Other physical parameters of the solution depend on the choice of R_out and ρ_max, as well as on the actual values of constants such as κ and G. Possible values of R_out and ρ_max are restricted by the assumptions that were made when deriving equations of the model. The assumption of the thin-gas approximation requires that the mean free path of the photon λ = m_p/(σρ) = (cκρ)^-1 should be comparable with the size of the entire configuration. The radial velocity |U| should be lower than the angular one rω. Finally, the temperature of the disk should be high enough so that the assumption that the gas is fully ionised is justified. The radial velocity U can be computed from Eq. (8) based on the density distribution and the assumed mass accretion rate Ṁ. Assuming the perfect gas approximation and a value of the mean molecular weight of the gas μ, we can also find the distribution of the gas temperature

$\mathit{T}\mathrm{=}\mathit{p\mu }{\mathit{m}}_{\mathrm{p}}\mathit{/}\mathrm{\left(}\mathit{\rho }{\mathit{k}}_{\mathrm{B}}\mathrm{\right)}\mathrm{=}\mathit{K\mu }{\mathit{m}}_{\mathrm{p}}{\mathit{\rho }}^{\mathrm{\Gamma }\mathrm{-}\mathrm{1}}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{,}$where k_B denotes the Boltzmann constant. For the fully ionised hydrogen one can take μ = 1/2. The maximal temperature in the disk can be computed as

${\mathit{T}}_{\max }\mathrm{=}\mathit{u\mu }{\mathit{m}}_{\mathrm{p}}\begin{array}{}\mathrm{˜}\\ \mathit{K}\end{array}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{.}$We take R_out = 10parsecs for both solutions, and assume ρ_max = 10^-18 g   cm^-3 and ρ_max = 10^-17 g   cm^-3 for solutions (a) and (b), respectively. This gives the central mass M_c = 2.9544 × 10⁷ M_⊙ for solution (a). The central mass for solution (b) is higher by a factor of 10. Thus, if solutions (a) and (b) were to be applied in an astrophysical context, they would serve as models of accretion disks around ultramassive galactic black holes.

Other physical quantities characterising solutions (a) and (b) are as follows: for solution (a): M_fluid = (3.18720 ± 0.00049) × 10⁶ M_⊙, L = (3.1101 ± 0.0092) × 10⁶ L_⊙, T_max = (3.73508 ± 0.00031)    ×    10⁴ K. For solution (b): M_fluid = (3.20779 ± 0.00056) × 10⁷ M_⊙, L = (2.84443 ± 0.00088) × 10⁹ L_⊙, T_max = (3.72160 ± 0.00036) × 10⁵ K. The error estimates given here were obtained by comparing solutions computed on numerical grids of different resolution (spanning the range of 2000 × 2000 to 4000 × 4000), with different numbers of the Legendre polynomials (l_max = 16 to l_max = 22), and different convergence level (ϵ = 10^-6 up to ϵ = 10^-8).

The estimate of the mean free path of a photon λ – evaluated for the maximal gas density – is of the order of 1 parsec for solution (a), which is roughly the maximal vertical height of the disk. For solution (b) this estimate is 10 times smaller, but it is clear from the shape of the disk displayed in Fig. 1 that the optical thicknes is smaller than one, i.e.,

$\mathit{\tau }\mathrm{=}\underset{\mathit{r}\mathrm{\in }\mathrm{\left[}{\mathit{R}}_{\mathrm{in}}\mathit{,}{\mathit{R}}_{\mathrm{out}}\mathrm{\right]}}{\sup }{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\mathrm{0}}}\frac{\mathit{\sigma \rho }\mathrm{\left(}\mathit{r,z}\mathrm{\right)}}{{\mathit{m}}_{\mathrm{p}}}\mathrm{d}\mathit{z}\lesssim \mathrm{1.}$Here z₀ = z₀(r) is the disk height at the cylinder of radius r (Mihalas & Weibel-Mihalas 1984).

Numerical correctness of our solutions can be tested using the following virial theorem. Define ${\mathit{E}}_{\mathrm{pot}}\mathrm{=}{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{3}}\mathit{x\rho }\mathrm{(}\mathrm{\Phi }\mathrm{+}{\mathrm{\Phi }}_{\mathrm{Kep}}\mathrm{)}\mathit{/}\mathrm{2}$, ${\mathit{E}}_{\mathrm{kin}}\mathrm{=}{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{3}}\mathit{x\rho }\mathrm{|}{U}{\mathrm{|}}^{\mathrm{2}}\mathit{/}\mathrm{2}$, ${\mathit{E}}_{\mathrm{therm}}\mathrm{=}{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{3}}\mathit{x}\mathrm{3}\mathit{p}\mathit{/}\mathrm{2}$, and $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{E}\end{array}\mathrm{=}\mathit{\kappa }{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{3}}\mathit{x\rho }{x}\mathrm{·}{j}$, where Φ_Kep =  −GM_c/R is the Keplerian potential of the point mass only (note that Φ_Kep enters the formula for E_pot twice – Φ is the total gravitational potential, i.e., the sum of Φ_Kep and the contribution from the accretion disk). It can be shown (Mach 2012) that for a stationary configuration consisting of a point mass M_c and the fluid satisfying the Euler Eq. (4) the following virial identity holds:

${\mathit{E}}_{\mathrm{pot}}\mathrm{+}\mathrm{2}{\mathit{E}}_{\mathrm{kin}}\mathrm{+}\mathrm{2}{\mathit{E}}_{\mathrm{therm}}\mathrm{+}\begin{array}{}\mathrm{˜}\\ \mathit{E}\end{array}\mathrm{=}\mathrm{0.}$It differs from a standard formulation of the virial theorem by the presence of the point mass term in E_pot, and the radiation coupling term $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{E}\end{array}$. The virial check of a numerical solution can be performed by computing ${\mathit{ϵ}}_{\mathrm{v}}\mathrm{=}\mathrm{|}{\mathit{E}}_{\mathrm{pot}}\mathrm{+}\mathrm{2}{\mathit{E}}_{\mathrm{kin}}\mathrm{+}\mathrm{2}{\mathit{E}}_{\mathrm{therm}}\mathrm{+}\begin{array}{}\mathrm{˜}\\ \mathit{E}\end{array}\mathrm{|}\mathit{/}\mathrm{\left|}{\mathit{E}}_{\mathrm{pot}}\mathrm{\right|}\mathit{.}$(40)We have obtained ϵ_v ≈ 10^-8 for all numerical solutions. All constituent terms appearing in (40) were always much greater than this value. The following quantities are roughly the same for both solutions (a) and (b): E_kin/|E_pot| ≈ 0.478, E_term/|E_pot| ≈ 2.22 × 10^-2, ⁽∫d³xρΦ/2⁾/E_pot ≈ 0.525, ⁽∫d³xρΦ_Kep/2⁾/E_pot ≈ 0.475. The ratio $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{E}\end{array}\mathit{/}\mathrm{\left|}{\mathit{E}}_{\mathrm{pot}}\mathrm{\right|}$ equals approximately 3.08 × 10^-6 for solution (a) and 2.81 × 10^-4 for solution (b). That confirms the validity of our results. Clearly, our solutions pass the virial test surprisingly well. That is not very common, but not exceptional either (see e.g. Axenov & Blinnikov 1994).

5. Summary

We have formulated a consistent model of a self-gravitating disk with steadily accreting matter. The radiation emitted by the disk interacts with the gas by Thompson scattering (thin-gas approximation). The most interesting feature of our solutions is that the accretion mass rate flux density is concentrated in the equatorial plane z = 0. The slow radial drift of gas is observed in the equatorial plane and the bulk of gas remains approximately stationary. It appears that the conservation laws of the energy and the momentum together with the assumption of approximate stationarity suffice to obtain the structure of the disk. The emissivity index of accreting matter can be deduced after solving equations of the model. The approximation of stationarity demands that the radial inflow speed of matter has to be negligible compared with the rotational velocity of the gas in the disk. The secular change in the mass of the central accreting object should also be negligible. These two assumptions have been verified post factum for the numerical solutions presented in this paper.

The mathematical description of the radiating disk reduces to a pair of elliptic partial differential equations. They can be solved iteratively in a way similar to that routinely used in the literature when finding self-gravitating equilibria of non-radiating gas (cf. Ostriker & Mark 1968; Eriguchi & Müller 1985; Nishida et al. 1992). In our case each iteration step consisted of finding new distributions of the density, the gravitational potential and the radiation potential. One of the main results of this paper is that this procedure numerically converges. Analytic results can be obtained in simplified cases. We investigated the influence of the emitted radiation onto the disk structure in the test fluid approximation. The interesting conclusion is that approximately stationary solutions do exist only when the mass accretion (and thus the luminosity) is not too large. This intuitively well-understood feature of solutions has been revealed also in our numerical analysis of heavy self-gravitating disks.

We assumed thin-gas approximation, slow radial flow of matter, and approximate stationarity. These conditions appear quite restrictive – in the parameter space of solutions there is only an island of parameters for which the above assumptions can be satisfied. They lead to solutions with essential features that agree with the conditions encountered in some AGNs. We found solutions with the central mass of the order of 10⁷ M_⊙ to 10⁸ M_⊙ and disk masses of the order of 10⁶ M_⊙ to 10⁷ M_⊙. The luminosity varies between 10⁶ L_⊙ and 10⁹ L_⊙. We assumed the quasi-Keplerian rotation law ω ~ r ^− 3/2, but the rotation appears in fact to be Keplerian in the examples described in this paper.

It is interesting that the luminosity in numerical examples has been much lower than the Eddington limit. This differs from the results of Paczyński & Abramowicz (1982) or Hashimoto et al. (1995), where supercritical luminosity has been discovered. This discrepancy has a physical explanation – we describe disks in the thin gas approximation, while the quoted works dealt with the disks in thermal equilibrium or with convection transport.

A future investigation of our model can be aimed in two directions: the study of stability of disks and the formulation of a general-relativistic version. Unlike for the standard models, the structure of our Newtonian disks is dependent on the accretion flux and radiation. We have discovered that Bondi-type, spherically symmetric solutions are stable also in the self-gravitating regime (Mach & Malec 2008). The Bondi accretion models are spherically symmetric in contrast to accreting disks, but they share with our model the property that their structure depends on the accretion, and that can have a stabilising effect also in the non-spherical case. Thus we expect that the solutions discussed in this paper are stable.

There are two extreme classes of general-relativistic radiating accretion disks. Disks characterised by the size of an inner

boundary exceeding 2GM_c/c² (quasi-Newtonian case), where M_c is the central mass, should have a structure similar to those presented in this paper. However, their stability properties can still be different owing to the influence of gravitational radiation. Disks that are inherently relativistic, with the inner boundary located close to the minimum stable surface (6GM_c/c² for the Schwarzschild solution), can be significantly brighter than their Newtonian counterparts – again, this expectation is based on our experience with Bondi-type solutions (Karkowski et al. 2006; Malec & Rembiasz 2010).

Acknowledgments

The research was carried out with the supercomputer “Deszno” purchased thanks to the financial support of the European Regional Development Fund in the framework of the Polish Innovation Economy Operational Program (contract No. POIG. 02.01.00-12-023/08).

References

All Figures

	Fig. 1 Density $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\rho }\end{array}$ obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane. The density is greyscale-coded.
In the text

Fig. 2 Part of the gravitational potential $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{\Phi }\end{array}$ due to the accretion disk only, i.e., $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Phi }\end{array}\mathrm{+}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{M}}_{\mathrm{c}}\end{array}\mathit{/}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{R}\end{array}$, obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane.

In the text

	Fig. 3 Radiation potential obtained for solution (a). The quantity ${\log }_{\mathrm{10}}\mathrm{\left|}\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{\Psi }\end{array}\mathrm{\right|}$ is plotted in greyscale. The graph shows a section of the upper hemisphere in a meridian plane.
In the text