Accretion and structure of radiating disks
M. Smoluchowski Institute of Physics, Jagiellonian University, Reymonta 4, 30059 Kraków, Poland
email: mach@th.if.uj.edu.pl
Received: 14 September 2010
Accepted: 19 March 2012
We studied a steadily accreting, geometrically thick disk model that selfconsistently takes into account selfgravitation of the polytropic gas, its interaction with the radiation and the mass accretion rate. The accreting mass is injected inward in the vicinity of the central z = 0 plane, where also radiation is assumed to be created. The rest of the disk remains approximately stationary. Only conservation laws are employed and the gasradiation interaction in the bulk of the disk is described in the thingas approximation. We demonstrate that this scheme is numerically viable and yields a structure of the bulk that is influenced by the radiation and (indirectly) by the prescribed mass accretion rate. The obtained disk configurations are typical for environments in active galactic nuclei (AGNs), with the central mass of the order of 10^{7} M_{⊙} to 10^{8} M_{⊙}, quasiKeplerian rotation curves, disk masses ranging from about 10^{6} M_{⊙} to 10^{7} M_{⊙}, and the luminosity ranging from 10^{6} L_{⊙} to 10^{9} L_{⊙}. These luminosities are much lower than the corresponding Eddington limit.
Key words: accretion, accretion disks / hydrodynamics / gravitation / radiative transfer
© 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; LyndenBell & 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 thickdisk equations, but assumed elements of thindisk approximation in determining the vertical structure. Paczyński & Wiita (1980) imposed some ingredients of thindisk 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 selfgravitating 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 thingas 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 selfcontained 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 Borntype 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 – thingas approximation, steady accretion and approximately stationary disks – appear quite stringent. We find disk configurations with the central mass of the order of 10^{7} 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 (pointdependent) 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
where σ denotes the total Thompson scattering cross section, ϕ describes the angular dependence, and
Here and 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 (1)Here κ is the scattering opacity divided by the speed of light c. For the fully ionised hydrogen κ = σ/(m_{p}c).
Introducing the radiative flux
and integrating Eq. (1) over the solid angle, we obtain ∇F_{ν} = ρj_{ν}, where we assumed that j_{ν} is independent of the direction . 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
so that (2)The density and the velocity of the fluid must obey the continuity equation (3)The momentum conservation – stationary Euler equations with the radiation term – is given by (4)where Φ is the gravitational potential. The last term describes the interaction of gas and radiation in the thingas 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: (5)where h is the specific gas enthalpy (dh = dp/ρ).
Finally, the gravitational potential is given by the Poisson equation (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 = ∫_{S}dS·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) (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 .
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 (8)Equation (3) implies that Ṁ depends only on the variable z, i.e., Ṁ = Ṁ(z).
The Euler equations read (9)For barotropes we have ∇p/ρ = ∇h, so that ∇(h + Φ) = −(U·∇)U + κj. The expression on the righthand side has a vanishing curl. Writing this explicitly, one finds that the consistency condition is ∂_{z}(ω^{2}r) − ∂_{z}(U∂_{r}U) = −κ(∂_{z}j^{r} − ∂_{r}j^{z}). We search for solutions with small U, and therefore we neglect the term ∂_{z}(U∂_{r}U). If in addition the rotation velocity ω depends only on the radius, i.e., ω = ω(r), then ∂_{z}j^{r} − ∂_{r}j^{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 (10)Introducing the centrifugal potential
and integrating Eq. (10), we obtain (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, (12)Combining Eqs. (11) and (12), one arrives at (13)We have to keep in mind that the righthand 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^{3}.
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 (14)where − GM_{c}/R is the contribution due to the central mass, V ⊂ R^{3} is the region occupied by the disk and denotes the Green function of the laplacian in the open space R^{3}, i.e., . 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 (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
(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 (16)where z_{0} is a constant. In this case the centrifugal potential is Φ_{c} = GM_{c}/. Other simple possibilities are the rigid rotation ω = const., the rotation law ω = v_{0}/r, where v_{0} = const. is the linear angular velocity, and ω = j_{0}/r^{2}, where j_{0} = const. is the specific angular momentum. We will refer to the last relation as to the jconst. rotation law.
The jconst. rotation is exceptional. In this case we have
and the radiation equation takes the form
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^{3}. Integrating by parts and employing the fact that ∂_{r}Ṁ = 0 one arrives at
The lefthand side is nonpositive, while the righthand side is nonnegative. Therefore . That is a torus with a jconst. rotation law is not emitting any radiation in r and zdirections. Since the j^{φ} component of the radiation flux vector vanishes identically for the jconst. rotation, we conclude that this rotation law is not compatible with the radiation. This is consistent with the picture emerging from the analysis of LyndenBell & 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 twosurface 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, (17)where we employed the condition U ≪ ωr. This agrees with the formula derived by LyndenBell & Pringle (1974) for the central star that is corotating 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 (LyndenBell & 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^{3}. 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 (18)Consider the rotation law given by Eq. (16). One can check that (19)solves Eq. (18), with the boundary given by two planes z = z_{0}. For small z_{0} one recovers the solution of Paczyński (1978):
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 zerothorder term h_{0} is given by Eq. (19), that is (20)and the h_{1} term is (21)Higher order terms h_{k} (k = 2,3,... ) are defined by (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^{3}x_{k} equals d^{3}x_{k} = dr_{k}dz_{k}dφ_{k}r_{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_{0} + h_{1} + δh, where δh = h_{2} + h_{3} + ... It is easy to check that all individual functions h_{1}, h_{2}, h_{3}, ...are nonpositive. 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 (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 selfgravity 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_{0} + h_{1} + δh ≡ 0. In the absence of accretion, i.e., when Ṁ = 0, we have h_{k} = 0 for k = 1,2,... Thus at z = z_{0} 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_{1} and δh are nonnegative. This means that the boundary is pushed inward to a location with a value z < z_{0}.
The directional derivative of the specific enthalpy h vanishes along the boundary: dh = 0, dh = ∂_{r}hdr + ∂_{z}hdz. Notice that and . This leads to the equation (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 wellknown formula .
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 deltalike 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^{2}ρ, 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_{0} enthalpy term according to Eq. (20). With this assumption one can obtain the cusplike 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, (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
in the region (r_{in},b). It is continuous everywhere, albeit it is not differentiable at r = b.
Let us define the elliptic function (here k = 1,2,... ). The term h_{1} is now given by (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
We can use estimate (23) to obtain the inequality (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
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_{0} + h_{1} ≈ 0; thus (keeping only the leading term of h_{1}) we have (28)for r > b and (29)for r ≤ b. Assuming that the outer part of the disk extends up to r_{out}, one obtains the value of F(30)where
One can check that (31)for small z_{0}. 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
is satisfied. Inequality (31) combined with Eq. (30) yields
Obviously, for thin disks the righthand side is much smaller than one; this implies the convergence of our approximation scheme. The consistency condition F ≪ ∫dzωr^{2}ρ, 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_{0} one can always achieve α ≪ 1, which gives the total luminosity close to the value corresponding to the Keplerian rotation curve. Formula (17) now yields
Other quantities, including the gas density, can be found from the Euler equations and the radiation transport equations.
4. Heavy selfgravitating disks
4.1. Numerical methods
Solutions corresponding to heavy selfgravitating 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 Kepleriantype 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 has the dimension of the potentials Φ, Φ_{c} and . It can be used to define following dimensionless quantities , , , , . In similar fashion we introduce and . Dimensionless spatial coordinates are defined as . The relevant equations of the model can be rewritten in the form \arraycolsep1.75ptwhere , is the laplacian with respect to the rescaled coordinates , and
The system of Eqs. (32) and (33) with , 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 selfconsistent field method. It has been used successfully by many authors and in many variants, including computation of generalrelativistic 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 Poissonlike, the numerical method of this paper follows the selfconsitent field pattern.
Knowing the density distribution, one can compute the gravitational potential from a rescaled version of Eq. (14), i.e., (35)Equation (35) is clearly not suited for a direct numerical implementation due to the singularity at present in the Green function . 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 we can write (36)where
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 can be solved in a similar fashion. Assuming that is known, we can obtain as (37)Knowing is necessary to compute , which in turn is needed to find . We demonstrate that in the case investigated here an iterative process can be applied that eventually yields an appropriate solution both for and . The expression for is valid in the interior of the disk, i.e., where . The integral in Eq. (37) has to be evaluated only over this region, hence is implicitly assumed to be known.
Conversely, once and are known, we can compute the density 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 in domains where the enthalpy given by Eq. (32) is positive; otherwise we will assume .
The strategy of solving Eqs. (32), (33) and (34) is as follows. We fix the value of the central mass , 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 , and construct an initial density distribution in the form of a toroid with assumed ratio R_{in}/R_{out} and the maximum density . 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 selfgravitating disk, we “switch on” the radiation. We fix the function 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 and , 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 and the multiplicative constant in the rotation law, provided that the potentials and are known. Expressing the rotation law as with denoting the appropriate multiplicative constant, we can write (38)The value of can be now obtained from Eq. (32) taken at R_{in} or R_{out} with . The constant 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 , where the index i numbers subsequent iterations. We continue to iterate until ϵ reaches the level of 10^{6}, or even 10^{8}. Convergence in and 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 obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane. The density is greyscalecoded. 

Open with DEXTER 
Fig. 2 Part of the gravitational potential due to the accretion disk only, i.e., , obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane. 

Open with DEXTER 
Fig. 3 Radiation potential obtained for solution (a). The quantity is plotted in greyscale. The graph shows a section of the upper hemisphere in a meridian plane. 

Open with DEXTER 
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 . We assume the quasiKeplerian rotation law . We refer to this rotation curve as “quasiKeplerian” because , although the actual difference between and appears small. Notice, for the subsequent discussion, that setting is equivalent to . We take κṀ distribution in the form (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 , gravitational potential of the disk , and the radiation potential , 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 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 thingas 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
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
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^{7} 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^{6} M_{⊙}, L = (3.1101 ± 0.0092) × 10^{6} L_{⊙}, T_{max} = (3.73508 ± 0.00031) × 10^{4} K. For solution (b): M_{fluid} = (3.20779 ± 0.00056) × 10^{7} M_{⊙}, L = (2.84443 ± 0.00088) × 10^{9} L_{⊙}, T_{max} = (3.72160 ± 0.00036) × 10^{5} 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.,
Here z_{0} = z_{0}(r) is the disk height at the cylinder of radius r (Mihalas & WeibelMihalas 1984).
Numerical correctness of our solutions can be tested using the following virial theorem. Define , , , and , 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:
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 . The virial check of a numerical solution can be performed by computing (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^{3}xρΦ/2^{)}/E_{pot} ≈ 0.525, ^{(}∫d^{3}xρΦ_{Kep}/2^{)}/E_{pot} ≈ 0.475. The ratio 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 selfgravitating disk with steadily accreting matter. The radiation emitted by the disk interacts with the gas by Thompson scattering (thingas 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 selfgravitating equilibria of nonradiating 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 wellunderstood feature of solutions has been revealed also in our numerical analysis of heavy selfgravitating disks.
We assumed thingas 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^{7} M_{⊙} to 10^{8} M_{⊙} and disk masses of the order of 10^{6} M_{⊙} to 10^{7} M_{⊙}. The luminosity varies between 10^{6} L_{⊙} and 10^{9} L_{⊙}. We assumed the quasiKeplerian 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 generalrelativistic version. Unlike for the standard models, the structure of our Newtonian disks is dependent on the accretion flux and radiation. We have discovered that Bonditype, spherically symmetric solutions are stable also in the selfgravitating 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 nonspherical case. Thus we expect that the solutions discussed in this paper are stable.
There are two extreme classes of generalrelativistic radiating accretion disks. Disks characterised by the size of an inner
boundary exceeding 2GM_{c}/c^{2} (quasiNewtonian 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^{2} for the Schwarzschild solution), can be significantly brighter than their Newtonian counterparts – again, this expectation is based on our experience with Bonditype 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.0012023/08).
References
 Abramowicz, M. A., Jaroszyński, M., & Sikora, M. 1978, A&A, 63, 221 [NASA ADS] (In the text)
 Abramowicz, M. A., Calvani, M., & Nobili, L. 1980, ApJ, 242, 772 [NASA ADS] [CrossRef] (In the text)
 Axenov, A. G., & Blinnikov, S. I. 1994, A&A, 290, 674 [NASA ADS] (In the text)
 Blinnikov, S. I. 1975, Sov. Astron., 19, 151 [NASA ADS] (In the text)
 Bondi, H. 1952, MNRAS, 112, 195 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Clement, M. J. 1974, ApJ, 194, 709 [NASA ADS] [CrossRef] (In the text)
 Eriguchi, Y. 1978, Publ. Astron. Soc. Japan, 30, 507 [NASA ADS] (In the text)
 Eriguchi, Y., & Müller, E. 1985, A&A, 146, 260 [NASA ADS] (In the text)
 Hashimoto, M., Eriguchi, Y., Arai, K., & Müller, E. 1993, A&A, 268, 131 [NASA ADS] (In the text)
 Hashimoto, M., Eriguchi, Y., & Müller, E. 1995, A&A, 297, 135 [NASA ADS] (In the text)
 Karkowski, J., Kinasiewicz, B., Mach, P., Malec, E., & Swierczyński, Z. 2006, Phys. Rev. D, 73, 021503(R) [NASA ADS] [CrossRef] (In the text)
 Komatsu, H., Eriguchi, Y., & Hachisu, I. 1989, MNRAS, 237, 355 [NASA ADS] (In the text)
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] (In the text)
 Mach, P. 2012, MNRAS, 422, 772 [NASA ADS] [CrossRef] (In the text)
 Mach, P., & Malec, E. 2008, Phys. Rev. D, 78, 124016 [NASA ADS] [CrossRef] (In the text)
 Malec, E., & Rembiasz, T. 2010, Phys. Rev. D, 82, 124005 [NASA ADS] [CrossRef] (In the text)
 Mihalas, D., & WeibelMihalas, B. 1984, Foundations of Radiation Hydrodynamics (Oxford University Press) (In the text)
 Müller, E., & Steinmetz, M. 1995, Comput. Phys. Commun. 89, 45 (In the text)
 Nishida, S., Eriguchi, Y., & Lanza, A. 1992, ApJ, 401, 618 [NASA ADS] [CrossRef] (In the text)
 Odrzywołek, A. 2003, MNRAS, 345, 497 [NASA ADS] [CrossRef] (In the text)
 Ostriker, J. P., & Mark, J.WK. 1968, ApJ, 151, 1075 [NASA ADS] [CrossRef] (In the text)
 Paczyński, B. 1978, Acta Astron., 28, 91 [NASA ADS] (In the text)
 Paczyński, B., & Abramowicz, M. A. 1982, ApJ, 253, 897 [NASA ADS] [CrossRef] (In the text)
 Paczyński, B., & Wiita, P. J. 1980, A&A, 88, 23 [NASA ADS] (In the text)
 Padmanabhan, T. 2000, Theoretical Astrophysics, Vol. I (Cambridge Univ. Press) (In the text)
 Papaloizu, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 [NASA ADS] [CrossRef] (In the text)
 Papaloizu, J. C. B., & Pringle, J. E. 1985, MNRAS, 213, 799 [NASA ADS] (In the text)
 Piran, T. 1978, ApJ, 221, 652 [NASA ADS] [CrossRef] (In the text)
 Pringle, J. E. 1976a, MNRAS, 177, 65 [NASA ADS] (In the text)
 Pringle, J. E. 1976b, ARA&A, 19, 137 [NASA ADS] [CrossRef] (In the text)
 Qian, L., Abramowicz, M. A., Fragile, P. C., et al. 2009, A&A, 498, 471 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
 Shakura, N. I., & Sunyaev, R. A. 1976, MNRAS, 175, 613 [NASA ADS] [CrossRef] (In the text)
 Sikora, M. 1981, MNRAS, 196, 257 [NASA ADS] [CrossRef] (In the text)
 Stoeckly, R. 1965, ApJ, 142, 208 [NASA ADS] [CrossRef] (In the text)
All Figures
Fig. 1 Density obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane. The density is greyscalecoded. 

Open with DEXTER  
In the text 
Fig. 2 Part of the gravitational potential due to the accretion disk only, i.e., , obtained for solution (a). The plot shows a section of the upper hemisphere in a meridian plane. 

Open with DEXTER  
In the text 
Fig. 3 Radiation potential obtained for solution (a). The quantity is plotted in greyscale. The graph shows a section of the upper hemisphere in a meridian plane. 

Open with DEXTER  
In the text 