Protoplanetary Disks as (Possibly) Viscous Disks

0.42 MB·

by  Roman R. Rafikov

Draftversion January23,2017
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
PROTOPLANETARY DISKS AS (POSSIBLY) VISCOUS DISKS
Roman R. Rafikov1,2
Draft version January 23, 2017
ABSTRACT
Protoplanetary disks are believed to evolve on Myr timescales in a diffusive (viscous) manner as
7 a result of angular momentum transport driven by internal stresses. Here we use a sample of 26
1 protoplanetary disks resolved by ALMA with measured (dust-based) masses and stellar accretion
0 ratestoderivethedimensionlessα-viscosityvaluesforindividualobjects,withthegoalofconstraining
2 the angular momentum transport mechanism. We find that the inferred values of α do not cluster
around a single value, but instead have a broad distribution extending from 10−4 to 0.04. Moreover,
n
they correlate with neither the global disk parameters (mass, size, surface density) nor the stellar
a
characteristics (mass, luminosity, radius). However, we do find a strong linear correlation between α
J
and the central mass accretion rate M˙ . This correlationis unlikely to result from the direct physical
0
2 effect of M˙ on disk viscosity on globalscales. Instead, we suggest that it is caused by the decoupling
of stellar M˙ from the global disk characteristics in one of the following ways. (1) The behavior (and

] range) of α is controlled by a yet unidentified parameter (e.g. ionization fraction, magnetic field
P
strength, or geometry), ultimately driving the variation of M˙ . (2) The central M˙ is decoupled from
E
theglobalviscousmassaccretionrateasaresultofaninstabilityormassaccumulation(orloss)inthe
.
h inner disk. (3) Perhaps the most intriguing possibility is that angular momentum in protoplanetary
p disks is transported non-viscously, e.g. via magnetohydrodynamic winds or spiral density waves.
-
o Subjectheadings: accretion,accretiondisks—protoplanetarydisks—planetsandsatellites: formation
r
t
s
a 1. INTRODUCTION temperature) and Ω ≡ (GM⋆/r3)1/2 (M⋆ is the mass of
[ Protoplanetary disks are thought to persist around the central star) being the local sound speed and Kep-
their parent stars for relatively short span of time. Ob- lerianangular frequency. Viscous models invariably pre-
2 dict that on long timescales (exceeding the viscous time
servations present a clear evolutionary picture, in which
v
of the initial, more compact disk, so that r grows be-
2 boththefractionofsystemsexhibitingIRexcess(Hillen- out
yond the initial disk radius) the central mass accretion
5 brand2005)andthemassaccretionrateontothecentral
3 star(Calvetet al.2000)decline ofMyr timescales. Both rate M˙ should be related to the total disk mass Md as
2 observational indicators are thought to be the clear sig-
M˙ ≈M /t , (3)
0 natures of the presence of the circumstellar disks. d ν
. Astrophysical accretion disks are believed to evolve
1 (uptoaconstantfactoroforderunity)witht evaluated
predominantlyunder the actionoftheir internalstresses ν
0 at r , see equation (1).
(Shakura&Sunyaev1973;Lynden-Bell&Pringle1974), out
7 The idea of the viscous evolution of the protoplane-
and protoplanetary disks are no exception to the rule.
1 tary disks, diffusive in character and characterized by
: For a long time evolutionary models of the protoplane- equations (1)-(3), has gained certain observational sup-
v tary disks have been developed assuming that the disks
port. In particular, Hartmann et al. (1998) and Calvet
Xi spread viscously, losing mass to the central star, while et al. (2000) found that the observed average properties
at the same time providing the birth site for planet for-
r ofprotoplanetarydiskscanbeexplainedbytheirviscous
a mation. Thecharacteristictime forthe disk evolutionin evolution with efficiency corresponding to α≈10−2.
these models is simply the viscous time t at the outer
ν Thevalueofthedimensionlessparameterαisbelieved
radius of the disk r ,
out to directly reflect the physics of the mechanism respon-
t ≈r2 /ν. (1) sible for the angular momentum transport in the disk.
ν out
In hot and well ionized accretion disks around compact
Here ν is the kinematic viscosity, which is conveniently objectstransportisgenerallythoughtto be mediatedby
parametrized using the α-prescription (Shakura & Sun- the magneto-rotationalinstability (MRI; Velikhov 1959;
yaev 1973) Chandrasekhar1960; Balbus & Hawley 1991). Situation
ν =αc2Ω−1, (2) is less clear in the cold and poorly ionized protoplan-
s etary disks, where the non-ideal magnetohydrodynamic
withα.1beingconstant,c ≡(kT/µ)1/2 (T isthedisk (MHD)effectsareknowntoweakenorevensuppressthe
s
transportdrivenbythe MRI(Turneretal.2014). Other
1Department of Applied Mathematics and Theoretical potentialcandidatessuchasgravitoturbulence(Gammie
Physics, Centre for Mathematical Sciences, University of 2001; Rafikov 2015), Rossby-wave instabilitiy (Lovelace
Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK;
et al. 1999), convective over-stability (Latter 2016), ver-
[email protected]
2Institute for Advanced Study, 1 Einstein Drive, Princeton tical shear instability (Urpin & Brandenburg 1998; Stoll
NJ08540 & Kley 2014),and so on, have been proposed to explain
2
the observedevolution ofthe protoplanetarydisks prop-
erties.
On the other hand, recent studies argue that the non-
MRI related transport mechanisms can hardly be re-
sponsible for the observed disk evolution on the Myr
timescales, because of the weakness of the purely hy-
drodynamic transport mechanisms (Stoll & Kley 2014;
Turner et al. 2014). Partly for that reason, the non-
diffusive angularmomentumtransportmechanismssuch
as MHD winds (Wardle & Koenigl 1993; Suzuki & Inut-
suka 2009; Bai & Stone 2013) or spiral shocks (Rafikov
2002, 2016) have been gaining popularity. The distinc-
tivefeatureofthesemechanismsisthattheydonotneed
toobeytheequations(1)-(3),thusresultinginadifferent
relation between M˙ and M .
d
The advent of ALMA made possible more precise and
focused efforts to understand protoplanetary disk evolu-
tion. Recent measurements of the continuum and CO
line emission for a large sample of protoplanetary disks
in Lupus by Ansdell et al. (2016) and Miotello et al.
(2016b), coupled with the most up-to-date determina-
tionsofthemassaccretionrateM˙ ontotheirparentstars
byAlcala´etal.(2014,2016),allowedManaraetal.(2016)
to identify a correlation between the disk mass M and
d
thecentralmassaccretionrate. Thediskmasseswerede-
rivedusing the dustmassesinferredfromthe continuum
sub-mm emission assuming fixed gas-to-dust ratio. This
correlationhas been interpretedby Manaraet al. (2016)
as providing evidence for the viscous character of the
protoplanetary disk evolution, in which the global disk
properties directly determine the mass accretion rate at
its center.
In this work we focus on a different diagnostics of the Fig.1.— (a)CentralaccretionrateM˙ and(b)disksizerout for
viscousdiskevolution. Usingasampleofprotoplanetary thesampleof objects usedinthis work(Table 1), plotted vs. the
disksdirectlyresolvedbyALMA,withmeasureddustand disk mass Md inferred from the continuum dust emission. While
gas masses (Ansdell et al. 2016), as well as the central we do not find strong evidence for a correlation between M˙ and
accretion rates (Alcala´ et al. 2014, 2016), we provide a Md in our sample [cf. Manara et al. (2016)], we do observe a
correlationbetweenrout andMd,withthelinearregressionshown
direct determination of the value of the α-parameter in asthebluedashedline. Reddotshowsobject2MASS-J16081497-
individual systems. Given that different mechanisms of 3857145, which is close to the brown dwarf regime. See text for
the angular momentum transportin disks predict differ- moredetails.
entvalues ofα, oureffortcanprovidedirectinformation
onthephysicalnatureoftheinternalstressesdrivingthe
disk evolution. Unlike other studies (Hartmann et al. bination yields
1998; Jones et al. 2012), in this work we (1) utilize in-
formation about the individual disk sizes provided by M˙ µ
ALMAand(2)donotuseinformationontheagesofthe α≈ Ωr2, (4)
M kT
d
parent stars, which are known to be very uncertain.
Our work is organized as follows. We describe our whereΩ,r, T areevaluatedattheouterdiskradiusr ,
out
methodology for inferring α in §2,and our observational and M is the disk mass enclosed within r .
d out
sample in §2.1. Our results, including correlations of α In equation (4) the values of r , M , M˙ , and M
out d ⋆
withdifferentcharacteristicsoftheobservedsystems,can come directly from observations. However, to obtain α
be found in §3. We provide extensive discussion of our we still need to make assumptions about the disk tem-
findings in §4, and summarize the results in §5. perature T(r ). We try three different thermodynamic
out
prescriptions in this work.
First, we simply assume that
2. METHODOLOGY T(rout)=20 K, (5)
In our analysis we will assume that, as a result of ex- forall disks inour sample. This prescriptionis the same
pansiondrivenby internalstresses,the presentday sizes asthatusedbyAnsdelletal.(2016)forderivingthedust
of the protoplanetary disks in Lupus exceed their initial massesM of the disks from their continuum sub-mm
dust
radii, set at the mass infall phase. Then the disk can fluxes, providing certain internal consistency.
be approximatelyconsideredas evolvingin a self-similar Second, we take the T(r ) to correspondto the tem-
out
fashion, and equations (1)-(3) should apply. Their com- perature of optically thin dust, directly illuminated by
3
the central star, in which case
L 1/4
⋆
T(r )= . (6)
out (cid:18)16πσr2 (cid:19)
out
Thisexpressionneglectsthedifferencebetweentheemis-
sion and absorption efficiencies of the grains. Stellar lu-
minosity L is known to us from observations.
⋆
Finally, we also use a prescription for the optically
thick,externallyirradiatedpassiveprotoplanetarydisks,
motivated by Chiang & Goldreich (1997), that reads
L⋆ 2/7 M⊙ 1/7 AU 3/7
T(r )=120 K . (7)
out
(cid:18)L⊙(cid:19) (cid:18)M⋆(cid:19) (cid:18)rout(cid:19)
This prescriptionpredicts T(r ) lower than in the case
out
(6).
We determine the full disk masses using the dust
massesM derivedfromthecontinuumsub-mmfluxes
dust
(Ansdell et al. 2016). To convert M to the full
dust
disk mass M we use a uniform gas-to-dust mass ratio
d
χ = 100. In certain cases (§3.2) we also use the infor-
mation on the gas masses coming from the 13CO and
C18O line measurements by ALMA. However, it should
be kept in mind that the disk masses inferred this way
Fig.2.— Distributionoftheinferredvaluesofα,computedvia
are believed to be systematically underestimated, often equation(4). Solidblack,dottedred,anddashedgreenhistograms
by more than an order of magnitude, as a result CO correspondtoT(rout)prescriptionsgivenbytheequations(5),(6),
freeze out on dust grains or sequestrationof carboninto and(7),correspondingly. Onecanseealargespreadinthevalues
ofα,coveringmorethantwoordersofmagnitude.
large bodies (Ansdell et al. 2016; Miotello et al. 2016a).
As a result, it is expected that Md = χMdust, which we and M˙ and stellar parameters from Alcala´ et al. (2014,
employ in this work, should provide a better estimate of
2016).
the disk mass.
Forsimplicity,inthisstudyweassociatetheouterdisk
radiusr with the semi-majoraxisobtainedin Ansdell
2.1. Observational sample out
et al. (2016) by simple Gaussian fitting of the resolved
Our approach to determining α via the equation (4) continuum intensity pattern. This alone may introduce
works only for disks that have non-trivialmeasurements asystematicduncertaintyinthedeterminationofr at
out
of M˙ , M , and r , as well as of L and M . For this theleveloftensofpercent. Evenmoreseriouserrormay
d out ⋆ ⋆
reason, we are interested only in resolved disks with sig- arisefromthepossibledifferencebetweentheradiiofthe
nificant detections of both the continuum dust emission gasanddustdisks,evidencefor whichhasbeen foundin
byALMAandthestellarmassaccretionrateM˙ viaspec- anumber of systems(Andrews etal.2012;Cleeves etal.
troscopy. 2016; Walsh et al. 2016). We discuss the impact of the
Ansdelletal.(2016)havecarriedoutALMAsurveyof rout uncertainty on our results in §4.1.
89 protoplanetary disks in Lupus star-forming complex We show some characteristics of our systems in Fig-
at ∼ 150−200 pc away from the Sun (age 3±2 Myr, ure 1. We single out one object — 2MASS-J16081497-
Alcala´ et al.2014). They directly resolvedmany sources 3857145—whichisclosetothebrowndwarfregimeand
and provided initial measurements of the dust and gas is different from the rest of the sample (shown as a red
masses for about 2/3 and 1/3 of their sample, corre- dot).
spondingly. Miotello et al. (2016b) carried out a more Thetoppanelshowsthat,unlikeManaraetal.(2016),
sophisticated analysis of this data set based on work of we do not observe a significant correlation between Md
Miotelloetal.(2016a),providingmoreaccuratedustand andM˙ (see Table 2 for the statistical parametersofcor-
gas mass measurements. At the same time, Alcala´ et al. relations shown in our plots: Pearson correlation coeffi-
(2014,2016)carriedoutX-shooterspectroscopyformany cient ρ, Spearman’s rank correlation coefficient, r , and
s
ofthese targets,derivingcentralmassaccretionratesM˙ p-value—probabilityofthenullhypothesisthatthetwo
based on the UV excesses. variables have zero correlation). This is most likely ex-
By examining the samples presented in these studies plained by the modest size of our sample compared to
weselected26objects,whichpossessresolveddisks with that of Manara et al. (2016).
well measured M and M˙ . Out of these disks 18 also On the other hand, we do find an appreciable correla-
dust
tionbetweenthedisksizeandmass,asFigure1bdemon-
have significant measurements of the gas mass based on
strates—moreextendeddiskstypicallyhavelargerdust
CO line measurements. Two disks — Sz84 and MYLup
masses. The best fit bisector regression (Isobe et al.
— fall in the transitional disk category (Alcala´ et al.
1990) describing this correlation is lgr = (3.05 ±
2016). The parameters of all 26 system are listed in out
Table 1. We adopt Mdust and disk sizes from Ansdell 0.16)+(0.6±0.06)lgM˙ (withroutandM˙ measuredinAU
etal.(2016),gasmassesMg fromMiotelloetal.(2016b), andM⊙ yr−1), but there is significantscatter aroundit.
4
Fig.3.— Values of α plotted against some global characteristics of the disk: (a) disk mass Md, (b) outer radius rout, (c) Md/ro2ut,
which is a proxy for the surface density at rout, and (d) (M⋆rout)1/2/Md, which is a combination of variables entering equation (4). No
statisticallysignificantcorrelationsbetween αandtheseglobalvariablesarefound(seeTable2forquantitative metrics).
This relationmay seem to suggestthat the values of the
disksurfacedensityattheouteredgeΣ(r )∝M /r2
out d out
should be roughly the same. This could raise a worry
thatr ,interpretedbyAnsdelletal.(2016)astheouter
out
extent of the disk, in fact corresponds to the detection
limit of ALMA. However, Figure 3c shows that this is
not the case, and that M /r2 for the observed sam-
d out
ple spans almost two ordersof magnitude, thanks to the
large scatter in Figure 1b.
3. RESULTS
InFigure2we showthe histogramsforthe values ofα
computed through equation (4) for different thermody-
namic assumptions, as shown on the panel. One can see
thatdifferentmethods ofcalculatingthe outerdisk tem-
perature do not result in large differences in the values
of α. Regardless of our assumptions, the distribution
of α does not seem to show complicated substructure,
roughly consistent with either being peaked (for T(r )
out
givenbyequations(6)and(7))orapproximatelyuniform
(for T(r )=20 K).
out
The most important feature of these distributions is
their broad range. Irrespective of the T(r ) prescrip-
out
tion, we find that in our sample of 26 disks the values of
αspanmorethantwoordersofmagnitude—from10−4
to0.04. Thisspreadishardlycompatiblewiththesimple
ideaofasingleangularmomentumtransportmechanism
setting the value of α, as one would then expect a nar-
rowly peaked distribution of α values. Nor could it be
several physical mechanisms operating in different sys-
tems (e.g. different instabilities driving the transport),
asthenonewouldexpecttoseemoresubstructureinthe
distribution of α.
Asweshowin§4.1,this spreadisunlikely tobe forced Fig.4.— Values of α computed assuming T(rout) given by (a)
equation(6)and(b)equation(7),plottedasafunctionofthedisk
bytheintrinsicscatterorobservationalerrorsinoursam- mass Md. Nocorrelations emerge here, in agreement with Figure
ple. Thus, we are left to hypothesize that there may be 3a.
some real physical reasons for this behavior of α in dif-
ferent systems, and we try to identify them next. the global evolution of the disk. For this reason it is
natural to seek possible connection of α with the global
3.1. Dependence of α on global disk properties variables characterizing disk on scales ∼r .
out
What we are calculating via equation (4), given the In Figure 3 we plot α computed for T(rout) = 20 K
observables, is the value of α at r , which determines versusthe disk mass M , its radialextentr , M /r2 ,
out d out d out
5
Fig.5.—Effectiveviscosityα,computedforT =20K(equation
5), as a function of the central mass accretion rate M˙. One can
clearly see a strong correlation between the two variables. Blue
line is the best fit to the data given by equation (8). Horizontal
cyan line is the dependence, around which the data points would
be expected to cluster if the angular momentum transport were
characterizedbyasinglevalueofα(takentobeequalto10−2 for
illustrative purposes). Such clustering is clearly not exhibited by
oursample,necessitatingmodificationstothesimplepictureofthe
viscousevolutionoftheprotoplanetary disks.
Fig.6.— Correlations between the central mass accretion rate
M˙ and α, computed using different assumptions about disk tem-
which characterizes the surface density at r , and
out perature: (a) T(r) given by equation (6), and (b) T(r) given by
(M⋆rout)1/2Md−1,whichappearsinequation(4)together equation(7). CorrelationbetweenM˙ andαclearlypersistsinboth
with M˙ . The errors on α were calculated quadratically cases.
fromtheuncertaintiesoftheobservables,asfollowsfrom line
equation (4). It is clear that these plots do not reveal
significant correlations of α with these global variables. lgα=(5.98±0.91)+(0.97±0.1)lgM˙ , (8)
This is also confirmedby the quantitative metrics of the
possiblerelationsbetweenthepairsofvariablesshownin (withM˙ measuredinM⊙ yr−1)whichisconsistentwith
a linear dependence. This relation links the broad dis-
Table 2.
tribution of α seen in Figure 2 with the spread of M˙ in
One may wonder whether this lack of correlationwith
our sample.
the global disk parameters is forced by our simple as-
sumptionaboutthethermalstateofthedisk,represented The α-M˙ correlation is robust with respect to our as-
by the equation (5). To address this issue, in Figure 4 sumptionsaboutT(rout),asfurtherdemonstratedinFig-
we show the analog of Figure 3a, i.e. α vs. M , but cal- ure 6. There we again observe that M˙ and the effective
d
culated for T(r ) given by equations (6) and (7). One viscosity parameter computed using equations (6) and
out
cansee that, again,no correlationis presentinthe data, (7)arestronglycorrelated,despite the differentassump-
implying that this result is robust with regard to our tions about disk thermodynamics.
assumptions about the disk temperature structure. To Correlation persists even if use the CO-based gas
summarize,we do not find any clear dependence of α on masses M , available for 18 objects in our sample, in-
g
the most obvious global characteristics of the disk. stead of the dust-based masses M when computing α.
d
This is illustrated in Figure 7. The spread in the val-
3.2. Dependence of α on M˙ ues of α measured this way is considerably larger than
in Figures 5-6, and the best fit line is significantly offset
Effective viscosity computed via equation (4) depends
from the relation (8), illustrating the problem with the
not only on the global disk characteristics, but also on
CO-baseddiskmasses(Ansdelletal.2016;Miotelloetal.
the central mass accretion rate M˙ . In Figure 5 we plot
2016a).
the effective viscosity α for T(rout)=20 K vs. M˙ . This If we were to take the CO-based masses at face value,
Figureclearlyrevealsastrongcorrelation(Pearsoncoef- wewouldconclude fromFigure 7 that, inthe framework
ficientρ(α,M˙ )=0.877)betweenαandM˙ . Simplelinear oftheviscousmodelforthediskevolutionbasedonequa-
bisectorregression(Isobe etal.1990)resultsinabestfit tions(1)-(3),somesystemsrequireα&1. Assuchvalues
6
is a certain causal relation between M˙ and the inferred
disk viscosity.
3.3. Dependence of α on stellar properties
Havingfoundcorrelationofα withM˙ , whichis a local
characteristic measured at the star, we also checked if α
could have some relation to other stellar parameters.
In Figure 9 we examine this possibility, finding no
significant correlations between α and the stellar mass,
luminosity, or radius. Weak correlations that may be
present in the full data set vanish when we remove
the brown dwarf-like object 2MASS-J16081497-3857145
(which has very distinct properties and strongly affects
covariances between the variables) from the sample.
Thislackofcorrelationisnotsurprisingfromthephysi-
calpointofview,asonemayexpectonlyrelativelyweak
effect of M (e.g. through local shear, proportional to
⋆
M1/2) or L (on which the disk temperature might de-
⋆ ⋆
pend) on the global disk properties.
4. DISCUSSION
Having established a close relation between the mass
accretion rate onto the central star M˙ and the inferred
Fig.7.— αcomputedusinggasmassesMginferredfromtheCO value of α on the global scale of the disk, we now seek
line observations, instead of the dust-based masses Md computed to understand the implications of this finding. When
using the continuum sub-mm emission(displayed inFigures 5-6), doing this, it is also important to keep in mind the lack
shown as a function of M˙. Despite the use of a different tracer of any significant correlations of α with other obvious
of the disk mass, the α-M˙ correlation (dahed line) is still present characteristics of the system, be it global (like M or
at high significance. Solidline corresponds to the correlation (8), d
r ) or local, stellar (e.g. M , or L ).
whichisclearlyoffsetfromthebestfitforCO-basedαandM˙. out ⋆ ⋆
There are different ways, in which such a correlation
could emerge. First, it may result from various system-
ofαareunlikely,thiscould,again,demonstratetheprob-
aticeffectsrelatedtothemeasurementoftheobservables
lem with the disk mass determinations based on the CO
(§4.1). Second, there may be a physical reason for the
line emission.
The existence of a tight correlation between α and M˙ correlation. Thiswouldbe thecaseif,forexample,some
processesrelated to accretionof gas onto the stellar sur-
is a non-trivial and rather unexpected result. Indeed,
faceareabletoinfluencethevalueofαglobally,onscales
if the angular momentum transport in the disk were ef-
∼ r (§4.2). Alternatively, the value of α-parameter
fected by a mechanism characterized by a unique value out
may depend on some yet unidentified property of the
of α, then Figure 5 would look very differently, with α
clustering around a well-defined value regardless of M˙ , protoplanetary disk, resulting in observed spread, and
and the slope of α(M˙ ) relation being close to zero, as giving rise to a variation of M˙ (§4.3). Third, the α-
illustratedbythe cyanlineinthis Figure. Theα=10−2 M˙ relation (8) may simply reflect the way, in which M˙
corresponding to this line is for illustrative purposes, al- enters the determination of α in equation (4), with M˙
thoughthis valuehasbeensuggestedbythe paststudies being, in fact, largely unrelated to the global disk char-
(Hartmann et al. 1998). In that case the variation of M˙ acteristics. This wouldbe the case ifthe centralM˙ were
wouldhavebeenexactlycompensatedbythevariationof decoupled from the global accretion rate set by the disk
Ωr2 /(M T). Figure 3d shows that the latter variable properties, e.g. as a result of some instability operating
out d
(proportionalto(M r )1/2M−1)exhibitsessentiallyno in the inner disk, or mass accumulation in a dead zone
⋆ out d (§4.4). Decoupling would also be natural if the angular
correlationwith α, unlike M˙ entering the expression (4)
momentum transport in protoplanetary disks does not
for α in an identical fashion. Thus, the very fact that
havea diffusive character(§4.5) and is not characterised
a strong α-M˙ correlationexists tells us something inter- by equations (1)-(3).
esting. We now examine each of these possibilities in detail.
Wealsoexploredthebehaviorofthedimensionalkine-
matic viscosity ν ≈ r2M˙ /M , which plays the role of a 4.1. Observational biases
d
diffusion coefficient for viscous spreading of the disk, see Ourcalculationofαinvolvesseveralobservables—M˙ ,
Figure 8. One possible advantage of using ν instead of
M , r — and we need to make sure that the origin of
α is that its determination does notinvolveassumptions d out
the α-M˙ correlation is not related to the possible sys-
about the thermodynamic properties of the disk. One
tematic biases in their measurement. We do this next
can see that ν is also strongly correlated with M˙ (al-
for each of these variables.
though the spread around the best fit line is larger than
in Figures 5-6), while at the same time being indepen- 4.1.1. Uncertainty in M˙
dentofeitherM orr . Thisagainsuggeststhatthere
d out
7
Fig.8.—Dimensionalkinematicviscosityν characterizingglobaldiskevolution, computedforT =20K,asafunctionofseveralglobal
parameters: (a)centralmassaccretionrateM˙,(b)diskmassMd,and(c)theouter diskradiusrout. Thereisaclearcorrelationbetween
ν andM˙,butnocorrelationbetweenν andeitherMd orrout.
Fig.9.— Effective viscosity α, computed for T =20 K, plotted versus stellar parameters: (a) stellar mass M⋆, (b) luminosity L⋆, and
(c)radiusR⋆. Noclearcorrelationsarepresentinthedata,especiallywhenthenear-browndwarfobject2MASS-J16081497-3857145 (red
dot)isnotincludedinthesample.
of stellar M˙ is a challenging task, which was accom- One way to do this might involve the unobserved por-
plishedinAlcala´etal.(2014,2016)bymeasuringtheUV tionoftheaccretionluminosity,whichcouldskewtheM˙
excess above the stellar photospheric emission. A vari- determinationinasystematicway. Indeed,itmaybethe
ety of factors,including the differences between the stel- casethatinmanysystemsmostoftheaccretionalenergy
lar evolution tracks computed by different groups, con- ≈GM M˙ /R is re-emittedinthe (highenergy)spectral
⋆ ⋆
tribute to the uncertainty in the subsequent derivation
regioninaccessibleto ground-basedinstruments. Inthat
of M˙ , which we conservatively adopted to be about 0.4 case the accretion luminosity measured from the ground
dex (Alcala´ et al. 2016). However, it is not easy to see would account for only a small fraction of the bolomet-
how they could enforce a systematic (and not random) ric accretion flux. If gas accretion is mediated by stellar
correlation(8). magnetosphere, which truncates the disk, then one may
8
expect (Calvet & Gullbring 1998) the discrepancy in M˙ on 12CO emission (Schwarz et al. 2016; Cleeves et al.
determination to correlate with the virial temperature 2016). This statement seems to hold in our sample too
T = (µ/k)GM /R of the gas striking the stellar sur- (Ansdell et al. 2016),based on the disk images obtained
vir ⋆ ⋆
face in free fall — the higher T would shift emission using different tracers.
vir
fromaccretionshocktoshorterwavelengthsandresultin
a more severe underestimate of M˙ . This deviation of M˙ 4.1.3. Uncertainty in Md
from its true value would then lead to an underestimate Finally, we discus the effect of the uncertainty in the
of α inferred through equation (4), perfectly correlated disk mass measurement. Miotello et al. (2016a) derived
with the biased estimate of M˙ . moreaccuratedustcontinuum-basedmassesofthe disks
To address this issue, in Figure 10 we plot M˙ vs. the from the sample of Ansdell et al. (2016) using detailed
virial temperature calculated using stellar parameters radiative transfer calculations of the thermal structure
from Table 1. One can see no correlation of a kind sug- of the disk (instead of assuming a single T = 20 K as
gested above,with systems having higher T not show- in Ansdell et al. 2016). They found that Ansdell et al.
vir
ing systematically lower values of M˙ . The two variables (2016)systematicallyoverestimateMd byaboutafactor
appearcompletelyuncorrelatedinoursample. Thissug- of2forMd .10−2M⊙. However,thisbiaswouldsimply
gests that the determination of M˙ does not suffer from uniformly shift our α-M˙ relation, without affecting its
the bias related to the unobserved accretion luminosity. scatter or slope. A similar effect would be causedby the
possibility of aninwarddrift ofsolids (see §4.1.2),which
4.1.2. Uncertainty in rout tends to decrease the gas-to-dust ratio χ in the disk re-
In this work we also implicitly assumed that r , ob- gion probed by the sub-mm continuum measurements.
out
tained in Ansdell et al. (2016) by fitting a Gaussian to However, such bias would just shift down (roughly uni-
the observed intensity pattern, is the true outer disk ra- formly) the disk mass enclosed within the dust disk ra-
dius, which encloses its full mass. One may worry that, dius,withoutbreakingtheα−M˙ correlationofincreasing
in fact, this radius corresponds to the sensitivity limit the spread of α.
of ALMA and in reality the disk extends beyond r , Moreover, equation (4) remains valid even if r and
out out
so that both r and M underestimate their true val- M do not characterize the full disk: as long as M (<
out d d d
ues. However, Figure 3c demonstrates that this is not r) accounts for the full disk mass enclosed within some
the case: the values of M /r2 proportional to the sur- radiusr, theirvalues canbe usedinsteadofM andr
d out d out
face brightness of the outer disk do not cluster around a for the determination of α via equation (4).
singlevalue(whichcouldbeinterpretedasthesensitivity Based on this discussion we conclude that observa-
limit of observations)but rather extend over almost two tional uncertainties and biases are unlikely causes of the
orders of magnitude. correlation(8)andcanhardlyaccountforthefullspread
A potentially more serious issue with r may arise in the inferred values of α seen in Figure 2.
out
in systems with different apparent sizes of the gas and
dustdisks. Evidenceforthisdiscrepancyhasbeenfound 4.2. M˙ setting α
recently in TW Hya (Andrews et al. 2012), IM Lup Anotherpossibilityforthe originofα-M˙ correlationis
(Cleevesetal.2016),HD97048(Walshetal.2016),with
that the central M˙ has a direct physical effect on α. It
the dust continuum emission being radially more cen-
is difficult to see how this connection can be realized in
trally concentrated by a factor of 2-3 than the gaseous
diskemittingin12COlines. Thishasbeeninterpretedas practice,sinceαissetbythediskphysicsonglobalscales,
the evidence for the radialinwarddrift of solids in these whileM˙ isalocalproperty,characterizingtheinnermost
disks(Birnstiel&Andrews2014),whichdecouplesradial region of the disk.
distributionsofthegasanddustsurfacedensities. Ifthis One possiblity forestablishing this connectionis if the
interpretationiscorrect,thenthedustmasseswouldstill accretion energy release at the stellar surface has a di-
properly reflect the full disk mass, but the size of the rectimpact onthe value of α onglobalscales. This may
main mass reservoir (gas disk) would be underestimated be the case if the value of α depends on the degree of
by a factorofseveral. Althoughthis issue shouldbe fur- ionization (as may be expected for the non-ideal MRI),
ther explored observationally for our Lupus sample, we and the accretional luminosity plays a major role in de-
believe that it is unlikely to affect our main conclusions termining the ionization balance in the outer disk. If
for the following reasons. that were the case, one would expect to see a correla-
First, the inferred α depends on r rather weakly, tion between the global α and the accretion energy flux
out
e.g. as α ∝ r1/2 if T = 20 K, see equation (4). Thus, a F =GM⋆M˙ /(4πR⋆ro2ut) at r =rout.
out Figure 11 demonstrates that such correlation does in-
possible underestimate of r by a factor of 2−3 would
out
deed exist. However, it shows larger scatter around the
not explain the broad distribution of the inferred values
bestfit line thanthe correlationin Figure5. This would
ofα. Second, it is notclear that the gasdisk sizes based
on12COmeasurementsrepresenttheradiiwheremostof not be expected if it were F rather than M˙ alone being
gas mass is concentrated (which is what the actual r therealculpritbehindtheα-M˙ correlation. Moreover,it
out
should correspond to). Because of the optical thickness is also unlikely that the spectral range used for inferring
ofthe12COlines,itisgenerallybelievedthattheCOiso- M˙ (longward of 310 nm, Alcala´ et al. 2016) dominates
topologuesarebettertracersofthegasmassdistribution the ionization balance of the disk. Nor is it clear that
than the 12CO molecule. And the sizes of regions emit- the accretionenergy releaseprovidesmajor contribution
ting in 13CO and C18O lines tend to be less discrepant to the flux of ionizing photons impinging on the disk
withthedustcontinuum-basedradiithantheonesbased (Glassgold et al. 2000). Furthermore, it is not obvious
9
Fig.10.— Stellar M˙ plotted vs. the virial temperature at the Fig. 11.— Effective viscosity plotted as a function of the
stellar surface Tvir. No obvious correlation is seen in the data, accretion energy flux at the outer radius of the disk, F =
demonstratingthelackofbiasesrelatedtotheunobservedfraction GM⋆M˙/(4πR⋆ro2ut).
oftheaccretionluminosity.
Coolingtimeofthedisk,whichdirectlydependsonthese
whythis physicalmechanismshouldgiverisetoanα(F)
globaldiskcharacteristics,alsocannotbethecontrolling
dependence with a slope so close to unity.
parameter. This likely excludes the vertical shear in-
Onefinalargumentagainstaphysicaleffectofthecen-
stability (Urpin & Brandenburg1998),which sensitively
tral M˙ on the global value of α is that the disk with depends on the local cooling time (Stoll & Kley 2014;
α(F) ∝ F0.8 ∝ r−1.6 dependence, suggested by the best Lin & Youdin 2015), from being a candidate for driving
fit in Figure 11, should have a rather unusual structure. the viscous evolution of the protoplanetary disks.
Indeed,insider thediskshouldconvergetoaconstant
out Atthesametime,thereisavarietyofpossiblecontrol-
M˙ structure, meaning that Σ ∝ M˙ /ν ∝ r0.1T−1, where lingparametersif transportinthe disk is effectedby the
weusedequation(2)forν andtookα∝r−1.6. Sincethe non-ideal MRI. They include (but are not limited to)
disk temperature T does not increase with r, this would ionization fraction on scales ∼ r (Jin 1996; Fleming
out
mean that Σ should be an increasing function of r. This et al. 2000; Bai & Stone 2011),strength of the magnetic
conclusion is hardly compatible with our understanding field inthe disk (Bai& Stone 2011),or its geometry (Si-
of the protoplanetary disk structure. mon et al. 2013b,a). All of these physical characteristics
Forallthesereasonswe donotfindthe directphysical aredifficulttodetermineobservationallyatthemoment.
effect of stellar M˙ on α to be a plausible explanation of Nevertheless, we believe that this way of producing α-
the α-M˙ correlation. M˙ correlationis more plausible than the one outlined in
§4.2.
4.3. α setting M˙
PhysicalconnectionbetweenαandM˙ mayalsoemerge 4.4. Decoupling of stellar M˙ from the global mass
accretion rate
inthe directionopposite to that consideredin §4.2,with
αdirectlyaffectingM˙ . Suchconnectionisrathernatural A correlation between the inferred α and M˙ would
inlightoftheequation(4). However,intheconventional also naturally emerge if the M˙ measuredthrough stellar
pictureαhasaunique value,whichisincompatible with accretional luminosity is, in fact, unrelated to the global
thedistributionshowninFigure2. Thus,intheobserved mass accretion rate Md/tν.
sample the broad range of α has to be caused by some In the standard viscous disk theory the two should be
additional environmental parameter, which controls an- equal, as demonstrated by the equation (3). However,if
gular momentum transport and allows α to vary over stellarM˙ issomehowdecoupledfromtheglobalaccretion
almost three orders of magnitude. And then α-M˙ cor- rate,thenequation(4)wouldnaturallyresultinastrong
relation would naturally emerge from equation (4), with linear correlation between the inferred α (unrelated to
the distribution of α directly translating into the broad the real global α) and M˙ . This would be true even if
range of the M˙ values. the real α set by the physics of the angular momentum
The hidden parameter controlling α cannot be one of transportonscales∼routtakesonauniquevalue. Errors
the globaldisk variables— M , r , globalsurfaceden- in measuring M˙ could lead to this situation, but they
d out
sity — as Figure 3 shows no correlationof α with them. would need to be very dramatic (potentially exceeding
10
twoordersofmagnitude),whichisunlikely,asweshowed observedrecentlyinseveralprotoplanetarydisksinscat-
in §4.1. tered light, e.g. in SAO 206462 (Garufi et al. 2013),
A decoupling between M˙ and M /t (by more than MWC 758 (Benisty et al. 2015), HD 100453 (Wagner
d ν
twoordersofmagnitude,toexplaintherangeofinferred et al. 2015), etc.
α) requires a modification of the simple picture of vis- Webelievethatinlightoftheperceiveddifficultyofthe
cous disk accretion. It can arise, for example, if some knownlocalturbulenttransportmechanismstodrivethe
instability operates in the inner regions of protoplane- protoplanetarydiskevolutiononMyrtimescales(Turner
tary disks, dramatically modulating local M˙ compared etal.2014),thenon-diffusivemechanismsfordrivingdisk
toitsglobalvaluesetonscales∼r . Thecharacteristic evolution should be considered very seriously. Our work
out
maythusprovidestrongindirectevidenceinfavorofthis
timescale for such variability shouldbe substantialfor it
possibility.
to have escaped detection until now. One may suspect
FUorandEXoroutbursts(Audardetal.2014)tobe the 4.6. Comparison with previous studies
known realizations of such an instability. However, one
would then expect the distribution of α to be bimodal, There has been a handful of studies trying to under-
standviscousevolutionoftheprotoplanetarydisksbased
with most disks being in quiescence and having low α,
on observational data. Using the mean properties of
and a small population of disks undergoing an outburst
a large sample of the protoplanetary disks Hartmann
and having high inferred α (Audard et al. 2014).
Another way of decoupling stellar M˙ from the global et al. (1998) have concluded that their effective viscos-
ity should be narrowly clustered around α ≈ 10−2. In
accretionrateisiftheviscousmassflowtowardsthestar
this workwe determine αfor individual objects, andfind
accumulatesinasubstantialmassreservoiratsomeinter-
a much broader distribution of α, extending over more
mediate radii, e.g. in a dead zone (Gammie 1996). This
thantwoordersofmagnitude(Figure2). Thisdifference
reservoir should be able to accumulate large amounts of
suggests that care should be taken when making infer-
mass,comparabletothetotaldiskmassatthestartofits
ences based on the averagedproperties of the sample.
evolution. This may be difficult to realize on timescales
Some studies of viscous evolution of the protoplane-
comparabletothedisklifetime (Myrs),necessitatingpe-
tary disks have tried to verify the equation (3) with t
riodicdepositionofmassfromthereservoirontothestar, ν
set equal to the age of the system t (Hartmann et al.
and making this scenario similar to the aforementioned ⋆
1998; Jones et al. 2012; Manara et al. 2016). Identifi-
instability in the inner disk.
cation of t with the age of the central star is a valid
Alternatively,gasreachingtheinnerregionsofthedisk ν
procedure as long as t exceeds the viscous time of the
maybelostinawind(photoevaporativeorMHD-driven) ⋆
disk at its initial radius r . This assumption is similar
orconsumedbyvigorouslyaccretingplanets. Ifthatwere 0
to our implicit assumption that the current disk size ex-
thecase,theninmanysystemsthemasslossrateshould
ceeds its initial size, r &r , so that viscous evolution
be matchingthe globalaccretionratesetonlargescales, out 0
enters the self-similar regime and the memory of initial
withonlyasmallamountofmassreachingthestar. The
conditionsgetserased. However,the problemwithusing
wind is likely to also affect the angular momentum bud-
this metric of disk evolutionis thatthe determinationof
getofthe disk inanon-trivialmanner,apossibility that
agesoftheyoungstarsisnotoriouslydifficult(Soderblom
we consider next.
et al. 2014).
4.5. Non-viscous evolution of the protoplanetary disks Despite this drawback Hartmann et al. (1998) and
One final, very intriguing possibility, is that the angu- Jones et al. (2012) have tried to verify that Md/M˙ ∼t⋆
larmomentumandmasstransportintheprotoplanetary (up to a constant factor of order unity) using observa-
diskshasa non-viscous(non-diffusive)character. Inthis tional data. Both studies found significant deviations
case equations (1)-(3) do not hold, and α-M˙ correlation (up to two orders of magnitude) from this simple rela-
tion. In Figure 12 we show the characteristic accretion
emergessimplyasaconsequenceofcalculatingαthrough
equation (2), with no real physical meaning for α. Also, time for the objects in our sample, plotted against M˙
stellarM˙ mayhavelittle to do with withthe globaldisk and Md. Dotted lines show the range of ages for the
Lupus objects in our sample (Alcala´ et al. 2016) (we do
parameters, although the work of Manara et al. (2016)
not attempt to use uncertain ages of the individual ob-
does showevidence for a correlationbetweenM andM˙
d jects). It is clear that accretion times of many objects
(which is not obvious in our sample of resolved disks).
fall outside this range, by more than an order of magni-
Such non-viscous transport may be effected in pro-
tude in some cases. This agrees with the conclusions of
toplanetary disks by the magnetically-controlled winds
Hartmann et al. (1998) and Jones et al. (2012).
(Blandford & Payne 1982; Konigl 1989). Outflows from
InFigure12aonecanalsoseeastronganti-correlation
the disks of YSOs are a well studied observational phe-
between M /M˙ and M˙ , with a much weaker statistical
nomenon (Frank et al. 2014). Recently self-consistent d
connection for the M (see Figure 12b). This demon-
launching of the magnetocentrifugal winds has been d
observed in simulations of magnetized accretion disks strates the key role of M˙ for the accretion time, just as
(Suzuki&Inutsuka2009;Bai&Stone2013),addingsup- we found for α in §3.2.
port to this possibility. Jones et al. (2012) were not able to account for the
Another potential driverof the non-diffusive evolution discrepancy between M /M˙ and t even using sophis-
d ⋆
of the protoplanetary disks could be the density waves ticated disk models including the effects of dead zones,
excited by massive perturbers, e.g. planets or stellar photoevaporation, planet formation, etc. Instead, they
companions (Goodman & Rafikov 2001; Rafikov 2002, concluded that it follows from the systematic errors in
2016; Dong et al. 2016). Global spiral waves have been the determinationofM . However,we findthat the sys-
d