Accretion disk dynamics
αviscosity in selfsimilar selfgravitating models
^{1}
Institut für Theoretische Physik und Astrophysik,
ChrisitanAlbrechtsUniversität zu Kiel,
Leibnizstraße 15,
24118
Kiel,
Germany
email: mmeissner@astrophysik.unikiel.de; tillense@astrophysik.unikiel.de; wjd@astrophysik.unikiel.de
^{2}
Steward Observatory, The University of Arizona,
933 N. Cherry Ave.,
Tucson, AZ
85721,
USA
Received: 30 July 2015
Accepted: 25 January 2016
Aims. We investigate the suitability of αviscosity in selfsimilar models for selfgravitating disks with a focus on active galactic nuclei (AGN) disks.
Methods. We use a selfsimilar approach to simplify the partial differential equations arising from the evolution equation, which are then solved using numerical standard procedures.
Results. We find a selfsimilar solution for the dynamical evolution of selfgravitating αdisks and derive the significant quantities. In the Keplerian part of the disk our model is consistent with standard stationary αdisk theory, and selfconsistent throughout the selfgravitating regime. Positive accretion rates throughout the disk demand a high degree of selfgravitation. Combined with the temporal decline of the accretion rate and its low amount, the model prohibits the growth of large central masses.
Conclusions. αviscosity cannot account for the evolution of the whole mass spectrum of supermassive black holes (SMBH) in AGN. However, considering the involved scales it seems suitable for modelling protoplanetary disks.
Key words: accretion, accretion disks / turbulence / hydrodynamics / methods: analytical
© ESO, 2016
1. Introduction
Over half a century ago, the rise of radio astronomy lead to the discovery of quasars. The extreme luminosity of those point like sources baffled astronomers at first. At the time, nuclear fusion, the wellknown energy source of stars, was deemed the energy source of the universe. However, it quickly became apparent that nuclear fusion could not account for the tremendous amounts of energy released. The solution was proposed by Zeldovich (1964), Salpeter (1964), and LyndenBell (1969): the release of gravitational energy from material falling from a rotating disk onto a massive central object at its center, i.e. an accretion disk. Since then, accretion disks have become the paradigm and a major topic of astrophysical research. The huge interest in these disks is fueled by the fact that they can explain phenomena ranging from the evolution of supermassive black holes (SMBH) in active galactic nuclei (AGN) to star and planet formation.
However, owing to the law of momentum conservation, a successful accretion process depends on a mechanism that directs angular momentum from the inner regions to the outer regions. The mechanism of choice is viscosity. Originally derived by Goldreich & Schubert (1967) for rotating stars^{1}, it is now generally agreed that molecular viscosity, because its corresponding viscous timescale is too long to accord for the observed quasar luminosities, is a negligible process. LyndenBell (1969) himself considered magnetic fields resulting from nonuniform rotation to be the best candidates for angular momentum transport, which Balbus & Hawley (1991) confirmed, although under the condition that charged particles and a small magnetic field acting as a seed exist in the accretion disk. However, if the velocity field in an accretion disk and the corresponding timescales are considered, extremely high Reynolds numbers are obtained (e.g. LyndenBell & Pringle 1974). In consequence, turbulence is considered to emerge. Turbulent flows and the resulting eddies offer a transport mechanism for angular momentum and mass (Weizsäcker 1948; and Lüst 1952) which is independent of the existence of magnetic fields.
Consequently, some kind of viscosity prescription ν is needed to describe accretion disks. The widely used α ansatz by Shakura & Sunyaev (1973; hereafter SS73) produces physically unreasonable results such as completely isothermal disks when applied to geometrically thin and stationary AGN disks (Duschl et al. 2000). Furthermore, Shlosman et al. (1990) found that the α ansatz is not efficient enough to fuel the observed rapid mass growth of AGN (Fan et al. 2003). Thus, Illenseer & Duschl (2015; hereafter ID15) neglected the α ansatz and used – among others – the more general and physically faithful β description (Duschl et al. 2000) when they developed their dynamical model for selfgravitating accretion disks.
However, since the viscosity prescription is left undetermined in the most general form of the selfsimilar model by ID15, we investigate the consequences of applying the α prescription in this paper. Accordingly, we derive a disk equation using αviscosity and apply the methodology of similarity solutions to arrive at a system of ordinary differential equations (ODEs) which describe the time evolution of the disk. These equations are solved using numerical standard procedures. In turn, these results are discussed, which will finally allow us to review the suitability of using the α ansatz to model AGN.
It is out of the scope of this paper to provide a detailed discussion of the numerical schemes and the concept of utilizing selfsimilarity to solve differential equations. For a short overview of the latter see Dresner (1998); for a comprehensive treatment we recommend Bluman et al. (2010).
2. Disk evolution equation
For a thin, axisymmetric disk which is in hydrostatic balance in the vertical direction, ID15 derived the following disk evolution equation (1)with angular velocity and kinematic viscosity ν(r,t). In combination with (2)which gives the local power law exponent of the rotation law, this secondorder nonlinear partial differential equation (PDE) characterizes the transport processes of angular velocity under the influence of selfgravity and viscous friction. With given boundary and initial conditions and a viscosity prescription this equation can be solved.
2.1. Viscosity prescription
Viscosity is important for the transport of angular momentum and the heating of the disk. A widely used prescription is achieved^{2} via the α parametrization of SS73^{3}, (3)which is related to the disk model via the rϕ component of the sheer stress tensor (SS73), (4)Here, Σ and Π are the respective vertically integrated quantities of density ρ and pressure p which – assuming vertical isothermy – are connected via (5)where c_{s} denotes the speed of sound; ℜ the gas constant; m_{mol} molar mass; γ the adiabatic coefficient, which can be assumed to be of the order of 1; and T the temperature^{4}. Because , the heating and cooling mechanisms of the disk will eventually be of interest.
By equating Eqs. (3)and (4), solving for ν, and using Eq. (5)with vertically integrated quantities, we obtain (6)which can now be inserted into Eq. (1): (7)Applying the definition of x (Eq. (2)) leads to the expression (8)Since the speed of sound is dependent on the radius, a relationship between c_{s} and the other parameters of the system is constructed to eliminate c_{s} from the evolution equation. This can be achieved via the energy equation. Presuming a vertically optical thin disk, the effective temperature can be assumed to be equal to the central temperature (T = T_{eff}), i.e. generated heat is immediately radiated locally and consequently the thermal timescale must be much shorter than the matter diffusion timescale (Shakura & Sunyaev 1976). Thus the viscous dissipation rate (SS73) can be equated to the effective temperature of a blackbody spectrum and – using the vertically integrated form of Eq. (5)– can be expressed in terms of Σ and the speed of sound: (9)With help of the relation (10)derived in ID15 and Eq. (5)Σ and T can be eliminated from Eq. (9), which gives (11)This result can now be used to eliminate the speed of sound from the evolution equation, which gives (12)After a short calculation using the definition of x and setting (13)the subsequent result is obtained: (14)To remove unnecessary complexities in the further derivations, dimensionless scales for length, mass, and time (, , and ) are used in the next sections. The dimensions in SI units of the basic quantities involved are as follows: (15)Now, a scaling relation for corresponding to (16)is defined which allows the dimensionless solutions to be rescaled. Except for equations containing Newton’s constant, which has to be set to unity, all expressions preserve their form^{5}. When the dimensionless problem has been solved, it is possible to return to physical quantities via the scaling transformations arising from Eq. (16). Two out of the three scales (mass, length, time) are sufficient to be able to calculate the third. Consequently, η becomes a dimensionless parameter . Its value depends on the dimensions of the actual problem via (17)Furthermore, substituting allows Eq. (14)to be rewritten, which yields the following partial differential equation: (18)Thus, is basically the viscous coupling parameter equivalent to the β used by ID15.
2.2. Selfsimilar solution
A first useful step to reduce the complexity of solving the disk evolution PDE (Eq. (14)) is to introduce new variables and rewrite the equation in those terms. In order to eliminate the roots of the radial coordinate, the new variable can now be defined (19)and plugged into Eq. (14). After a short calculation (20)is obtained. Obviously, x has to be transformed (as defined in Eq. (2)) as well: (21)
2.3. Scaling transformation
The next step in applying the similarity method is to determine the group invariants. To this end, a oneparameter scaling transformation with group parameter λ and family parameters a,b,c is used: (22)Inserting the primed variables in Eq. (21)gives (23)which allows x to remain unchanged when repeating the operation for Eq. (20): (24)Consequently, it follows immediately that transformations of the evolution equation (20)are invariant if and only if a = b and c remains a constant. After a short calculation using these parameters, it becomes apparent that with the following relations hold: (25)The set of group invariants (26)where ψ = ϖ^{2}Ω^{3} allows the original PDE (Eq. (20)) to be rewritten as (27)where f(x) collects all terms depending on x and f′ is its derivative The remaining transformation of x (Eq. (21)) yields (30)Now, a set of two firstorder ordinary differential equations can be derived by applying the group invariants. Starting with the last Eq. (30), the first step is to convert to . After two subsequent transformations, it is found that . Consequently, the second and last step is to solve for : (31)To obtain the second ODE, thus arriving at an expression for , an analogous transformation of Eq. (27)has to be performed, which results in (32)Contrary to the case ID15 faced, the ODEs are not coupled, i.e. we can solve Eq. (32)independently of Eq. (31).
2.4. Auxiliary conditions
To arrive at a welldefined problem, it is necessary to specify auxiliary conditions to solve the initialboundaryvalue problem which the PDE (20)poses (e.g. Ince 1956). The PDE describes the dynamical development of Ω(ϖ,τ), i.e. the evolution of the angular velocity field of an accretion disk. Obviously, the spatial domain of physical relevance is given by 0 <ϖ ≤ ∞. Consequently, two boundary values and one initial condition Ω_{0}(ϖ) at time τ = 0 are needed.
Furthermore, these conditions have to be consistent with the properties of the physical system which is modeled to arrive at meaningful results. The model demands a ubiquitously positive surface density Σ and enclosed mass M. In addition, ∂_{r}M ≥ 0 is required. Since radial momentum transport denoted by (33)holds for this model (ID15), the auxiliary conditions must also comply with (34)in the relevant domains of the parameters (ibid.).
2.4.1. Initial conditions
To arrive at an initial condition for Ω_{0}, it is useful to take a look at Eq. (26), which describes the group invariants. After solving for Ω one can take the limit and arrive at the following relation: (35)Since Ω_{0} should not vanish (Eq. (33)), it follows that (36)which consequently leads to the initial condition for selfsimilar solutions being a power law of radius (37)where the second expression follows from the definition of ϖ in Eq. (19).
As ID15 point out, rotation laws of the kind of Ω ∝ r^{μ} cause infinite centrifugal forces for r → ∞, if the exponent μ exceeds . Additionally, the monopole approximation used in the model breaks down when μ passes −1 and approaches . However, for up to the error still appears acceptable (ibid.). After identifying μ with , this imposes an upper limit on κ. Equations (21)and (34)with (37)yield (38)as a lower limit for the exponent of the power law. Therefore, (39)holds, which translates to a domain of κ of (40)
2.5. Boundary conditions
The dependence^{6} of ξ on ϖ and τ expressed in Eq. (26)yields the following: (41)It is easy to see that the outer boundary condition coalesces with the initial condition. This reduction of auxiliary conditions is required by demanding selfsimilar solutions in terms of the independent variable ξ (Ames 1965).
At the inner rim of the disk, three reasonable boundary conditions exist: increasing, decreasing, and constant torque. The viscous torque is given by (42)which, with the help of Eqs. (4), (6), (10), (11), (13), and (28)can be rewritten^{7} and expressed in terms of the similarity variables to yield (43)The last step before being able to calculate physical meaningful solutions which correspond to specific physical situation is to specify a value for α, κ, τ_{0}, the central mass M_{⋆} at some time τ_{0}, and torque at the inner rim at a specific time τ_{0}. Expressing M_{⋆} and in terms of the group variants yields with integration constants ζ_{1} and ζ_{2}. To arrive at these results, an analysis of the phase plane and critical point is necessary. Since this analysis is, in principle, analogous to the analysis presented in ID15, it was moved to Appendix A: phase plane and critical point.
3. Results
All plotted solutions in this and the following section were obtained with the program lsode (Hindmarsh 1983) using the BDF scheme. The value of ξ_{0} can be chosen relatively arbitrarily, as long as the corresponding errors occurring during the calculation of the initial conditions via the linearized solutions do not exceed the inherent numerical errors. Before elaborating on the results in more detail, we note that ξ is based on both time and radius. Thus, each diagram can be read in two ways, i.e. for a fixed time the diagram shows the evolution of x or y according to radius, while for a fixed radius the diagram shows the time evolution of the respective dependent variable.
Fig. 1 Solutions of x for different values of κ. 

Open with DEXTER 
Fig. 2 Solutions of y for different values of κ. 

Open with DEXTER 
3.1. Similarity solutions
Figures 1 and 2 show solutions for x and y (Eqs. (32)and (31)) obtained with , and varying values of κ.
The results depicted in Fig. 1 confirm the results from the last section and the Appendix. As predicted, x approaches −1.5 for small values of ξ for any value of κ and grows to a somewhat higher finite limit for large values of ξ depending on the value of κ. The limit of x for large ξ depends on κ (see Eq. (39)): (46)Using this result, it is possible to find an approximation of y for large radii: (47)Furthermore, the influence of the saddle point on the xaxis located at x = −1.25 is manifest in the form of a maximum for solutions with κ ≲ −0.25, corresponding to the position of the saddle at x = −1.25. For solutions with κ ≥ −0.25, the maximum becomes invisible because the limits of x exceed −1.25. Consequently, these solutions do not approach the saddle point close enough to be influenced by it.
A numerical solution for κ = −0.7 could not be computed owing to the function’s stiffness in this regime. However, as mentioned above, the solution for this case is a horizontal line at x = −1.5 corresponding to a fully evolved system with Keplerian rotation and no temporal evolution or radial dependence.
Taking a more physical perspective, the results are also satisfactory. Assuming constant time and a radial dependency of x, one arrives at the following picture: at the inner rim, self gravity can be neglected and rotation is consequently Keplerian. With growing radius the influence of selfgravity grows, leading to a flatter rotation law. Looking at the alternative picture of constant radius and time evolution the development of the disks becomes apparent. At a certain radius where rotation was initially governed by selfgravity, the rotation law shifts towards Keplerian rotation as a result of the mass transfer towards the inner parts of the disk through the accretion process. In addition, Fig. 3 shows that using the thin disk assumption and the slow accretion limit is justified because v_{ϕ} is always highly supersonic, while the ratio is ≪1.
In addition, from the discussion above it becomes clear what κ actually represents: the similarity parameter describes the mass distribution within the disk.
Figure 2 also confirms the results from the last section concerning Eq. (31). For small values of ξ all solutions are proportional to while for large values of ξ and no torque the value of κ determines the slope of the linear term in Eq. (A.4). Consequently, the value of κ determines the value of y for large ξ.
Fig. 3 Solutions of and with α = 0.01 for κ = 0.2. 

Open with DEXTER 
Fig. 4 Solutions of x for different torques T_{⋆} applied at the inner rim. 

Open with DEXTER 
3.2. Influence of the torque
All the solutions we have presented so far were obtained with no torque acting on the inner rim of the disk. In Fig. 4, the solutions were obtained using τ_{0} = 1,M_{⋆}(τ_{0}) = 1,ln ξ_{0} = −25, and the given . However, the values in Figs. 9 and 10 depend explicitly on the value of η and thus on the value of α, which was set to 0.01. The value was chosen since it lies well within the known range of α (0.1 to 0.001, King et al. 2007).
A further mandatory ingredient is to calculate a numerical value for η. According to Eq. (13), η is dependent on the molecular weight. Assuming a gas composition of 75% helium and 25% hydrogen in a rather cold environment^{8}, according to Kippenhahn et al. (2012) the mean molecular weight can be set to ≈, which yields .
Using AGN scales , , and α = 0.01, a value of 1.5 × 10^{6} is obtained for η. This value is significantly lower than the value of β = 10^{3} used by ID15 and within the parameter range of α one can only arrive at η = 3 × 10^{5}. Consequently, one has to choose different scales to obtain η within the range of β. Downsizing both scales until a satisfactory value is reached leads to the scales of protoplanetary disks, while downsizing only one quantity leads either to a very spread out disk or to a very small and massive disk. Both of the latter scenarios correspond to Keplerian rotation and hence little evolution. To conclude, as Shlosman et al. (1990) point out, it seems that the α prescription is not able to describe the evolution of AGN faithfully. However, in Sect. 4 we investigate further into the efficiency of the α prescription for AGN scales.
According to Eq. (43), the torque is increasing (positive κ), decreasing (negative κ), or constant (κ = 0). In all three cases, the transition from the Keplerian to the selfgravitating regime starts at significantly smaller radii with respect to earlier times than in the solution obtained with zero torque. A significant difference to the notorque solution is that the solutions now have a third, intermediate step which serves as a maximum for the negative κ and as a saddle for positive κ.
3.3. Comparison of viscosity prescriptions
To be able to compare our results with the results from ID15, it is necessary to calculate r from ξ because the definitions of the similarity variable differ. Since the definitions of y also differ substantially, it is only reasonable to compare the different dependencies of x on the radial coordinate r. To avoid the influence of a viscosity parameter, only solutions with zero torque are compared since they are still independent of a viscosity parameter for all viscosity prescriptions. Thus, Fig. 5 shows solutions for x with initial conditions , and the same or the corresponding values^{9} of the similarity parameter κ for four different viscosity prescriptions. The relation between the κ introduced in this paper and the κ used by ID15 (hereafter ) is
Results for βviscosity (Duschl et al. 1998, (DSB)), RZ viscosity (Richard & Zahn 1999), and LP viscosity (Lin & Pringle 1987) were taken from ID15. In general, the results are similar, i.e. all solutions are basically step functions with an abrupt transition marking the change from a Keplerian to a selfgravitating rotation regime. All in all, the α prescription seems to take a middle position between LP viscosity, to which it is most closely related, and β with respect to RZ viscosity.
Fig. 5 Solutions of x for different viscosity prescriptions. 

Open with DEXTER 
4. Discussion
Using the definition of ϖ, x, and the invariants in Eqs. (19), (21), and (26), it is possible to express the defining physical quantities of the disk such as angular velocity Ω, central mass M, and surface density Σ. Applying the solutions of x and y, the temporal and spatial development of these sizes can now be investigated. The value of Ω can be directly derived via solving the definition of y for Ω and substituting ϖ: (48)From the auxiliary conditions, Eq. (33)can be used to obtain (49)In addition, the auxiliary conditions also yield the torque at the inner rim of the disk expressed in terms of the similarity variables in Eq. (43). For the surface density, it is simply necessary to solve the relation 2πGΣ = rΩ^{2}(2x + 3) derived in ID15 for Σ and insert the group invariants, which yields^{10}(50)From these three basic quantities, the accretion rate Ṁ and radial velocity v_{r} can be derived quite easily. Since , the accretion rate is given by (51)The same rationale was used to derive an expression for the radial velocity. Again, starting with a formula derived in ID15, the expressions derived above are inserted to arrive at the following relation: (52)All the necessary ingredients to find an expression for the vertically integrated viscous dissipation rate Q_{vis} based on the similarity variables are now available and can be written as (53)It is a wellknown and important result for Keplerian disks that in the radial direction Q_{vis} is proportional to the product of central mass and accretion rate over the third power of radius and not explicitly dependent on the viscosity prescription (SS73). To investigate whether this also holds for our selfsimilar model, we have to investigate the equation close to the critical point located at x = −1.5, which is equivalent to small values of ξ and thus small radii. Applying the linearization of x(ξ) for the notorque condition in Eq. (A.4)to and transforming the equation back to r yields (54)Taking the limit for yields (55)This is a remarkable result since it confirms the validity of the model in the Keplerian regime. Furthermore, it is in concordance with the results presented in ID15.
In the same fashion one can calculate the asymptotic behaviour of the other physical quantities for small and large radii and for initial conditions with and without torque supplied at the inner rim of the disk. The results of these calculations are shown in Table 1. For solutions containing κ, an equivalent solution with is given for the sake of comparison. For Ω and M, the results are identical with the solutions ID15 obtained for DSB, RZ, and LP viscosity. The same is valid for Q_{vis} and Ṁ in the limit r → 0. The results differ for the other quantities, but generally speaking point in the same direction. However, the radial velocity at the outer rim is independent of the similarity parameter, which poses an interesting deviation.
Fig. 6 Solutions of Σ(r) for different times and torques . 

Open with DEXTER 
Power law exponents of the radial coordinate for small and large radii.
Fig. 7 Solutions of Ω(r) for different times and torques . 

Open with DEXTER 
Power law exponents of the time dependence for small radii.
The solutions displayed in Fig. 6 and the following figures are given in nondimensional units and are based on solutions obtained with τ_{0} = 1,M_{⋆}(τ_{0}) = 1,ln ξ_{0} = −25, and κ = 0.2. The value of κ corresponds to which allows the solutions to be easily compared with those of ID15. Solutions with a torque provided at the inner rim are dependent on α and were obtained with α = 0.01 and η = 1.5 × 10^{6}, congruous with an AGN disk (see Sect. 3.2).
The results depicted in Fig. 6 fit very well with the results noted in Table 1. All solutions have a kink where the transition between the Keplerian and the selfgravitating regime occurs which moves to regions further outward as time progresses, a behaviour already described by Mineshige & Umemura (1996, 1997) for their selfsimilar α disk solutions. We note that the surface density decreases through time and that the transition point wanders to larger radii. This is a sensible result because the surface density has to decrease because the disk is losing mass due to accretion. Furthermore, since mass is moving inwards, the point where the mass of the disk becomes relevant for its gravitational potential moves outwards. In addition, the temporal development is consistent with the values given in Table 2. Moreover, it becomes apparent that a torque provided at the inner rim of the disk only affects this part of the disk. Such a torque allows a generally higher surface density but leads to a steeper decline in the Keplerian regime. The effect that the torque provides at the inner rim does not affect the proportionality in the selfgravitating regime, which is consistent throughout all quantities (Figs. 6 to 9), because the approximation of x for large radii (Eq. (46)) does not contain a factor which differentiates between notorque and finite torque. After a comparison of the figure to the one given in ID15, a qualitatively identical behaviour becomes evident. Moreover, the solutions in the Keplerian regime agree with the wellknown results for standard stationary α disks for the outer region (SS73). This result is not surprising since the assumption that T_{eff} = T_{c} made in Eq. (9)is only true for the rather cold and thus optically thin outer parts of the disk.
The results concerning Ω are depicted in Fig. 7. Again, the transition point between regimes moves outward with time for the same reason as in the previous case, but contrary to the case of Σ – consistently with the results from Table 2 – Ω grows with time. The value of Ω grows at a given radius because the central mass grows through accretion. Thus, to compensate centrifugal forces, Ω has to grow. In addition, the results are also qualitatively and quantitatively consistent with those of ID15. However, whether a torque is applied at the inner rim or not only seems to affect the position of the transition point.
Fig. 8 Solutions of M(r) for different times and torques . 

Open with DEXTER 
Fig. 9 Solutions of Ṁ(r) with α = 0.01 for different times and torques . 

Open with DEXTER 
The evolution of the contained mass visible in Fig. 8 is as expected from the literature (e.g. Filipov 1988). The central mass continually grows and the effect of a torque applied at the inner rim is that of an offset in the inner rim direction. Similar to the results obtained by ID15, the case of the Kelperian valued similarity parameter belongs to a solution with constant disk mass. Although we cannot acquire numerical solutions for this case, the behaviour and dependencies are well known (e.g. LyndenBell & Pringle 1974).
While at first glance the results for accretion onto the central object depicted in Fig. 9 also seem to hold no surprises, they are in fact quite different from those obtained by ID15. Depending on the similarity parameter, they report solutions with increasing, decreasing, and steady accretion rate. However, from Table 2 and the range of κ (Eq. (39)) it can be inferred that for αviscosity there are only solutions with decreasing accretion rate Ṁ, i.e. there are no quasistationary selfgravitating solutions. On the one hand, this makes it easy to check whether the results – when scaled to AGN dimensions – exceed the Eddington limit, but it is also questionable whether this model will be able to describe the full mass range (Fan et al. 2003) of SMBHs in the early universe (Duschl & Strittmatter 2011).
Fig. 10 Solutions of Ṁ(r) with α = 0.01 for different values of κ. 

Open with DEXTER 
However, our results for Ω, Σ, and M are consistent with the results of Mineshige & Umemura (1996), i.e. we find the same radial power laws (see Table 1) if we set κ = 0.2. Moreover, from Table 2 it can be inferred that the temporal development of our accretion rate is consistent with Mineshige & Umemura (1997). Furthermore, we find the same power laws for Keplerian rotation for Ω and v_{r} in Table 1 as SS73 in their cold region, thus confirming that our more general model reproduces the results for standard stationary α disks.
The similarity parameter κ has an influence on both, the temporal development of the accretion rate (Table 2) and, as one infers from Fig. 10, its total amount. Although, as Fig. 10 indicates and Table 1 predicts, the degree of selfgravity only affects the accretion process in the outer regions of the disk. Large values of κ, i.e. corresponding to selfgravity, yield the highest accretion rates and the slowest temporal decrease. Thus, the result of ID15 that “objects embedded in selfgravitating disks with flatter rotation laws grow faster than those embedded in nearly Keplerian disks” is confirmed. Moreover, our behaviour of Ṁ is qualitatively identical with the one reported for βviscosity (ID15); depending on the value of the similarity parameter, Ṁ is either positive throughout the disk or – at a certain radius – becomes negative and converges towards zero, i.e. the selfgravitating part of the disk does not accrete material at all. For the α ansatz, the value of κ separating the two regimes is κ = 0 respectively , while for the β ansatz it is . Thus, α disks require an even more selfgravitating scenario than the β ansatz to describe disks which accrete material at all radii.
Application to AGN disks
Finally, we return to physical dimensions to investigate whether the model complies with the scales known from observations. Of the three scales involved in Eq. (16), two can now be specified and the third calculated. With these three basic scales all the other scales (angular velocity, surface density, etc.) can then be calculated. For an AGN example we choose , , and α = 0.01. The consequently arising scales are listed in Table 3.
Scales for an AGN example with α = 0.01.
These scales can now be applied to the figures presented in this section. Applying^{11} the scales to Figs. 9 and 10 shows that the accretion rate lies at 4.5 × 10^{6} × 2 × 10^{5}M_{⊙} yr^{1} = 0.9 M_{⊙} yr^{1}. For black holes ≲10^{7.7}M_{⊙} this value far exceeds the Eddington limit (Eddington 1921) (assuming 10% accretion effiency). However, for the late accretion phase of a relatively heavy black hole, this is a fitting result.
Nevertheless, owing to the temporal decline of the accretion rate (Table 2) this is not sufficient to form SMBHs exceeding 10^{9}M_{⊙} in less than 10^{9} yr (Duschl & Strittmatter 2011) – even for the most selfgravitating systems – as can be seen from the following calculation: Considering that the latest observations (Wu et al. 2015) suggest even more rapid growth of the mass of the central black hole, a lack of roughly an order of magnitude of mass is discouraging.
5. Conclusions
In addition to demonstrating the validity of the model in the Keplerian and in the selfgravitating case, the discussion in the previous section quite impressively points out the problem of the α prescription: it is too inefficient to faithfully describe the evolution of AGN containing SMBHs, a result in agreement with Shlosman et al. (1990), Mineshige & Umemura (1997), and Duschl et al. (2000). Since we utilize a dynamic model, this result is not an artefact of the constraint that the accretion rate is constant. Furthermore, in contrast to Duschl et al. (2000), who have shown that αviscosity produces physically nonsensical results, we developed a physically consistent model which fails to produce the observed central masses within the time those observations suggest (Fan et al. 2003; Wu et al. 2015). Hence, our results weigh even more against the application of the α prescription to describe AGN.
However, from a simple review of scales for which the α ansatz might produce sensible results, we find that those scales correspond to protoplanetary disks. As Laughlin & Bodenheimer (1994) mention, αviscosity has already quite successfully been applied to model such disks, although mostly in
steady models. Here, because it is dynamic, our model poses an interesting alternative for future research. It is a possible vantage point for future research and development since the admittedly rather simple treatment of the thermodynamic structure of the disk would have to be modified. Furthermore, the Eddington limit could be incorporated directly into the model.
Finally, all this yields one further conclusion. The viscosity prescription does not describe the actual viscous process, but parameterizes it, ultimately based on temperature. Obviously, our understanding of the physical processes governing the temperature within an AGN disk are too poor to give the correct relation. However, since that relation seems to be a better assumption for protoplanetary disks than for AGN disks, this hints at hitherto unknown or unconsidered physics in AGN.
See e.g. Kato et al. (2008) for an application to accretion disks.
For a modern treatment of the standard α theory, see e.g. Kato et al. (2008) or Frank et al. (2002).
Owing to the definition of ϖ in Eq. (19), the behaviour of ϖ corresponds to the behaviour of r.
Through the definition of our scales in Eq. ((16)), we also account for the numerical value of the gravitational constant G.
References
 Ames, W. F. 1965, Nonlinear Differential Equations in Engineering (New York: Academic Press) [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Bluman, G. W., Cheviakov, A. F., & Anco, S. C. 2010, Applications of symmetry methods to partial differential equations (Springer) [Google Scholar]
 Dresner, L. 1998, Applications of Lie’s theory of ordinary and partial differential equations (CRC Press) [Google Scholar]
 Duschl, W. J., & Strittmatter, P. A. 2011, MNRAS, 413, 1495 [NASA ADS] [CrossRef] [Google Scholar]
 Duschl, W. J., Strittmatter, P. A., & Biermann, P. L. 1998, BAAS, 30, 917 [NASA ADS] [Google Scholar]
 Duschl, W. J., Strittmatter, P. A., & Biermann, P. L. 2000, A&A, 357, 1123 [NASA ADS] [Google Scholar]
 Eddington, A. S. 1921, Zeitschrift für Physik A Hadrons and Nuclei, 7, 351 [Google Scholar]
 Fan, X., Strauss, M. A., Schneider, D. P., et al. 2003, AJ, 125, 1649 [NASA ADS] [CrossRef] [Google Scholar]
 Frank, J., King, A., & Raine, D. 2002, Accretion power in astrophysics (Cambridge University Press) [Google Scholar]
 Frommer, M. 1928, Die Integralkurven einer gewöhnlichen Differentialgleichung erster Ordnung in der Umgebung rationaler Unbestimmtheitsstellen (Springer) [Google Scholar]
 Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [NASA ADS] [CrossRef] [Google Scholar]
 Hindmarsh, A. C. 1983, IMACS transactions on scientific computation, 1, 55 [Google Scholar]
 Illenseer, T. F., & Duschl, W. J. 2015, MNRAS, 450, 691 [NASA ADS] [CrossRef] [Google Scholar]
 Ince, E. L. 1956, Ordinary differential equations (Dover, New York) [Google Scholar]
 Kato, S., Fukue, J., & Mineshige, S. 2008, BlackHole Accretion Disks – Towards a New Paradigm (Kyoto University Press), 1 [Google Scholar]
 King, A., Pringle, J., & Livio, M. 2007, MNRAS, 376, 1740 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar structure and evolution (Springer) [Google Scholar]
 Filipov, L., Shakura, N. I., & Liubarokii, Iu 1988, Adv. Space Res., 8, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Laughlin, G., & Bodenheimer, P. 1994, ApJ, 436, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Pringle, J. E. 1987, MNRAS, 225, 607 [NASA ADS] [CrossRef] [Google Scholar]
 Lüst, R. 1952, Z. Naturforsch. Teil A, 7, 87 [NASA ADS] [Google Scholar]
 LyndenBell, D. 1969, Nature, 223, 690 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Mineshige, S., & Umemura, M. 1996, ApJ, 469, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Mineshige, S., & Umemura, M. 1997, ApJ, 480, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, D., & Zahn, J.P. 1999, A&A, 347, 734 [Google Scholar]
 Salpeter, E. E. 1964, ApJ, 140, 796 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shakura, N., & Sunyaev, R. 1976, MNRAS, 175, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Shlosman, I., Begelman, M. C., & Frank, J. 1990, Nature, 345, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Tufillaro, N. B., Abbott, T., & Reilly, J. 1992, An experimental approach to nonlinear dynamics and chaos (Redwood City, CA: AddisonWesley) [Google Scholar]
 Weizsäcker, C. F. 1948, Z. Naturforsch. A, 3, 524 [NASA ADS] [Google Scholar]
 Wu, X.B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Zeldovich, Y. B. 1964, Doklady Akademii Nauk SSSR, 155, 67 [Google Scholar]
Appendix A: Phase plane and critical point
After having determined the auxiliary conditions in Sect. 2.4 and taking into consideration that Eq. (32)is independent of y, it is now time to have a first look at the phase plane of the equation in Fig. A.1 to find out which points need further inquiry. Time and radius are always positive. Thus, according to our definition of ξ (see Eq. (26)), ξ is positive at any time. Since x is defined as the local power law exponent of Ω (Eq. (2)), its range is limited to negative values and must be ≥−1.5 (the limit of Keplerian rotation) and ≤−0.75 (where the model breaks down). In general, the particulars of the phase plane hinge on the value of κ. However, the general features remain the same and thus the value of κ is set to 0.2, which is well within the parameter range and is equivalent to the value of κ = −1 used by ID15 – thereby often allowing convenient comparison.
The two dashed black vertical lines located at x = −1.5 and x = −0.3 indicate the values of x at which the function f (Eq. (28)) becomes zero. Furthermore, the line at x = −1.5 shows the solution for κ = −0.7, which describes an already completely developed and thus completely Keplerian system. This is easily explained if one remembers that from the definition of ξ in Eq. (26)and the discussion of the auxiliary conditions one knows that a small ξ corresponds to a small radius or late times. After the passing of sufficient viscous time scales, everything will have been accreted and consequently, lacking any selfgravity, a Keplerian rotation law is valid throughout the whole disk. In a similar fashion, the rotation law has to become Keplerian for small radii since the gravitational potential there is equivalent to a point mass potential.
Furthermore, the lines where the nominator respectively denominator of become zero and three critical points located on the xaxis at x = −1.5,x = −1.25, and x = 0 are visible. Owing to the constraints of the model, only the critical points located at x = −1.5,x = −1.25 are potentially relevant for further analysis.
In addition, Fig. A.1 shows four numerically obtained (ode, Tufillaro et al. 1992, RungeKutta scheme) solutions with varying initial condition and the line which probably defines the separatrix. Since a small ξ corresponds to a small radius, physically sensible solutions must approach x = −1.5 for small ξ and, depending on the value of κ, a somewhat larger value of x for large ξ since the rotation law should reflect the growing selfgravity at the outer parts of the disk. From the four solutions depicted in Fig. A.1 only one solution appears to meet these requirements, i.e. the solution which seemingly enters the critical point located at (− 1.5  0). In order to be able to deliberately generate physically meaningful solutions, one has to be able to determine precise initial conditions. To find these conditions, the next step is to analyse the precise nature of the critical point.
By the same token a more sophisticated analysis of the system around the critical point located at x = −1.25 can be neglected because a physically sensible solution will not be in its vicinity. Classifying it as a saddle point (e.g. Ince 1956) is sufficient.
Fig. A.1 Phaseplane of for κ = 0.2 

Open with DEXTER 
Critical point
Although the phase plane presented above already shows satisfying results, one still has to determine the behaviour of the system in the vicinity of the critical point located at (− 1.5  0) to be able to generate solutions corresponding to specific physical situations. Via standard methodology (e.g. Ince 1956), the critical point can easily be classified qualitatively as an unstable node. Solutions which enter the critical point approach it from an area between the xaxis and the yet to be determined separatrix, i.e. the two characteristic directions Frommer (1928). However, further quantitative results are necessary to arrive at sensible initial conditions. In order to achieve these conditions, one can linearize the problem in the vicinity of the critical point.
To begin with, let us consider Eq. (32). Owing to its highly nonlinear dependence on variable x, which cannot be approximated, the equation is transformed to (A.1)and x is now considered a function of f. Now, one can use a Taylor series to approximate f^{3}(x) at x = −1.5, solve for x, and utilize f = 0 at the critical point to arrive at (A.2)The consequent simplification by neglecting terms of higher order in Eq. (A.1)leads to a formula which can actually be integrated analytically and solved for f which gives (A.3)where ζ_{1} denotes the integration constant yet to be determined. To arrive at an equation including x, one simply has to use Eq. (A.2)to convert f back to x: (A.4)This solution shows that is the dominant term close to the critical point where ξ → 0 as long as ζ_{1} ≠ 0. For ζ_{1} = 0, only the linear term remains and denotes the separatrix. All solutions for ζ_{1} ≠ 0 will eventually converge to the solution of the linearized equation.
To solve Eq. (31)in the vicinity of the critical point, one can drop the ξ dependent term in Eq. (A.4), and insert the result (x = −1.5) into the aforementioned equation, which can now be easily integrated to obtain (A.5)What remains to be done is to define the integration constants ζ_{1} and ζ_{2}.
For ζ_{2}, one can use that Eq. (33)gives a relation to the central mass for a certain point in time if r → 0: (A.6)If we express this with the help of the group invariants in Eq. (26)and the definition of y in Eq. (A.5), we obtain (A.7)Thus, if one specifies κ and the central mass M_{⋆τ0} at some time τ_{0} one can calculate ζ_{2}.
To compute the second integration constant ζ_{1}, the viscous torque at the inner boundary of the disk will be of interest. There are two sensible possibilities for any time 0 <t< ∞, i.e. vanishing and finite torque. To analyze the limit at the inner rim, it is helpful to express y in Eq. (43)via x. In consequence, one has to eliminate y and subsequently ξ. The first step is easily done using Eq. (A.5): (A.8)For the second step, one has to solve Eq. (A.4)for ξ with ξ ≪ 1, which gives two solutions; one for ζ_{1} = 0 and one for ζ_{1} ≠ 0: From Eq. (A.8), one can infer that a vanishing torque at the inner boundary demands that (A.11)If one plugs in the solution from (A.10)and f as defined in Eq. (28), one can see that this is indeed the case: (A.12)Hence, one can conclude that ζ_{1} = 0 denotes a special notorque solution and defines the separatrix.
The behaviour of ξ in the proximity of the critical point derived in this section can indeed be seen in Fig. A.2, where six solutions of Eq. (32)and the separatrix are depicted. Three solutions were obtained using ode (Tufillaro et al. 1992) with initial conditions above the separatrix and three solutions with initial conditions below the separatrix. The slope of the solutions below the separatrix is in concordance with the slope predicted in Eq. (A.9)and the slope of the separatrix corresponds with the one given by Eq. (A.10), thereby confirming the results acquired with the linearized version of Eq. (32).
Fig. A.2 Logarithmic phaseplane of for κ = 0.2. 

Open with DEXTER 
Additionally, one can conclude that exactly one solution, i.e. the notorque solution, enters the critical point along the singular characteristic direction defined by the separatrix and that infinitely many solutions enter the critical point along the other multiple characteristic direction, i.e. the xaxis (Frommer 1928). The initial conditions to obtain those solutions which enter the critical point along the xaxis are determined by the torque acting on the inner rim of the disk (Eq. (A.9)). To obtain these conditions, one has to investigate the limit of for ξ → 0 and ζ_{1} ≠ 0: (A.13)Thus, using Eq. (A.7), ζ_{1} can be calculated if one prescribes a torque at the inner rim at a specific time τ_{0} via (A.14)
All Tables
All Figures
Fig. 1 Solutions of x for different values of κ. 

Open with DEXTER  
In the text 
Fig. 2 Solutions of y for different values of κ. 

Open with DEXTER  
In the text 
Fig. 3 Solutions of and with α = 0.01 for κ = 0.2. 

Open with DEXTER  
In the text 
Fig. 4 Solutions of x for different torques T_{⋆} applied at the inner rim. 

Open with DEXTER  
In the text 
Fig. 5 Solutions of x for different viscosity prescriptions. 

Open with DEXTER  
In the text 
Fig. 6 Solutions of Σ(r) for different times and torques . 

Open with DEXTER  
In the text 
Fig. 7 Solutions of Ω(r) for different times and torques . 

Open with DEXTER  
In the text 
Fig. 8 Solutions of M(r) for different times and torques . 

Open with DEXTER  
In the text 
Fig. 9 Solutions of Ṁ(r) with α = 0.01 for different times and torques . 

Open with DEXTER  
In the text 
Fig. 10 Solutions of Ṁ(r) with α = 0.01 for different values of κ. 

Open with DEXTER  
In the text 
Fig. A.1 Phaseplane of for κ = 0.2 

Open with DEXTER  
In the text 
Fig. A.2 Logarithmic phaseplane of for κ = 0.2. 

Open with DEXTER  
In the text 