aa r X i v : . [ a s t r o - ph ] J a n RIVISTA DEL NUOVO CIMENTO
Vol. ?, N. ? ? Self-gravitating accretion discs
Giuseppe Lodato ( ) ( ) Department of Physics and Astronomy, University of Leicester, University Road, LE1 7RH,Leicester, UK (ricevuto ?)
Summary. —I review recent progresses in the dynamics and the evolution of self-gravitating ac-cretion discs. Accretion discs are a fundamental component of several astrophysicalsystems on very diverse scales, and can be found around supermassive black holes inActive Galactic Nuclei (AGN), and also in our Galaxy around stellar mass compactobjects and around young stars. Notwithstanding the specific differences arisingfrom such diversity in physical extent, all these systems share a common featurewhere a central object is fed from the accretion disc, due to the effect of turbulenceand disc instabilities, which are able to remove the angular momentum from the gasand allow its accretion. In recent years, it has become increasingly apparent thatthe gravitational field produced by the disc itself (the disc’s self-gravity) is an im-portant ingredient in the models, especially in the context of protostellar discs andof AGN discs. Indeed, it appears that in many cases (and especially in the colderouter parts of the disc) the development of gravitational instabilities can be one ofthe main agents in the redistribution of angular momentum. In some cases, theinstability can be strong enough to lead to the formation of gravitationally boundclumps within the disc, and thus to determine the disc fragmentation. As a result,progress in our understanding of the dynamics of self-gravitating discs is essentialto understand the processes that lead to the feeding of both young stars and ofsupermassive black holes in AGN. At the same time, understanding the fragmenta-tion conditions is important to determine under which conditions AGN discs wouldfragment and form stars and whether protostellar discs might form giant gaseousplanets through disc fragmentation.PACS – Hydrodynamics.PACS – Accretion and accretion disks.
1. – INTRODUCTION
Disc-like or flattened geometries are very common in astrophysics, from the large scaleof spiral galaxies down to the small scales of Saturn’s rings. In both these examples, thediscs are characterized by a prominent structure either in the form of a spiral structure(in the case of galaxy discs) or in the form of gaps and spirals on small scales (in the c (cid:13) Societ`a Italiana di Fisica GIUSEPPE LODATO case of Saturn’s rings). The dynamics underlying the development of such structuresis determined by the propagation of density waves, and the important role of the disc’sself-gravity in their development has been clearly recognized (see, for example, [1]). Inthe above examples the system is either collisionless , such as in the case of spiral galax-ies, where the dominant component is constituted of stars, or particulate , such as in thecase of Saturn’s rings, where the dominant components (mainly rocks and pebbles) doundergo inelastic collisions, but the system cannot be simply described in terms of hy-drodynamics. In the last thirty years increasing attention has been given to fluid discs,where the dynamically active component is gaseous. Here, dissipative effects, associatedwith friction or viscosity, can significantly affect the dynamics of the disc. Through dis-sipation the fluid elements in the disc can lose their energy or, more fundamentally, theirangular momentum, as I describe in more detail below, and can fall towards the bottomof the potential well, hence accreting on to a central gravitating body. Such a system,where a disc feeds a central object through accretion, under the effect of viscous forces,is called an accretion disc .Accretion discs are found in very different contexts and over a wide range of physicalscales. On the largest scale, they are one of the main physical ingredient to power thecentral engine of Active Galactic Nuclei (AGN), through accretion on to a supermassiveblack hole. The mass of the central black hole ranges from 10 M ⊙ to 10 M ⊙ and theaccretion discs can extend out to large distances, of the order of about 1 pc, where rotatinggaseous discs have often been detected through water maser emission [2-4]. On thegalactic scale, accretion discs around compact objects, such as white dwarfs, neutron starsand stellar mass black holes, are often found in binary systems, where the compact object(with a mass of the order of 1 M ⊙ ) is fed by material outflowing from the companion.Historically, this has been one of the first contexts where accretion disc theory has foundwidespread application (a detailed review with an emphasis on accretion in galacticbinary system and AGN can be found in the textbook by Frank, King and Raine [5]).Another class of objects where accretion discs play an important role are Young StellarObjects (YSO). Here, the central accreting object is a young protostar, which receivesmost of its mass from a surrounding protostellar discs. In this case, protostellar discs playalso another key role, as the site where planet formation takes place. Understanding thedynamics of the disc and the development of the various instabilities that might occur init, and possibly lead to turbulence, is therefore essential if one wants to understand someproperties of our own Solar system, as well as those of extra-solar planets. A thoroughdescription of the current observational state of protostellar discs can be found, forexample, in the recent Protostars and Planets V book [6].As for the case of spiral galaxies, also for accretion discs the effects of the disc self-gravity can be very important. Indeed, also in this case the development of gravitationalinstabilities can lead to the formation of grand design spiral structures, which deeplyaffect the structure and the dynamics of the disc. The analysis of such instabilities natu-rally reveals several similarities between the case of spiral galaxies and that of accretiondiscs, but also some important differences. First of all, as mentioned above, accretiondiscs are fluid and, as such, they are intrinsically dissipative. Indeed, the evolution ofgravitationally unstable accretion discs, as I will discuss below, is strongly dependent onthe gas cooling rate. On the other hand, gravitationally unstable discs of stars, in theabsence of any gas, have a natural tendency to dynamically heat up. Thermal energy, ordisordered kinetic energy, plays a key role in stabilizing self-gravitating discs. However,while for stellar discs this energy is mostly in the form of disordered stellar motions,with the possibility of an anisotropic velocity dispersion, for gaseous discs the disordered ELF-GRAVITATING ACCRETION DISCS kinetic energy content is mostly in the form of thermal energy. Furthermore, for gaseousdiscs it is expected that the role of resonances is much reduced with respect to a collision-less system. Finally, while for galaxy discs the disc mass can be a substantial fraction ofthe total mass of the system and therefore give a large contribution to the gravitationalpotential, in the case of accretion discs the disc mass is often much smaller than the massof the central object (with some important exceptions, described below).Historically, accretion disc theory has concentrated on the non self-gravitating case,the effects of self-gravity having only been discussed occasionally [7-10]. In the last tenyears, on the other hand, the important role of the disc self-gravity has been clearly rec-ognized, partly due to improved observations, which have shown that in several observedsystems (on all scales, from AGN to protostars) the disc mass can be high enough to havea dynamical role, and partly due to the increased computational resources, which haveallowed a detailed numerical investigation of the development of gravitational instabilitiesin the non-linear regime. However, the early papers mentioned above already reveal themain directions around which research in this field would have later developed. Currently,most of the interest revolves around the three issues of self-regulation [7], disc fragmen-tation [8] and angular momentum transport induced by gravitational instabilities [9, 10].This again reveals a difference with respect to the galaxy discs case, where most of theattention has been traditionally given to the morphology induced by self-gravity.In this paper, I present an overview of the dynamics of self-gravitating accretion discs.Accretion disc theory is a very vast topic and has been treated extensively elsewhere[5, 11], so I will not repeat it here. I will rather summarize the most salient generalfeatures, with particular emphasis on those features which are most affected by the discself-gravity. Similarly, a detailed description of all the different astrophysical systemswhere self-gravitating accretion discs play a role (protostellar discs, AGN,...) is obviouslybeyond the scope of this review. The structure of the paper is as follows. In section Idescribe the basic equations that determine the evolution of accretion discs and the maineffects of the disc self-gravity. In section I discuss the development of gravitationalinstabilities, and present some models which incorporate in a simple way the salientphysics behind their development. In section I review the current state of numericalsimulations of gravitationally unstable gaseous discs. In section I focus on the transportproperties induced by gravitational instabilities, while in section I describe the relatedissue of disc fragmentation. Finally, in the last two sections I will present some examplesof observed systems where the concepts described in this paper find a natural application,and in particular, in section I will discuss the impact of self-gravity in the dynamics ofprotostellar discs with some application related to planet formation. In section I focuson the impact of self-gravity in AGN discs also showing how self-gravitating accretiondiscs might have played an important role at high red-shifts, by allowing the formationof the seeds of supermassive black holes.
2. – SELF-GRAVITATING ACCRETION DISC DYNAMICS2 .1.
The thin disc approximation . – One of the main properties of accretion discs isthat they are often thin. This means that the typical lengthscale in the vertical direction,the disc thickness H , is much smaller that the radial distance from the central object R .For example, AGN discs are generally very thin, with aspect ratios H/R ≈ . − . H/R ≈ .
1. As we willdiscuss below, this has a significant impact on the conditions under which the disc is self-gravitating, and hence this intrinsic difference between the AGN case and the protostellar
GIUSEPPE LODATO case will be reflected in a significantly different behaviour of these two kind of systemsin the self-gravitating regime.The thin disc condition allows us to treat the disc, to a first approximation, as in-finitesimally thin, and to introduce a small quantity , the aspect ratio
H/R ≪
1. Thisimplies that most of the equations we are going to use can be integrated in the verticaldirection, and rather than dealing with quantities per unit volume (such as the density ρ ), we will deal instead with quantities per unit surface (such as the surface densityΣ = R ∞−∞ ρ d z ). When “volume” quantities are needed (for example the viscosity ν ),these will generally be understood as vertically averaged.The fundamental ordering of lengthscales H/R ≪ .4) that this orderingis equivalent to requiring that the sound speed c s is much smaller than the rotationalvelocity v φ . It is useful to anticipate another important relation between the relevantvelocities, derived by requiring that accretion takes place on a long timescale. This con-dition implies that the radial velocity v R should be smaller than both the sound speedand the rotational speed. We can therefore summarize these relations as: v R ≪ c s ≪ v φ . (1) Finally, it is worth adding a word of caution on the applicability of the thin disccondition. Indeed, in some cases the disc might not be very thin. This happens forexample, in the hotter inner regions of AGN discs. Protostellar discs, as mentionedabove, are generally relatively thick. One should therefore be cautious when using thethin disc approximation, especially when considering quantities averaged in the verticaldirection. .2. Basic equations . – The evolution of accretion discs is determined by the basicequations of viscous fluid dynamics: the continuity equation and Navier-Stokes equations.Consider an accretion disc rotating around a central object of mass M . Given thegeometry of the problem, we adopt cylindrical polar coordinates centred on the centralobject and we assume that to first order the disc is axisymmetric, so that no quantitydepends on the azimuthal angle φ . Consider then a disc with surface density Σ( R, t ).In cylindrical coordinates, the continuity equation (integrated in the vertical direction,following the thin disc approximation) reads ∂ Σ ∂t + 1 R ∂∂R ( R Σ v R ) = 0 . (2) To determine the velocity v , we use the momentum equation. For an inviscid fluidthis takes the form of Euler’s equation. For the moment, we consider the full three-dimensional (i.e. not vertically integrated) equation ρ (cid:20) ∂ v ∂t + ( v · ∇ ) v (cid:21) = − ∇ P − ρ ∇ Φ . (3)On the right-hand side there are the various forces acting on the fluid. First of all, wehave the term describing pressure forces, where P is the pressure and ρ is the density.We then have gravity, represented by the last term on the right-hand side, where Φ isthe gravitational potential. This in general includes both the contribution of the central ELF-GRAVITATING ACCRETION DISCS point mass and that of the disc. A detailed description of these different contributions isprovided in the next section. For an infinitesimally thin disc, we can obtain the potentialassociated with the gas surface density Σ (the ‘self-gravity’) from the solution to Poisson’sequation: ∇ Φ sg = 4 πG Σ δ ( z ) , (4)where δ ( z ) is the Dirac δ -function, indicating that the disc density is confined to themidplane.For an accretion disc, however, Euler equation is not sufficient, because it does notinclude the important terms due to viscosity. Including viscous forces, the relevantmomentum equation is the Navier-Stokes equation. This reads ρ (cid:20) ∂ v ∂t + ( v · ∇ ) v (cid:21) = − ( ∇ P − ∇ · σ ) − ρ ∇ Φ . (5)The second term on the right hand side in equation (5) contains the stress tensor σ anddescribes the effect of viscous forces. This term plays a very important role in accretiondisc dynamics. The nature of the disc viscosity and of the stress tensor σ has beentraditionally one of the main unsolved theoretical issues in accretion disc dynamics. Inits simplest form, it can be assumed to be given by classical shear viscosity, so that theonly non-vanishing component of σ in a circular shearing flow is the Rφ component,proportional to the rate of strain R Ω ′ (where Ω ′ = dΩ / d R ): σ Rφ = ηR dΩd R , (6)where Ω = v φ /R is the angular velocity and we have introduced the shear viscositycoefficient η . Equation (5) has three components in the radial, vertical and azimuthaldirections, respectively, and each one defines some important properties of accretion discs.In the following three sections we will examine in turn each of these three equations. .3. Radial equilibrium: centrifugal balance . – Let us first consider the radial compo-nent of equation (5). Here, the ordering of velocities described above turns out to beparticularly useful. In particular, on the left-hand side, the first term, that contains theradial velocity v R is negligible with respect to the second term, which (when using theappropriate expression for the differential operators in cylindrical coordinates) gives riseto the centrifugal term v φ /R . On the right-hand side, the viscous term vanishes and thepressure term can be obtained using the equation of state. If the gas is barotropic (i.e.if pressure only depends on density), the sound speed is simply defined as: c = d P d ρ . (7)We then have: 1 ρ ∂P∂R = c ρ ∂ρ∂R ∼ c R . (8)
GIUSEPPE LODATO
Since c s ≪ v φ , also this term is second order with respect to the leading term. The onlyterm on the right hand side able to balance the centrifugal force is therefore gravity. Wethen have, to first order: v φ R ≃ ∂ Φ ∂R . (9)For accretion discs the gravitational field is generally dominated by the central compactobject, but in some cases the disc self-gravity can also produce a sizable effect. Thecentral object contribution to the field is simply given by ∂ Φ ∂R = GMR (10)Thus, in the limit where we neglect self-gravity, the disc rotation is Keplerian, with v φ = v K ≡ p GM/R and where the angular velocity is given by Ω = Ω K ≡ p GM/R .The disc contribution to the gravitational field depends on the matter distributionthrough Poisson’s equation. For an infinitesimally thin disc, the relation between Σand the disc gravitational field can be written, for example, using the complete ellipticintegrals of the first kind, K and E∂ Φ disc ∂R ( R, z ) = GR Z ∞ (cid:20) K ( k ) − (cid:18) k − k (cid:19) (cid:18) R ′ R − RR ′ + zRR ′ (cid:19) E ( k ) (cid:21) r R ′ R k Σ( R ′ )d R ′ , (11)where k = 4 RR ′ / [( R + R ′ ) + z ] [12, 13]. The above equation in general describesa complicated relation between Σ and the gravitational field, which can be obtainednumerically. However, there are a few cases where the relation is particularly simple andcan be derived analytically. One interesting case is the one described by Mestel [14],where Σ = Σ R /R and the disc extends out to infinity. In this case the gravitationalfield at the disc midplane takes a particularly simple form: ∂ Φ disc ∂R ( R ) = 2 πG Σ( R ) . (12)In the extreme limit, where the central object contribution in negligible and a Mesteldisc dominates the gravitational field [15], the rotation curve is thus flat: v φ = R ∂ Φ disc ∂R ( R ) = 2 πG Σ( R ) R = 2 πG Σ R , (13)and the angular velocity Ω ∝ /R . Interestingly, as described in more detail in section .3 below, it turns out that a Mestel disc is the natural solution to the accretion discequations in a steady state, in the limit of strongly self-gravitating discs [15]. Now, whilein actual accretion discs the central object mass is not expected to be negligible, it isinteresting to note that in many cases observations of protostellar discs imply surfacedensities Σ ∝ /R [16]).Having obtained the rotation curve in the two limiting cases of negligible disc self-gravity and negligible central object gravity, we can now ask under which conditions dowe expect the different contributions to dominate. Comparing Eq. (10) to Eq. (12), ELF-GRAVITATING ACCRETION DISCS we see that in order for the disc contribution to be comparable to the central object werequire M disc ( R ) ≈ π Σ( R ) R ≈ M, (14)that is, we require disc masses of the order of the central object mass. This is confirmed bydetailed models of self-regulated discs [13]. In most cases, the accretion disc mass is muchlower than the central object and so the Keplerian approximation is generally appropriate.However, there are a few important exceptions where non-Keplerian rotation has beenobserved and can be interpreted in terms of relatively massive discs. This occurs forexample for some protostellar discs around massive stars [17, 18] and also in some AGNdiscs, such as in the nucleus of the active galaxy NGC 1068 [3, 19]. We defer a moredetailed description of such cases to section and below.The above derivation is only valid to first approximation, since we have neglectedthe pressure forces in the radial direction. This is generally a good approximation, butin some cases the corrections due to pressure effects might play an important role. Inparticular, this is the case when one considers the dynamics of small solid bodies withinthe disc (a very important component involved in the process of planet formation, asdiscussed below in section .5). The solids are not subject to pressure and would thusmove on exactly Keplerian orbits, therefore generating a small velocity difference withrespect to the gas (in a process somewhat analogous to the so called asymmetric drift instellar dynamics), which in turn determines a very fast migration of the solids [20-22].Let us then calculate this correction. The radial equation of motion, including pressureterms is: v φ R = 1 ρ ∂P∂R + v R , (15)where we have indicated with v = R (dΦ / d R ) the circular velocity in the absence ofpressure. If we introduce the sound speed c s and make the simple assumption that thedensity ρ has a power law dependence on radius, with power law index − β (so that ρ ∝ R − β ), we obtain: v φ = v circ " − β (cid:18) c s v circ (cid:19) / . (16)In cases where the disc is hot, in the sense that the sound speed is non negligible withrespect to v circ , then pressure effects will give a sizable contribution to the rotation curve. .4. Vertical structure: hydrostatic equilibrium . – We now consider the vertical com-ponent of equation (5). Here, since the velocity in the vertical direction is very small (thedisc being confined to the equatorial plane), we can neglect the left-hand side of the equa-tion altogether. The viscous force vanishes as well, since the only non-zero componentof the stress is in the Rφ direction. We are therefore left with just two terms to balance:gravitational force and pressure force in the vertical direction. Such a situation, wheregravity is balanced by pressure, is called “hydrostatic balance”. The equation reads:1 ρ ∂P∂z = − ∂ Φ ∂z . (17) GIUSEPPE LODATO
We again first consider the case of non self-gravitating discs. In this case, the verticalcomponent of gravity is simply: − GMr zr ∼ − GMR zR , (18)where r is the spherical radius and the approximation is valid for small z . We can thenrewrite Equation (17) as: c ρ ∂ρ∂z = − GM zR . (19)The solution to this equation is straightforward if the sound speed is independent of z .The vertical density profile in this case turns out to be a Gaussian: ρ ( z ) = ρ exp (cid:20) − GM z R c (cid:21) = ρ exp (cid:20) − z H (cid:21) , (20)where ρ is the midplane density and we have introduced a typical vertical scale-lengthin the non-self-gravitating case H nsg : H nsg = c s Ω K . (21)It is also worth introducing a slightly modified scale-length( ) H ⋆ = p π/ H nsg , so thatthe surface density Σ = R ∞−∞ ρ d z is related to the midplane density by Σ = 2 ρ H ⋆ . Notethe simple relation between the thickness, the sound speed and the angular velocity. Thedisc aspect ratio is then: H nsg R = c s v K , (22)which then demonstrates, as anticipated above, that the requirement that the disc be thinis equivalent to requiring that the disc rotation is highly supersonic, i.e. that v K ≫ c s .Let us now consider the case of a disc dominated by self-gravity. In this case, wesimply have to replace the gravitational force on the right hand side of equation (19)with the force produced by a slab of gas with surface density Σ. If the slab is radiallyhomogeneous we have c ρ ∂ρ∂z = − πG Σ( z ) , (23)where Σ( z ) = R z − z ρ ( z ′ )d z ′ . The solution of the hydrostatic balance in this case is moredifficult but can be done analytically in the case of constant sound speed (such case is ( ) The small difference between H nsg and H ⋆ is generally neglected in many papers, and thetwo expressions are both commonly used to indicate the disc thickness in the non-self-gravitatingcase. ELF-GRAVITATING ACCRETION DISCS generally referred to as the “self-gravitating isothermal slab”) [23]. The density profilein this case is not Gaussian, but is given by: ρ ( z ) = ρ ( z/H sg ) , (24)where the thickness in the self-gravitating case is: H sg = c πG Σ , (25)where in this case the surface density and the midplane density are related by Σ = 2 ρ H sg .Now, we can ask the same question that we asked above for the radial contributionto the gravitational field. How massive has to be the disc, in order for self-gravityto significantly affect its vertical structure? Note that the vertical component of thegravitational field produced by the disc is of the order of 2 πG Σ ∼ GM disc /R . On theother hand, the vertical component of the star’s gravitational field is smaller than theradial ( GM/R ) by a factor H nsg /R (cf. Eq. (18)) and thus [7] the vertical structure ofthe disc is affected by self-gravity already when the disc mass is of the order of: M disc M ≈ H nsg R ≪ . (26)Thus, for thin discs the vertical structure of the disc is affected by the disc self-gravity,even when the disc mass is much smaller than the central object mass. Another way oflooking at the same question is to compare the two expressions for the thickness obtainedabove: H sg H nsg = c s Ω K πG Σ . (27)The two expressions become comparable when the dimensionless parameter on the righthand side is of order unity. It is then easy to see that this occurs when M disc /M ≈ H nsg /R . The parameter on the right (to within factors of order unity) is nothing else thanthe well known Q parameter that controls the development of gravitational instabilitiesin the disc [24, 25]. These will be discussed at length starting from section below.It is also possible to solve the vertical hydrostatic balance equation in the combinedcase where both the central object and the disc contribute to the gravitational field. Inthis case, we have: 1 ρ ∂ρ∂z = − πG Σ( z ) + Ω zc . (28)This equation does not have an analytical solution and has to be solved numerically.However, an approximate but accurate solution can be obtained also in this case [13],with the thickness given by H = H sg (cid:18) H ⋆ H sg (cid:19) s (cid:18) H sg H ⋆ (cid:19) − , (29) GIUSEPPE LODATO
Fig. 1. – Disc thickness from on Eq. (29) as a function of H sg /H ⋆ . The two dashed lines indicatethe limiting cases H = H sg and H = H ⋆ . which in this form holds for Keplerian rotation. Further correction terms for non-Keplerian cases can be easily incorporated in the above expression [13]. Note that Eq.(29) correctly reproduces the limiting cases of fully self-gravitating discs, H sg ≪ H ⋆ ,where H ≈ H sg and that of non-self-gravitating disc, H ⋆ ≪ H sg , where H ≈ H ⋆ . Inthe following we will simply use H to indicate the disc thickness. Figure 1 shows thethickness calculated from Eq. (29) as a function of H sg /H ⋆ , along with the two limitingcases H = H sg and H = H ⋆ . Note that when H sg ≈ H ⋆ (which, as I discuss below, is acommon case for self-gravitating discs), we have H ≈ . H sg = 0 . H ⋆ . .5. Angular momentum conservation: the issue of ‘anomalous’ viscosity . – Let usfinally discuss the last component of Navier-Stokes equation, in which the viscous termplays an important role. First, we integrate Eq. (5) in the vertical direction:Σ (cid:18) ∂ v ∂t + ( v · ∇ ) v (cid:19) = − ( ∇ P − ∇ · T ) − ρ ∇ Φ , (30)where T is the vertical integral of the stress tensor (Eq. (6)), whose only non-vanishingcomponent is T Rφ = ν Σ R dΩd R . (31)
ELF-GRAVITATING ACCRETION DISCS In Eq. (31) the vertically averaged kinematic viscosity ν is defined through ν Σ = Z ∞−∞ η d z. (32)Combining the φ component of Eq. (30) with continuity equation (Eq. (2)) andafter a little algebra, we obtain the expression for angular momentum conservation for aviscous disc: ∂∂t (Σ Rv φ ) + 1 R ∂∂R ( Rv R Σ Rv φ ) = 1 R ∂∂R ( ν Σ R Ω ′ ) , (33)where Ω ′ = dΩ / d R . The physical interpretation of the above expression is readily ap-parent. The left hand side is the Lagrangian derivative of the angular momentum perunit mass Σ Rv φ , while the right hand side is the torque exerted by viscous forces.With the help of the continuity equation (Eq. (2)), we can obtain the radial velocityfrom equation (33): v R = 1 R Σ( R Ω) ′ ∂∂R ( ν Σ R Ω ′ ) , (34)which can be inserted back in equation (2) to finally give: ∂ Σ ∂t = − R ∂∂R (cid:20) R Ω) ′ ∂∂R (cid:0) ν Σ R Ω ′ (cid:1)(cid:21) , (35)which is the master equation that describes the evolution of an accretion disc. In general,in Eq. (35), ν can depend on both radius and Σ, and even Ω, for self-gravitating discs,may be a function of Σ, thus making Eq. (35) a non-linear evolutionary equation.However, if Ω and ν are independent of Σ, then Eq. (35) is a diffusion equation and asimple analysis of its structure shows that the disc surface density evolves on a timescaleof the order of t visc = R ν . (36)In particular, in the important case of Keplerian rotation, the diffusion equation has thefollowing form: ∂ Σ ∂t = 3 R ∂∂R (cid:20) R / ∂∂R (cid:16) ν Σ R / (cid:17)(cid:21) . (37)The above discussion shows the extremely important role of viscosity in the evolutionof accretion discs. However, unfortunately, the ultimate physical origin of viscosity is alsothe major unsolved problem in accretion disc theory. As it was soon realized, standardkinetic viscosity due to collisions between gas molecules is far too low to account for thetransport of angular momentum needed in observed accretion discs. In order to see this,let us consider the viscous timescale introduced above, t ν = R /ν . It is convenient toexpress this quantity in units of the dynamical timescale t dyn = Ω − : t ν t dyn = R Ω ν . (38) GIUSEPPE LODATO
Note that the ratio above is nothing else than the Reynolds number Re of the flow. Stan-dard, collisional viscosity can be expressed as the product between the typical randomvelocity of molecules (that will be of order of the sound speed c s ) times the collisionalmean free path λ = 1 / ( nσ coll ), where n is the number density of the gas and σ coll is thecollisional cross-section. We then have: λ = 1 nσ coll = µm p ρσ coll = (cid:18) µm p Σ σ coll (cid:19) H, (39)where m p ≈ − g is the proton mass, µ is the average molecular weight (we can takeit to be ≈
2, for simplicity) and where ρ and Σ are the usual volume and surface densityof the disc, respectively. Now, substituting ν = λc s in equation (38), we get: t ν t dyn = (cid:18) Σ σ coll µm p (cid:19) (cid:18) HR (cid:19) − . (40)To give some ideas of the numbers involved, let us assume that the collisional cross sectionis simply of the order of the size of an hydrogen molecule σ coll ≈ − cm . In orderto give a rough estimate of the disc surface density let us consider a typical protostellardisc, and assume (to be conservative) that the disc mass is ≈ . M ⊙ and that thedisc size is ≈
50 AU, which gives us Σ ≈ . M ⊙ / (50AU) ≈ . Inserting thisnumbers in Eq. (40), and also assuming that H/R ≈ .
1, we get t ν /t dyn ≈ . Sincethe dynamical timescale for protostellar discs is of the order of a few years, the aboveestimates would lead to the conclusion that in order to accrete a disc of mass 0 . M ⊙ from a distance of 50 AU, it would take much longer than the Hubble time! Clearly, themagnitude of viscosity must be much larger than the simple collisional one estimatedabove. Similar estimates can also be done for discs in other contexts, such as AGN andgalactic binaries, where the gas is expected to be ionized and the transport coefficientwill be mostly determined by plasma processes [5]However, the very small magnitude of collisional viscosity offers a possible way outfrom the apparent paradox. Indeed, as remarked above, the ratio of the viscous timescaleto the dynamical one is also equal to the Reynolds number of the flow. The fact thatcollisional viscosity gives such a large estimate of this ratio also implies that the Reynoldsnumber of the accretion flow is correspondingly large. Now, it is well known that forhigh Reynolds numbers a flow is subject to the development of turbulence and we shouldtherefore expect that the flow in an accretion disc should be highly turbulent. In theseconditions, viscosity can be much higher because in this case angular momentum isexchanged not through collisions of individual gas molecules, but by the mixing of fluidelements moving around in the disc due to turbulence. The typical length-scale of suchmotion can be several orders of magnitude larger than the collisional mean free pathand consequently transport becomes much more effective. The question is then how tocharacterize such an ‘anomalous’ viscosity. For decades, modelers have relied on a verysimple and successful parameterization of viscosity in terms of an unknown dimensionlessparameter, called α [26]. The argument is very simple and draws from the fact that thestress tensor has the physical dimension of a pressure, that is a density times the squareof a velocity. The simplest assumption is then to take the stress tensor to be justproportional to the (vertically integrated) pressure Σ c : T Rφ = d ln Ωd ln R α Σ c , (41) ELF-GRAVITATING ACCRETION DISCS where d ln Ω / d ln R ( ∼ -3/2 for a quasi Keplerian rotation curve) is just a number (toremind us that the viscous stress is proportional to the rate of strain), and where α isthe proportionality factor between the stress tensor and the pressure.Another way of expressing the α -prescription is by considering the kinematical vis-cosity ν . Indeed, Eq. (41) is equivalent to: ν = αc s H. (42)This form of the α -prescription offers a simple way to put some constraints on α . Themagnitude of turbulent viscosity is given roughly by ν ∼ ˆ vl , where l is the typical size ofthe largest eddies in the turbulent pattern ad where ˆ v is the typical turbulent velocity.Now, it is unlikely that the turbulence will be highly supersonic, since otherwise it wouldbe easily dissipated through shocks and we thus have ˆ v . c s . An upper limit to thesize of the largest eddies l is obviously given by the disc thickness H (if we consideran isotropic turbulence). These two upper limits, taken together, clearly imply that α <
1. Typically quoted values for α , derived from comparing theoretical evolutionarymodels to observations range from α ≈ .
01 [27] for protostellar discs to α ≈ . α ≈ − [29, 30].Note that the arguments above do not necessarily imply that the stress tensor isstrictly proportional to the gas pressure. Rather they merely indicate that we mayexpect the ratio of the two to be smaller than unity. While in many cases the parameter α is taken to be a constant, it obviously does not need to be so. In this respect, the α -prescription is not a theory of viscosity in accretion discs, since it only replaces anunknown parameter, ν , with another unknown parameter, α . By using the α -prescriptionwe simply measure the stress in units of the local pressure, and we realize that in theseunits the stress should not be larger than unity, if it is driven by local turbulence.Modern accretion disc theory is thus based on the assumption that ‘anomalous’ trans-port of angular momentum is provided by turbulent viscosity. Clearly, in order to main-tain the turbulence for long enough, energy needs to be continuously injected at thelargest scales. This is usually associated with the development of instabilities in the disc.Hydrodynamic, magneto-hydrodynamic and gravitational instabilities have all been con-sidered in this respect and in recent years several numerical simulations have shown thatmany of the instabilities considered can indeed lead to sustained turbulence. MHD in-stabilities are often considered to be the most likely cause of ‘anomalous’ viscosity. Adetailed discussion of these processes can be found elsewhere [31]. In the context of thispaper, I will put more emphasis on the role of gravitational instabilities, and the trans-port of angular momentum and energy associated with gravitational instabilities will bediscussed in detail below in section . Gravitational instabilities are the most likely al-ternative to MHD instabilities in providing a source of transport in accretion discs, andtheir role is now well established, especially in the context of protostellar discs [32-34].Before moving on, it is worth pointing out that the term “turbulence” is sometimemisused in astronomy, often indicating a simply fluctuating velocity field, which doesnot necessarily display the typical turbulent cascade from large to small scales. While insome cases simulations of hydrodynamic instabilities do produce turbulent vortices, oftenthe structure consists of larger scale coherent disturbances. This is particularly true forthe case of gravitational instabilities, that develop in the form of spiral structures andfor which the term “gravo-turbulent”, sometimes used, is probably inappropriate. GIUSEPPE LODATO .6. Simple Keplerian disc models . – Let us see what kind of evolution is associatedwith the equation of motion described above in a few simple cases. Let us first considerthe case of Keplerian rotation, described by Eq. (37). A particularly simple but instruc-tive configuration is when the viscosity ν is a constant, independent of radius. In thiscase Eq. (37) can be solved analytically [35]. We take as initial condition for the surfacedensity an infinitesimally thin ring of mass m , whose shape is a δ function centered atsome radius R : Σ( R, t = 0) = m πR δ ( R − R )(43)The solution to equation (37) in this case is:Σ( x, τ ) = mπR x − / τ e − (1+ x ) /τ I / (cid:18) xτ (cid:19) , (44)where I / is a modified Bessel function of the first kind, x = R/R and τ = 12 νt/R .The evolution of the surface density is shown in the left panel of Fig. 2, where the differentlines refer (from top to bottom) to different times τ = 0.01, 0.05, 0.1 and 0.15. Thisfigure illustrates quite clearly why this example is called “the spreading ring”. We see herethat, indeed, under the action of viscous forces the disc, rather than merely accreting,spreads both inwards ad outwards. This outward spreading is needed in order to conserveangular momentum. Viscosity simply acts on the fluid to redistribute angular momentumbetween different annuli transporting it from the inside out, the total angular momentumbeing conserved. In this way, as the inner parts of the disc lose their angular momentumand accrete, the outer parts take up the excess angular momentum and move outward. Adetailed analysis of the analytic solution [5,11] shows that the transition between inwardand outward motion occurs at a radius of order R tr ∼ tν/R ∼ R ( t/t ν ) and is thereforean increasing function of time. Therefore, at late times, only the outermost parts of thedisc move outwards, and most of the mass is eventually accreted. For t → ∞ , all themass in the ring is accreted and the angular momentum is transported to infinitely largeradii by a negligibly small amount of mass.Another important class of solution [35, 36] is found in cases where the viscosity hasa simple power-law dependence on radius, ν ∝ R b . In this case, it is possible to finda self-similar analytical solution, where the initial density profile is such that accretionproceeds almost steadily out to a typical radius R at which the density is exponentiallytruncated. An example of the evolution in this case is shown in the right-hand panel ofFig. 2, for the case b = 1. Also in this case, as for the spreading ring, the disc spreadssignificantly outwards, even if most of the mass is eventually accreted. The transitionbetween inward and outward moving portion of the disc for the self-similar solutionsoccurs at a transition radius R tr ∼ R ( t/t ν ( R )) / (2 − b ) and therefore increases with time(for the particular case where b = 1 described above, we have that R tr increases linearlywith time). As the disc empties out, the accretion rate drops steadily. .7. Steady-state solutions . – In a steady state, the continuity equation takes the form:˙ M = − πRv R Σ , (45)where the constant ˙ M is the mass accretion rate and the signs have been chosen in sucha way that when v R is negative (i.e. directed inwards) the accretion rate is positive. ELF-GRAVITATING ACCRETION DISCS Fig. 2. – Left panel: evolution of the surface density of a Keplerian spreading ring, under theeffect of a constant viscosity ν . The lines refer (from top to bottom) to τ =0.01, 0.05, 0.1 and0.15 (see text for the definition of τ ). Right panel: evolution of a truncated power-law surfacedensity, in the case where ν ∝ R . Units are arbitrary. Analogously, angular momentum conservation becomes:˙ M Ω R − πν Σ | Ω ′ | R = ˙ J (46)where ˙ J is the constant net flux of angular momentum, and is determined by two con-tributions: the first term on the left hand side, which indicates the angular momentumadvected with the accretion process, and the second term, which indicates the outwardflux produced by viscous torques.˙ J is sometimes determined by applying the so-called “no torque” condition, accordingto which at the disc inner radius R in the angular velocity profile flattens (due to thepresence of a boundary layer where the disc connects to the central object) so that itsgradient and therefore the viscous torque vanish. For Keplerian rotation, this implies˙ J = ˙ M (Ω R ) in = ˙ M √ GM R in . Inserting this in equation (46), we obtain:3 πν Σ = ˙ M − r R in R ! . (47)At large radii, R ≫ R in , where the effects of a finite ˙ J are negligible, the surfacedensity and the viscosity satisfy the following simple relation:˙ M = (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) πν Σ , (48)that is, surface density and viscosity are inversely proportional to each other. An in-teresting consequence of the above relation arises when we use the α -prescription in the GIUSEPPE LODATO self-gravitating case, where H = c /πG Σ. Indeed, in this case Eq. (48) becomes:˙ M = 2 α c G (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) , (49)which shows that for a self-gravitating disc the accretion rate ˙ M is only dependent on α and on the sound speed c s . In particular, a self-gravitating disc with a constant α in asteady state is approximately isothermal. .8. Temperature profile and spectral energy distribution . – Let us conclude the analysisof the main properties of accretion discs by considering the energetics associated with theaccretion process. One of the most important consequences of accretion is the release ofgravitational potential energy as the accreting matter falls into the potential well. Thepower per unit surface D ( R ) dissipated by viscous stresses is given by: D ( R ) = ν Σ( R Ω ′ ) = (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) ν ΣΩ ≈ (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) ˙ M π Ω , (50)where the last equality holds approximately at large radii, where we can use Eq. (48). Thedissipation rate is therefore proportional to Ω , which can be a relatively steep decreasingfunction of radius. For example, if the potential is dominated by the central point massand the rotation curve is Keplerian, we have D ( R ) ∝ R − . In the other extreme conditionof a purely self-gravitating disc, with a flat rotation curve, the decrease is slightly moregentle, with D ( R ) ∝ R − . In fact, whenever the disc contributes to the rotation curve,it results in a gentler decrease of Ω at large radii, therefore producing a relatively largerenergy dissipation, especially in the outer disc.The above discussion provides a way to estimate to first approximation what is thespectral energy distribution (SED) produced by an actively accreting disc. Let us assumethat the energy dissipated through viscosity is released as blackbody emission at thedisc surface (which would occur if the disc is sufficiently opaque). Then the effectivetemperature T eff of the emitted spectrum is simply given by 2 σ B T , where σ B is Stefan-Boltzmann constant and where the factor 2 is introduced because the energy is emittedby the two sides of the disc. This then leads to T eff ( R ) = ˙ M πσ SB (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) Ω ! / = GM ˙ M πσ SB R ! / , (51)where the last equality holds for a Keplerian rotation curve. Each annulus of the discthus emits as a blackbody at a different temperature. The full spectrum of the disc isthen obtained by superimposing the different blackbody spectra from different annuli. If F λ is the flux emitted at wavelength λ , we thus have F λ = cos id Z R out R in πRB λ [ T s ( R )]d R, (52)where d is the distance between the disc and the observer, i is the inclination and B λ ( T )is the Planck function. The shape of the SED is a ‘stretched’ blackbody spectrum,extending over a much wider wavelength range than a simple Planck spectrum, from the ELF-GRAVITATING ACCRETION DISCS -2 0 2 4-4-202 Fig. 3. – SED of a Keplerian optically thick disc. At short wavelength the SED approachesthe Wien’s tail corresponding to the temperature in the inner disc T in . At large wavelengthsthe SED approaches the Rayleigh-Jeans tail corresponding to the outer disc temperature T out .At intermediate wavelengths the SED is a power law with index depending on the temperaturescaling with radius. For the Keplerian disc shown here the power law index is n = − /
3. Thetwo dashed lines show for comparison two pure blackbody spectra corresponding to T out and T in . Units are arbitrary. shortest wavelength λ min ≈ hc/kT in , associated to the temperature in the inner disc, tothe longest λ max ≈ hc/kT out , corresponding to the outer disc temperature. Fig. 3 showsa typical SED produced by a Keplerian disc, with the two dashed lines correspondingto blackbody spectra at T = T min and T = T max . The short wavelength part of theSED approaches the Wien’s tail appropriate to the inner disc temperature T in . Thelong wavelength SED approaches the Rayleigh-Jeans tail appropriate to the outer disctemperature T out . At intermediate wavelengths, it can be shown [37, 38] that, if thetemperature profile is a power law with respect to radius, with power law index − q ,then the SED λF λ is also a power law with respect to wavelength, with power law index n = 2 /q −
4. We thus get the typical Keplerian accretion disc spectrum, where q = 3 / n = − /
3. For a flat rotation curve disc, we have q = 1 /
2, which then leads GIUSEPPE LODATO to n = 0 and the SED is flat at intermediate wavelengths. This reflects the fact that fora flatter rotation curve the emission at large distances is enhanced, thus enhancing thelong-wavelength contribution to the SED.It is interesting to notice that in many cases the observed SED from accretion discs(for example, in the case of protostellar discs) is indeed relatively flat, indicating a flattemperature profile. However, in most cases, the disc mass is not high enough to producesubstantial deviations from Keplerian rotation in order to reproduce such flat SED,except possibly for some interesting cases (and in particular for the so called FU Orionisobjects), which will be discussed in more detail in Section .1 below. The observationalevidence for a relatively hot outer disc is generally interpreted in terms of an additionalsource of heating, coming for example from irradiation from the central star or from theinner disc. A full treatment of the energy balance in an accretion disc requires solvingthe radiative transfer equation through the disc, and therefore depends on the detailedoptical properties of the gas. This has been done under a variety of conditions in severalpapers [39, 40].Finally it is worth pointing out the relationship that is established between the viscos-ity parameter and the cooling time in thermal equilibrium. Indeed, in equilibrium (andassuming that there is no other source of heating for the disc) the heating rate providedby viscosity has to balance the cooling rate. The cooling rate per unit surface can besimply written as C = Σ c γ ( γ − τ cool , (53)where γ is the ratio of specific heats and where we have introduced a cooling timescale τ cool . Equating this expression to the heating rate due to viscosity (Eq. (50)) and usingthe α -prescription we get α = (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) − γ ( γ − τ cool . (54)
3. – GRAVITATIONAL INSTABILITIES AND SELF-REGULATION
In the sections above, we have described the equilibrium configuration of accretiondiscs and their secular evolution due to viscosity. We have described what is the effect ofthe disc self-gravity on such equilibrium state. In particular, we have seen that for smallbut non negligible total disc masses M disc /M ≈ H/R , self-gravity affects the disc verticalstructure, while the rotation curve is only affected by self-gravity when the disc massbecomes comparable with M . However, one of the most important effects associated withself-gravity is the possibility of developing gravitational instabilities. The physical originof the instability is essentially related to the standard Jeans instability for a homogeneousfluid, where large scale disturbances cannot be stabilized by pressure gradients, thusdetermining a typical wavelength (the so-called Jeans wavelength) associated with theinstability. In a rotating disc, as I will show here below, the additional effect of rotationcan also stabilize these large scale perturbations.As I have already remarked above, gravitational instabilities might be particularlyimportant in providing a significant source of transport of angular momentum throughthe disc. Under certain conditions, the non-linear outcome of the instability can be rela-tively violent, leading to the fragmentation of the disc into gravitationally bound clumps. ELF-GRAVITATING ACCRETION DISCS Most of the recent research on self-gravitating discs has concentrated on exploring theconditions under which the instability takes the form of a long-lived spiral structure,leading to a steady transport of angular momentum, or alternatively determines the discfragmentation. In this section and in the following ones, I discuss in detail such issues. .1. Development of gravitational instabilities . – The dynamical processes underlyingthe development and the maintenance of gravitational instabilities in discs have beendiscussed in several papers and reviews, mostly in view of an application to the discsof spiral galaxies [12, 25]. On the other hand, many processes relevant to the dynamicsof stellar discs are also applicable to the case of a fluid disc. The starting point ofsuch analysis is the determination of the propagation properties of density waves ina self-gravitating disc, through a standard linear perturbation analysis of the relevantdynamical equations, i.e. the continuity equation (Eq. (2)), the momentum equationfor an inviscid disc (Eq. (3)), supplemented by the equation of state (Eq. (7)), torelate pressure and density, and by Poisson’s equation (Eq. (4)) in order to connectthe disc density to the gravitational potential. The difficulty here arises from the long-range nature of gravitational force, so that the gravitational potential at a given locationdepends on the whole distribution of the gas density Σ, rather than on its local value.A local analysis can nonetheless be carried out in the WKB approximation, for tightlywound perturbations, for which the pitch angle m/ ( kR ) ≪
1, where k is the radialwavenumber and m is the azimuthal wavenumber (corresponding to the number of armsof a spiral disturbance) . In this approximantion, the dispersion relation is:( ω − m Ω) = c k − πG Σ | k | + κ , (55)where ω is the frequency of the perturbation and the epicyclic frequency κ is given by κ = 2Ω R d(Ω R )d R . (56)Note that the epicyclic frequency κ is generally of the order of the angular velocity Ω,the exact proportionality depending on the rotation curve. For example, for Keplerianrotation, we have κ = Ω, while for a flat rotation curve, we have κ = √ m = 0. Clearly, whenever ω is positive,a given perturbation simply propagates in a wave-like way. On the other hand, if ω isnegative, an exponentially growing instability arises. We thus see from the right-hand sideof Eq. (55) that the second term, which is negative and is associated with the disc self-gravity, is a potentially destabilizing term, while the first, associated with pressure forces,and the third, related to rotation, are stabilizing terms. In particular, while pressurestabilizes the disc at small wavelengths (large k ), rotation is more effective at stabilizinglarge scale disturbances (small k ). Equation (55) is a simple quadratic expression forthe wavelength k and it can be easily seen that for axisymmetric perturbations the righthand side is positive definite (i.e. ω is positive and the perturbation is stable) if: Q = c s κπG Σ > . (57)The dimensionless parameter Q is therefore essential in determining the local axisymmet-ric stability of self-gravitating accretion discs (a similar parameter also determines the GIUSEPPE LODATO stability properties of stellar discs, [24]). We have already encountered a parameter verysimilar to Q , when discussing the relative importance of the central object and of thedisc in determining the disc thickness (Eq. 27). We had seen that when c s Ω /πG Σ ≈ m does not enter explicitly in the dispersion relation (except inthe Doppler shifted frequency) so that both axisymmetric perturbations ( m = 0) andnon-axisymmetric ones ( m = 0) share the same instability properties.For discs that are gravitationally unstable, that is in cases where Q . k for which the right-hand side of Eq. (55) attains a minimum. This isgiven by k = πG Σ c = 1 H sg . (58)We thus see that the most unstable wavelength is of the order of the disc thickness. Thecondition for the WKB approximation to hold is then kR = R/H sg ≫ M disc ≪ M (cf. Eq. (26)).It is well known (for example, from early N -body numerical simulations [41, 42])that massive discs are subject to violent large scale (long-wavelength) non-axisymmetricinstabilities, even in cases where the parameter Q is above unity, and therefore the disc ispredicted to be stable in the WKB approximation. Indeed, a local dispersion relation canalso be obtained in the case of relatively long radial wavelengths in the non-axisymmetriccase [43]. This relation is more complicated than the simple quadratic form describedabove and depends on a new parameter, in addition to Q , defined as: J = m πG Σ Rκ κ (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) / . (59)Being proportional to Σ, J is a measure of how massive the disc is. Indeed, to give anidea of the numbers involved, consider a quasi-Keplerian disc, where κ ≈ Ω and whereΣ ∝ R − . In this case we have J = √ mM disc ( R ) /M . If J ≪ m ), thedispersion relation reduces to the standard quadratic form. On the other hand, for J ≈ m disturbances can be violently unstable. Interestingly, asI will discuss below, this change of behaviour between light and massive discs has beenobserved numerically in a recent survey of the development of gravitational instabilitiesin fluid discs [45, 46].To summarize, in the tightly wound approximation we expect a gaseous disc to be-come gravitationally unstable when Q <
1. Non-axisymmetric disturbances, with largeazimuthal wavenumber m can be unstable even for relatively larger value of Q , eventhough if Q rises above ≈
2, then even non-axisymmetric disturbances are stabilized.Finally, it is worth noting that taking into account finite-thickness effects in the tightlywound approximation leads to a stabilization of the disc, so that the marginal stabilityvalue of Q is reduced below unity [47]. ELF-GRAVITATING ACCRETION DISCS .2. Self-regulation . – Having established that massive discs are linearly unstable, weshould now ask what is the non-linear evolution of the instability and whether it can leadto a sustained transport of angular momentum. The stability criterion, Eq. (57) offersa natural way to describe this. Indeed, we should note that the stability parameter Q is proportional to the sound speed (and hence to temperature), so that colder discs aremore unstable. Now, let us consider a disc which is initially hot, so that Q ≫ Q ≈
1. At this stage, the disc developsa gravitational instability in the form of a spiral structure. Compression and shocksinduced by the instability lead to an efficient energy dissipation and as a result the discwill heat up. In turn, the stability condition works like a kind of ‘thermostat’, so thatheating turns on only if the disc is colder than a given temperature. If the ‘thermostat’works, we should then expect the instability to self-regulate in such a way that the discis always kept close to marginal stability, and thus a self-gravitating disc will evolve to astate where Q ≈ ¯ Q , where ¯ Q is a constant of order unity. Here again, it is reassuring tosee that numerical simulations [45] reproduce closely this behaviour.Note here an important difference between a stellar and a fluid disc. In order forthe above self-regulation process to work, the presence of a dissipative component isessential. While this is obviously the case for a gaseous disc, a pure stellar disc lacks thisdissipation channel. For a galaxy disc, therefore, self-regulation requires the presence ofa gaseous component in addition to the dominant stellar one [48]. .3. Self-regulated accretion discs models . – The self-regulation mechanism describedabove can be easily introduced in steady state models of self-gravitating accretion discs[15]. Indeed, for a self-gravitating disc, as shown in Section .7 at large radii we have:˙ M = 2 αc G (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) , (60)where we have also used the α -prescription. If we supplement this equation by theself-regulation condition: Q = c s κπG Σ = ¯ Q ≈ , (61)we thus easily see that in the case of a disc dominated potential (that is, if the centralstar mass is negligible) a Mestel disc with Σ ∝ R − provides an analytical self-similarsolution to the model. In particular, we have: c s = G ˙ M α ! / , (62) v c = 2 √ Q G ˙ M α ! / , (63) 2 πG Σ R = 8¯ Q G ˙ M α ! / . (64) GIUSEPPE LODATO
The solution above has remarkably simple properties, with a flat profile of v c , c s and Q and where the surface density has a simple power-law dependence on radius. However,it is clearly approximate and only applies to large radii, where the effects of ˙ J andespecially of the presence of a central point mass M can be neglected. In particular, inrealistic accretion discs the latter condition is difficult to achieve and we should thereforeextend the analysis to include a central mass in the model. In particular, we expect thetransition from an inner Keplerian disc to an outer disc, where the rotation curve wouldtend to flatten, to occur at a distance of the order of R s = 2 GM (cid:18) ¯ Q (cid:19) G ˙ M α ! − / . (65)In this way we can use the simple model above to check the influence of the disc self-gravity on the rotation curve of some observed systems. Let us then estimate the tran-sition radius R s for a number of different cases. We first consider a typical protostellardisc, for which ˙ M ≈ − M ⊙ / yr, the central star mass is of the order of 1 M ⊙ and wherewe assume α = 0 .
01. Inserting these numbers in Eq. (65), we get R s ≈ .1 for more details) the accretion rate can reach ˙ M ≈ − M ⊙ / yr. Insuch cases, R s is reduced to ≈ −
20 AU, so that the outer parts of the disc mightshown some deviation from Keplerian rotation. Finally, let us consider a typical disc inan Active Galactic Nucleus, for which ˙ M ≈ . M ⊙ / yr and M ≈ M ⊙ . In this case,we get R s ≈ below. Additionally, it is worth mentioning that numerical simulations whichtry to reproduce the formation of protostellar discs [50], tend to be characterized byfairly massive discs in their earliest stages, with properties that resemble closely the onesdescribed semi-analytically by these models.Self-regulated disc models also tend to be relatively hotter in their outer regions, com-pared to their non-self-gravitating counterparts. This is due both to the self-regulationmechanism, and the effective heating associated with it, and by the fact that the dischas a flatter rotation curve in its outer parts, as discussed in Section .8. As a con-sequence, self-regulated SEDs tend to be flatter than the standard disc SED, shown inFig. 3. and often display a secondary “bump” at long wavelengths [49]. An exampleSED is shown in Fig. 4, where the different curves refer to different choices of the discouter radius (see [49] for details). These models were originally proposed in order toexplain some observations of protostellar discs, and will be discussed more in detail insection .1 below. Interestingly, very similar SEDs are also produced in numerical sim- ELF-GRAVITATING ACCRETION DISCS Fig. 4. – Left: examples of rotation curves for self-regulated disc models in the presence ofa central point mass. In the inner region the rotation curve is Keplerian, while it becomesprogressively flatter in the outer disc, asymptotically reaching the self-similar configuration.More details can be found in [13]. Right: example SEDs for self-regulated discs. At largeradii self-regulated discs are hotter than their non-self-gravitating counterparts, resulting in asecondary bump in the SED at long wavelengths [49], which becomes more pronounced if thedisc outer radius is increased. ulations of self-gravitating protostellar discs [51], and similar semi-analytic models forself-gravitating disc SEDs have later been proposed also for AGN discs [52].
4. – NUMERICAL SIMULATIONS OF SELF-GRAVITATING ACCRETIONDISCS
In the sections above we have discussed the development of gravitational instabilitiesessentially in terms of linear perturbation analysis. We have already noted that the non-linear evolution of such instabilities can result in a feedback loop which self-regulates thedisc and keeps it close to marginal stability. It is then possible, as we have seen, to developanalytical models where such complex behaviour is introduced in a simple way in the basicsteady state disc equations. Clearly, in order to check whether self-regulation is indeedestablished once the instability becomes non-linear, and to investigate secular effects,like the transport of angular momentum induced by the spiral structure, it is importantto perform numerical simulations of the disc well into the non-linear regime. Manysimulations of this kind were performed in the early ’70s in the context of collisionlessstellar systems by means of N -body techniques [41, 42, 53, 54].Here, we are interested in the description of a dissipative, gaseous disc and thus theseearly N -body simulations were not be able to catch some aspects of the disc hydrodynam-ics. In the last fifteen years complex, three-dimensional hydrodynamic simulations havebecome an essential tool to study accretion discs dynamics and have provided an impor-tant ‘experimental’ counterpart to analytical investigations, such as those described inthe previous sections. In this respect, a set of simulations should be considered as a ‘nu- GIUSEPPE LODATO merical experiment’ where our models are checked against measurements of a controlledsystem, often unavailable from astronomical observations. Simulations have evolved sig-nificantly throughout the years, partly due to the increased computational power (whichhas enable to reach higher resolution and accuracy), and partly to the development ofprogressively more refined algorithms to include new physical processes in the models.A key role in the evolution of self-gravitating accretion discs is played by the coolingprocess. As discussed above, it is the interplay between cooling and internal heatingdue to gravitational instability that allows to reach a self-regulated condition. As will bediscussed in more detail in section below, the cooling rate also determines the efficiencywith which the gravitational instability redistributes angular momentum and allows ac-cretion. Most of the effort in this area of research is therefore aimed at introducinginto the codes progressively more realistic implementations of the cooling processes inaccretion discs. In this section we provide an overview of the development of numeri-cal simulations of accretion discs in the last few years. Specific results concerning thetransport induced by gravitational instability and the possibility of disc fragmentationare discussed in Sections and below.Some early attempts at simulating self-gravitating accretion discs were made by usingcollisionless N -body methods [55]. However, the earliest hydrodynamical simulations ofself-gravitating accretion discs were performed later [56] in the context of the forma-tion of protostellar discs. Such simulations employed a very simplified thermodynamicaldescription of the gas, that was considered to be ‘locally isothermal’ in the sense thatevery parcel of fluid was assumed to keep its internal energy throughout the simulation.Several other early simulations were performed using different equations of state, such asa polytropic equation [57,58] and were used to describe the development of gravitationalinstability, but the effects of the choice of thermal state of the gas were generally over-looked. More emphasis on thermodynamical aspects was placed later [59, 60], when acomparison between models evolved isentropically and models which enforced an isother-mal behaviour revealed that gravitational instabilities developed much more violently inan isothermal configuration, leading to a very fast redistribution of angular momentumand in extreme cases to the disruption of the disc. A posteriori, this behaviour is simplyunderstood based on the arguments described in Section .3 above. The assumptionof isothermality prevents the development of the instability to feed back thermal en-ergy in the equilibrium state, therefore not allowing the non-linear stabilization of thedisturbance and the saturation of the gravitationally unstable modes.All the simulations described so far start from an initial configuration which is ex-pected to be highly unstable, with Q values of the order of unity or below. It is thereforenot surprising that as soon as the simulation is set to evolve, strong disturbances sud-denly appear and grow out of control in the case of isothermal simulations. The nextquestion is therefore what would be the response of a disc which is initially stable andis driven to an unstable configuration through cooling. Clearly, in order to consider thiscase, one has to move away from either an isothermal or a polytropic equation of stateand consider instead the full evolution of the energy equation. This was initially doneby considering a simple cooling prescription [45,46,61-68] where the cooling rate is givenby: d u d t (cid:12)(cid:12)(cid:12)(cid:12) cool = − uτ cool , (66)where u is the internal energy of the fluid and τ cool is the cooling timescale, which is set as ELF-GRAVITATING ACCRETION DISCS a free parameter in the simulations. It is important to stress that the above descriptionis not meant to reproduce any specific cooling law, but is just a convenient way (a ‘toy’model) of exploring the role of the cooling timescale in the outcome of the gravitationalinstability. Different investigators use different prescriptions for the specific form of thecooling time, which in some cases is taken to be a constant while in other cases is takento be proportional to the dynamical timescale, so that Ω τ cool is kept constant. Despitethe different prescriptions for the cooling time and the different type of codes used (localshearing sheet models [61], global Eulerian grid-based models [66, 67], global Lagrangianparticle-based models [45,46,62-65,68]) the results appear to converge into an essentiallycoherent picture, at least for low mass discs. It is now well established that the longterm behaviour depends on the relative values of the cooling and dynamical timescale.If the cooling timescale is larger than a few dynamical timescales, an initially stable(large Q ) disc cools down until Q becomes of the order of unity. At this stage, thedisc becomes gravitationally unstable and develops a spiral structure which provides aheating source, through compressional heating and shock dissipation, able to balance theexternally imposed cooling. Once in thermal equilibrium, the disc is characterized byan approximately constant value of Q very close to marginal stability. In such a statea spiral structure persists in the disc, to provide the required heating. Therefore theself-regulation mechanism described above determines the disc structure and evolution.Figure 5 shows the result of one such simulations, where in this case τ cool Ω = 7 . M disc = 0 . M [45]. In the left panel the colour plot shows the disc surfacedensity, in which a spiral structure is clearly seen. The right panel shows the azimuthallyand vertically averaged value of Q as a function of radius. The disc in this case extendsfrom R = 0 .
25 to R = 25 in code units. It is then seen that far from the boundaries(where the density drops and Q correspondingly grows) the disc is self-regulated, with Q ≈ M disc /M = 0 .
05, 0.1 and 0.25. This allows a comparison of the morphology andof the dynamics of the spiral structure developed in the different cases. The upper panel ofFig. 6 shows a snapshot of the density structure three simulations once they have reacheda self-regulated state, while the lower panel shows the power associated with Fouriercomponents of different m . Several interesting features can be noted. First of all, notethat the typical radial wavelength of the instability increases as the disc mass increases.This can be easily understood even from the simple quadratic dispersion relation, Eq.(55), which predicts that the typical wavelength of the instability is proportional to thedisc thickness H , Eq. (58). Now, as we know, self-regulated discs have H/R ≈ M disc /M and therefore a more massive disc is also thicker and the resulting spiral structure tendsto be more open. Also, note that while for the low mass case the different Fouriermodes have similar amplitudes, as the disc mass increases, modes with lower m tend tobecome more important. Also this behaviour can be understood in terms of the localdispersion relations described above. Indeed, when the disc mass is low, the parameter J ≈ mM disc /M is also small (except for very large m which are not easily excited)and we can use the quadratic dispersion relation which predicts a similar behaviour forall m . On the other hand, for larger disc masses, J can be of order unity already forsmall m which then tend to be much more unstable, as confirmed in the simulations.If we increase the mass ratio further, going up to M disc ∼ M [46] this trend continuesand strong m = 2 disturbances become violently unstable. In this high mass case,self-regulation cannot be established and the disc shows a recurrent pattern of highly GIUSEPPE LODATO
Fig. 5. – Numerical simulation of the development of a gravitational instability in a cooled discwith τ cool Ω = 7 . M disc = 0 . M [45]. Left: Surface density of thedisc, showing the prominent spiral structure. Right: radial profile of the azimuthal and verticalaverage of Q . The disc is self-regulated over a wide radial range. The three curves refer to threedifferent times during the evolution, showing that the disc has reached thermal equilibrium. temporally variable instabilities, as shown in Fig. 7. During such episodes, the torquesinduced by self-gravity can be very large, and cause a relatively fast redistribution of gaswithin the disc.The behaviour described above changes when the cooling time is decreased to smallervalues [61]. In this case, the disc does not reach a quasi-steady self-regulated state butrather fragments into several bound objects. Figure 8 show the results of a simulationvery similar to the one displayed in Fig. 5, but where the cooling time is decreased to τ cool = 3Ω − [65]. This result can be understood in the following way, by adopting alocal approach to describe the instability. In a gravitationally unstable disc, the typicalgrowth timescale of unstable perturbations is of the order of the dynamical timescaleΩ − . The non-linear stabilization of the perturbation only works if the heat generatedby compression and shocks is not removed efficiently from the disc through cooling. Sincethe perturbation grows on the dynamical timescale, if we want to avoid fragmentation, werequire that cooling acts on a longer timescale. Note that the requirement that the coolingtimescale be shorter than the dynamical timescale in order to result in fragmentationhas been known for several years, even outside the context of disc instability [69, 70].The exact value of the threshold for fragmentation does depend slightly on the specificnumerical set-up and ranges from τ cool = 3Ω − to τ cool = 6Ω − [61,65,71]. An additionaldiscussion on the dependence of the threshold on the ratio of specific heats and on itsinterpretation is provided in Section below.More recent simulations have further improved the thermodynamical description ofthe disc, including a more realistic cooling prescription based on the opacity propertiesof the gas. In the simplest case, one can implement such ‘realistic’ cooling in local ELF-GRAVITATING ACCRETION DISCS Fig. 6. – Top panels: Surface density structure for a set of self-gravitating disc simulations with τ cool Ω = 7 .
5. From left to right the mass ratio is M disc /M = 0 .
05, 0.1 and 0.25. Lower panels:corresponding mode analysis. The amplitude of various modes with different m is shown. simulations of infinitesimally thin discs, where the vertical radiative transfer is simplytreated with a one-zone approximation [72]. More realistic, three-dimensional simulationsincluding radiative transfer are being developed in the last few years [73] and start toprovide their first results [51, 74, 75].
5. – TRANSPORT INDUCED BY GRAVITATIONAL INSTABILITIES
The discussion in the last two sessions has shown that (a) accretion discs can besubject to gravitational instabilities whenever the parameter Q (which describes linearstability in the tightly wound approximation) is of order unity and (b) under appropriateconditions the non-linear evolution of the instability tends to self-regulate in such a waythat the disc is kept very close to marginal stability ( Q ≈
1) over a wide radial range.We have also seen that, in the limit of Keplerian rotation, the condition Q ≈ M disc /M ≈ H/R , so that if the thickness is small, even a relatively lightdisc can be gravitationally unstable.We now turn to the question of the transport of energy and of angular momentuminduced by the instability. In a standard viscous accretion disc, both energy and angularmomentum are transported through the disc by viscosity, and are therefore described bya single transport coefficient. As we will see, this might not be the case for self-gravitatingaccretion discs, and we will then have to determine if and under which conditions can GIUSEPPE LODATO
Fig. 7. – Evolution of the surface density of the disc during the development of a transientepisode of spiral activity in a simulation of a massive M disc = 0 . M accretion disc. The foursnapshots refer to four different times during the simulation, t = 0 (upper left), t = 1 . t = 2 . t = 3 (lower right), where the times are given in units of the orbitalperiod at the disc outer edge. From [46]. such a self-gravitating disc be described in terms of a single ‘viscosity’ coefficient. That aspiral density wave transports angular momentum was early recognized in the context ofstellar dynamics [76, 77]. In the context of accretion discs, where most of the theoreticaldevelopment has been achieved through the use of an ad hoc (cf. the α -prescription)anomalous viscosity, the key question is understanding to what extent can this form oftransport be described in terms of a ‘viscous’ process. In other words, the α -prescriptionfor viscosity has always been used with the (often tacit) assumption that it actuallycorresponds to the transport induced by some kind of disc instability. It is therefore ELF-GRAVITATING ACCRETION DISCS Fig. 8. – Numerical simulation of a self-gravitating disc with M disc = 0 . M and τ cool = 3Ω − .Once unstable, the disc breaks up into numerous gravitationally bound clumps. natural to ask whether the development of gravitational instability can indeed be thislong sought cause of angular momentum transport.A possible way to proceed is to simply assume that the transport of angular momen-tum takes the same form of a viscous diffusion and to introduce a Q -dependent viscosityparameter α [9, 10], such that α sg = ξ (cid:18) ¯ Q Q − (cid:19) Q < ¯ Q Q > ¯ Q (67)Here ¯ Q is the value of Q at which the disc becomes unstable to non axisymmetric per-turbations and ξ is a parameter to measure the strength of the induced torques. Theabove formulation is useful in practical cases, for example, when one wants to incorpo-rate in a simple way the self-regulation mechanism in simple time-dependent models ofself-gravitating discs. However, it lacks one important feature elucidated from the nu- GIUSEPPE LODATO merical simulations described above. In this picture, α sg only depends on the local valueof Q and not on the cooling timescale τ cool , that we have seen controls so efficiently thedevelopment of the instability. In particular, for self-regulated discs, we expect Q ≈ ¯ Q and the formula above would then produce a negligibly small α sg , while we know thata finite amplitude spiral structure is present in self-regulated discs and indeed it is thisspiral structure that provides the heating to balance the imposed cooling rate.If we maintain the assumption that the transport induced by self-gravity can be de-scribed in terms of viscosity, a prescription for α sg consistent with the apparent behaviourof the numerical simulations described above can be obtained by the following argument,which is inspired from the self-regulation mechanism. We have seen that the behaviourof the disc is strongly dependent on the cooling timescale. In particular, for a givencooling timescale, the disc eventually reaches a steady thermal equilibrium state, wherethe amplitude of the spiral instability saturates at a value such that heating balances theimposed cooling. If we now decrease the cooling time, the disc will cool down and becomemore unstable, and the saturation amplitude of the instability will grow, thus providingan enhanced transport and dissipation, to match the new cooling configuration. A simpleprescription can then be provided from the requirement of thermal equilibrium (cf. Eq.(54)) α sg = (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) − γ ( γ − τ cool . (68)Note that the above equation is simply a restatement of the energy balance equation,where the viscous heating term on the left-hand side is balanced by the cooling term onthe right-hand side. A relation like the one described above necessarily holds wheneverself-gravity provides the necessary heating to reach thermal equilibrium, and if transportcan be described viscously. The practical inconvenient of Eq. (68) is due to the fact thatin general the cooling time in a given configuration is not known a priori , and it requiresa careful evaluation of the cooling processes in order to determine the amount of stressin a specific case.We still need to verify whether the condition for the validity of Eq. (68) describedabove (i.e. that the transport induced by self-gravity is viscous) does hold. It is herethat numerical simulations can provide the missing piece of information. Indeed, fromthe results of the simulation one is able to directly measure the relevant torques andstresses and compare the results with the models described above. The ‘viscosity’ can beestimated from the simulations in two different ways. First, one can compare the secularevolution of the surface density as seen in the simulation with the one predicted forexample from Eq. (35) and fit the value of the viscosity to the observed evolution. Thishas been done [56, 57] for isothermal or polytropic equations of state and the evolutionappears to be reproduced reasonably well with values of α sg in the range of 0.01-0.03.It should be noted that these simulations had relatively massive discs, with M disc ≈ M .However, such estimates cannot be directly compared to what should be expected basedon Eq. (68), since these simulations, being isothermal, do not have a well defined coolingtime. In some cases [57], the simulations require the adoption of a spatially varying α sg , prompting the authors to speculate that this is a consequence of a ‘failing’ of the α -prescription due to the global nature of transport in their models. Now, while globaltransport might well operate in the massive ( M disc ≈ M ) discs considered in thosesimulations, we note here that, as already remarked in Section .5 above, the simple factthat α is not constant does not call into question neither the use of the α -prescription, ELF-GRAVITATING ACCRETION DISCS Fig. 9. – Upper panel: gravitational (solid line) and Reynolds (dotted line) contributions to thestress as a function of radius for the simulation shown in Fig. 5. Lower panel: comparisonbetween the total (Reynolds plus gravitational) stress (solid line) and the expected value fromEq. (68). nor the use of a local viscosity or a diffusion equation to describe the disc evolution.A second possibility to estimate viscosity is to compute directly the torque exertedby the non axisymmetric disturbances developing in self-gravitating discs. As we willsee, this will lead to more insight on the issue of whether transport can be describedin terms of local quantities or rather global effects invalidate such a treatment. Thestarting point for this analysis is Euler equation for a dissipationless fluid, Eq. (3).Consider the case where the fluid velocity is made up of two contributions, one from a(secularly varying) average component, that we call v , and one from a rapidly varyingfluctuating component, that we call u . In a similar way, the self-gravitating potential canbe divided in a quasi-steady component and a fluctuating one. It can be shown (see, forexample, [31,78]) that in cylindrical coordinates, the vertically intergated Euler equationtakes the following form: ∂∂t (Σ Rv φ ) + 1 R ∂∂R ( Rv R Σ Rv φ ) = − X i R ∂∂R (cid:16) R Σ h u ( i ) R u ( i ) φ i (cid:17) . (69) GIUSEPPE LODATO
In Eq. (69) the angle brackets indicate a vertical and azimuthal average and the sum-mation extends to the various fluctuating fields that contribute to the stress. In the caseconsidered here, the two relevant fields are the velocity field and the gravitational field.The stress associated with the fluctuating velocity is usually called the “Reynolds” stressand is given by T Re Rφ = − Σ h u R u φ i , (70)while the term associated with self-gravity is given by T g Rφ = − Σ h u g R u g φ i , (71)where u g = g / √ πGρ , and g = − ∇ Φ is the gravitational field [77]. Equation (69)shows that in the presence of a non-zero correlation between the radial and azimuthalcomponents of any of the fluctuating fields described above, angular momentum can beremoved from the mean flow. Indeed, comparing Eq. (69) to Eq. (33), it is readilyseen that such correlations play exactly the same role as a viscous stress tensor for whatconcerns angular momentum conservation. Of course, we can also use the α -prescription,that is, we can measure the stress in units of the local vertical integrated pressure: α sg = (cid:12)(cid:12)(cid:12)(cid:12) d ln Ωd ln R (cid:12)(cid:12)(cid:12)(cid:12) − h u R u φ i + h u g R u g φ i c . (72)The above formula is simply a measure of the torque provided by non-axisymmetricstructures in units of the local pressure and obviously it can always be used to charac-terized the angular momentum transport due to self-gravity. The question of whetherthis can be interpreted as a ‘local’ viscosity or not requires a consideration of two furtherissues: (a) Is the value of α sg determined locally or does it depend on the global discconfiguration? (b) Does this process dissipate energy as a viscous term would do, thatis according to Eq. (50)?Both questions above can be answered using the results of numerical simulations.Indeed, we can simply measure, based on the simulations, α sg from Eq. (72). As wehave seen above, in many cases the disc reaches a self-regulated state in which thermalequilibrium holds. We can then compare the torque computed from Eq. (72) with theexpectations based on Eq. (68), that is derived under the assumption that the coolingterm is balanced by local viscous heating. This exercise has been done in a number ofpapers [45,51,61] and the general results is that the two expressions agree relatively wellover a large radial range. Note that the simulations reported in [61] were set up to belocal (i.e. they only described a small patch in the disc) and there is therefore no surprisethat possible global effects did not show up. On the other hand, the simulations of [45]and those of [51] are global models of the disc (which in these two cases is taken to be ofrelatively low mass M disc /M ≪
1) and global effects could in principle be present. Fig.9 shows an example of such calculation, where the stress computed from Eq. (72), (solidline in the lower panel) is compared to the expectation from Eq. (68) (dotted line) forthe simulation shown in Fig. 5. An even clearer example is shown in Fig. 13 of [51],where the same comparison is made, again showing a remarkable agreement over a largeportion of the disc. Note that the simulations of [45] and those of [51] employ significantlydifferent cooling prescriptions, in that while in [45] τ cool Ω is set to be a constant, resultingin an approximately constant value of α sg (cf. Eq. (68)), in [51] τ cool is either constant ELF-GRAVITATING ACCRETION DISCS or self-consistently computed based on the solution of the radiative transfer equation,therefore resulting in a non constant value of α sg . In all cases, however, energy appearsto be dissipated according to the expectations from a viscous model. Clearly, such acomparison is only meaningful when the disc is in thermal equilibrium. As discussedabove, for very massive discs [46], where M disc . M , a quasi-steady thermal equilibriumis generally not reached and this kind of comparison cannot be made.It is also possible to estimate whether the value of α sg is determined locally or not. Forexample one can calculate the contribution to the gravitational stress at a given referencelocation coming from a region of finite size around it [45]. It can then be seen that thedominant contribution to the stress arise from a neighbouring region of size ≈ H . Asimilar result has also been obtained using a different method for local simulations [61].The local approximation thus breaks down if 10 H & R , or H/R ≈ M disc /M & .
1. Noteonce again the different behaviour of light and massive discs in this respect.We have thus seen that, at least in the case of low mass discs, the transport of energyand angular momentum induced by gravitational instabilities can be well reproduced interms of a local, ‘viscous’ model. Is this result to be expected? Low mass discs can sup-port a relatively large spectrum of density waves, with different m and with wavenumbersof the order of 1 /H . In the presence of effective cooling, many of such waves can be easilyexcited. One possibility is that such waves interact with one another, undergo shocksand dissipation, and therefore do not propagate over large radial distances, especiallyif their typical wavelength is small. In such a case, a local description of the transportinduced by the waves would appear to be appropriate. Alternatively, another possibleoutcome is that a relatively small number of waves can propagate significantly and es-tablish a pattern of global modes, which persist over a few orbits and extend over a largeradial range. Such persistent, global modes have been observed in some hydrodynamicalsimulations of accretion discs, even in the case of low mass discs [51], but even in thiscase the calculated stress tensor agrees with the value required by energy balance in aviscous disc.The energy flux transported by self-gravitating density waves can also be calculatedanalytically [76]. It can be shown [79] that in general the energy flux (contrary to whathappens for a viscous process) is not proportional to the stress tensor, and that someother, ‘anomalous’ terms are present, which would preclude a local treatment. Theseterms can be evaluated within a WKB analysis and it is shown [79] that they vanishat corotation (where the pattern speed of the wave Ω P = Ω( R )), their contributionincreasing as we move away from corotation and becoming significant if | Ω P − Ω( R ) | ≈ Ω( R ). Let us consider the case of an approximately Keplerian disc. In this case, for awave generated close to corotation such global terms become important only if the wavehas travelled a radial distance of the order of ∆ R/R ≈ / ≈ .
58. Now, it is plausibleto assume that the minimum distance travelled by a density wave is no less than its radialwavelength ≈ πH . A necessary (but not sufficient!) condition for global effects to benegligible is thus H/R . ∆ R/ πR ≈ .
25. Once again, we see that thin and thick discs(or, equivalently, light and massive discs) behave in a significantly different way. Whilefor thick and massive discs a local approach appears inadequate, in the case of lighterand thinner discs it is possible to have either a local or a global behaviour, dependingon the coherence of the spiral modes involved. The numerical evidences described aboveappear to imply that a local approach is generally valid for thin discs, and indeed evenin the cases where a global pattern is observed, its radial extent is not large enough toprevent the use of a local model [51].Let us summarize what we have discussed so far. The behaviour of self-gravitating GIUSEPPE LODATO discs depends on three dimensionless parameters: the parameter Q (which measures howhot the disc is), the parameter M disc /M (which measures how massive the disc is) andthe parameter Ω τ cool (which measures how fast does the disc cool). If Q is significantlylarger than unity the disc is gravitationally stable, while if Q approaches unity it becomesunstable. For an approximately Keplerian rotation the condition Q ≈ M disc /M ≈ H/R . The behaviour of a gravitationally unstable disc depends on M disc /M and Ω τ cool . In particular, if Ω τ cool . M disc /M . .
25, thedisc settles down in a quasi-steady self-regulated state, the instability takes the form ofa tightly wound spiral structure and the transport of energy and angular momentuminduced by it are usually well described within a local ‘ α -like’ approach. On the otherhand more massive discs, with M disc /M & .
25, cannot be described within the localapproximation and global transport is expected to occur. In this case, a quasi-steadyself-regulated state is not generally observed in the simulations (see Fig. 7), which showinstead a series of recurrent strong episodes of instability, during which the gravitationaltorques can be very large ( α sg ≈ .
6. – THE RELATION BETWEEN FRAGMENTATION AND TRANSPORT
The above discussion on the transport induced by self-gravity in a light disc (i.e. forthe case where a local description in terms of an α viscosity is valid) allows us to givea more detailed look at the issue of fragmentation. Let us then concentrate on light,quasi Keplerian discs. If these discs are self-gravitating, they settle down in a quasi-steady self-regulated state, with Q ≈ ¯ Q , in which they transport angular momentumwith a torque strength that can be expressed in terms of α through Eq. (68). However,from the discussion of Section we know that if the cooling rate is smaller than agiven threshold τ cool Ω < β crit ≈ −
6, the disc undergoes fragmentation. Therefore,through Eq. (68), a minimum value of Ω τ cool to prevent fragmentation corresponds toa maximum value of α sg . It is then possible to interpret the fragmentation threshold ina different way, in terms of how much angular momentum can be transported throughself-gravity in a steady state. This interpretation is further reinforced by numericalsimulations which employ different values of the ratio of specific heats γ [65]. It appearsthat as γ decreases, discs are more susceptible to fragmentation. For example, while for γ = 5 /
3, the threshold for fragmentation is at β crit ≈ .
5, it increases to β crit ≈ . γ = 7 /
5. The results of this analysis is shown in Fig. 10, where the three curves show therelation between gravitationally induced viscosity α and the cooling time τ cool in thermalequilibrium, expressed in Eq. (68), for three different values of γ = 2 (solid line), γ = 5 / γ = 7 / τ cool at which fragmentation is observed in [65], while the triangle refers to the 2Dsimulations with γ = 2 presented in [61]. It is then interesting to note that fragmentationoccurs whenever the value of α needed to be provided by self-gravity in order to reachthermal balance exceeds a maximum value of order α max ≈ .
06. For relatively longcooling times, the strength of the torques induced by self-gravity does not need to bevery large in order to provide the required heating to reach equilibrium. For shortercooling times, the saturation amplitude of the spiral structure increases until it becomestoo large to be sustained by the disc, which then undergoes fragmentation. The torqueprovided by the spiral structure at this maximum sustainable amplitude corresponds to
ELF-GRAVITATING ACCRETION DISCS Fig. 10. – The relation between gravitationally induced viscosity α and the cooling time τ cool in thermal equilibrium, expressed in Eq. (68), for three different values of γ = 2 (solid line), γ = 5 / γ = 7 / τ cool at which fragmentation is observed in [65], while the triangle refers tothe 2D simulations presented in [61]. It appears that fragmentation occurs whenever the valueof α needed to be provided by self-gravity in order to reach thermal balanceexceeds a maximumvalue of order α max ≈ . α max ≈ . α discussed above refers toa quasi-steady state, in which the disc stays in thermal equilibrium and the relevantphysical quantities do not change significantly on time scales shorter than the thermaltimescale. This does not happen for massive discs (with masses comparable to that of thecentral object), that can generate transient strong spiral episodes, with correspondinglylarge values of the stress α , which, however, do not last for longer than one dynamicaltimescale.If the fragmentation threshold is ultimately due to an inability of the disc to redis-tribute angular momentum on short timescales for a prolonged period of time, then it ispossible to envisage a further possible route to fragmentation, in cases where the disc,even if the cooling timescale is long, is fed at a high accretion rate. This possibility,although it may play an important role in several contexts, has not yet been tested nu- GIUSEPPE LODATO merically. In particular, the maximum accretion rate sustainable at any given radius inthe disc is only dependent on the local temperature and is given by ˙ M max = 2 α max c /G .
7. – APPLICATION TO OBSERVED SYSTEMS: STAR AND PLANETFORMATION
The discussion of the sections above has focussed on the basic physical processesthat determine the evolution of self-gravitating accretion discs, without referring to anyspecific astrophysical system where such mechanisms can be at work. We have thusseen what are the basic parameters involved in the determination of the disc propertiesand evolution. We now turn to an examination of a few concrete cases and we will tryto connect the theoretical picture described above to the observed properties of thesesystems. .1. Protostellar discs . – The formation of protostellar discs is the natural outcomeof the process of star formation. Stars are supposed to form from dense condensationswithin turbulent molecular clouds, called molecular cloud cores (see [80] for a recentreview). The typical specific angular momentum of one such core is of the order of10 cm / sec and is therefore much larger than the specific angular momentum of a star(for comparison, the Sun has a specific angular momentum of the order of 10 cm / sec).The gas that is going to form a star therefore has to get rid of most of its angularmomentum. The collapse of a molecular cloud core is thus bound to form initially aquite extended flattened object, a disc, with sizes of the order of tens to hundreds ofAU. Angular momentum redistribution inside the disc allows accretion to proceed and tofeed the forming protostar. Clearly, during this process the mass balance in the star-discsystem moves from a disc-dominated initial stage, to a star-dominated evolved stage,after the main infall phase has finished. The effects of the disc self-gravity are thengoing to be more important for younger object than for evolved ones. This evolution hasan observational counterpart in a frequently used classification [81], which identifies asClass 0 and Class I the youngest objects, which are still embedded into a dusty gaseousenvelope. During this phase the disc is still actively accreting from the envelope andthe disc mass is expected to be very large. Unfortunately, during this phase most of theactivity which occurs in the disc is enshrouded by the envelope and is difficult to observe,with the exception of a few cases. Class II objects have cleared out most of the envelopeand the protostar-disc system is now visible. Such objects are often called Classical TTauri stars. Since the main infall phase has ended, disc masses in the T Tauri phase areexpected to be smaller than in the Class 0 and Class I case. Finally, Class III objects arethe most evolved class of protostars, where the gaseous disc has almost entirely dispersed,leaving a remnant ‘debris disc’ mostly formed of rocks and pebbles (akin to the Zodiacaldust in our Solar System).The properties of protostellar discs are obviously a strong function of the stellar mass.Most of the data available at the moment refers to discs around solar type stars. However,in recent years a significant amount of data has started to become available for discs atthe extremes of the stellar mass distribution, that is for brown dwarfs at the low massend, and for massive O and B type stars at the upper end of the spectrum. In the lattercase, in particular, there is strong evidence pointing to the presence of very massive discs,and these objects are therefore an interesting laboratory to test the models described inthe previous sections. We will therefore discuss both the well known case of solar typeprotostars and the case of massive O-B protostars. ELF-GRAVITATING ACCRETION DISCS .2. Solar mass young stars . – Disc masses in the T Tauri phase are generally esti-mated from sub-mm emission, since at these wavelengths dust emission is usually opti-cally thin and we can thus easily convert the observed flux into a mass, if we know whatis the dust opacity [16]. Recent surveys of protostellar discs around solar mass starsin the Taurus-Auriga complex [82] indicate disc masses between 10 − M ⊙ and 10 − M ⊙ ,with a median value of 5 10 − M ⊙ . In another case, [83] the results for a smaller sample(possibly biased to larger masses) indicate a median mass one order of magnitude largerand in a few cases discs as massive as 0 . M ⊙ have been found. Similar values are alsoreported for the Orion region, where in a few cases masses as large as 0 . M ⊙ have alsobeen reported [84]. Since the stellar mass for these samples ranges from 0.1 to 1 M ⊙ , themass ratios between the disc and the star is between 10 − and a few times 10 − . Thethickness of protostellar discs can be estimated from measurements of the temperature inthe disc, and it turns out that in general H/R ≈ .
1. Bearing in mind that gravitationalinstabilities develop when M disc /M ⋆ ∼ H/R , we see that objects in the upper end of themass distribution for T Tauri stars are expected to be gravitationally unstable.It should be noted, however, that such measurements suffer from very large systematicerrors, mostly due to uncertainties in the dust opacity. Indeed, in many cases there isevidence for relatively evolved dust grains, so that the dust size distribution can besignificantly different than in the interstellar medium [85,86]. If dust grains have evolvedto relatively large sizes, their opacity would be reduced and the inferred disc mass wouldcorrespondingly increase. Additionally, it is important to stress that such measurementsare sensitive only to the dust content of the disc, which is expected to be a very smallfraction of the total gas mass (in the interstellar medium, for example, the dust to gasratio is of the order of 0.01). The dust mass measurements are then rescaled to obtain thegas mass, by adopting the standard interstellar medium scaling factor, the applicabilityof which is highly uncertain. In summary, it is likely that most of the above estimatesare actually underestimating the real disc masses [32]. We thus see that even in manyT Tauri objects gravitational instabilities might contribute substantially to the angularmomentum transport. The importance of self-gravity becomes even higher when oneconsider that T Tauri stars are expected to have already accreted much of their massfrom the disc. If by this rather evolved stage self-gravity is still important, it must havedominated the dynamics of younger systems, such as Class I objects.Accretion rates are generally measured from the strength of emission lines emitted asthe gas reaches the star in the so called ‘magneto-spheric accretion’, when the innermostparts of the disc are truncated by the stellar magnetic field, and the gas falls almostfreely onto the star along magnetic field lines, or from the veiling of photospheric linesdue to accretion shocks [87, 88]. It should be noted here that these measurements referto the accretion rate onto the star, which does not need to correspond to the accretionrate at larger distance, if significant sinks of matter are present, such as a disc wind(which is generally very small) or, more likely, as a massive planet, whose presence mightact like a dam for the accretion flow [89]. In any case, typical values range between10 − − − M ⊙ / yr. With such relatively low values of the accretion rate, the expecteddisc luminosity is relatively small and the energy balance in the disc is determined mostlyby irradiation from the central star [39]. While the detailed balance of course might varyfrom object to object, it is generally hard to look for accretion signatures in the broadbandSED of T Tauri stars.What about the role of self-gravity? Theoretical models of the evolution of protostellardiscs clearly indicate that in the early phases of star formation (in the Class I phase) thedisc is massive enough to be self-gravitating and the transport should be dominated by GIUSEPPE LODATO self-gravity [32-34]. Such effects are expected to become less important in the subsequentT Tauri phase.From the observational point of view, an intriguing possibility would be to directlyobserve a spiral morphology from high angular resolution observations. In the nearfuture, the Atacama Large Millimeter Array (ALMA), might in principle detect suchdetailed structures in the dust continuum. We have seen in section above that forlow mass discs (such as those that might be expected around T Tauri stars), the spiralstructure is tightly wound and characterized by high- m and short wavelengths (see Fig.6). This can be difficult to observe even at the high resolution achievable with ALMA.On the other hand, more open features, typical of more massive discs, might be detectedmore easily. In a couple of cases such extended spiral structure has been observed, inthe case of AB Aur (an evolved intermediate mass protostar) [90, 91] and in a case of ayoung Class 0 object [92]. The latter case is particularly interesting, since estimates ofthe disc mass indicate that the disc makes up more than 30% of the system. On the otherhand, in the case of the more evolved AB Aur, the disc mass is very low and estimatesof the disc surface density and temperature indicate Q ≈
10, which is too high to arguein favour of a gravitational instability. .3. FU Orionis objects . – A remarkable group of solar mass young stars are the so-called FU Orionis objects. FU Orionis objects are a small class of protostellar systemsundergoing large outbursts, during which their luminosity increases by as much as threeorders of magnitude. It is currently believed that the origin of the outburst lies in asudden enhancement of the accretion rate in the disc, that can rise to values of theorder of 10 − M ⊙ / yr. Unfortunately only a very small number of such systems have beenobserved and in only three cases (the prototypical FU Ori, V1057 Cyg and V1515 Cyg)we have a detailed knowledge of the light curve over a long period of time [93], whichenables us to have some insight on their evolution and therefore on the dynamics of theoutburst.FU Orionis objects are generally surrounded by dusty envelopes [94], and substantialinfall from the envelope, at a rate of approximately 10 − M ⊙ / yr is often required bydetailed modeling of the outburst [29, 30]. These characteristics tend to identify FUOrionis objects as transient episodes of enhanced activity in an otherwise quiescent ClassI object.The disc origin of the large luminosities observed in these systems has been con-firmed recently by high-resolution interferometric observations [95], which indicate thatemission comes from an extended (a few AU) source, like a disc. Additionally, it hasbeen shown [96, 97] that the spectral energy distribution of FU Orionis objects is to firstapproximation consistent with the one expected from a standard steady state accretiondisc (as sketched, for example, in Fig. 3). Indeed, this is essentially the only clear casewhere a protostellar disc displays the typical emission feature of an active disc. Detailedmodeling [97] shows that the typical accretion rates are of the order of 10 − M ⊙ / yr. Theagreement however only holds at relatively short wavelengths, in the optical and nearinfrared, while in the mid infrared a substantial excess emission with respect to the ex-pected disc spectrum is observed. Given the youth of the system, we expect that the discmass in these cases should be relatively large and it is therefore plausible to associatethe excess emission at mid-infrared wavelengths with the effect of self-gravity. Indeed,self-regulated discs are hotter in their outer parts, due to the combined effect of theself-regulation process and of the flattening of the rotation curve when the disc mass islarge. As a result, the SED tends to show some excess emission at long wavelengths (see ELF-GRAVITATING ACCRETION DISCS
12 13 14 1533343536 0 20 40 600100200300400500 self-gravitatingKeplerian 0 20 40 6005101520 self-gravitatingKeplerian
Fig. 11. – Modeling of FU Ori SED with a self-regulated self-gravitating disc. Left panel:observed (data points with error bars) and model SED (solid line). Middle panel: temperatureprofile implied by the self-reulated model (solid line) compared to a non-self-gravitating modelwith the same parameters (dashed line). Right panel: Rotation curve of the model (solid line)compared to a Keplerian curve (dashed line).
Fig. 4), consistent with observations of FU Orionis objects. Detailed modeling [49] hasshown how it is possible to fit accurately the long-wavelength SED of the best studiedFU Orionis objects, based on steady-state models of self-regulated discs (note that whileclearly these outbursting objects are evolving, the evolution timescale is longer than thedynamical timescale, so that a steady-state model is reasonably appropriate). An exam-ple of such modeling is shown in Fig. 11 and refers to the prototypical case of FU Ori(full details can be found in [49]). These models do however require fairly massive discs,with M disc ≈ M ⋆ . Such models predict deviations from Keplerian rotation, that mightbe detected from observations of the line profiles emitted by these systems [98]Traditional models to explain the cause of the outburst assume that this is due to athermal instability in the disc, caused by the steep change in disc opacity as hydrogenbecomes ionized [29, 99]. Such models have been extensively tested against observationsand reproduce the broad features of the light curve. However, agreement is only obtainedfor rather extreme values of the parameters involved (for example, the α viscosity param-eter is required to be as low as 10 − − − ) and by introducing some ad hoc triggeringevent [30, 100].An alternative outburst model is that the increase of mass accretion rate is due tothe development of a gravitational instability. As we have seen above, SED modeling isconsistent with the presence of fairly massive discs. Now, when the disc mass is of theorder of the central object mass, we are in the regime where we expect the developmentof large scale global gravitational instabilities (cf. Section ). Numerical simulations ofthis regime [46, 101] show that accretion can be highly variable, and indeed character-ized by periods of quiescence followed by recurrent episodes of enhanced accretion. Agravitational origin for FU Orionis outburst has been proposed [102] (or a combinationof gravitational and MHD instabilities, [103]), even though a detailed comparison of themodel light curves with observations is still lacking. .4. Massive protostars . – In the previous two subsections, we have concentrated onthe properties of relatively low mass young stars, i.e. stars with a mass of the order ofone solar mass. Stars of higher mass behave in a significantly different way. One of the GIUSEPPE LODATO main reasons behind this different behaviour lies in the fact that the typical contractiontime of the forming protostar (that is roughly the time it takes for a forming star to ignitenuclear reactions in its center) is much shorter for high mass stars than for low mass stars.In particular, for stars of mass larger than ≈ M ⊙ , the star reaches the main sequencewhile still accreting matter from the surrounding protostellar core [104]. Hence, while forlow mass stars we can observe optically visible pre-main-sequence stars as T Tauri (ClassII) stars after the main infall phase has terminated, which are therefore surrounded bya relatively low mass disc, this is not the case for stars of mass larger than 8 M ⊙ . Inthe latter case, the protostar is always going to be found during the main infall phase,and we thus expect that, if discs are present also here, they should be characterized bymuch larger masses than in T Tauri stars, and their properties should reflect more thoseof younger protostars (Class 0 or Class I).Clearly, observations of high mass protostars and of their circumstellar environment(including the possible presence of discs) are difficult, due to the fact that they are stillembedded within their parent core. Nevertheless, observations in the last few years havebeen able to reveal a significant number of them, and there is now a substantial evidencethat protostellar discs are present also for high mass stars (see [105] for a recent review).It should be mentioned however that the evidence for the presence of discs is mostlycircumstantial. In most cases, the ‘disc’ is found as a structure that lies perpendicular toa large scale outflow, which is common around young stars, and is easy to detect. Suchstructures are either seen in infrared continuum observations as dark silhouette discsagainst a bright background [106], or through observations of thermal [17, 18] and non-thermal (usually maser) emission lines [107-110], which have the advantage of providinginformation about the velocity structure of the circumstellar environment, so that a disccan be associated with a velocity gradient perpendicular to the outflow axis, usuallyinterpreted as a sign of rotation.The general result is consistent with the expectations outlined above. The discs areindeed found to be very massive, with masses comparable to, and in some cases evenlarger than, that of their central star. Although the data are still very uncertain, thereis also some evidence for non-Keplerian rotation, as expected when the disc mass ishigh [18]. Observed structures can be broadly divided in two groups: some are veryextended (with sizes of ≈ AU), very massive (with masses much larger than theircentral stars) and the luminosity of the central object is generally consistent with an O-type protostar (mass larger than 20 M ⊙ ). The second group of objects is characterized bysmaller sizes ( ∼ AU), smaller mass ( M disc . M ), and are generally associated withB-type stars (with masses between 8 and 20 M ⊙ ). Objects in the first group are so largethat their typical orbital timescale is of the order of the lifetime of the structure. Suchobjects are therefore unlikely to be in centrifugal equilibrium and is therefore difficultto identify them as accretion discs (even though a smaller size disc can in principle behidden within these structure). On the other hand, objects in the second group may beregarded as proper discs.Let us now consider what accretion rates are associated with the discs in the secondgroup of objects. Here again, observations are very difficult. A well constrained parame-ter is the outflow rate, which is produced in the inner parts of the disc (at a distance . − − − M ⊙ / yr. Clearly, the accretionrate through the disc should be at least of this order of magnitude to feed the outflow.Can the disc sustain such large accretion rates? Let us first consider the case wheretransport can be simply described in terms of a local, viscous phenomenon. As we haveseen above, this is the expectation when the transport is induced by a gravitational insta- ELF-GRAVITATING ACCRETION DISCS bility in low mass discs. We have also seen that in this case a quasi-steady, self-regulatedspiral structure does not provide a torque larger than α ≈ .
06 (when measured in termsof the α -prescription) without undergoing fragmentation. It is then easy to estimatethe accretion rate expected under these conditions, from ˙ M ≈ αc /G , where c s is thesound speed. Typical temperatures in the outer disc for these objects are of the order of170K (this is the case of IRAS 20126+4104, the best studied among these objects [18]),which correspond to a maximum accretion rate of ˙ M ≈ − M ⊙ / yr, well below therate needed to steadily feed the outflow from the inner disc. On the other hand, we haveto keep in mind that the picture where gravitational instabilities provide a local, viscoussource of angular momentum transport is only appropriate for relatively low mass discs.Here, the disc masses are much higher and hence the situation is more similar to thecase shown in Fig. 7, than to the one shown in Fig. 6. We should thus expect thedevelopment of a series of recurrent global spiral structures, which can temporarily (fora period of the order of the outer dynamical timescale, which in this case is of the orderof 10 years) transport angular momentum much more efficiently than a local processwould do. As discussed above, during such episodes the stress tensor induced by thespiral structure can be up to a factor ten larger than in a quasi-steady self-regulatedcase, and we can thus reconcile in this way the observationally implied accretion rateswith the rates expected from theoretical arguments. .5. Planet formation . – The process of planet formation is intimately linked to theprocess of star formation. The interest in this subject has grown significantly in the lastten years, and in particular after the discovery of the first planetary system around asolar type star [111], outside our own Solar System. To date, hundreds of extra-solarplanets have been discovered, and the list is ever-growing. Dedicated space missions,such as Darwin and the Terrestrial Planet Finder, are going to be launched in the nearfuture. This large set of data allow us to compare theoretical models of planet formationwith a statistical sample, rather than to the only one example (the Solar system) thatwas available until very recently. Our theories about planet formation are going to besignificantly improved and modified as new data become available.The properties of extra-solar planets appear to be very different from the ones ofthe planets in our Solar System, although in many cases they reflect the limitationsintrinsic to the detection methods used. In fact, most extra-solar planets are discoveredthrough the ‘radial velocity’ technique, based on the measurement of the radial velocityinduced by the presence of the planet on the motion of the host star. This methodnaturally favours the detection of planets of high mass and large eccentricities, whichwould induce the largest radial velocity. Also, planets at large separations, which havevery long periods, cannot be detected this way. An additional observational limitationis given by the fact that the radial velocity provides information about M p sin i , where M p is the planetary mass and i is the inclination of the planetary orbit. If (as it oftenhappens) the system is not eclipsing, we do not have information on the inclination andthus the planetary mass is only known to within this potentially large unknown factor.In any case, the inferred planetary masses range from 5 Earth masses, in the case ofthe ‘super-Earths’ recently discovered [112], up to 20 Jupiter masses, much larger thanthe most massive planets in the Solar System. It should be noted that many of thesemassive planets are found at very small separations from their host stars, smaller than0.1 AU, where the locally available mass is too small to form the planet, thus suggestingthat during the course of the planet formation substantial radial migration can happen.Additionally, while in the Solar System the orbital eccentricity is generally low, below GIUSEPPE LODATO M ⊕ , inthe case of Jupiter the upper limit is 14 M ⊕ and models with no core at all are also notexcluded.Currently, there are two models for the formation of planets. The first (and themore widely accepted) is the so called ‘core accretion’ model, while the second (whichwas usually disregarded, but has recently gained progressively more attention) is the socalled ‘gravitational instability’ model. Clearly, in the context of the present paper, thesecond model appears to be more interesting. However, as we shall see, gravitationalinstabilities may play an important role even in the core accretion model. While athorough review of planet formation theories is clearly beyond the scope of this paper, Ihere briefly summarize the main features of the two models. .5.1. The core accretion model . In the core accretion model [114] (see also [115]and [116] for two recent reviews), planets are supposed to form sequentially out of thesolid component present in protostellar discs. The first step in planet formation is thegrowth of interstellar dust through inelastic collisions up to a size of tens of centimeters.This growth is seen in observations of protostellar discs, whose opacities often indicatesthe presence of large dust grains [85], and where high-resolution spectroscopy in theinfrared indicate the presence of structured silicate grains [117, 118] (incidentally, thelarge uncertainties in protostellar disc masses mentioned in section .1s are essentiallydue to the unknown level of dust growth). For sizes larger than one meter, the collisionalagglomeration of dust is not effective. The growth in this mass range, up to the formationof kilometer-sized planetesimals, is still largely unknown (but see below). Once km-sized planetesimals have formed they essentially decouple from the gas and keep growingthrough their mutual interaction, enhanced by gravitational focussing. In this way a fewplanetary cores emerge from the background, and if the core reaches a size of the orderof 5 − M ⊕ it is able to further accrete a large gaseous envelope, hence becoming a giantplanet, like Jupiter. An important attractive of this model is that it is able to produceboth small, Earth-like, rocky planets, and the large Jupiter-like gaseous planets. Notethat within this model, Jupiter is bound to have a massive rocky core, the evidence infavour of which, as mentioned above, is not robust. Apart from this, two other problemsaffect this model for planet formation. First of all, it is a relatively slow process. Detailedmodels [114] require 10 million years in order to form Jupiter, which is just within theestimates of the lifetime of protostellar discs. A second important problem is related tothe formation of planetesimals. It is in this respect that the disc self-gravity might playa fundamental role in planet formation and I thus discuss it here in some detail.The issue arises when one considers the motion of meter sized solid particles withina gaseous disc. Even in the case where the disc contribution to the radial gravitationalfield is negligible, the gas rotation can be non-Keplerian, due to the effects of pressuregradients. As discussed in section .3, pressure gradients provide additional support tothe gas against gravity, so that the rotation curve can be non-Keplerian to some extent(Eqs. (15) and (16)). The degree of non-Keplerianity in this case is determined by ELF-GRAVITATING ACCRETION DISCS -2 0 2 423456 R=1 auR=3 auR=5 auR=10 au Fig. 12. – Solid particles migration timescale as a function of particle size for a typical protostellardisc configuration. The disc is here assumed to have a smooth pressure distribution, uniformlydecreasing outward as a power law. The four lines refer to four different radial locations in thedisc. From [21]. ( c s /v φ ) ≈ ( H/R ) . The sign of this contribution, depends on the pressure gradient,so that if the pressure decreases with increasing radius, the rotation is sub-Keplerian.Solid particles, on the other hand, do not feel the effect of pressure and therefore theirrotation curve to first approximation is Keplerian. This thus determines a velocity dif-ference between the solids and the gas (note that this effect is analogous to the so called‘asymmetric drift’ which occurs in galaxy dynamics, cf. [119]). If the solids and the gasare coupled via a drag force (for example, the usual Stokes or Epstein drag, which areusually dependent on the solid particle size) this may have some sizable effect on thesolid’s motion. Let us consider a smooth disc with pressure decreasing outwards [20],in this case the solids would move faster than the gas and the drag force therefore actsin such a way to remove angular momentum from them and thus determines an inwardradial migration. The magnitude of such radial motion and the timescale for radial mi-gration depends on the particle size. In particular, for very small sizes the solids are sostrongly coupled to the gas that they essentially follow the gas motion and their velocitydifference with respect to the gas is negligible. For very large sizes, on the other hand, thecoupling becomes very weak, and the solids move in essentially Keplerian orbits, with avery small inward velocity. The most interesting behaviour occurs for intermediate sizes,for which the solids are coupled to the gas, but not so strongly to completely follow the GIUSEPPE LODATO , Fig. 13. – Surface density of the solid component in a self-gravitating gaseous discs, where thegas and the solids are coupled through drag. Left panel: solid size equal to 10 meters. Rightpanel: solid size equal to 50 cm. While in the 10 meter case the solid density is very similarto the corresponding gas density, in the 50 cm case the solid concentrate very efficiently intofilamentary structures at the peaks of the spiral arms. From [21]. gas streamlines. The migration timescale as a function of particle size is shown in Fig.12 for a typical smooth protostellar disc configuration at different radial distances fromthe star. It can be seen that for sufficiently small and for sufficiently large particle sizesthe migration timescale is very long. For meter-sized objects however, the effect becomesstronger and the migration time can be as small as a few thousand years, a very shorttimescale with respect to the typical solid growth time at these size. This then poses asevere problem for the core accretion model, in that as soon as the dust grows in size andreaches the meter-size region, it suffers from this very fast radial migration and wouldessentially disappear into the innermost regions of the disc before being able to growfurther and form planetesimals.The above behaviour occurs for a smooth disc, in which the pressure is a monotonicfunction of radius. If the disc is characterized by a long-lasting spiral structure, inwhich the density (and hence the pressure) presents local maxima, the situation can besignificantly different [120, 121]. Indeed, if the pressure increases with radius, the gasmotion is locally super-Keplerian and the direction of the drag force is reverted, leadingto a fast outward migration. In turn, in the presence of a local pressure maximum, thesolid tend to migrate and concentrate to the maximum. The timescale for this processcan be even faster than the migration timescale. First of all, recall that the typical sizeof a spiral perturbation due to the disc self-gravity is of the order of H . This meansthat the pressure gradient is enhanced with respect to a smooth disc by a factor R/H .Secondly, in order to migrate to the peak of the spiral structure, the solids have to movea distance
H/R shorter than in a smooth disc. We thus expect this migration to the peakto occur on a timescale (
H/R ) shorter than the timescale shown in Fig. 12. Simulationsof a two component self-gravitating disc, where both a gaseous and a solid component ELF-GRAVITATING ACCRETION DISCS are present have been run [21, 22] and essentially confirm the scenario outlined above.Fig. 13 shows a comparison of the density structure of solids of different sizes, i.e. 10meters (left) and 50 cm (right) [21]. Note that in these simulations the parameter Ω τ cool was taken to be much larger than unity in order to prevent the gaseous disc to fragment.In the 10 meter case, the solid surface density essentially follows that of the gas, whilein the 50 cm case, the solid concentrate in a filamentary structure at the peak of thespiral arm. This strong concentration can have important effects on the growth rateof solids in this size range. First of all, the enhanced density favours and speeds upcollisional growth. Additionally (and possibly more importantly), the solid componentmight become gravitationally unstable, both because of the high density and becausethe velocity dispersion of solid particles is strongly reduced as the particles get ‘trapped’within the spiral. Simulations of the same setup, and including the self-gravity of thesolid component, have shown that meter sized solids in the filaments easily collapse andform gravitationally bound larger objects [22]. Unfortunately, the mass resolution ofsuch simulations is too small to determine whether the resulting solid fragments wouldbe of the right mass range for planetesimals. .5.2. The gravitational instability model . In the gravitational instability modelgaseous giant planets are simply formed from the fragmentation of a self-gravitatinggaseous disc [122-125]. As we have seen in section above, this naturally follows fromthe development of gravitational instabilities, if the cooling rate in the disc is sufficientlylarge , such that Ω τ cool <
1. If this happens, then gravitationally bound objects areable to form rapidly, within a few orbital timescales. With respect to the core accretionmodel, the gravitational instability model has the advantage of being a fast process,hence obviating to one of the difficulties associated with core accretion. A downside of itis that it can only account for the gas giant planets and would not explain the formationof terrestrial ones. Additionally, in this model gas giants would not have a solid core.Now, while for Jupiter this might also be the case, a solid core is thought to be presentin Saturn. In a few cases, extra-solar planets have also been detected through eclipsesof their parent stars [126], which allows to put constraints on their size and hence ontheir inner density, which can then give hints on the presence of a solid core. In oneextreme case, HD 149026b, which has a total mass of 0 . M Jupiter , it is estimated thatthe core would be as massive as 70 M ⊕ [127], which would then exclude the gravitationalinstability model (but would also be very difficult for the core accretion one!).What would be the typical mass of a planet formed by gravitational instability? Asimple estimate can be obtained by considering that, in a thin disc, the most unstablewavelength is λ = 2 πH . The typical mass of a fragment would then be of the order of M frag ≈ Σ λ . This can be simply rewritten in terms of the mass of the star M and theaspect ratio H/R , as M frag ≈ πQ (cid:18) HR (cid:19) M. (73)Now, a gravitationally unstable disc has Q ≈
1, the typical aspect ratio of protostellardiscs is
H/R ≈ .
1, and if we take the central star mass M = M ⊙ , we get M frag ≈ M Jupiter . It then appears that this model would only form the most massive observedplanets.However, most of the debate over the gravitational instability model for planet forma-tion concentrates on whether the actual cooling time in protostellar discs can be short GIUSEPPE LODATO enough to induce fragmentation. This has been estimated analytically [128] and theresults is that the cooling time is short enough to induce fragmentation only at verylarge distances from the central star, of the order of 50-100 AU. On the other hand, ithas also been argued that if the disc becomes convectively unstable, then cooling mightbecome more effective, producing the correct conditions to induce fragmentation ( [129],but see [130] for a rebuttal of this hypothesis). Recent simulations that employ sophis-ticated (but still uncertain) radiative cooling prescriptions [51, 74] indicate that indeedfragmentation of protoplanetary discs only occurs at very large distances from the centralstar.Even if the cooling time is long, it has been argued that the possible interaction witha companion or the close passage of a nearby star might trigger disc fragmentation bytidally increasing the local density. Early numerical isothermal simulations did indeedshow such a behaviour [131-133]. However, it should be noticed that, as mentioned above,the response of an isothermal disc to a gravitational instability is bound to be relativelyviolent, since the stabilizing feedback loop associated with self-regulation cannot happenunder isothermal conditions. More recently, the same issue has also been investigatedwith a non-isothermal equation of state and the results appear to be different. It hasbeen found that both in the case of a disc in a binary system [134] and in the case ofa stellar fly-by [135], the effect of tidal heating is much more effective in stabilizing thedisc than the de-stabilizing effect of the associated density enhancement, so that even inthis case the ultimate criterion for disc fragmentation is that the cooling time be shorterthan the dynamical time (see [136] for an alternative view).However, not all hope is lost for the gravitational instability model. As we have seen,this mechanism is indeed able to form large gaseous planets (or small brown dwarfs)at large radial distances from the host star, of the order of 100 AU. Such a new plan-etary population is still not excluded from observations. Note that the radial velocitytechnique requires tracking the stellar velocity for a period of the order of the planetaryorbital period in order to detect a planet. Now, planets at very large distances have verylong periods and might therefore require more time in order to be detected. Addition-ally, in the future such a different planet population might become detectable by othertechniques, such as direct imaging.In this respect, a very interesting case is that of the system 2MASS 1207. It consistsof a brown dwarf primary, with a mass of 25 M Jupiter , and a planetary mass companion,with a mass of 5 M Jupiter , orbiting at a distance of 55 AU from the central brown dwarf[137] (recent revision of the distance to the system has provided slightly different orbitalparameters [138], with revised masses of 21 and 3 M Jupiter for the primary and thesecondary and a semi-major axis of 41 AU, but the following arguments remain essentiallyunchanged). The companion, 2MASS 1207b, is the first planetary mass object to be everimaged. While clearly the absolute mass of the companion is in the planetary regime,it is questionable whether this system should be really considered as a planetary systemaround a brown dwarf or rather as an extremely low mass binary, since the mass ratioand the separation are perfectly consistent with the statistical properties of binaries oflarger mass [139]. In any event, what is relevant here is to understand whether sucha planetary mass object formed through the core accretion scenario or rather throughgravitational instability. It has been shown [139] (see also [140]) that given the largeseparation of the system and the very low mass of the primary, it is extremely unlikelyto form this object via core accretion. On the other hand, this is a good example wherethe gravitational instability scenario is expected to work, given the large mass ratio andthe large separation. Whether such an object is a rare case, and especially whether cases
ELF-GRAVITATING ACCRETION DISCS like this also occur around larger mass primaries, is still unknown but future observationsmight provide interesting constraints.
8. – APPLICATION TO OBSERVED SYSTEMS: AGN
After describing the impact of the disc self-gravity on the dynamics of accretion discsat the smallest scales, that is on the AU scale of protostellar and planet forming discs,we now turn to explore how it may affect the structure and dynamics of accretion discsat the largest scales, that is at the parsec scales of discs feeding supermassive black holesin Active Galactic Nuclei. .1. Active Galactic Nuclei . – .1.1. Feeding the hole . It is well known that accretion discs feed the growth ofsupermassive black holes in the nuclei of galaxies. Here the typical scale of the systemis obviously much larger than in the case of protostellar discs. The central supermassiveblack hole has a mass that ranges from 10 to 10 M ⊙ . The radial extent of the disc ismore uncertain. While the high energy emission from the AGN comes from the hotterinnermost parts of the disc (from a few to a few tens Schwarzschild radii), there is alsosubstantial evidence that the disc can extend much further out. One important pieceof information in this context comes from the observation of water maser emission fromAGN [2-4, 141], which reveals the presence of a thin disc of gas in almost Keplerianrotation out to a distance of the order of 1 pc, which corresponds to 10 R • for a 10 M ⊙ black hole, where R • is the Schwarzschild radius. In some cases the kinematics of thesemaser spots presents some interesting features that will be discussed in more detail below.A very important difference between AGN discs and protostellar discs is that AGN discsare generally much thinner. Detailed models of the structure of the non-self-gravitatinginner disc [142] indicate that the aspect ratio is of the order of H/R ≈ − , althoughin the self-gravitating regime it can be larger, of the order of a few times 10 − [52].This has important consequences for the issue of the role of the disc self-gravity. Aswe have often mentioned above, a simple way to check the importance of the disc self-gravity is to compare the mass ratio M disc /M with the aspect ratio H/R . Now, given thesmall aspect ratio of AGN discs, a very light disc, with M disc /M ≈ − can already besubject to gravitational instabilities. It is then relatively easy to calculate the distanceat which we would expect the cumulative disc mass to be high enough to make the discunstable [19, 143]. It tuns out that this distance is of the order of 10 R • , or 10 − pc fora 10 M ⊙ black hole.There is still some debate about the radial extent of the disc feeding the central blackhole. Some models [144] postulate that the disc outer radius corresponds to the distancewhere it becomes self-gravitating, i.e. ≈ . − . −
100 pc [145]. Probably reality falls somewhere in between. In fact,on the one hand observations of water maser emission clearly indicate the presence of arotating gaseous disc at radii of the order of 1 pc, i.e. well within the self-gravitatingportion of the disc. On the other hand it might be difficult to reconcile with observationsdiscs which extend much further beyond this point. Let us see this in more details.The main difference between between AGN discs and protostellar discs with respect totheir behaviour in the self-gravitating regime lies in the different cooling rate, which forAGN is typically very short. This is usually seen by noting that the heating rate neededto keep the disc marginally stable at Q ≈ GIUSEPPE LODATO a viscous accretion disc, with reasonable values of α [143,146]. Put it in other terms, thisessentially means that the cooling timescale in AGN is much shorter than the dynamicaltimescale. If we now consider the results of the simulations described above in Section , we would thus conclude that the disc would fragment and form stars [147, 148]. Thiswould be disastrous in terms of AGN feeding because star formation would essentiallyremove most of the fuel for accretion and turn it into stars. Indeed, such arguments arethe main arguments in favour of AGN feeding from a small scale disc, but leave us withthe problem of explaining the presence of large maser emitting discs.One possible way out of this apparent paradox can be found if energy and angularmomentum are transported by global spiral waves in the disc, rather than by a local α -like process. As we have seen above, this usually happens if the disc mass is a sizablefraction of the central object mass. In this case, as already discussed, the energy balanceshould include some extra “global” terms [79], arising from wave transport of energy,that might provide the required energy to prevent fragmentation in the outer disc. Inthis picture, a density wave might remove free rotational energy from the inner disc,but rather than dissipating it locally (as would a standard viscous process do), it mightcarry it a long way out along the wave, and release it at large radii, where the wave isdissipated ( ). As seen above, this requires the presence of quite massive discs, and theevolution in this case is generally highly variable, with episodes of strong accretion andblack hole feeding followed by more quiescent periods where the accretion rate is small.Such a time variable accretion model has also been sometimes proposed [149]. In someway, this solution is similar to the proposed solution for the high mass accretion ratesobserved around high mass young stars above (Section .4).Alternatively, even if fragmentation does occur, it is still unclear whether accretionwould be inhibited altogether, or whether we might still have accretion with only arelatively small amount of mass ending up in stars. Numerical simulations [150] show thatwhile it is true that for small cooling times the disc undergoes fragmentation, if a smallamount of energy is released from the forming stars (either from accretion luminosityfrom the stars or from the ignition of nuclear burning), this is sufficient to preventfurther fragmentation (see also [151]). Depending on the parameters (and in particular ifthe cooling time is close to the fragmentation threshold) the lifetime of the gaseous disccan be very long. In this case, the missing energy to keep the disc self-regulated wouldbe provided by star formation feedback rather than by global energy transport, as in themodel above.In both models described above, the disc is then allowed to extend outside ∼ R • ,that is the radius where it becomes self-gravitating, provided that some extra-heatingsource is able to maintain it close to marginal stability in a self-regulated way. We havealready noticed that self-regulated discs are thus characterized by a secondary bumpin the broadband Spectral Energy Distribution, which can be prominent at infraredwavelengths. As shown in Fig. 4, the luminosity in this long wavelength bump is anincreasing function of the location of the outer radius of the disc [49]. This fact canthen be used to put constraints on the disc outer edge [52]. It turns out that if the discextends beyond ≈ ( ) Note that even a standard viscous process does transport energy between different regionsof the disc, but this does not affect the arguments above, where we assume that in a wave-dominated transport the relative contributions of dissipation and transport of energy are signif-icantly different from a viscous one. ELF-GRAVITATING ACCRETION DISCS by thermal pressure, rather than by turbolent motions of gas clumps and clouds withinthe disc) then the infrared bump would become too luminous, exceeding the typicalobserved SED of Active Galactic Nuclei. We thus see that, while an approximately self-regulated accretion disc can be present well beyond 0.01 pc, thus potentially explainingthe presence of dense gaseous discs to produce the water maser emission, it is unlikelyto extend far beyond 1 pc. .1.2. Star formation in the Galactic Center . In the context delineated above,an important role is played by the evidence that has been gathered in the last fewyears, which points to the presence of a large number of young stars very close to thesupermassive black hole at the center of our own Milky Way [152, 153]. In particular,most of these stars appear to belong to two distinct stellar discs orbiting at roughly thesame distance to the black hole, i.e. at a distance of 0.05-0.5 pc [154, 155]. The mostlikely explanation for the origin of these stars in that they formed in situ and in particularfrom the fragmentation of a self-gravitating accretion disc [154, 156]. Such observationsthus fit naturally in the context described above, since we know that at parsec distancesan AGN accretion disc would be self-gravitating and its cooling time is expected to beshort enough to induce fragmentation. The conditions in the Galactic Center might betypical of other galaxies, where a nuclear starburst can be a result of the very samemechanism [151, 157]. .1.3. Kinematics of water masers . In the discussion above, we have often referredto the detection of water maser emission from the inner parts of AGN. Indeed, suchobservations are a valuable tool to study the disc at the radial distances where self-gravityis expected to be important. In the section above, we used water maser observationssimply to indicate the mere presence of discs at radii of the order of 0.1-1 pc. On theother hand they also offer a way to probe the disc geometry (for example, the presenceof a warp [158-160]) and, most importantly, a way to probe its kinematic. In fact, watermasers are usually observed to arise from a series of distinct “maser spots”, at differentradial distances from the nucleus, and from the Doppler shift of the maser emission lineit is then possible to reconstruct the rotation curve. In some cases, for example for NGC4258, the resulting rotation curve is very close to Keplerian [2], and it thus allows a veryprecise determination of the mass of the central BH, which for NGC 4258 is 3 . M ⊙ .However, in many other cases the rotation curve, while still displaying a smoothdeclining profile, as would be expected for a rotating disc, does not follow exactly Kepler’slaw. This is for example the case of NGC 1068 [3, 141], of the Circinus galaxy [160] andof NGC 3079 [161]. In particular, for the case of NGC 1068 the maser data are consistentwith a circular velocity v φ ∝ r − . [141]. Given the discussion above, which shows thatat a scale of a fraction of a parsec, where the maser spots are detected, the disc canbe self-gravitating, it is then tempting to attribute such (often small) deviation fromKeplerian rotation to the contribution of the disc self-gravity.A detailed fit to the circular velocity traced by water masers in NGC 1068 with amodel which incorporates both the gravitational field of the black hole and that of thedisc has been performed [19] (see also [162]), by using the self-regulated disc modelsdescribed in Section .3. The results of such modeling are shown in Fig. 14, whichillustrates the observed rotation curve and the fit with a self-regulated model [19]. Asseen in Fig. 14, the observed rotation curve of the maser spots can be divided in twodifferent regions: 1) at small impact parameter the masers show an apparently risingrotation curve; 2) starting from a radial distance ≈ . GIUSEPPE LODATO
Fig. 14. – The data points show the observed rotation curve measured from water maser emissionin the AGN NGC 1068. The solid line is the best -fit (non-Keplerian) model to the data, obtainedfrom a self-regulated, self-gravitating disc model (see Section .3). In this model both the discand the black hole contribute to the radial gravitational field, and we have M disc ≈ M =(8 . ± .
3) 10 M ⊙ . The dashed line is a fit fit obtained from a Keplerian model, where the massof the black hole is M ≈ (1 . ± .
02 10 M ⊙ . curve declines, following a sub-Keplerian profile.The natural interpretation of the declining part of the rotation curve is that it arisesfrom material that moves parallel to our line of sight (i.e. that lies on a disc diameterperpendicular to the line of sight). The best argument in favor of this interpretation isthat maser amplification is largest for material that lies close to the line of the nodes. Onthe other hand, the rising part of the observed “rotation curve” is thought to originatefrom one quarter of the disc at the inner maser disc radius, so that the rising curve isan effect of velocity projection along the line of sight (see also [2]). According to thisinterpretation, the inner radius of the maser disc is located at ≈ . ELF-GRAVITATING ACCRETION DISCS disc radius is at ≈ M and the disc surface density Σ, which in turn (see Section .3) givesus information on ˙ M /α . The resulting black hole mass is M = (8 . ± .
3) 10 M ⊙ and thedisc mass is approximately equal to the black hole mass. From the required disc surfacedensity it is then possible to obtain ˙ M = (28 . ± . αM ⊙ / yr. The mass accretion rate˙ M can be estimated for example from the bolometric luminosity as ˙ M ≈ . M ⊙ / yr,and we thus obtain α ≈ . − , which is of the right order of magnitude as would beexpected from the transport induced by gravitational instabilities.The best-fit curve resulting from the self-gravitating disc model is not a power-law.However, as mentioned above the data are consistent with a rotation curve of the form V ∝ r − . [141]. Indeed, computing the quantity d ln v φ / d ln R for the best fit model, itis easily seen that it ranges from − .
35 at the inner disc edge of the disc to − .
30 at theouter disc edge.The dashed line in Fig. 14 shows the result of a fit done with a Keplerian disc model.The fit is significantly worse than the non-Keplerian one (a detailed statistical analysiscan be found in [19]). The resulting best-fit value of the black hole mass in the Keplerianfit is M ≈ (1 . ± .
02) 10 M ⊙ . Note that the total mass (disc + black hole) of the self-gravitating best fit model is roughly the same as the black hole mass of the Keplerian fit.Therefore, a non-self-gravitating fit, which attributes all the mass to the central object,gives the correct value for the total mass, but fails to provide the correct slope of therotation curve. .2. Supermassive black hole formation . – In the section above we have concentratedon supermassive black holes in the local Universe, i.e. in the Milky Way and in nearbyAGN. We now turn our attention to very high redshifts, where we expect the emergenceof the seeds of the first supermassive black holes. We will then see that even at suchearly epochs, phenomena related to the development of gravitational instabilities mightplay a key role in determining the properties of the black hole population.Observations of AGN at high redshifts, up to z ∼ M ⊙ were already in place when the Universe wasonly 10 years old. This clearly requires that the black hole growth occurred at very highrates, with an average of 1 M ⊙ / yr. Such a rapid early growth poses serious challenges tomodels of their formation. Some models [165-167] assume that the seeds of supermassiveblack holes are the remnants of the zero-metallicity first stars (the so-called PopulationIII stars), that are expected to be relatively massive [168, 169] and thus produce blackholes with a mass of up to 100 M ⊙ . However, unless the efficiency of conversion of mat-ter into energy through the accretion process is very low, it is impossible to grow theseeds to the required masses by z ∼ ǫ = L/ ˙ M c exceeds ≈ . L is the accretion luminosity and c is the speed of light), the Eddingtonlimit does not allow the large accretion rates needed to grow the seeds fast enough tobecome bright AGN by z ∼ GIUSEPPE LODATO small. Once the black hole mass has exceeded 10 − M ⊙ , further growth at ≈ M ⊙ / yrcan proceed at sub-Eddington rates.The efficiency is in turn dependent on the spin of the black hole, with high spinproducing very large efficiencies ǫ ∼ .
5. Accretion of matter naturally tends to spin upthe hole [166] and hence to increase the efficiency, thus exceeding the Eddington limitfor relatively low ˙ M and preventing a fast growth of the hole. While recent calculations[171, 172] show that it is possible to keep the hole spin low if the growth occurs throughseveral small randomly oriented accretion episodes [173], we still have to face the issueof how to produce the high infall rates required.Alternative models propose the direct formation of more massive seeds with masses ofabout 10 M ⊙ directly out of the collapse of dense gas [174-180]. The key limiting factorfor these models is the disposal of the angular momentum. Recently, it has been proposed[179-182] that large scale gravitational instabilities developing during the growth of pre-galactic discs is the missing ingredient, able to funnel the required amount of gas intothe center of the galaxy.According to such models, the formation of the seeds of supermassive black holesoccurs at a redshift z ∼ −
15, when the intergalactic medium had not been yet enrichedby metals forming in the first stars. As a consequence, the chemical composition of thegas at this early epoch is essentially primordial, i.e. the gas is mostly hydrogen andhelium. The cooling properties of this gas are therefore relatively simple. In particular,in the absence of molecular hydrogen, the main coolant is provided by atomic hydrogen,for which the cooling timescale becomes extremely long for temperatures smaller than ∼ K, and we thus expect the gas to reach thermal equilibrium at a temperature T gas of the order of 10 K.Now, consider a dark matter halo (modeled, for simplicity, as a truncated singu-lar isothermal sphere) of mass M halo and circular velocity V h , extending out to r h = GM halo /V . We also assume that the halo contains a gas mass M gas = m d M halo , where m d is of the order of the universal baryonic fraction, ≈ .
1, whose angular momentumis J gas = j d J , where j d ∼ m d . The angular momentum of the dark matter halo J isexpressed in terms of its spin parameter λ = J | E | / /GM / , where E is its total energy.The probability distribution of the spin parameter of dark matter halos can be obtainedfrom cosmological N-body simulations [183] and is well described by a log-normal distri-bution peaking at λ = 0 . T vir ∝ V is larger than the gas temperature T gas , the gas collapses and forms a rotationally supported disc, with circular velocity V h ,determined by the gravitational field of the halo. For low values of the spin parameter λ the resulting disc can be compact and dense. In this case, during the infall of gas ontothe disc, its density rises until the stability parameter Q becomes of the order of unity.At this point, the disc starts developing a gravitational instability, which as we have seenabove is able to efficiently redistribute angular momentum and allow accretion. Furtherinfall of gas does not cause the density to rise much further, but rather it promotes anincreasingly high accretion rate into the center. This process goes on until infall is overand the disc has attained a surface density low enough to be marginally gravitationallystable, i.e. with Q = ¯ Q . It is then possible to calculate what fraction of the infallingmass needs to be transported into the center to make the disc marginally stable, as a ELF-GRAVITATING ACCRETION DISCS Fig. 15. – Mass available for the formation of the seed of a supermassive black hole in the centerof pre-galactic discs as a function of the mass of the parent dark matter halo. The plots refer tothe following choice of parameters: ¯ Q = 2, T gas = 4000K, m d = j d = 0 . λ = 0 .
01 (solid line), λ = 0 .
015 (long-dashed line), λ = 0 .
02 (short-dashed line). The red curve shows the thresholdfor fragmentation from Equation (75), with α c = 0 .
06. Halos on the right of the red line giverise to fragmenting discs. function of the main parameters involved. In this way, we get [181, 182]: M BH = m d M halo − s λm d ¯ Q (cid:18) j d m d (cid:19) (cid:18) T gas T vir (cid:19) / , (74)where we have suggestively called M BH the accreted mass, since this mass is the totalmass available for the formation of the black hole seed in the center.However, for large halo mass, the internal torques needed to redistribute the excessbaryonic mass become too large to be sustained by the disc, which might then undergofragmentation. We have seen in the previous sections that the maximum torque that canbe delivered by a quasi-steady self-regulated disc is of the order of α c ≈ .
06. Since theinfall rate of gas from the halo is proportional to T / , we expect fragmentation whenthe virial temperature exceeds a critical value T max , given by (see [181] for details): T max T gas > (cid:18) α c m d
11 + M BH /m d M halo (cid:19) / , (75) GIUSEPPE LODATO
Fig. 16. – Left panel: Mass function of seed black holes predicted by the model based on Eqs.(74) and (75). The black solid line refers to z = 10, while the red line refers to z = 20. Thelong-dashed line shows the effect of reducing ¯ Q from 2 to 1.5. The short dashed line showsthe effect of not including the possibility of fragmentation (more details can be found in [182]).Right panel: the integrated density of black holes predicted from a merger tree evolution of theblack hole seed population. The three solid line refer to three different choices of the parameters(see [184] for details). The dashed area indicate the observationally permitted region fromestimates of the black hole density at low redshifts. At high redshift ( z &
6) the density reflectsthe one attained after the seed formation phase, while the rise at lower redshifts indicates thegrowth through AGN activity.
Although it is possible, as often mentioned in this paper, that accretion proceeds evenfor larger values of α in a highly time-variable way when the disc mass is large, and itis also possible that accretion proceeds even in a fragmenting disc, we make here theconservative assumption that all halos that violate Eq, (75) do fragment and do notaccrete. Fig. 15 illustrates the relationship between halo mass and black hole massbased on Equation (74) for three different values of the spin parameter λ . The red line inFig. 15 correspond to Eq. (75), so that halos on the right of the red line are expected tofragment. We can thus see that the typical mass fed into the center of such pre-galacticdisc is of the order of 10 M ⊙ up to 10 M ⊙ . The typical accretion rates during thisearly epochs is of the order of 10 − M ⊙ / yr [181]. If such high masses are assembled asseeds of supermassive black holes at redshift 10 −
15, it is then easy to grow throughEddington-limited accretion to 10 M ⊙ by z = 6, as required by observations.Equation (74) provides a powerful link between the properties of dark matter haloesand the mass of massive seed black holes that can grow within them. As shown, theamount of mass that will be concentrated in the central regions of these pre-galacticdiscs depends only on halo properties (such as the spin parameter λ and the fraction ofbaryonic mass that collapses to the disc m d ), on the ratio between gas temperature andhalo virial temperature, and on the threshold value of Q , which has a very small range ofvariation around ¯ Q ≈
1. This simple model has been used to calculate several properties
ELF-GRAVITATING ACCRETION DISCS of the black hole population at high redshift. In particular, from the distribution ofhalo masses and angular momentum it is straightforward to derive the mass functionof the supermassive black hole seeds [182], which turns out to be strongly peaked ataround 10 M ⊙ , as shown in the left panel of Fig. 16. Furthermore, it is also possible toinclude such a simple prescription within evolutionary models that track the propertiesof the black hole population along cosmic time, such as merger tree models [184] (suchmodels are usually based on the so called hierarchical model for the growth of structures,where the structures observed in the present day Universe are the result of the subsequentmerger of smaller structures at early times). It is then interesting to see that the evolutionof such a primordial seed population can naturally account for the current estimates of thedensity of black holes at low redshift (Fig. 16, right panel). In addition, an importantand testable prediction of such models is that dwarf galaxies, that did not have anyprogenitor massive enough to seed a black hole, should not host a supermassive blackhole. In particular if the velocity dispersion of the galaxy is below ∼
9. – CONCLUSIONS AND CHALLENGES FOR THE FUTURE
As we have seen, self-gravitating accretion discs are an important aspect of the model-ing of several systems, and thus play a key role in modern astrophysics. The influence ofthe disc self-gravity extends from the very small scales, associated with the formation ofplanetary systems, up to the large scales, associated with the formation and the feedingof supermassive black holes in the nuclei of galaxies. While clearly the specific propertiesof systems associated with such hugely different scales are going to differ substantially,there is yet a unifying framework related to the scale-free nature of the gravitationalforce.Because of self-gravity, accretion discs can develop gravitational instabilities. Theseare going to deeply affect the structure and the evolution of the accretion disc. Firstly,we have seen that the instability usually takes the form of a regular spiral structure,that can efficiently redistribute the disc angular momentum and thus allow the accretionprocess to take place. In this respect the role of the disc self-gravity is complementary tothat of other instabilities, such as magneto-hydrodynamic (MHD) instabilities [31, 78].Indeed, while MHD instabilities are likely to be the most important source of angularmomentum transport in the inner and hotter parts of the disc, it is not clear whether theyare able to provide the necessary torques in colder environments, where the ionizationlevel of the gas is low and the coupling to the magnetic field is weaker. This is, forexample, the case of the outer parts of protostellar discs [185]. On the other hand,such cold parts of the discs are naturally subject to gravitational instabilities. It isthen possible to envisage a situation where gravitational instabilities are responsible forbringing the accreting material from large distances (which might be of the order of tensof AU for protostellar discs, and a fraction of a parsec for AGN discs) into the inner disc,where MHD instabilities take over and are responsible for the ultimate deposition of thematter onto the central object. In this respect, it would be important to develop modelswhich include both the effects of MHD induced turbulence and that of gravitationalinstabilities. Semi-analytical models along these lines have been discussed [34], but it isstill unclear what would be the interplay between magnetic and gravitational instabilitiesin the general case [186, 187].A second aspect related to gravitational instabilities is the natural tendency associatedwith self-gravity to collapse and produce gravitationally bound objects. Disc fragmenta- GIUSEPPE LODATO tion is a possible outcome of the onset of gravitational instabilities. As discussed above,it is now well established that the conditions leading to fragmentation are essentiallyrelated to the disc thermodynamics and to the cooling rate. From the modeling pointof view, fragmentation has a twofold aspect. On the one hand, if one is interested inproducing smaller objects within the disc, it offers a natural way to accomplish the task.On the other hand, there is the danger that most of the disc mass would end up inthe fragments, with very little being accreted onto the central object. In this respect,fragmentation is a dangerous outcome, which might inhibit accretion.We have seen that the biggest difference between protostellar discs and AGN discswith respect to self-gravity lies indeed in their thermal properties, and in particular inthe fact that AGN discs are expected to be cold (in the sense that
H/R ≪
1) and witha short cooling timescale (in the sense that τ cool Ω ≪ H/R ≈ .
1) and with a long cooling timescale ( τ cool Ω ≫
1, at least within 100 AU).The behaviour of these two kind of systems in this respect is then significantly different.Most protostellar discs are not expected to fragment, except perhaps in their outermostregions. This is the main problem affecting models of planet formation which rely ondirect disc fragmentation. At the other end of the scales that we have considered, wesee that AGN discs instead are likely to fragment. This is good news in view of formingyoung massive stars in the vicinity of the central supermassive black hole, as are observed,for example, close to the black hole in our Galactic Center. On the other hand, we havediscussed above the problems that might arise from disc fragmentation in this context,related to the feeding of the black hole.The last 10-15 years have witnessed a significant improvement of our understandingof the evolution of self-gravitating accretion discs. This has been made possible throughthe interplay between simple analytical models and the results of detailed numericalsimulations of the disc hydrodynamics. The big advances in computing power in thelast decade have made it possible to run such numerical simulations with unprecedentedresolution, and progressively including the effects of several physical processes in moredetail. In this paper, emphasis has been put on the important role of the gas coolingrate in determining the outcome of the instability. This has been often introduced ina simplified way in the simulations, which has led to a valuable investigation of thedifferent regimes, in a somewhat academic way. The theoretical challenge in this respectis to introduce in the numerical codes a more realistic cooling function, in order tounderstand the behaviour of actual systems. For example, most of the current debate onthe applicability of the fragmentation scenario for planet formation is ultimately due toour ignorance of the actual thermodynamics at work in protostellar discs. The difficultyhere is to treat numerically the radiative transfer within very optically thick media. Onlyvery recently have numerical codes implemented radiation physics and, although stillpreliminary, some results are starting to emerge [51, 74]. Codes which couple radiativetransfer with the disc hydrodynamics are also desired to include the possibly importanteffect of irradiation from the central object, which has mostly been neglected in thesimulations run so far (with some exceptions [188]).If we look at the problem of star formation within AGN accretion discs, here theissue is ultimately whether accretion can proceed at significant rates even if the discfragments. In other words, it is still unclear how much gas needs to be turned into starsbefore fragmentation is quenched. Here, it would be important to include in the modelsthe energy feedback arising from the forming stars. Even in this context then, the keyingredient to be added into numerical codes is a proper treatment of radiative transfer,coupled with hydrodynamics.
ELF-GRAVITATING ACCRETION DISCS ∗ ∗ ∗ ACKNOLEDGMENTS. I would like to dedicate this paper to the memory of EduardoDelgado-Donate, who contributed to the development of some of the ideas presentedhere. I acknowledge several interesting discussions with Giuseppe Bertin, Cathie Clarke,Andrew King, Sergei Nayakshin and Jim Pringle. I would also like to thank the variouscollaborators with whom I have worked on these issues along the years, and in particularPhil Armitage, Priya Natarajan and Ken Rice.
REFERENCES[1]
Bertin G. and
Lin C. C. , Spiral Structure in Galaxies: a Density Wave Theory (MITPress, Cambridge) 1996.[2]
Miyoshi M. et al. , Nature , (1995) 127.[3] Greenhill L. J. and
Gwinn C. R. , Astrophysics & Space Science , (1997) 261.[4] Kondratko P. T., Greenhill L. J. and
Moran J. M. , ApJ , (2006) 136.[5] Frank J., King A. and
Raine D. , Accretion Power in Astrophysics (CambridgeUniversity Press, Cambridge) 2002.[6]
Reipurth B., Jewitt D. and
Keil K. (Eds.),
Protostars and Planets V (ArizonaUniversity Press) 2007.[7]
Paczy´nski B. , Acta Astronomica , (1978) 91.[8] Kolykhalov P. I. and
Sunyaev R. A. , Soviet Astronomy Letters , (1979) 180.[9] Lin D. N. C. and
Pringle J. E. , MNRAS , (1987) 607.[10] Lin D. N. C. and
Pringle J. E. , ApJ , (1990) 515.[11] Pringle J. E. , ARA&A , (1981) 137.[12] Binney J. and
Tremaine S. , Galactic dynamics (Princeton, NJ, Princeton UniversityPress, 1987, 747 p.) 1987.[13]
Bertin G. and
Lodato G. , A&A , (1999) 694.[14] Mestel L. , MNRAS , (1963) 553.[15] Bertin G. , ApJ , (1997) L71.[16] Beckwith S. V. W., Sargent A. I., Chini R. S. and
Guesten R. , AJ , (1990)924.[17] Cesaroni R. et al. , ApJ , (1994) 137.[18] Cesaroni R., Neri R., Olmi L., Testi L., Walmsley C. M. and
Hofner P. , A&A , (2005) 1039.[19] Lodato G. and
Bertin G. , A&A , (2003) 517.[20] Weidenschilling S. , MNRAS , (1977) 57.[21] Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J. and
Bonnell I. A. , MNRAS , (2004) 543.[22] Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J. and
Bonnell I. A. , MNRAS , (2006) L9.[23] Spitzer L. , ApJ , (1942) 329.[24] Toomre A. , ApJ , (1964) 1217.[25] Bertin G. , Dynamics of Galaxies (Cambridge University Press, Cambridge) 2000.[26]
Shakura N. I. and
Sunyaev R. A. , A&A , (1973) 337.[27] Hartmann L., Calvet N., Gullbring E. and
D’Alessio P. , ApJ , (1998) 385.[28] Lasota J. P. , New Astronomy Reviews , (2001) 449.[29] Bell K. R. and
Lin D. N. C. , ApJ , (1994) 987.[30] Lodato G. and
Clarke C. J. , MNRAS , (2004) 841.[31] Balbus S. A. and
Hawley J. F. , Reviews of Modern Physics , (1998) 1.[32] Hartmann L., D’Alessio P., Calvet N. and
Muzerolle J. , ApJ , (2006) 484.[33] Clarke C. J. , MNRAS , (2007) 1350. GIUSEPPE LODATO [34]
Kratter K. M., Matzner C. D. and
Krumholz M. R. , ArXiv e-prints , (2007) .[35] Lynden-Bell D. and
Pringle J. E. , MNRAS , (1974) 603.[36] Hartmann L. , Accretion Processes in Star Formation (Cambridge University Press,Cambridge) 1998.[37]
Lynden-Bell D. , Nature , (1969) 690.[38] Adams F. C., Lada C. J. and
Shu F. H. , ApJ , (1988) 865.[39] Chiang E. P. and
Goldreich P. , ApJ , (1997) 368.[40] Dullemond C. P., Hollenbach D., Kamp I. and
D’Alessio P. , Protostars and PlanetsV , (2007) 555.[41]
Hohl F. , ApJ , (1971) 343.[42] Ostriker J. P. and
Peebles P. J. E. , ApJ , (1973) 467.[43] Lau Y. Y. and
Bertin G. , ApJ , (1978) 508.[44] Bertin G., Lin C. C., Lowe S. A. and
Thurstans R. P. , ApJ , (1989) 104.[45] Lodato G. and
Rice W. K. M. , MNRAS , (2004) 630.[46] Lodato G. and
Rice W. K. M. , MNRAS , (2005) 1489.[47] Vandervoort P. O. , ApJ , (1970) 87.[48] Bertin G. and
Romeo A. B. , A&A , (1988) 105.[49] Lodato G. and
Bertin G. , A&A , (2001) 455.[50] Vorobyov E. I. and
Basu S. , MNRAS , (2007) 855.[51]
Boley A. C., Mej´ıa A. C., Durisen R. H., Cai K., Pickett M. K. and
D’AlessioP. , ApJ , (2006) 517.[52] Sirko E. and
Goodman J. , MNRAS , (2003) 501.[53] Miller R. H., Prendergast K. H. and
Quirk W. J. , ApJ , (1970) 903.[54] Hohl F. , ApJ , (1973) 353.[55] Anthony D. M. and
Carlberg R. G. , ApJ , (1988) 637.[56] Laughlin G. and
Bodenheimer P. , ApJ , (1994) 335.[57] Laughlin G. and
R`o˙zyczka M. , ApJ , (1996) 279.[58] Laughlin G., Korchagin V. and
Adams F. C. , ApJ , (1998) 945.[59] Pickett B. K. et al. , ApJ , (1998) 468.[60] Pickett B. K. et al. , ApJ , (2000) 1034.[61] Gammie C. F. , ApJ , (2001) 174.[62] Rice W. K. M., Armitage P. J., Bate M. R. and
Bonnell I. A. , MNRAS , (2003) 227.[63] Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., Jeffers S. V. and
Vine S. G. , MNRAS , (2003) L36.[64] Rice W. K. M., Armitage P. J., Bate M. R. and
Bonnell I. A. , MNRAS , (2003) 1025.[65] Rice W. K. M., Lodato G. and
Armitage P. J. , MNRAS , (2005) L56.[66] Pickett B. K., Mej´ıa A. C., Durisen R. H., Cassen P. M., Berry D. K. and
LinkR. P. , ApJ , (2003) 1060.[67] Mejia A. C., Durisen R. H., Pickett M. K. and
Cai K. , ApJ , (2005) 1098.[68] Mayer L., Quinn T., Wadsley J. and
Stadel J. , ApJ , (2004) 1045.[69] Rees M. J. , MNRAS , (1976) 483.[70] Silk J. , ApJ , (1977) 152.[71] Clarke C., Harper-Clark E. and
Lodato G. , MNRAS , (2007) 1543.[72] Johnson B. M. and
Gammie C. F. , ApJ , (2003) 131.[73] Whitehouse S. C., Bate M. R. and
Monaghan J. J. , MNRAS , (2005) 1367.[74] Mayer L., Lufkin G., Quinn T. and
Wadsley J. , ApJ , (2007) L77.[75] Stamatellos D., Whitworth A. P., Bisbas T. and
Goodwin S. , A&A , (2007)37.[76] Shu F. H. , ApJ , (1970) 99.[77] Lynden-Bell D. and
Kalnajs A. J. , MNRAS , (1972) 1.[78] Balbus S. A. , ARA&A , (2003) 555.[79] Balbus S. A. and
Papaloizou J. C. B. , ApJ , (1999) 650. ELF-GRAVITATING ACCRETION DISCS [80] McKee C. F. and
Ostriker E. C. , ARA&A , (2007) 565.[81] Lada C. J. and
Wilking B. A. , ApJ , (1984) 610.[82] Andrews S. M. and
Williams J. P. , ApJ , (2005) 1134.[83] Andrews S. M. and
Williams J. P. , ApJ , (2007) 705.[84] Eisner J. A. and
Carpenter J. M. , ApJ , (2006) 1162.[85] Testi L. et al. , ApJ , (2001) 1087.[86] Natta A., Testi L., Neri R., Shepherd D. S. and
Wilner D. J. , A&A , (2004)179.[87] Gullbring E., Hartmann L., Brice˜no C. and
Calvet N. , ApJ , (1998) 323.[88] Gullbring E., Calvet N., Muzerolle J. and
Hartmann L. , ApJ , (2000) 927.[89] Artymowicz P. and
Lubow S. H. , ApJ , (1996) 77L.[90] Lin S.-Y. et al. , ApJ , (2006) 1297.[91] Pi´etu V., Guilloteau S. and
Dutrey A. , A&A , (2005) 945.[92] Rodriguez L. F., Loinard L., D’Alessio P., Wilner D. and
P.T.P. H. , ApJ , (2005) L133.[93] Clarke C., Lodato G., Melnikov S. Y. and
Ibrahimov M. A. , MNRAS , (2005)942.[94] Kenyon S. J. and
Hartmann L. , ApJ , (1991) 664.[95] Millan-Gabet R. et al. , ApJ , (2006) 547.[96] Hartmann L. and
Kenyon S. J. , ApJ , (1985) 462.[97] Kenyon S. J., Hartmann L. and
Hewett R. , ApJ , (1988) 231.[98] Lodato G. and
Bertin G. , A&A , (2003) 1015.[99] Bell K. R., Lin D. N. C., Hartmann L. W. and
Kenyon S. J. , ApJ , (1995)376.[100] Clarke C. J., Lin D. N. C. and
Pringle J. E. , MNRAS , (1990) 439.[101] Vorobyov E. I. and
Basu S. , ApJ , (2005) L137.[102] Vorobyov E. I. and
Basu S. , ApJ , (2006) 956.[103] Armitage P. J., Livio M. and
Pringle J. E. , MNRAS , (2001) 705.[104] Palla F. and
Stahler S. W. , ApJ , (1993) 288.[105] Cesaroni R., Galli D., Lodato G., Walmsley C. M. and
Zhang Q. , in proc. of
Protostars and Planets V , edited by
Reipurth B., Jewitt D. and
Keil K.
Chini R. et al. , Nature , (2004) 155.[107] Torrelles J. M., G´omez J. F., Garay G., Rodr´ıguez L. F., Curiel S., CohenR. J. and
Ho P. T. P. , ApJ , (1998) 262.[108] Greenhill L. J., Gezari D. Y., Danchi W. C., Najita J., Monnier J. D. and
Tuthill P. G. , ApJ , (2004) L57.[109] Patel N. A., Curiel S., Sridharan T. K., Zhang Q., Hunter T. R., Ho P. T. P.,Torrelles J. M., Moran J. M., G´omez J. F. and
Anglada G. , Nature , (2005)109.[110] Edris K. A., Fuller G. A., Cohen R. J. and
Etoka S. , A&A , (2005) 213.[111] Mayor M. and
Queloz D. , Nature , (1995) 355.[112] Udry S., Bonfils X., Delfosse X., Forveille T., Mayor M., Perrier C., BouchyF., Lovis C., Pepe F., Queloz D. and
Bertaux J.-L. , A&A , (2007) L43.[113] Guillot T. , Planetary and Space Sciences , (1999) 10.[114] Pollack J. B. et al. , Icarus , (1996) 62.[115] Lissauer J. J. and
Stevenson D. J. , in proc. of
Protostars and Planets V , edited by
Reipurth B., Jewitt D. and
Keil K.
Armitage P. J. , ArXiv Astrophysics e-prints , (2007) .[117] van Boekel R. et al. , A&A , (2005) 189.[118] Apai D., Pascucci I., Bouwman J., Natta A., Henning T. and
Dullemond C. P. , Science , (2005) 834.[119] Bertin G. and
Cava A. , A&A , (2006) 333.[120] Haghighipour N. and
Boss A. P. , ApJ , (2003) 996.[121] Haghighipour N. and
Boss A. P. , ApJ , (2003) 1301.[122] Boss A. P. , ApJ , (2000) L101. GIUSEPPE LODATO [123]
Boss A. P. , ApJ , (2002) 462.[124] Mayer L., Quinn T., Wadsley J. and
Stadel J. , Science , (2002) 1756.[125] Durisen R. H., Boss A. P., Mayer L., Nelson A. F., Quinn T. and
Rice W. K. M. ,in proc. of
Protostars and Planets V , edited by
Reipurth B., Jewitt D. and
Keil K.
Charbonneau D., Brown T. M., Burrows A. and
Laughlin G. , in proc. of
Protostarsand Planets V , edited by
Reipurth B., Jewitt D. and
Keil K.
Fortney J. J., Saumon D., Marley M. S., Lodders K. and
Freedman R. S. , ApJ , (2006) 495.[128] Rafikov R. , ApJ , (2005) 69.[129] Boss A. P. , ApJ , (2004) 456.[130] Rafikov R. R. , ApJ , (2007) 642.[131] Boffin H. M. J., Watkins S. J., Bhattal A. S., Francis N. and
Whitworth A. P. , MNRAS , (1998) 1189.[132] Watkins S. J., Bhattal A. S., Boffin H. M. J., Francis N. and
Whitworth A. P. , MNRAS , (1998) 1205.[133] Watkins S. J., Bhattal A. S., Boffin H. M. J., Francis N. and
Whitworth A. P. , MNRAS , (1998) 1214.[134] Mayer L., Wadsley J., Quinn T. and
Stadel J. , MNRAS , (2005) 641.[135] Lodato G., Meru F., Clarke C. J. and
Rice W. K. M. , MNRAS , (2007) 590.[136] Boss A. P. , ApJ , (2006) 1148.[137] Chauvin G., Lagrange A.-M., Dumas C., Zuckerman B., Mouillet D., Song I.,Beuzit J.-L. and
Lowrance P. , A&A , (2004) L29.[138] Mamajek E. E. , ApJ , (2005) 138.[139] Lodato G., Delgado-Donate E. and
Clarke C. J. , MNRAS , (2005) L91.[140] Payne M. J. and
Lodato G. , MNRAS , (2007) 1597.[141] Greenhill L. J., Gwinn C. R., Antonucci R. and
Barvainis R. , ApJ , (1996)L21.[142] Collin-Souffrin S. and
Dumont A. M. , A&A , (1990) 292.[143] Goodman J. , MNRAS , (2003) 937.[144] King A. R. and
Pringle J. E. , MNRAS , (2007) L25.[145] Thompson T. A., Quataert E. and
Murray N. , ApJ , (2005) 167.[146] Bertin G. and
Lodato G. , A&A , (2001) 342.[147] Shlosman I., Frank J. and
Begelman M. C. , Nature , (1989) 45.[148] Collin S. and
Zahn J. P. , A&A , (1999) 433.[149] Collin S. and
Zahn J.-P. , ArXiv e-prints , (2007) .[150] Nayakshin S., Cuadra J. and
Springel V. , MNRAS , (2007) 21.[151] Nayakshin S. , MNRAS , (2006) 143.[152] Genzel R. et al. , ApJ , (2003) 812.[153] Ghez A. M. et al. , ApJ , (2005) 744.[154] Levin Y. and
Beloborodov A. M. , ApJ , (2003) L33.[155] Paumard T., Genzel R., Martins F., Nayakshin S., Beloborodov A. M., LevinY., Trippe S., Eisenhauer F., Ott T., Gillessen S., Abuter R., Cuadra J.,Alexander T. and
Sternberg A. , ApJ , (2006) 1011.[156] Nayakshin S. and
Sunyaev R. , MNRAS , (2005) L23.[157] Levin Y. , MNRAS , (2007) 515.[158] Herrnstein J. R., Greenhill L. J. and
Moran J. M. , ApJ , (1996) L17.[159] Papaloizou J. C. B., Terquem C. and
Lin D. N. C. , ApJ , (1998) 212.[160] Greenhill L. J., Kondratko P. T., Lovell J. E. J., Kuiper T. B. H., MoranJ. M., Jauncey D. L. and
Baines G. P. , ApJ , (2003) L11.[161] Kondratko P. T., Greenhill L. J. and
Moran J. M. , ApJ , (2005) 618.[162] Hur´e J.-M. , A&A , (2002) L21.[163] Fan X. et al. , AJ , (2004) 515.[164] Fan X. et al. , AJ , (2006) 1203.[165] Haiman Z. and
Loeb A. , ApJ , (1998) 505. ELF-GRAVITATING ACCRETION DISCS [166] Volonteri M., Madau P., Quataert E. and
Rees M. J. , ApJ , (2005) 69.[167] Wyithe J. S. B. and
Loeb A. , ApJ , (2005) 910.[168] Abel T., Bryan G. L. and
Norman M. L. , ApJ , (2000) 39.[169] Bromm V., Coppi P. S. and
Larson R. B. , ApJ , (2002) 23.[170] Volonteri M. and
Rees M. J. , ApJ , (2005) 624.[171] King A. R., Lubow S. H., Ogilvie G. I. and
Pringle J. E. , MNRAS , (2005) 49.[172] Lodato G. and
Pringle J. E. , MNRAS , (2006) 1196.[173] King A. R. and
Pringle J. E. , MNRAS , (2006) L90.[174] Haehnelt M. G. and
Rees M. J. , MNRAS , (1993) 168.[175] Umemura M., Loeb A. and
Turner E. L. , ApJ , (1993) 459.[176] Loeb A. and
Rasio F. A. , ApJ , (1994) 52.[177] Eisenstein D. J. and
Loeb A. , ApJ , (1995) 11.[178] Bromm V. and
Loeb A. , ApJ , (2003) 34.[179] Koushiappas S. M., Bullock J. S. and
Dekel A. , MNRAS , (2004) 292.[180] Begelman M. C., Volonteri M. and
Rees M. J. , MNRAS , (2006) 289.[181] Lodato G. and
Natarajan P. , MNRAS , (2006) 1813.[182] Lodato G. and
Natarajan P. , MNRAS , (2007) L64.[183] Warren M. S., Quinn P. J., Salmon J. K. and
Zurek W. H. , ApJ , (1992) 405.[184] Volonteri M., Lodato G. and
Natarajan P. , ArXiv e-prints , (2007) .[185] Fromang S., Terquem C. and
Balbus S. A. , MNRAS , (2002) 18.[186] Fromang S., Balbus S. A. and
De Villiers J.-P. , ApJ , (2004) 357.[187] Fromang S., Balbus S. A., Terquem C. and
De Villiers J.-P. , ApJ , (2004)364.[188] Cai K., Durisen R. H., Michael S., Boley A. C., Mej´ıa A. C., Pickett M. K. and
D’Alessio P. , ApJ ,636