Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed1February2008 (MNLATEXstylefilev1.4) Episodic accretion in magnetically layered protoplanetary discs ⋆ Philip J. Armitage1 , Mario Livio2 and J.E. Pringle2,3 1Max-Planck-Institut fu¨rAstrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany 2Space Telescope Science Institute, 3700 San MartinDrive, Baltimore, MD21218, USA 3Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK 1 0 1February2008 0 2 n ABSTRACT a We study protoplanetary disc evolution assuming that angular momentum transport J isdrivenbygravitationalinstabilityatlargeradii,andmagnetohydrodynamic(MHD) 6 turbulence in the hot inner regions. At radii of the order of 1 AU such discs develop 1 a magnetically layered structure, with accretion occurring in an ionized surface layer overlyingquiescentgasthatistoocooltosustainMHDturbulence.Weshowthatlay- 1 ereddiscsaresubjecttoalimitcycleinstability,inwhichaccretionontotheprotostar v occurs in ∼ 104 yr bursts with M˙ ∼ 10−5M⊙yr−1, separated by quiescent intervals 3 5 lasting∼105yrwhereM˙ ≈10−8M⊙yr−1.Suchburstscouldleadtorepeatedepisodes 2 ofstrongmassoutflow inYoung Stellar Objects.The transitionto this episodic mode 1 of accretion occurs at an early epoch (t ≪ 1 Myr), and the model therefore predicts 0 that many young pre-main-sequence stars should have low rates of accretion through 1 the inner disc. At ages of a few Myr, the discs are up to an order of magnitude more 0 massive than the minimum mass solar nebula, with most of the mass locked up in / the quiescent layer of the disc at r ∼1 AU. The predicted rate of low mass planetary h p migration is reduced at the outer edge of the layered disc, which could lead to an - enhanced probability of giant planet formation at radii of 1 – 3 AU. o r Key words: accretion,accretiondiscs—MHD—stars:pre-main-sequence—stars: t s formation — solar system: formation — planets and satellites a : v i X 1 INTRODUCTION at larger radii of r ∼1 AU, where the temperature is typi- r a cally a few hundred K, magnetic field instabilities are sup- Thestructureandevolutionofprotoplanetarydiscsdepend pressedbythelowionizationfraction(Matsumoto&Tajima upontherateat whichgascansheditsangularmomentum 1995; Gammie 1996; Gammie & Menou 1998; Livio 1999; and thereby flow inwards. Two widely applicable physical Wardle1999; Sano&Miyama 1999; Sanoet al. 2000). This mechanisms are known to lead to the required outward an- leads (Gammie 1996) to the formation of a layered disc gular momentum transport. If the gas is coupled to a mag- structure, in which the gas near the disc midplane is cold, neticfield,instabilities that inevitablyarise indifferentially shielded from ionizing high energy radiation, and quiescent rotatingdiscs(Balbus&Hawley1991;Chandrasekhar1961; (non-turbulent).Turbulence and accretion occurs only in a Velikhov 1959) lead to turbulence and angular momentum thinsurfacelayerthatisionizedbycosmicrays.Movingstill transport (Stoneet al. 1996; Brandenburget al. 1995; for a further outwards the entire thickness of the disc again be- reviewseee.g.Hawley&Balbus1999).Ifthediscismassive comeviscous,eitherattheradiuswherethesurfacedensity enough, gravitational instability leads to additional trans- issmallenoughforcosmicraystopenetratetothemidplane, port (Toomre 1964; Laughlin & Bodenheimer 1994; Nelson orwheretheonsetofdiscself-gravityprovidesanalternative et al. 1998; Pickett et al. 2000). non-magnetic source of angular momentum transport. Applying these findings to the construction of proto- Thepredictionsofastaticlayereddiscmodelfortheac- planetarydiscmodelsleadstothestructureshownschemat- cretionrateandspectralenergydistributionofTTauristars ically in Fig. 1 (after Gammie 1996). In the inner disc, were discussed by Gammie (1996), and are broadly consis- MHD turbulence transports angular momentum. However, tentwithobservations(e.g.withtheaccretionrateforClas- sicalTTauristarsmeasuredbyGullbringetal.1998).Inthis ⋆ Presentaddress:School of Physicsand Astronomy,University paper we consider the evolution of the layered disc, which ofStAndrews,NorthHaugh,StAndrewsKY169SS cannot be in a steady state (Gammie 1996, 1999; Stepinski (cid:13)c 0000RAS 2 P.J. Armitage, M. Livio & J.E. Pringle n κ0 a b Tmax /K Inner disc Layered disc Self−gravitating disc 1 1.0×10−4 0 2.1 132 2 3.0×100 0 -0.01 170 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 3 1.0×10−2 0 1.1 377 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 4 5.0×104 0 -1.5 389 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 5 1.0×10−1 0 0.7 579 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) 6 2.0×1015 0 -5.2 681 Quiescent layer 7 2.0×10−2 0 0.8 964 8 2.0×1081 1 -24 1572 T>Tcrit T<Tcrit T<Tcrit 9 1.0×10−8 2/3 3 3728 r~0.1 AU r~2 AU 10 1.0×10−36 1/3 10 10330 11 1.5×1020 1 −5/2 45060 Figure1.Illustrationoftheradialstructureofthediscinthelay- 12 0.348 0 0 — ereddiscmodel. Atsmallradii,the central temperature exceeds Tcrit,thetemperatureabovewhichMHDturbulenceprovidesan efficientsourceofangularmomentumtransport.Intheintermedi- Table1.Opacityregimesinascendingorderoftemperature,fit- ate,layeredregion,thecentraltemperatureistoolowforthedisc tedbyanalyticfunctionsoftheformκ=κ0ρaTb.Wehaveused tosupportMHDturbulence. Accretionoccurs throughasurface fits providedbyBell& Lin(1994), as modifiedforlow tempera- layerwhichiskeptionizedbycosmicrays,whilethediscnearthe tures byBellet al.(1997). Themaximumtemperature Tmax for midplaneisquiescent. Atlargerradii,theentirethickness ofthe each regime is quoted for a typical disc density of 10−9 gcm−3 disc is again viscous, with angular momentum transport being (wherespecificationofthedensityisnecessary). drivenbyself-gravity. Here c is the disc specific heat, which for temperatures 1999), andexaminetheimplications for theoutflowhistory p T ∼ 103 K is given by c ≃ 2.7R/µ, where R is the gas of young stars and for the predicted disc mass. The most c p significantchangestothediscstructureoccurattheradiiof constant and µ = 2.3 is the mean molecular weight. Q+ representslocalheatingduetoviscousdissipation,givenby, greatestinterestforplanetformation(Reyes-Ruiz&Stepin- ski1995), andwediscusstheimplications forthemigration Q+= 9 νΣΩ2 (5) oflowmassplanets,andfortheeccentricityofmassiveplan- (cid:16)8(cid:17) ets interacting with thedisc. if theentire disc is viscous and 9 Q+= 8 νΣlayerΩ2 (6) (cid:16) (cid:17) 2 LAYERED PROTOPLANETARY DISC EVOLUTION otherwise. For Q−, the local cooling rate, we assume that eachannulusofthediscradiates asablackbodyat temper- 2.1 Equations atureT , so that e Describing the evolution of the surface density Σ(r,t) and Q− =σTe4, (7) midplanetemperatureT (r,t)ofalayereddiscrequiresonly c where σ is the Stefan-Boltzmann constant. Finally, we in- minor modifications to the usual time-dependent equations clude an advective term in the energy equation, which de- forthinaccretiondiscs.Wedenotethesurfacedensityofthe pendson thevertically averaged radial velocity, ‘active’ (viscous) disc by Σ . If, a 3 ∂ Tc >Tcrit (1) vr =−Σr1/2∂r νΣr1/2 , (8) (cid:0) (cid:1) or, and theradial temperature gradient. Σ<2 Σ (2) layer thenthediscisviscousthroughoutitsthicknessandΣ =Σ. 2.2 Vertical structure a Otherwise only the surface layers are viscous and Σ = a Completingthemodelrequiresspecificationofboththevis- 2Σ . The values of these parameters are determined by layer cosity ν and the vertical structure, which sets the relation therequirementthat thedisc besufficiently ionized tosup- between the central temperature T and the surface tem- port MHD turbulence (Gammie 1996). We adopt Tcrit = peratureT .Weadopt thesimplest,cvertically averaged ap- 800 K, and Σ = 102 gcm−2. For a Keplerian disc, the e layer proach, for which, angular velocity is Ω= GM∗/r3, where M∗ is the stellar 3 mass. The surface densitpy evolution is then described by, T4= τT4, (9) c 8 e ∂Σ 3 ∂ ∂ ∂t = r∂r r1/2∂r νΣar1/2 +Σ˙(r,t), (3) where τ = (Σ/2)κ is the optical depth for a given opacity h (cid:0) (cid:1)i κ(ρ ,T ).Whenanannulusmakesthetransition tothelay- where ν is the kinematic viscosity and Σ˙(r,t) is the rate of c c ered state, we crudely account for this by replacing (Σ/2) change of the surface density dueto infall onto thedisc. in theexpression for τ by Σ . Notethat this meansthat layer For the energy equation, we adopt a simplified form of wedonotattempttotreattheverticalstructureduringthe that used byCannizzo (1993), transition consistently. ∂Tc = 2(Q+−Q−) −v ∂Tc. (4) Analytic expressions for low temperature Rosseland ∂t c Σ r ∂r meanopacitiesaregivenbyBelletal.(1997).Thebehaviour p (cid:13)c 0000RAS,MNRAS000,000–000 Evolution of layered protoplanetary discs 3 of the disc depends primarily on the opacity near the tran- ‘α’, νgrav =αgravc2s/Ω, with sition temperature Tcrit, for which thefit is, Q2 κ=2×10−2Tc0.8 cm2 g−1. (10) αgrav = 0.01(cid:18) Qcr2it −1(cid:19) if Q<Qcrit ThefulllistofopacitiesusedisquotedinTable1.Thesefits αgrav = 0 otherwise. (13) havebeentakenfrom Bell&Lin (1994), withthemodifica- Operationally, we solve equations (3) and (4) using an ex- tions for low temperatures quoted in Bell et al. (1997). plicit finite difference method, and include the self-gravity asanadditionalsweepoverthegridusingνgrav andthefull surface density Σ in the diffusive part of equation (3), and 2.3 Viscosity as an additional heating term in theenergy equation. Weadopt an alpha prescription (Shakura& Sunyaev1973) for the viscosity, c2 3 RESULTS ν =α s, (11) Ω ForournumericalcalculationsweadoptM∗ =M⊙,aninner where cs = RTc/µ is the midplane sound speed and α a disc radius rin = 3.5×1011 cm (5 R⊙), and an outer disc dimensionlespsparametermeasuringtheefficiencyofangular radius rout = 6.0×1014 cm (40 AU). The computational momentum transport. Numerical simulations of MHD tur- grid of 120 radial mesh points is uniform in a scaled radial bulence in discs suggest that α ≈ 10−2 (Stone et al. 1996; variableX ∝r1/2.Azero torqueboundarycondition isim- Brandenburget al. 1995). posedatrin,whileatrout wepreventoutflowbysettingthe radial velocity, v , to zero. When mass is added to the disc r (representing infall of further material from the molecular 2.4 Self-gravity cloud),we takeΣ˙ to beaGaussian centeredat 10 AUwith a width of 1 AU. This is exterior to the layered section of An immediate consequence of the layered disc structure is the disc that we are interested in studying, so the details that the accretion rate through the surface layer is a fixed, of how mass is added are unimportant here. We add mass increasing function of radius (Gammie 1996). Inevitably, assuming that it has the same specific angular momentum therefore, material builds up in the quiescent layer over and temperature as thelocal disc material. time, until the surface density becomes high enough that self-gravity of the disc becomes important. The condition for this is that the Toomre Q parameter (Toomre 1964), 3.1 Steady infall models given by, To demonstrate theevolution of thelayered disc model, we Q= csΩ , (12) first take the infall rate onto the outer disc M˙infall to be πGΣ a constant. We take α = 10−2, as suggested by numerical is less than some critical value Qcrit, which is of the order simulations of MHD disc turbulence, and allow the disc to of unity. Here, we adopt Qcrit = 2. When Q < Qcrit, non- evolve from an initially low mass (10−2 M⊙) until either a axisymmetric instabilities set in and act to transport angu- steady state or a limit cycle has developed. lar momentum outwards. In doing so, the surface density Figure 2 shows the accretion rate onto the star for a evolvesinsuchawaythatthestabilityofthediscgenerally range of mass infall rates between 5×10−7 M⊙yr−1 and increases, thereby self-regulating the violence of the insta- 3×10−6 M⊙yr−1.Forthehighestinfallrate,thediscisable bility. totransport gas steadily ontothestar. Theinnerregion, in In a local approximation, the efficiency of angular mo- which Tc > Tcrit, is able to be supplied by an outer disc in mentum transport due to self-gravity can be determined which Q<Qcrit and theangular momentumis transported by the requirement that viscous heating balances radiative byself-gravity.Qistypicallyaround1.5inthisregion,with losses (Gammie 1999). This approach, in which the cool- thetransition occurring at a radius of ≈3 AU. ingtimedirectlydeterminesα,ismostappropriatefordiscs For lower infall rates such steady accretion is not pos- whose thermal balance is determined by viscous heating. sible.Afteraninitialphaseinwhichthediscmassincreases Even in this limit, global effects may be important (Balbus steadily, a limit cycle is obtained in which outbursts, with & Papaloizou 1999). A global redistribution of surface den- a peak accretion rate M˙ ≈ 10−5 M⊙yr−1, alternate with sity, driven by self-gravity, can also occur even if the disc quiescent intervals where M˙ ≈10−8 M⊙yr−1. Both thedu- temperature is set externally (for example by irradiation). ration of the outbursts, and the recurrence time, vary with This is the limit that has been studied extensively using therateatwhichgasissuppliedtotheouterdisc.However, isothermal or polytropic simulations (e.g. Laughlin & Bo- to order of magnitude, for M˙infall <∼ 10−6 M⊙yr−1 we find denheimer 1994). For protoplanetary discs, models suggest toutburst ∼104 yr, and trecur ∼105 yr. that both viscous heatingand irradiation can beimportant The mechanism for these outbursts is shown in Figure at different epochs and radii (Bell et al. 1997). 3,which plotsthecentraltemperatureduringthequiescent Given the above complications, it is evident that at- intervalleadinguptoanoutburst.Following theendofone tempts to include theeffects of self-gravity into one dimen- outburst,thereisanextendedlayeredregion ofthediscbe- sional disc models can only be approximate. We adopt the tween ≈ 0.3 AU and ≈ 3 AU. This region is stable against form suggested by Lin & Pringle (1987, 1990), and use a gravitational instability, and supports only the small qui- variantofequation(11)alongwithaneffectivegravitational escent accretion rate. Gas flowing inwards from the more (cid:13)c 0000RAS,MNRAS000,000–000 4 P.J. Armitage, M. Livio & J.E. Pringle Figure 2. Protostellar accretion for models with a constant Figure 3.Evolutionofthedisccentraltemperature duringqui- rate of infall onto the outer disc. From top down, the pan- escence. An outburst is triggered when the central temperature els show models with M˙infall = 3×10−6 M⊙yr−1, M˙infall = inthelayeredregionfirstexceedsTcrit(dashedline),thetemper- 1.5×10−6M⊙yr−1,andM˙infall=5×10−7M⊙yr−1respectively. ature at which the disc is well enough ionized to support MHD Allmodelsareshownafterinitialtransientshavedecayed.Steady turbulence. The curves show temperature profiles plotted at in- accretion occurs forsufficiently highaccretion rates through the tervals of 104 yr, with the earliest timeslice having the lowest outer disc. At lower accretion rates, the mass flow onto the star temperature at 2 AU. The final timeslice, shown as the upper, isstronglytime-dependent. boldcurve,immediatelyprecedes thestartofanoutburst. activeself-gravitatingregionofthediscthusaccumulatesin tively large disc masses are required for the limit cycle to thequiescentregion,andthetransitionradiuswherethedisc operate. In our models, Mdisc fluctuates between 0.2 and becomes self-gravitating moves inward. At the same time, 0.3 M⊙. During an outburst, only a small fraction of the the greater dissipation rate in the self-gravitating disc cre- disc mass — around 20% — is accreted. In more realistic ates alocal peakin thedisctemperature, which also moves evolutionary models, in which the mass infall rate declines inwards. A new outburst is triggered when this peak first withtime,thediscwillthereforebeabletoproduceahand- exceedsTcrit,therebysatisfyingtheconditionforMHDtur- ful of outbursts, with a few hundredths of a solar mass of bulence to be restarted. Once the outburst begins thermal gas accreted during each, before settling into a permanent energy is rapidly advected inwards, triggering a transition quiescent state. of theentire inner disc into a hot turbulent state with high accretion rate. This continues until the reservoir of accu- mulatedmasshasbeendepleted,whereuponacoolingwave 3.2 Relation to FU Orionis events sweepsthroughthediscandreturnsittothequiescent,lay- Themoststrikingvariabilityobservedinpre-main-sequence ered, state. stars occurs in FU Orionis events (e.g. Kenyon 1995; Hart- The derivedoutburstdurations areconsistent with the mann&Kenyon1996),whicharelargeamplitudeoutbursts expected viscous timescale in protoplanetary discs at radii in the system luminosity originating in the accretion disc. of a few AU. For an outburst triggered at a radius of ap- Thestatisticsoftheseeventsaresubjecttoconsiderableun- proximately 2 AU, the viscous timescale once the disc is in certainty,butthepeakaccretion rateduringoutburstsisof outburst,given by the order of 10−4 M⊙yr−1, while the duration is generally r2 1 h −2 thought to bearound 102 yr. tν ≃ ν = αΩ r , (14) Themost popularexplanation for FUOrionis eventsis (cid:16) (cid:17) in terms of a disc thermal instability (e.g. Bell & Lin 1994; where h is thedisc scale height, is approximately, Kley &Lin 1999, andreferences therein),akin to thatused α −1 r 3/2 h/r −2 tomodeldwarfnovae(e.g.Cannizzo1993).Intheprotostel- tν ≈2×104(cid:16)10−2(cid:17) (cid:16)2 AU(cid:17) (cid:18)0.05(cid:19) yr. (15) larcase,athermalinstabilitycanonlyoperateatextremely smalldiscradii,typicallyatlessthan0.1AU.Matchingthe This is only a rough estimate, but it is of the correct order observed timescales then requires small values of α (10−3 of magnitude to match thenumerical model. to10−4),withcorrespondingly largevaluesforthediscsur- The involvement of disc self-gravity means that rela- face density. As noted by Gammie (1999), it would be at- (cid:13)c 0000RAS,MNRAS000,000–000 Evolution of layered protoplanetary discs 5 tractive to dispense with the need for a thermal instabil- ity by invoking an alternative limit cycle of a layered disc. Thismightbepossibleifthetransitiontotheoutburststate was triggered at the extreme inner edge of the layered re- gion. Whetherthishappensdependson wherein thequies- centlayermasspreferentiallyaccumulates,andtheoutcome is therefore sensitive to the detailed, and rather uncertain, physicsofthelayereddisc.Inourmodel,inwhichmatteris beingadded to thelayered region primarily at large radius, thetriggeringoccursattheouteredgeofthelayeredregion, and the timescales are inconsistent with those of FU Ori- onis events. We note, however, that the accretion rates of (1−10)×10−5M⊙yr−1 thatweobtainduringouroutbursts fall into the thermally unstable regime. The long outbursts obtained here could then feed shorter FU Orionis event of thesort calculated by Bell & Lin (1994). 3.3 Protostellar accretion history Continuousreplenishmentofthediscmassfrominfallisnot a realistic model for protostellar accretion. To explore how the accretion rate onto the star may evolve with time, we adoptasimplemodelfortheinfallofgasontothedisc.The detailsaremodel-dependent(seee.g.Larson1969;Shu1977; Figure 4. Accretion rate and disc mass as a function of time. Basu 1998), butat early times theinfall rate isexpectedto The solid curves are for layered disc models in which the mass be of the order of M˙ ∼ c3/G, where c is the sound infall s s infallrateontotheouterdiscdeclinesexponentially.Forcompar- speed in the collapsing cloud. For cloud temperatures T ∼ isonwiththese,thedashedcurvesshowtheevolutionofamodel 10K,thisimpliesaninfallrateof∼10−5 M⊙yr−1.Wetake without a layered structure. This model has identical parame- an initial infall rate of 2×10−5 M⊙yr−1, and assume that ters,exceptthatthediscisassumedtobeviscousatallradiiand this declines exponentially on the free fall timescale, t = centraltemperatures. ff (3π/32Gρ )1/2,whereρ istheclouddensity.Wetake cloud cloud t =105yr.Otherinputparametersareasdescribedearlier. ff timescales of 104 yr or greater. These longer timescales are With theseparameters, Fig. 4 showstheaccretion rate characteristicofprocessesoccurringatlargerdiscradiithan forthelayereddiscmodel.Forcomparison,wealsoshowthe thermalinstabilities,andareconsistentwithinstabilities,of resultsfromaviscousdiscmodelwhichhasidenticalparam- the kind discussed here, originating in the layered disc re- eters, but which has no layered region. This control model gion. would be appropriate, for example, if angular momentum Theoutburst accretion rateinthelayereddiscmodelis transportwasdrivenbyapurelyhydrodynamicmechanism several orders of magnitude higher than that in theequiva- whose efficiency was independent of the disc temperature. lentfullyviscousdiscmodelatt∼106yr.However,thetime Wecomputethecontrol model byusing thesame code and averagedaccretionrateissubstantiallylower,sothatatlate infall history, butwith Tcrit set tozero. times the disc mass in the layered model substantially ex- Initially, when the infall rate is higher than around 2×10−6 M⊙yr−1, both models support an accretion rate ceedsthatintheviscousmodel.WefindthatMdisc>∼0.1M⊙ after 1 – 2 Myr. By this time the mass of the fully viscous onto the star which tracks that with which matter is be- ing added to the disc. Subsequently, accretion through the dischasdroppedtoonly10−2 M⊙.AsshowninFig.5,most of this extra mass is tied up in the quiescent layer at radii layered disc occurs in outbursts with properties similar to of theorder of 1 AU,where thesurface density exceeds the thosedescribedearlier,butwithincreasingintervalsbetween minimum mass solar nebula estimate (Hayashi, Nakazawa events.Forthismodel,thelast outburstbefore thediscbe- & Nakagawa 1985) byabout two orders of magnitude. came permanently quiescent occurred after 2 Myr. The ac- cretionrateontothestarinthecontrolmodel,ontheother hand, declines smoothly with time. Similar results are ob- 3.4 Compatibility with observations of tained using other functional forms for the mass infall rate protoplanetary discs as a function of time. Thetimescalesofvariabilitypredictedbythemodelare Adiscmassofafew tenthsofasolar massatearly timesis much longer than any direct observational record. Indirect consistent with the upperend of the disc mass distribution evidence for long timescale variability of young stars, how- inferred from mm-wavelength observations (Beckwith et al. ever,isprovidedbystudiesofprotostellaroutflows,sincethe 1990;Osterloh&Beckwith1995).However,thelayereddisc rate of mass outflow via jets is widely believed to track the modelpredictsthatasubstantialmass–atleastafewhun- accretionrate.Forexample,observationsbyReipurth,Bally dredths of a solar mass – remains in the disc throughout & Devine (1997) suggest that large amplitude variations in most of the typical Classical T Tauri disc lifetime (Strom the outflow rate occur not only on the 102 yr timescales 1995). Although recent measurements of the mass of H2 in characteristicofFUOrionisevents,butalsoonmuchlonger debris discs indicate that several Jupiter masses of gas can (cid:13)c 0000RAS,MNRAS000,000–000 6 P.J. Armitage, M. Livio & J.E. Pringle Classical T Tauri stars (Gullbring et al. 1998). There have alsobeensuggestionsthattheobservationallyinferredaccre- tion rate declines with increasing stellar age (Hartmann et al.1998).Althoughcurrentlythescatterinmeasuredaccre- tion rates at agiven assumed age is very large, and theage estimatesthemselvesaresubjecttosubstantialuncertainties (Tout, Livio & Bonnell 1999), this observation poses diffi- cultiesfor thesimplest layered discmodel(Stepinski1998). We note, however, that changes in the quiescent accretion ratewithtimewouldbeexpectedinmorecompletemodels, for example if the opacity from dust decreased with time dueto grain growth. Including heating of thedisc from the changingstellarradiationwouldalsoleadtoadeclineinthe quiescent accretion rate. 4 PLANET FORMATION The environment which a layered disc presents for planet formation differsin severalrespectsfrom bothpreviousvis- cous disc models (e.g. Lin & Pringle 1990), and from static modelssuchastheminimummasssolarnebula.Atradiibe- tween a few tenths of an AU, and several AU, the gas near thediscmidplaneisalmostalwaysquiescent,unlesstheset- Figure 5. Cumulative disc mass as a function of (upper panel) tlingofcmsizedparticlesitselfdrivesadditionalinstabilities surfacedensityΣand(lowerpanel)radius.Thesolidcurveshows the layered disc model at t= 1.5×106 yr. For comparison, the (Goldreich&Ward1973;Cuzzi,Dobrovolskis&Champney dashed curve shows results for the fully viscous model at the 1993). The low viscosity of the disc in this region would same time. The vertical dashed line in the upper panel shows reduce the rate at which massive planets migrate inwards, the approximate surface density below which the disc would be and may allow planet-disc interactions to excite significant optically thin at λ = 1 mm. Essentially all the additional mass eccentricity (Papaloizou, Nelson & Masset 2001). present in the layered disc is locked up in highly optically thick Otherdifferencesincludethediscatlargeradiiremain- regions at small radii. The amount of mass at r > 10 AU is ing mildly unstable to gravitational instability for a rela- comparableinboththelayeredandthefullyviscousmodel. tivelylongtime–aroundaMyr–becausethereducedtime- averaged accretion rate leads to a larger disc mass at late times when compared to viscous disc models. Outbursts, persist beyond 10 Myr (Thi et al. 2001), our predicted disc which persist for a substantial fraction of the disc lifetime, masses areinapparentcontradiction withmm observations mean that much of the inner disc is expected to be subject ofmanysystemswithdiscmassesestimatedat10−3 M⊙ or to rapid heating and cooling episodes. However, we find no smaller. evidencethatthesecoolingwavesdrivethediscintoastate To check whether this is a serious problem for the lay- where gravitational collapse (Q <∼ 1) would occur. ered disc model, we plot in Fig. 5 the cumulative distribu- The survival of solid bodies in protoplanetary discs is tion of disc mass as a function of surface density. Most of limited by the rate at which they migrate inwards relative the ‘additional’ mass in the layered disc, compared to the tothegas. Migration can berapidbothforcm-sizedbodies equivalent fully viscous control model, is locked up in the (e.g. Godon & Livio 1999, and references therein), and is quiescentlayerat small radiiandhigh surface density.This particularlyproblematicforlowmassplanets,forwhichmi- material is highly optically thick, and thus would not con- gration occurs due to the influence of gravitational torques tribute to the observed emission at mm wavelengths. For (Goldreich & Tremaine 1979). In standard disc models, the a range of dust models, the opacity of dust at λ = 1 mm migration timescale for Earth mass planets at a few AU is thought (Men´shchikov & Henning 1997; Beckwith et al. can be as short as 104−105 yr, which leaves little time to 1990) tobeintherange0.3cm2g−1 <∼κ<∼10cm2g−1 (note assemblethecoresofgiantplanetsbeforetheputativebuild- that this is per gram of dust). For a gas to dust ratio of ing blocks are consumed by the star. The rate of migration 100, this implies that an optical depth of unity at 1 mm is depends only very weakly upon the surface density profile attained forΣ∼102 gcm−2.Atthissurfacedensitythresh- (Ward 1997), but is sensitive to the gradient of the central old,themassesofoptically thingasinthelayeredandnon- temperature,andcanbehaltedifthereexistsaregionofthe layered disc models are very similar. The masses of proto- disc where T ∝ r. During the quiescent phase, the models c planetarydiscsinferredfrommm-wavelengthmeasurements we have discussed here possess such a non-monotonic cen- couldthussubstantiallyunderestimatethetruediscmasses. tral temperature profile at radii of 1-3 AU (Fig. 3), leading The accretion rates inferred for Classical T Tauri stars tothepossibility ofanenhancedprobabilityofgiantplanet canalsoprovideatest.Inthesimplestversionofthelayered formationatthoseradii(Papaloizou&Terquem1999).Cur- disc model, the accretion rate during quiescence does not rent radial velocity surveys are sensitive to massive planets changeoverthelifetimeoftheprotostellarphase.Thisrateis at radii r <∼ 3 AU, i.e. within this region. If migration is consistentwiththeorderofmagnitudeinferredforaccreting as rapid as currently suspected, a consequence of the lay- (cid:13)c 0000RAS,MNRAS000,000–000 Evolution of layered protoplanetary discs 7 ered disc model would probably be fewer planets at larger ACKNOWLEDGMENTS radii than would be expected from an extrapolation from PJAthankstheInstituteofAstronomyfortheirhospitality thecurrent data. during the course of part of this work, and Charles Gam- mie, John Papaloizou and Henk Spruit for helpful discus- sions. This work was partially supported by the Training and Mobility of Researchers (TMR) program of the Euro- pean Commission. ML acknowledges support from NASA Grant NAG5-6857. JEP is grateful to STScI for continuing 5 SUMMARY support undertheirVisitors Program. Inthispaper,wehavepresentedmodelsfortheevolutionof magneticallylayeredprotoplanetarydiscs.Givenourpresent understandingofangularmomentumtransportindiscs,this REFERENCES model represents the best guess for the structure of proto- BalbusS.A.,HawleyJ.F.,1991,ApJ,376,214 planetarydiscsatradiioftheorderof1AU,wherethegasis BalbusS.A.,PapaloizouJ.C.B.,1999, ApJ,521,650 coolandpoorlycoupledtothemagneticfield.Theevolution- BasuS.,1998, ApJ,509,229 ary models presented here suggest a number of important BeckwithS.V.W.,SargentA.I.,ChiniR.S.,GuestenR.,1990,AJ, differences with non-layered viscous disc models. 99,924 Bell K.R., Cassen P.M., Klahr H.H., Henning Th., 1997, ApJ, (i) The disc cannot be in a steady state for outer disc 486,372 accretion rates in the range 10−8 M⊙yr−1 <∼ M˙ <∼ 2 × BellK.R.,LinD.N.C.,1994,ApJ,427,987 10−6 M⊙yr−1. A limit cycle is obtained, in which heating Brandenburg A., Nordlund A., Stein R.F., Torkelsson U., 1995, whenthelayeredregion ofthediscbecomesself-gravitating ApJ,446,741 periodically restarts MHD turbulence, leading to outbursts CannizzoJ.K.,1993,ApJ,419,318 Chandrasekhar S., 1961, in Hydrodynamic and Hydromagnetic ofaccretionontothestar.Allcalculationstodateoflayered Stability,OxfordUniversityPress,Oxford,p.384 discs obtain strongly episodic accretion, despite variations CuzziJ.N.,DobrovolskisA.R.,ChampneyJ.M.,1993,Icarus,106, in the physical processes included in the models (Gammie 102 1999; Stepinski1999). Thisis themost robust prediction of GammieC.F.,1996,ApJ,457,355 themodel. Gammie C.F., 1999, in Astrophysical Disks, eds. Sellwood & (ii) The duration of outbursts depends on where in the Goodman,ASPConf.Series,160,p.122 quiescentdisctheyaretriggered.Inourmodel,theoutbursts GammieC.F.,MenouK.,1998,ApJ,492,75 are long, with duration ∼104yr.They could driverepeated GodonP.,LivioM.,1999,ApJ,523,350 strongepisodesofmassoutflowfromtheinnerdisc,resulting GoldreichP.,TremaineS.,1979, ApJ,233,857 in ‘pulsing’ of theobserved jets. GoldreichP.,WardW.R.,1973, ApJ,183,1051 Gullbring E., Hartmann L., Briceno C., Calvet N., 1998, ApJ, (iii) The time-averaged accretion rate is reduced due to 494,323 thebottleneckcreatedbythelayeredregionofthedisc.Asa Hartmann L., Calvet N., Gullbring E., D’AlessioP., 1998, ApJ, result,thediscmassatlatetimesislarge,typically≈0.1M⊙ 495,385 at 1-2 Myr. Most of this mass is locked up in the quiescent HartmannL.,KenyonS.J.,1996,ARA&A,34,207 layer of thedisc at small radii. HawleyJ.F.,BalbusS.A.,1999,inAstrophysicalDisks,eds.Sell- (iv) The central temperature is strongly modified, and wood&Goodman,ASPConf.Series,160,p.108 often increases with radius,nearthetransition between the HayashiC.,NakazawaK.,NakagawaY.,1985,inProtostarsand layereddiscandtheouterself-gravitatingregion.Thiscould Planets II, eds D.C. Black & M.S. Matthews, University of slow the otherwise rapid rate of low mass planetary migra- ArizonaPress,Tuscon,p.1100 tion at radii r∼1−2 AU. Kenyon S.J., 1995, Rev. Mex. Astron. Astrophys. Conf. Ser., 1, 237 (v) A low viscosity in the layered region of the disc af- KleyW.,LinD.N.C.,1999, ApJ,518,833 fects the migration rate and eccentricity evolution of mas- LarsonR.B.,1969,MNRAS,145,271 sive planets at the radii currently probed by radial velocity LaughlinG.,BodenheimerP.,1994,ApJ,436,335 surveys. LinD.N.C.,PringleJ.E.,1987, MNRAS,225,607 LinD.N.C.,PringleJ.E.,1990, ApJ,358,515 Current observations provide only very limited constraints LivioM.,1999,inAstrophysicalDisks,eds.Sellwood&Goodman, on the properties of protoplanetary discs at the radii of ASPConf.Series,160,p.33 greatest interest for planet formation. Asa result, there re- MarcyG.W.,ButlerR.P.,2000, PASP,112,137 mainssubstantialuncertaintyinourknowledgeofhowthese MatsumotoR.,TajimaT.,1995,ApJ,445,767 discs evolve. Models, such as this one, that attempt to in- Men´shchikov A.B.,HenningTh.,1997,A&A,318,879 clude more realistic disc physics, can lead to very different Nelson A.F., Benz W., Adams F.C., Arnett D., 1998, ApJ, 502, 342 environments for planet formation and migration than the OsterlohM.,BeckwithS.V.W.,1995,ApJ,439,288 highly simplified models usually considered. The most ob- PapaloizouJ.C.B.,NelsonR.P.,MassetF.,2001,A&A,inpress vious observational predictions of the model are that some PapaloizouJ.C.B.,Terquem C.,1999, ApJ,521,823 very young stars should be accreting at the low ‘quiescent’ Pickett B.K.,CassenP.,DurisenR.H.,LinkR.,2000, ApJ, 529, rate of the order of 10−8 M⊙yr−1, and that high accretion 1034 rateoutburstsshouldcontinue,albeitwithlesserfrequency, ReipurthB.,BallyJ.,DevineD.,1997, AJ,114,2708 duringasubstantialfractionoftheClassical TTauriphase. Reyes-RuizM.,StepinskiT.F.,1995,ApJ,438,750 (cid:13)c 0000RAS,MNRAS000,000–000 8 P.J. Armitage, M. Livio & J.E. Pringle SanoT.,MiyamaS.M.,1999,ApJ,515,776 Sano T., Miyama S.M., Umebayashi T., Nakano T., 2000, ApJ, 543,486 ShakuraN.I.,SunyaevR.A.,1973, A&A,24,337 ShuF.H.,1977, ApJ,214,488 StepinskiT.F.,1998, ApJ,507,361 Stepinski T.F., 1999, 30th Annual Lunar and Planetary Science Conference,abstractno.1205 StoneJ.M.,HawleyJ.F.,GammieC.F.,BalbusS.A.,1996,ApJ, 463,656 Strom S.E., 1995, Rev. Mex. Astron. Astrophys. Conf. Ser., 1, 317 ThiW.F.etal.,2001, Nature,409,60 ToomreA.,1964, ApJ,139,1217 ToutC.A.,LivioM.,BonnellI.A.,1999,MNRAS,310,360 VelikhovE.P.,1959, SovietJETP,36,995 WardW.R.,1997,Icarus,126,261 WardleM.,1999, MNRAS,307,849 (cid:13)c 0000RAS,MNRAS000,000–000