Table Of ContentMon.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