EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 588 (April 2016) A&A, 588 (2016) A22 Full HTML
Free Access
Issue
A&A
Volume 588, April 2016
Article Number A22
Number of page(s) 11
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201527092
Published online 11 March 2016

© 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 well-known 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 Lynden-Bell (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 super-massive 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 stars1, 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. Lynden-Bell (1969) himself considered magnetic fields resulting from non-uniform 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. Lynden-Bell & 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 self-gravitating accretion disks.

However, since the viscosity prescription is left undetermined in the most general form of the self-similar 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 self-similarity 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 −r4Ω2∂Ω∂t=∂∂r(νr3Ω3x(2x+3))\begin{eqnarray} \label{original} -r^4\Omega^2\frac{\partial\Omega}{\partial t}=\frac{\partial}{\partial r}\left (\nu r^3 \Omega^3 x(2x+3) \right ) \end{eqnarray}(1)with angular velocity Ω=vϕr\hbox{$\Omega=\frac{v_{\varphi}}{r}$} and kinematic viscosity ν(r,t). In combination with x=rΩ∂Ω∂r=∂lnΩ∂lnr,\begin{eqnarray} \label{x} x=\frac{r}{\Omega} \frac{\partial \Omega}{\partial r} =\frac{\partial \ln \Omega}{\partial \ln r}, \end{eqnarray}(2)which gives the local power law exponent of the rotation law, this second-order non-linear partial differential equation (PDE) characterizes the transport processes of angular velocity under the influence of self-gravity 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 achieved2 via the α parametrization of SS733, Trϕ=\begin{eqnarray} T_{r \varphi}& =& -\alpha \Pi \label{alt_alpha} , \end{eqnarray}(3)which is related to the disk model via the rϕ component of the sheer stress tensor (SS73), Trϕ=\begin{eqnarray} T_{r \varphi}&=&\nu \Sigma r \frac{{\rm d}\Omega}{{\rm d}r}\cdot \label{2} \end{eqnarray}(4)Here, Σ and Π are the respective vertically integrated quantities of density ρ and pressure p which – assuming vertical isothermy – are connected via cs=γpρ=γℜTmmol=ΠΣ,\begin{eqnarray} \label{cs_t} c_{\rm s}=\sqrt{\frac{\gamma p}{\rho}} = \sqrt{\frac{\gamma \Re T}{m_{\text{mol}}}} =\sqrt{\frac{\Pi}{\Sigma}}, \end{eqnarray}(5)where cs denotes the speed of sound; ℜ the gas constant; mmol molar mass; γ the adiabatic coefficient, which can be assumed to be of the order of 1; and T the temperature4. Because cs∝T\hbox{$ c_{\rm s}\propto\sqrt{T} $}, 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 ν=−αcs2r ∂rΩ,\begin{eqnarray} \nu=-\frac{\alpha c_{\rm s}^{2}}{ r \,\partial_r \Omega}, \label{nu} \end{eqnarray}(6)which can now be inserted into Eq. (1): −r4Ω2∂Ω∂t=−∂∂r(αcs2r ∂rΩr3Ω3x(2x+3)).\begin{eqnarray} \label{4} -r^4\Omega^2\frac{\partial \Omega}{\partial t}=-\frac{\partial}{\partial r}\left (\frac{\alpha c_{\rm s}^{2}}{ r \,\partial_r \Omega} r^3 \Omega^3 x(2x+3) \right ). \end{eqnarray}(7)Applying the definition of x (Eq. (2)) leads to the expression r4Ω2∂Ω∂t=α∂∂r(cs2r3Ω2(2x+3)).\begin{eqnarray} r^4\Omega^2 \frac{\partial \Omega}{\partial t}=\alpha \frac{\partial}{\partial r}\left (c_{\rm s}^2 r^3 \Omega^2(2x+3) \right ). \end{eqnarray}(8)Since the speed of sound is dependent on the radius, a relationship between cs and the other parameters of the system is constructed to eliminate cs 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 = Teff), 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 black-body spectrum and – using the vertically integrated form of Eq. (5)– can be expressed in terms of Σ and the speed of sound: 2σT4=\begin{eqnarray} 2\sigma T^4&=&-\alpha \Pi r \frac{\partial \Omega}{\partial r}=-\alpha c_{\rm s}^2 \Sigma r \frac{\partial \Omega}{\partial r} \cdot \label{t_eff} \end{eqnarray}(9)With help of the relation Σ=rΩ2(2x+3)2πG\begin{eqnarray} \Sigma=\frac{r \Omega^2(2x+3)}{2\pi G} \label{sig} \end{eqnarray}(10)derived in ID15 and Eq. (5)Σ and T can be eliminated from Eq. (9), which gives cs2=−(αγℜ44σπGmmol4(r2Ω2(2x+3)∂rΩ))13.\begin{eqnarray} c_{\rm s}^2=-\left( \frac{\alpha \gamma \Re^4}{4 \sigma \pi G m_{\text{mol}}^4} \left ( r^2 \Omega^2 (2x+3) \partial_r \Omega \right )\right )^\frac{1}{3}. \label{cs} \end{eqnarray}(11)This result can now be used to eliminate the speed of sound from the evolution equation, which gives r4Ω2∂Ω∂t=α∂∂r(−(αγℜ44σπGmmol4(∂rΩ))13r113Ω83(2x+3)43).\begin{eqnarray} r^4\Omega^2\frac{\partial \Omega}{\partial t}&=\alpha \frac{\partial}{\partial r}\left (-\left( \frac{\alpha \gamma \Re^4}{4 \sigma \pi G m_{\text{mol}}^4} \left ( \partial_r \Omega \right )\right )^\frac{1}{3}r^{\frac{11}{3}} \Omega^{\frac{8}{3}} (2x+3)^{\frac{4}{3}} \right ). \end{eqnarray}(12)After a short calculation using the definition of x and setting η=(α4γℜ44σπGmmol4)13,\begin{eqnarray} \eta=\left (\frac{\alpha^4 \gamma \Re^4}{4 \sigma \pi G m_{\text{mol}}^4} \right )^{\frac{1}{3}}, \label{gamma} \end{eqnarray}(13)the subsequent result is obtained: r4Ω2∂Ω∂t=−η∂∂r(r103Ω3x13(2x+3)43).\begin{eqnarray} \label{6} r^4 \Omega^2 \frac{\partial \Omega}{\partial t}=-\eta \frac{\partial}{\partial r}\left ( r^{\frac{10}{3}} \Omega^3 x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}} \right ). \end{eqnarray}(14)To remove unnecessary complexities in the further derivations, dimensionless scales for length, mass, and time (˜r\hbox{$\tilde r$}, ˜M\hbox{$\tilde M$}, and ˜t\hbox{$\tilde t$}) are used in the next sections. The dimensions in SI units of the basic quantities involved are as follows: [Ω]=1s,[t]=s,[r]=m,[η]=m53s·\begin{eqnarray} [\Omega]=\frac{1}{{s}}, \; [t]={s}, \; [r]={m},\; [\eta]=\frac{{m}^{\frac{5}{3}}}{{s}}\cdot \end{eqnarray}(15)Now, a scaling relation for ˜t\hbox{$\tilde t$} corresponding to t̂=r̂3GM̂\begin{eqnarray} \hat t = \sqrt{\frac{\hat r^3}{G \hat M} } \label{dimless} \end{eqnarray}(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 form5. 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 ˜η\hbox{$\tilde \eta$}. Its value depends on the dimensions of the actual problem via ˜η=ηt̂r̂53·\begin{eqnarray} \tilde \eta = \eta \frac{\hat t}{\hat r^{\frac{5}{3}}}\cdot \end{eqnarray}(17)Furthermore, substituting τ=5˜η˜t\hbox{$\tau=5\tilde \eta \tilde t$} allows Eq. (14)to be rewritten, which yields the following partial differential equation: 5r4˜Ω2˜∂˜Ω∂τ=−∂∂˜r(r103˜Ω3˜x13(2x+3)43).\begin{eqnarray} 5 \tilde r^4 \tilde \Omega^2 \frac{\partial \tilde \Omega}{\partial \tau }=- \frac{\partial}{\partial \tilde r}\left ( \tilde r^{\frac{10}{3}} \tilde \Omega^3 x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}} \right ). \end{eqnarray}(18)Thus, ˜η\hbox{$\tilde \eta$} is basically the viscous coupling parameter equivalent to the β used by ID15.

2.2. Self-similar 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 ϖ=r53\begin{eqnarray} \label{defdelta} \varpi &=& r^{\frac{5}{3}} \end{eqnarray}(19)and plugged into Eq. (14). After a short calculation 3ϖ2Ω2∂Ω∂τ=−∂∂ϖ(ϖ2Ω3x13(2x+3)43)\begin{eqnarray} 3 \varpi^2 \Omega^2 \frac{\partial \Omega}{\partial \tau}&=- \frac{\partial}{\partial \varpi} \left ( \varpi^2 \Omega^3 x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}} \right ) \label{transevo} \end{eqnarray}(20)is obtained. Obviously, x has to be transformed (as defined in Eq. (2)) as well: x=(dlnϖdlnr)(∂lnΩ∂lnϖ)=53∂lnΩ∂lnϖ·\begin{eqnarray} \label{transx} x=\left (\frac{{\rm d} \ln \varpi}{{\rm d} \ln r}\right ) \left ( \frac{\partial \ln \Omega}{\partial \ln \varpi} \right) = \frac{5}{3}\frac{\partial \ln \Omega}{\partial \ln \varpi}\cdot \end{eqnarray}(21)

2.3. Scaling transformation

The next step in applying the similarity method is to determine the group invariants. To this end, a one-parameter scaling transformation with group parameter λ and family parameters a,b,c is used: ϖ′=λaϖ,τ′=λbτ,Ω′=λcΩ.\begin{eqnarray} \varpi'=\lambda^a\varpi, \; \tau'=\lambda^b \tau, \; \Omega'=\lambda^c \Omega. \end{eqnarray}(22)Inserting the primed variables in Eq. (21)gives x=53∂lnΩ∂lnϖ=53ϖΩ∂Ω∂ϖ=53λaλcϖ′Ω′λcλa∂Ω′∂ϖ′=53∂lnΩ′∂lnϖ′,\begin{eqnarray} x=\frac{5}{3}\frac{\partial \ln \Omega}{\partial \ln \varpi}=\frac{5}{3}\frac{\varpi}{\Omega}\frac{\partial \Omega}{\partial \varpi}=\frac{5}{3}\frac{\lambda^a}{\lambda^c}\frac{\varpi'}{\Omega'}\frac{\lambda^c}{\lambda^a}\frac{\partial \Omega'}{\partial \varpi'}=\frac{5}{3}\frac{\partial \ln \Omega'}{\partial \ln \varpi'}, \end{eqnarray}(23)which allows x to remain unchanged when repeating the operation for Eq. (20): 35ϖ′2Ω′2λbλ2aλ3c∂Ω′∂τ′=−λaλ2aλ3c∂∂ϖ′(ϖ′2Ω′3x13(2x+3)43).\begin{eqnarray} \frac{3}{5} \varpi'^2 \Omega'^2 \frac{\lambda^b}{\lambda^{2a} \lambda^{3c}}\frac{\partial \Omega'}{\partial \tau'}&=- \frac{\lambda^a}{\lambda^{2a} \lambda^{3c}}\frac{\partial}{\partial \varpi'} \left ( \varpi'^2 \Omega'^3 x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}} \right ). \end{eqnarray}(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 κ=ca\hbox{$\kappa = \frac{c}{a}$} the following relations hold: ϖ′τ′=ϖτ,Ω′τ′−κ=Ωτ−κ.\begin{eqnarray} \frac{\varpi'}{\tau'} = \frac{\varpi}{\tau}, \; \; \; \; \; \; \Omega'\tau'^{-\kappa} = \Omega \tau^{-\kappa}. \end{eqnarray}(25)The set of group invariants ξ=ϖτ,y(ξ)=ψτ−κ,\begin{eqnarray} \label{invariants} \xi = \frac{\varpi}{\tau}, \; \; \; \; \; \; y(\xi)=\psi \tau^{-\kappa}, \end{eqnarray}(26)where ψ = ϖ2Ω3 allows the original PDE (Eq. (20)) to be rewritten as ∂ψ∂τ=\begin{eqnarray} \frac{\partial \psi}{\partial \tau}&=&-\frac{\partial}{\partial \varpi} \left ( \psi f(x) \right ), \label{transpde} \end{eqnarray}(27)where f(x) collects all terms depending on x and f′ is its derivative f(x)=dfdx=(2x+3)13(10x+3)3x23·\begin{eqnarray} f(x)&=& x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}} \label{f} \\ \frac{{\rm d}f}{{\rm d}x}&=\frac{(2 x+3)^{\frac{1}{3}}(10 x+3)}{3 x^{\frac{2}{3}}} \cdot \end{eqnarray}The remaining transformation of x (Eq. (21)) yields x=59(∂lnψ∂lnϖ−2).\begin{eqnarray} x=\frac{5}{9} \left ( \frac{\partial \ln \psi}{\partial \ln \varpi}-2 \right ) . \label{lastx} \end{eqnarray}(30)Now, a set of two first-order ordinary differential equations can be derived by applying the group invariants. Starting with the last Eq. (30), the first step is to convert ∂lnψ∂lnϖ\hbox{$ \frac{\partial \ln \psi}{\partial \ln \varpi}$} to ∂lny∂lnξ\hbox{$ \frac{\partial \ln y}{\partial \ln \xi}$}. After two subsequent transformations, it is found that ∂lnψ∂lnϖ=∂lny∂lnξ\hbox{$ \frac{\partial \ln \psi}{\partial \ln \varpi} = \frac{\partial \ln y}{\partial \ln \xi}$}. Consequently, the second and last step is to solve for dydξ\hbox{$ \frac{{\rm d} y}{{\rm d} \xi}$}: dydξ=(9x+10)y5ξ·\begin{eqnarray} \frac{{\rm d} y}{{\rm d} \xi} = \frac{(9x+10)y}{5\xi}\cdot \label{ode1} \end{eqnarray}(31)To obtain the second ODE, thus arriving at an expression for dxdξ\hbox{$ \frac{{\rm d} x}{{\rm d} \xi}$}, an analogous transformation of Eq. (27)has to be performed, which results in dxdξ=−(f−ξ)(95x+2)+κξξf′·\begin{eqnarray} \frac{{\rm d} x}{{\rm d} \xi}=-\frac{\left ( f-\xi \right)\left( \frac{9}{5}x+2\right)+\kappa \xi}{\xi f'} \cdot \label{dx} \end{eqnarray}(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 well-defined problem, it is necessary to specify auxiliary conditions to solve the initial-boundary-value 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, ∂rM ≥ 0 is required. Since radial momentum transport denoted by r3Ω2=GM(r)>0→Ω≠0\begin{eqnarray} r^3\Omega^2=GM(r)>0 \: \: \: \rightarrow \: \: \: \Omega \not= 0 \label{aux_mass} \end{eqnarray}(33)holds for this model (ID15), the auxiliary conditions must also comply with ∂ ln r3Ω2∂ ln r=2x+3≥0→x≥−32\begin{eqnarray} \label{aux} \frac{\partial \, \ln \, r^3\Omega^2}{\partial \, \ln \, r}=2x+3\geq0 \: \: \: \rightarrow \: \: \: x \geq -\frac{3}{2} \end{eqnarray}(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: Ω0(ϖ)=ϖ(κ−23)limξ→∞(y13ξ−κ3).\begin{eqnarray} \Omega_0(\varpi)=\varpi^{(\frac{\kappa-2}{3})}\lim_{\xi \rightarrow \infty}(y^{\frac{1}{3}}\xi^{\frac{-\kappa}{3}}). \end{eqnarray}(35)Since Ω0 should not vanish (Eq. (33)), it follows that y(ξ)∝ξκforξ→∞,\begin{eqnarray} y(\xi) \propto \xi^{\kappa} \: \: \: \text{for} \: \: \: \xi \rightarrow \infty, \end{eqnarray}(36)which consequently leads to the initial condition for self-similar solutions being a power law of radius Ω0(r)∝r53(κ−23),\begin{eqnarray} \label{omega0} \Omega_0(r) &&\propto r^{\frac{5}{3}(\frac{\kappa-2}{3})}, \end{eqnarray}(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 −12\hbox{$-\frac{1}{2}$}. Additionally, the monopole approximation used in the model breaks down when μ passes −1 and approaches −12\hbox{$-\frac{1}{2}$}. However, for up to μ=−34\hbox{$\mu = -\frac{3}{4}$} the error still appears acceptable (ibid.). After identifying μ with 53(κ−23)\hbox{$\frac{5}{3}(\frac{\kappa-2}{3})$}, this imposes an upper limit on κ. Equations (21)and (34)with (37)yield −32≤x≤53(κ−23)\begin{eqnarray} -\frac{3}{2}\leq x \leq \frac{5}{3}\left(\frac{\kappa-2}{3}\right) \end{eqnarray}(38)as a lower limit for the exponent of the power law. Therefore, −32≤53(κ−23)≤−34\begin{eqnarray} -\frac{3}{2}\leq\frac{5}{3}\left(\frac{\kappa-2}{3}\right) \leq -\frac{3}{4} \label{kappa_range} \end{eqnarray}(39)holds, which translates to a domain of κ of −710≤κ≤1320·\begin{eqnarray} -\frac{7}{10} \leq \kappa \leq \frac{13}{20}\cdot \end{eqnarray}(40)

2.5. Boundary conditions

The dependence6 of ξ on ϖ and τ expressed in Eq. (26)yields the following: ξ→0⇔{\begin{eqnarray} \xi \rightarrow 0 &\Leftrightarrow \begin{cases} \tau \rightarrow \infty \: \: \: &\text{for any fixed, positive} \, \varpi \\ \varpi \rightarrow 0 \: \: \: &\text{for any fixed, positive} \, \tau. \end{cases} \end{eqnarray}(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 self-similar 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 𝒢(r,t)=2πr2νΣrdΩdr,\begin{eqnarray} \mathcal{G}(r,t)=2\pi r^2 \nu \Sigma r \frac{{\rm d}\Omega}{{\rm d}r}, \end{eqnarray}(42)which, with the help of Eqs. (4), (6), (10), (11), (13), and (28)can be rewritten7 and expressed in terms of the similarity variables to yield 𝒢(r,t)=\begin{eqnarray} \mathcal{G}(r,t)&=&\eta y \tau^{\kappa} f(x) \label{torque}. \end{eqnarray}(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 \hbox{$\mathcal{G}_{\star}$} at the inner rim at a specific time τ0. Expressing M⋆ and \hbox{$\mathcal{G}_{\star}$} in terms of the group variants yields M⋆(τ)=τ23κ+715ζ223,ζ2=M⋆τ032τ0κ+710,𝒢⋆(τ)=−ητκζ2ζ1,ζ1=𝒢⋆τ0τ0710−ηM⋆τ032\begin{eqnarray} M_{\star}(\tau)&=&\tau^{\frac{2}{3}\kappa +\frac{7}{15}}\zeta_2^{\frac{2}{3}}, \; \; \; \; \; \; \zeta_2= \frac{M_{\star \tau_0}^{\frac{3}{2}}}{{\tau_{0}^{\kappa +\frac{7}{10}}}}, \\ \mathcal{G}_{\star}(\tau)&=&-\eta \tau^{\kappa}\zeta_2\zeta_1, \; \; \; \; \; \; \zeta_1=\frac{\mathcal{G}_{\star \tau_{0}}\tau_{0}^{\frac{7}{10}}}{-\eta M_{\star \tau_{0}}^{\frac{3}{2}}} \end{eqnarray}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.

thumbnail Fig. 1

Solutions of x for different values of κ.

thumbnail Fig. 2

Solutions of y for different values of κ.

3.1. Similarity solutions

Figures 1 and 2 show solutions for x and y (Eqs. (32)and (31)) obtained with \hbox{$\tau_0 = 1, M_\star(\tau_0) = 1, \ln\, \xi_0 = -8, \mathcal{G}_\star = 0$}, 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)): limξ→∞x=53(κ−23)·\begin{eqnarray} \lim_{\xi \rightarrow \infty}x=\frac{5}{3}\left(\frac{\kappa-2}{3}\right) \cdot \label{x_limit} \end{eqnarray}(46)Using this result, it is possible to find an approximation of y for large radii: y=ζ2,∞ξκ.\begin{eqnarray} y&=&\zeta_{2,\infty}\xi^{\kappa}. \end{eqnarray}(47)Furthermore, the influence of the saddle point on the x-axis 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 self-gravity 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 self-gravity, 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 hr\hbox{$\frac{h}{r}$} 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 ξ-710\hbox{$\xi^{\frac{-7}{10}}$} 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 ξ.

thumbnail Fig. 3

Solutions of vϕc2\hbox{$\frac{v_{\varphi}}{c_{\text{2}}}$} and hr\hbox{$ \frac{h}{r}$} with α = 0.01 for κ = 0.2.

thumbnail Fig. 4

Solutions of x for different torques T⋆ applied at the inner rim.

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 \hbox{$\mathcal{G}_\star$}. 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 environment8, according to Kippenhahn et al. (2012) the mean molecular weight can be set to ≈0.002kgmol\hbox{$ 0.002 \frac{\text{kg}}{\text{mol}}$}, which yields η=α43×1.845×1010 m53s\hbox{$\eta = \alpha^{\frac{4}{3}} \times 1.845 \times 10^{10} \, \frac{\text{m}^{\frac{5}{3}}}{\text{s}}$}.

Using AGN scales \hbox{$\hat M = 10^3 \, M_{\odot}$}, \hbox{$\hat r = 1 \text{AU} $}, 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 self-gravitating regime starts at significantly smaller radii with respect to earlier times than in the solution obtained with zero torque. A significant difference to the no-torque 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 \hbox{$\tau_0 = 1, M_\star(\tau_0) = 1, \ln\, \xi_0 = -8, \mathcal{G}_\star = 0$}, and the same or the corresponding values9 of the similarity parameter κ for four different viscosity prescriptions. The relation between the κ introduced in this paper and the κ used by ID15 (hereafter ˜κ\hbox{$\tilde \kappa$}) is 95˜κ+2=κ.\hbox{$\frac{9}{5}\tilde \kappa + 2 = \kappa .$}

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 self-gravitating 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.

thumbnail Fig. 5

Solutions of x for different viscosity prescriptions.

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 ϖ: Ω=y13ξ-23τκ3−23.\begin{eqnarray} \Omega = y^{\frac{1}{3}} \xi^{\frac{-2}{3}} \tau^{\frac{\kappa}{3}-\frac{2}{3}}. \label{Omega_invariants} \end{eqnarray}(48)From the auxiliary conditions, Eq. (33)can be used to obtain M=y23ξ715τ2κ3+715.\begin{eqnarray} M = y^{\frac{2}{3}} \xi^{\frac{7}{15}} \tau^{\frac{2\kappa}{3}+\frac{7}{15}}. \label{M_invariants} \end{eqnarray}(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 yields10Σ=τ2κ3−1115ξ-1115y23(2x+3)2π·\begin{eqnarray} \Sigma = \frac{\tau^{\frac{2\kappa}{3}-\frac{11}{15}} \xi^{\frac{-11}{15}} y^{\frac{2}{3}} (2x+3)}{2\pi}\cdot \label{Sigma_invariants} \end{eqnarray}(50)From these three basic quantities, the accretion rate Ṁ and radial velocity vr can be derived quite easily. Since ∂τ∂t=5η\hbox{$\frac{\partial \tau}{\partial t} = 5\eta$}, the accretion rate is given by Ṁ=5ηMτ(−65x−43+23κ).\begin{eqnarray} \dot M = 5 \eta \frac{M}{\tau} \left ( -\frac{6}{5}x - \frac{4}{3}+\frac{2}{3}\kappa \right ). \label{mdot_invariants} \end{eqnarray}(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: vr=−5ηξ35τ-25(2κ3−43−6x5)(2x+3)·\begin{eqnarray} v_{\text{r}} = -\frac {5 \eta \xi^{\frac{3}{5}} \tau^{\frac{-2}{5}} \left ( \frac{2\kappa}{3}-\frac{4}{3} - \frac{6x}{5} \right ) }{(2x+3)}\cdot \label{v_r_invariants} \end{eqnarray}(52)All the necessary ingredients to find an expression for the vertically integrated viscous dissipation rate Qvis based on the similarity variables are now available and can be written as Qvis=ṀM10πτ95ξ145(2x+3)43x43(2κ3−43−6x5)·\begin{eqnarray} Q_{\text{vis}} = \frac{ \dot M M }{10 \pi \tau^{\frac{9}{5}} \xi^{\frac{14}{5}}} \frac{(2x+3)^{\frac{4}{3}} x^{\frac{4}{3}}}{ \left ( \frac{2\kappa}{3}-\frac{4}{3} - \frac{6x}{5} \right )}\cdot \end{eqnarray}(53)It is a well-known and important result for Keplerian disks that in the radial direction Qvis 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 self-similar 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 no-torque condition in Eq. (A.4)to (2x+3)43\hbox{$(2x+3)^{\frac{4}{3}}$} and transforming the equation back to r yields Qvis≈ṀM10πr3x43(2κ3−43−6x5)·\begin{eqnarray} Q_{\text{vis}} \approx \frac{ \dot M M }{10 \pi r^3} \frac{x^{\frac{4}{3}}}{ \left ( \frac{2\kappa}{3}-\frac{4}{3} - \frac{6x}{5} \right )}\cdot \end{eqnarray}(54)Taking the limit for x→−32\hbox{$x \rightarrow -\frac{3}{2}$} yields Qvis≈ṀM10πr3limx→−32x43(2κ3−43−6x5)=943321π(10κ+7)ṀMr3·\begin{eqnarray} Q_{\text{vis}} &\approx& \frac{ \dot M M }{10 \pi r^3} \lim_ {x \rightarrow -\frac{3}{2}} \frac{x^{\frac{4}{3}}}{ \left ( \frac{2\kappa}{3}-\frac{4}{3} - \frac{6x}{5} \right )} \notag \\ &=& \frac{9}{4} \sqrt[3]{\frac{3}{2}} \frac{ 1 }{ \pi \left (10\kappa + 7 \right )} \frac{\dot M M }{r^3 }\cdot \end{eqnarray}(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 ˜κ\hbox{$\tilde \kappa$} 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 Qvis 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.

thumbnail Fig. 6

Solutions of Σ(r) for different times and torques \hbox{$\mathcal{G}_\star$}.

Table 1

Power law exponents of the radial coordinate for small and large radii.

thumbnail Fig. 7

Solutions of Ω(r) for different times and torques \hbox{$\mathcal{G}_\star$}.

Table 2

Power law exponents of the time dependence for small radii.

The solutions displayed in Fig. 6 and the following figures are given in non-dimensional units and are based on solutions obtained with τ0 = 1,M⋆(τ0) = 1,ln ξ0 = −25, and κ = 0.2. The value of κ corresponds to ˜κ=−1\hbox{$\tilde \kappa = -1$} 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 self-gravitating regime occurs which moves to regions further outward as time progresses, a behaviour already described by Mineshige & Umemura (1996, 1997) for their self-similar α 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 self-gravitating 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 no-torque 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 well-known results for standard stationary α disks for the outer region (SS73). This result is not surprising since the assumption that Teff = Tc 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.

thumbnail Fig. 8

Solutions of M(r) for different times and torques \hbox{$\mathcal{G}_\star$}.

thumbnail Fig. 9

Solutions of Ṁ(r) with α = 0.01 for different times and torques \hbox{$\mathcal{G}_\star$}.

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. Lynden-Bell & 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 quasi-stationary self-gravitating 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).

thumbnail Fig. 10

Solutions of Ṁ(r) with α = 0.01 for different values of κ.

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 vr 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 self-gravity only affects the accretion process in the outer regions of the disk. Large values of κ, i.e. corresponding to self-gravity, yield the highest accretion rates and the slowest temporal decrease. Thus, the result of ID15 that “objects embedded in self-gravitating 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 self-gravitating part of the disk does not accrete material at all. For the α ansatz, the value of κ separating the two regimes is κ = 0 respectively ˜κ=−109\hbox{$ \tilde \kappa = -\frac{10}{9}$}, while for the β ansatz it is ˜κ=−54\hbox{$\tilde \kappa = -\frac{5}{4}$}. Thus, α disks require an even more self-gravitating 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 \hbox{$\hat M = 10^3 \, M_{\odot}$}, \hbox{$\hat r = 1 ~\text{AU} $}, and α = 0.01. The consequently arising scales are listed in Table 3.

Table 3

Scales for an AGN example with α = 0.01.

These scales can now be applied to the figures presented in this section. Applying11 the scales to Figs. 9 and 10 shows that the accretion rate lies at 4.5 × 10-6 × 2 × 105M⊙ yr-1 = 0.9 M⊙ yr-1. For black holes ≲107.7M⊙ 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 109M⊙ in less than 109 yr (Duschl & Strittmatter 2011) – even for the most self-gravitating systems – as can be seen from the following calculation: ∫t0=107.7yrt1=109yr0.9 M⊙ (t×10-3yr)-0.1dt=2.3×108M⊙2.3×108M⊙+107.7M⊙=2.8×108M⊙.\begin{eqnarray} \int_{t_0=10^{7.7} \text{yr}}^{t_1=10^9~\text{yr}} 0.9\, M_{\odot} \,(t \times 10^{-3}~\text{yr})^{-0.1} {\rm d}t &= 2.3 \times 10^8~ M_{\odot} \\ 2.3 \times 10^8 ~M_{\odot} + 10^{7.7}~ M_{\odot} &= 2.8 \times 10^8~ M_{\odot} . \end{eqnarray}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 self-gravitating 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.

1

See e.g. Kato et al. (2008) for an application to accretion disks.

2

This is just a variant of the more common form ν = αcsH, where H is the scale height of the disk.

3

For a modern treatment of the standard α theory, see e.g. Kato et al. (2008) or Frank et al. (2002).

4

T, ρ, and p denote the corresponding quantities in the equatorial plane.

5

To provide an easily legible type, the dimensionless equations in the next sections are written without tildes.

6

Owing to the definition of ϖ in Eq. (19), the behaviour of ϖ corresponds to the behaviour of r.

7

The gravitational constant G is set to unity.

8

In a cold environment a neutral gas can be assumed.

9

This is easily verified if one considers that all solutions have the same limit for large radii.

10

The gravitational constant G is set to unity.

11

Through the definition of our scales in Eq. ((16)), we also account for the numerical value of the gravitational constant G.

References

  1. Ames, W. F. 1965, Nonlinear Differential Equations in Engineering (New York: Academic Press) [Google Scholar]
  2. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bluman, G. W., Cheviakov, A. F., & Anco, S. C. 2010, Applications of symmetry methods to partial differential equations (Springer) [Google Scholar]
  4. Dresner, L. 1998, Applications of Lie’s theory of ordinary and partial differential equations (CRC Press) [Google Scholar]
  5. Duschl, W. J., & Strittmatter, P. A. 2011, MNRAS, 413, 1495 [NASA ADS] [CrossRef] [Google Scholar]
  6. Duschl, W. J., Strittmatter, P. A., & Biermann, P. L. 1998, BAAS, 30, 917 [NASA ADS] [Google Scholar]
  7. Duschl, W. J., Strittmatter, P. A., & Biermann, P. L. 2000, A&A, 357, 1123 [NASA ADS] [Google Scholar]
  8. Eddington, A. S. 1921, Zeitschrift für Physik A Hadrons and Nuclei, 7, 351 [Google Scholar]
  9. Fan, X., Strauss, M. A., Schneider, D. P., et al. 2003, AJ, 125, 1649 [NASA ADS] [CrossRef] [Google Scholar]
  10. Frank, J., King, A., & Raine, D. 2002, Accretion power in astrophysics (Cambridge University Press) [Google Scholar]
  11. Frommer, M. 1928, Die Integralkurven einer gewöhnlichen Differentialgleichung erster Ordnung in der Umgebung rationaler Unbestimmtheitsstellen (Springer) [Google Scholar]
  12. Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [NASA ADS] [CrossRef] [Google Scholar]
  13. Hindmarsh, A. C. 1983, IMACS transactions on scientific computation, 1, 55 [Google Scholar]
  14. Illenseer, T. F., & Duschl, W. J. 2015, MNRAS, 450, 691 [NASA ADS] [CrossRef] [Google Scholar]
  15. Ince, E. L. 1956, Ordinary differential equations (Dover, New York) [Google Scholar]
  16. Kato, S., Fukue, J., & Mineshige, S. 2008, Black-Hole Accretion Disks – Towards a New Paradigm (Kyoto University Press), 1 [Google Scholar]
  17. King, A., Pringle, J., & Livio, M. 2007, MNRAS, 376, 1740 [NASA ADS] [CrossRef] [Google Scholar]
  18. Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar structure and evolution (Springer) [Google Scholar]
  19. Filipov, L., Shakura, N. I., & Liubarokii, Iu 1988, Adv. Space Res., 8, 163 [NASA ADS] [CrossRef] [Google Scholar]
  20. Laughlin, G., & Bodenheimer, P. 1994, ApJ, 436, 335 [NASA ADS] [CrossRef] [Google Scholar]
  21. Lin, D. N. C., & Pringle, J. E. 1987, MNRAS, 225, 607 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lüst, R. 1952, Z. Naturforsch. Teil A, 7, 87 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lynden-Bell, D. 1969, Nature, 223, 690 [NASA ADS] [CrossRef] [Google Scholar]
  24. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
  25. Mineshige, S., & Umemura, M. 1996, ApJ, 469, L49 [NASA ADS] [CrossRef] [Google Scholar]
  26. Mineshige, S., & Umemura, M. 1997, ApJ, 480, 167 [NASA ADS] [CrossRef] [Google Scholar]
  27. Richard, D., & Zahn, J.-P. 1999, A&A, 347, 734 [Google Scholar]
  28. Salpeter, E. E. 1964, ApJ, 140, 796 [CrossRef] [Google Scholar]
  29. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  30. Shakura, N., & Sunyaev, R. 1976, MNRAS, 175, 613 [NASA ADS] [CrossRef] [Google Scholar]
  31. Shlosman, I., Begelman, M. C., & Frank, J. 1990, Nature, 345, 679 [NASA ADS] [CrossRef] [Google Scholar]
  32. Tufillaro, N. B., Abbott, T., & Reilly, J. 1992, An experimental approach to nonlinear dynamics and chaos (Redwood City, CA: Addison-Wesley) [Google Scholar]
  33. Weizsäcker, C. F. 1948, Z. Naturforsch. A, 3, 524 [NASA ADS] [Google Scholar]
  34. Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  35. 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 self-gravity, 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 dxdξ\hbox{$\frac{{\rm d}x}{{\rm d}\xi}$} become zero and three critical points located on the x-axis 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, Runge-Kutta 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 self-gravity 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.

thumbnail Fig. A.1

Phaseplane of dxdξ\hbox{$\frac{{\rm d} x}{{\rm d} \xi}$} for κ = 0.2

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 x-axis 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 non-linear dependence on variable x, which cannot be approximated, the equation is transformed to dfdξ=−(f−ξ)(95x(f)+2)+κξξ\appendix \setcounter{section}{1} \begin{eqnarray} \frac{{\rm d} f}{{\rm d} \xi}=-\frac{\left ( f-\xi \right)\left( \frac{9}{5}x(f)+2\right)+\kappa \xi}{\xi} \label{df} \end{eqnarray}(A.1)and x is now considered a function of f. Now, one can use a Taylor series to approximate f3(x) at x = −1.5, solve for x, and utilize f = 0 at the critical point to arrive at x≈−32+12(32)−14(−f)34.\appendix \setcounter{section}{1} \begin{eqnarray} x\approx-\frac{3}{2}+\frac{1}{2}\left( \frac{3}{2}\right)^{-\frac{1}{4}}(-f)^{\frac{3}{4}}. \label{taylor} \end{eqnarray}(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 f=−(103κ+73)ξ+ζ1ξ710(for ξ≪1),\appendix \setcounter{section}{1} \begin{eqnarray} f=-\left (\frac{10}{3}\kappa +\frac{7}{3} \right )\xi+\zeta_1 \xi^{\frac{7}{10}} \; \; \; \; \; \; (\text{for\,} \xi \ll 1), \end{eqnarray}(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: x=−32+12(32)−14{(103κ+73)ξ−ζ1ξ710}34(for ξ≪1).\appendix \setcounter{section}{1} \begin{eqnarray} x=-\frac{3}{2}+\frac{1}{2}\left( \frac{3}{2}\right)^{-\frac{1}{4}}\!\left \{ \left (\frac{10}{3}\kappa +\frac{7}{3} \right )\xi-\zeta_1 \xi^{\frac{7}{10}} \right \}^{\frac{3}{4}} \; \; (\text{for\,} \xi \ll 1). \label{x(xi)} \end{eqnarray}(A.4)This solution shows that ζ1ξ710\hbox{$\zeta_1 \xi^{\frac{7}{10}}$} 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 y=\appendix \setcounter{section}{1} \begin{eqnarray} y&=&\zeta_2\xi^{-\frac{7}{10}}. \label{y(xi)} \end{eqnarray}(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: M⋆(τ)=limr→0r3Ω2.\appendix \setcounter{section}{1} \begin{eqnarray} M_{\star}(\tau)=\lim_{r \rightarrow 0}r^3\Omega^2. \end{eqnarray}(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 M⋆(τ)=τ23κ+715ζ223,ζ2=M⋆τ032τ0κ+710·\appendix \setcounter{section}{1} \begin{eqnarray} M_{\star}(\tau)=\tau^{\frac{2}{3}\kappa +\frac{7}{15}}\zeta_2^{\frac{2}{3}}, \; \; \; \; \; \; \zeta_2= \frac{M_{\star \tau_0}^{\frac{3}{2}}}{{\tau_{0}^{\kappa +\frac{7}{10}}}} \cdot \label{ceta2} \end{eqnarray}(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): 𝒢=\appendix \setcounter{section}{1} \begin{eqnarray} \mathcal{G}&=& \eta \tau^{\kappa}\zeta_2\xi^{-\frac{7}{10}} f(x). \label{G_trans} \end{eqnarray}(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: ζ1≠0ζ1=0\appendix \setcounter{section}{1} \begin{eqnarray} \zeta_1 \neq 0 &&\; \; \rightarrow \; \; \xi = \left ( \frac{3}{2}\right )^{\frac{10}{21}} \frac{\left (2x+3 \right )^{\frac{40}{21}}}{\zeta_1^{\frac{10}{7}}} \label{xi_x_0} \\ \zeta_1 = 0 &&\; \; \rightarrow \; \; \xi = \left ( \frac{3}{2}\right )^{\frac{1}{3}} \frac{\left (2x+3 \right )^{\frac{4}{3}}}{\left (\frac{10}{3}\kappa +\frac{7}{3}\right )} \cdot \label{xi_x_1} \end{eqnarray}From Eq. (A.8), one can infer that a vanishing torque at the inner boundary demands that limξ→0ξ−710f(x)=limx→−32ξ−710f(x)=0.\appendix \setcounter{section}{1} \begin{eqnarray} \lim_{\xi \rightarrow 0} \xi^{-\frac{7}{10}} f(x) = \lim_{x \rightarrow -\frac{3}{2}} \xi^{-\frac{7}{10}} f(x) = 0 . \end{eqnarray}(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: limx→−32((32)13(2x+3)43(103κ+73))−710x13(2x+3)43=  (32)−730(103κ+73)710limx→−32x13(2x+3)25=0.\appendix \setcounter{section}{1} \begin{eqnarray} &&\lim_{x \rightarrow -\frac{3}{2}} \left (\left ( \frac{3}{2}\right )^{\frac{1}{3}} \frac{\left (2x+3 \right )^{\frac{4}{3}}}{\left (\frac{10}{3}\kappa +\frac{7}{3}\right )}\right )^{-\frac{7}{10}} x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}} = \notag \\ &&\qquad\left ( \frac{3}{2}\right )^{-\frac{7}{30}} \left (\frac{10}{3}\kappa +\frac{7}{3}\right )^{\frac{7}{10}} \lim_{x \rightarrow -\frac{3}{2}} x^{\frac{1}{3}} (2x+3)^{\frac{2}{5}} = 0. \end{eqnarray}(A.12)Hence, one can conclude that ζ1 = 0 denotes a special no-torque 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).

thumbnail Fig. A.2

Logarithmic phaseplane of dxdξ\hbox{$\frac{{\rm d} x}{{\rm d} \xi}$} for κ = 0.2.

Additionally, one can conclude that exactly one solution, i.e. the no-torque 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 x-axis (Frommer 1928). The initial conditions to obtain those solutions which enter the critical point along the x-axis 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 \hbox{$\mathcal{G}$} for ξ → 0 and ζ1 ≠ 0: limx→−32((32)1021(2x+3)4021ζ1107)−710x13(2x+3)43=  (32)−13ζ1limx→−32x13(2x+3)43(2x+3)43=−ζ1.\appendix \setcounter{section}{1} \begin{eqnarray} &&\lim_{x \rightarrow -\frac{3}{2}} \left (\left ( \frac{3}{2}\right )^{\frac{10}{21}} \frac{\left (2x+3 \right )^{\frac{40}{21}}}{\zeta_1^{\frac{10}{7}}}\right ) ^{-\frac{7}{10}} x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}}= \notag \\ &&\qquad\left ( \frac{3}{2}\right )^{-\frac{1}{3}} \zeta_1 \lim_{x \rightarrow -\frac{3}{2}} \frac{x^{\frac{1}{3}} (2x+3)^{\frac{4}{3}}}{\left (2x+3 \right )^{\frac{4}{3}}} = - \zeta_1. \end{eqnarray}(A.13)Thus, using Eq. (A.7), ζ1 can be calculated if one prescribes a torque \hbox{$\mathcal{G}_{\star}$} at the inner rim at a specific time τ0 via 𝒢⋆(τ)=−ητκζ2ζ1,ζ1=𝒢⋆τ0τ0710−ηM⋆τ032·\appendix \setcounter{section}{1} \begin{eqnarray} \mathcal{G}_{\star}(\tau)=-\eta \tau^{\kappa}\zeta_2\zeta_1, \; \; \; \; \; \; \zeta_1=\frac{\mathcal{G}_{\star \tau_{0}}\tau_{0}^{\frac{7}{10}}}{-\eta M_{\star \tau_{0}}^{\frac{3}{2}}}\cdot \end{eqnarray}(A.14)

All Tables

Table 1

Power law exponents of the radial coordinate for small and large radii.

In the text
Table 2

Power law exponents of the time dependence for small radii.

In the text
Table 3

Scales for an AGN example with α = 0.01.

In the text

All Figures

thumbnail Fig. 1

Solutions of x for different values of κ.
In the text
thumbnail Fig. 2

Solutions of y for different values of κ.
In the text
thumbnail Fig. 3

Solutions of vϕc2\hbox{$\frac{v_{\varphi}}{c_{\text{2}}}$} and hr\hbox{$ \frac{h}{r}$} with α = 0.01 for κ = 0.2.
In the text
thumbnail Fig. 4

Solutions of x for different torques T⋆ applied at the inner rim.
In the text
thumbnail Fig. 5

Solutions of x for different viscosity prescriptions.
In the text
thumbnail Fig. 6

Solutions of Σ(r) for different times and torques \hbox{$\mathcal{G}_\star$}.
In the text
thumbnail Fig. 7

Solutions of Ω(r) for different times and torques \hbox{$\mathcal{G}_\star$}.
In the text
thumbnail Fig. 8

Solutions of M(r) for different times and torques \hbox{$\mathcal{G}_\star$}.
In the text
thumbnail Fig. 9

Solutions of Ṁ(r) with α = 0.01 for different times and torques \hbox{$\mathcal{G}_\star$}.
In the text
thumbnail Fig. 10

Solutions of Ṁ(r) with α = 0.01 for different values of κ.
In the text
thumbnail Fig. A.1

Phaseplane of dxdξ\hbox{$\frac{{\rm d} x}{{\rm d} \xi}$} for κ = 0.2
In the text
thumbnail Fig. A.2

Logarithmic phaseplane of dxdξ\hbox{$\frac{{\rm d} x}{{\rm d} \xi}$} for κ = 0.2.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.