Protoplanetary Disks as (Possibly) Viscous Disks
Protoplanetary disks are believed to evolve on Myr timescales in a diffusive (viscous) manner as a result of angular momentum transport driven by internal stresses. Here we use a sample of 26 protoplanetary disks resolved by ALMA with measured (dust-based) masses and stellar accretion rates to derive the dimensionless -viscosity values for individual objects, with the goal of constraining 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 to . Moreover, they correlate with neither the global disk parameters (mass, size, surface density) nor the stellar characteristics (mass, luminosity, radius). However, we do find a strong linear correlation between and the central mass accretion rate . This correlation is unlikely to result from the direct physical effect of on disk viscosity on global scales. Instead, we suggest that it is caused by the decoupling of stellar 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 strength, or geometry), ultimately driving the variation of . (2) The central is decoupled from the global viscous mass accretion rate as a result of an instability or mass accumulation (or loss) in the inner disk. (3) Perhaps the most intriguing possibility is that angular momentum in protoplanetary disks is transported non-viscously, e.g. via magnetohydrodynamic winds or spiral density waves.
Subject headings:accretion, accretion disks — protoplanetary disks — planets and satellites: formation
Protoplanetary disks are thought to persist around their parent stars for relatively short span of time. Observations present a clear evolutionary picture, in which both the fraction of systems exhibiting IR excess (Hillenbrand, 2005) and the mass accretion rate onto the central star (Calvet et al., 2000) decline of Myr timescales. Both observational indicators are thought to be the clear signatures of the presence of the circumstellar disks.
Astrophysical accretion disks are believed to evolve predominantly under the action of their internal stresses (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974), and protoplanetary disks are no exception to the rule. For a long time evolutionary models of the protoplanetary disks have been developed assuming that the disks spread viscously, losing mass to the central star, while at the same time providing the birth site for planet formation. The characteristic time for the disk evolution in these models is simply the viscous time at the outer radius of the disk ,
Here is the kinematic viscosity, which is conveniently parametrized using the -prescription (Shakura & Sunyaev, 1973)
with being constant, ( is the disk temperature) and ( is the mass of the central star) being the local sound speed and Keplerian angular frequency. Viscous models invariably predict that on long timescales (exceeding the viscous time of the initial, more compact disk, so that grows beyond the initial disk radius) the central mass accretion rate should be related to the total disk mass as
(up to a constant factor of order unity) with evaluated at , see equation (1).
The idea of the viscous evolution of the protoplanetary disks, diffusive in character and characterized by equations (1)-(3), has gained certain observational support. In particular, Hartmann et al. (1998) and Calvet et al. (2000) found that the observed average properties of protoplanetary disks can be explained by their viscous evolution with efficiency corresponding to .
The value of the dimensionless parameter is believed to directly reflect the physics of the mechanism responsible for the angular momentum transport in the disk. In hot and well ionized accretion disks around compact objects transport is generally thought to be mediated by the magneto-rotational instability (MRI; Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991). Situation is less clear in the cold and poorly ionized protoplanetary disks, where the non-ideal magnetohydrodynamic (MHD) effects are known to weaken or even suppress the transport driven by the MRI (Turner et al., 2014). Other potential candidates such as gravitoturbulence (Gammie, 2001; Rafikov, 2015), Rossby-wave instabilitiy (Lovelace et al., 1999), convective over-stability (Latter, 2016), vertical shear instability (Urpin & Brandenburg, 1998; Stoll & Kley, 2014), and so on, have been proposed to explain the observed evolution of the protoplanetary disks properties.
On the other hand, recent studies argue that the non-MRI related transport mechanisms can hardly be responsible for the observed disk evolution on the Myr timescales, because of the weakness of the purely hydrodynamic transport mechanisms (Stoll & Kley, 2014; Turner et al., 2014). Partly for that reason, the non-diffusive angular momentum transport mechanisms such as MHD winds (Wardle & Koenigl, 1993; Suzuki & Inutsuka, 2009; Bai & Stone, 2013) or spiral shocks (Rafikov, 2002, 2016) have been gaining popularity. The distinctive feature of these mechanisms is that they do not need to obey the equations (1)-(3), thus resulting in a different relation between and .
The advent of ALMA made possible more precise and focused efforts to understand protoplanetary disk evolution. 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 determinations of the mass accretion rate onto their parent stars by Alcalá et al. (2014, 2016), allowed Manara et al. (2016) to identify a correlation between the disk mass and the central mass accretion rate. The disk masses were derived using the dust masses inferred from the continuum sub-mm emission assuming fixed gas-to-dust ratio. This correlation has been interpreted by Manara et 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 viscous disk evolution. Using a sample of protoplanetary disks directly resolved by ALMA, with measured dust and gas masses (Ansdell et al., 2016), as well as the central accretion rates (Alcalá et al., 2014, 2016), we provide a direct determination of the value of the -parameter in individual systems. Given that different mechanisms of the angular momentum transport in disks predict different values of , our effort can provide direct information on the physical nature of the internal stresses driving the disk evolution. Unlike other studies (Hartmann et al., 1998; Jones et al., 2012), in this work we (1) utilize information about the individual disk sizes provided by ALMA and (2) do not use information on the ages of the parent stars, which are known to be very uncertain.
Our work is organized as follows. We describe our methodology for inferring in §2, and our observational sample in §2.1. Our results, including correlations of with different characteristics of the observed systems, can be found in §3. We provide extensive discussion of our findings in §4, and summarize the results in §5.
In our analysis we will assume that, as a result of expansion driven by internal stresses, the present day sizes of the protoplanetary disks in Lupus exceed their initial radii, set at the mass infall phase. Then the disk can be approximately considered as evolving in a self-similar fashion, and equations (1)-(3) should apply. Their combination yields
where , , are evaluated at the outer disk radius , and is the disk mass enclosed within .
In equation (4) the values of , , , and come directly from observations. However, to obtain we still need to make assumptions about the disk temperature . We try three different thermodynamic prescriptions in this work.
First, we simply assume that
for all disks in our sample. This prescription is the same as that used by Ansdell et al. (2016) for deriving the dust masses of the disks from their continuum sub-mm fluxes, providing certain internal consistency.
Second, we take the to correspond to the temperature of optically thin dust, directly illuminated by the central star, in which case
This expression neglects the difference between the emission and absorption efficiencies of the grains. Stellar luminosity is known to us from observations.
Finally, we also use a prescription for the optically thick, externally irradiated passive protoplanetary disks, motivated by Chiang & Goldreich (1997), that reads
This prescription predicts lower than in the case (6).
We determine the full disk masses using the dust masses derived from the continuum sub-mm fluxes (Ansdell et al., 2016). To convert to the full disk mass we use a uniform gas-to-dust mass ratio . In certain cases (§3.2) we also use the information on the gas masses coming from the CO and CO line measurements by ALMA. However, it should be kept in mind that the disk masses inferred this way are believed to be systematically underestimated, often by more than an order of magnitude, as a result CO freeze out on dust grains or sequestration of carbon into large bodies (Ansdell et al., 2016; Miotello et al., 2016a). As a result, it is expected that , which we employ in this work, should provide a better estimate of the disk mass.
2.1. Observational sample
Our approach to determining via the equation (4) works only for disks that have non-trivial measurements of , , and , as well as of and . For this reason, we are interested only in resolved disks with significant detections of both the continuum dust emission by ALMA and the stellar mass accretion rate via spectroscopy.
Ansdell et al. (2016) have carried out ALMA survey of 89 protoplanetary disks in Lupus star-forming complex at pc away from the Sun (age Myr, Alcalá et al. 2014). They directly resolved many sources and provided initial measurements of the dust and gas masses for about 2/3 and 1/3 of their sample, correspondingly. Miotello et al. (2016b) carried out a more sophisticated analysis of this data set based on work of Miotello et al. (2016a), providing more accurate dust and gas mass measurements. At the same time, Alcalá et al. (2014, 2016) carried out X-shooter spectroscopy for many of these targets, deriving central mass accretion rates based on the UV excesses.
By examining the samples presented in these studies we selected 26 objects, which possess resolved disks with well measured and . Out of these disks 18 also have significant measurements of the gas mass based on CO line measurements. Two disks — Sz84 and MYLup — fall in the transitional disk category (Alcalá et al., 2016). The parameters of all 26 system are listed in Table 1. We adopt and disk sizes from Ansdell et al. (2016), gas masses from Miotello et al. (2016b), and and stellar parameters from Alcalá et al. (2014, 2016).
For simplicity, in this study we associate the outer disk radius with the semi-major axis obtained in Ansdell et al. (2016) by simple Gaussian fitting of the resolved continuum intensity pattern. This alone may introduce a systematic duncertainty in the determination of at the level of tens of per cent. Even more serious error may arise from the possible difference between the radii of the gas and dust disks, evidence for which has been found in a number of systems (Andrews et al., 2012; Cleeves et al., 2016; Walsh et al., 2016). We discuss the impact of the uncertainty on our results in §4.1.
We show some characteristics of our systems in Figure 1. We single out one object — 2MASS-J16081497-3857145 — which is close to the brown dwarf regime and is different from the rest of the sample (shown as a red dot).
The top panel shows that, unlike Manara et al. (2016), we do not observe a significant correlation between and (see Table 2 for the statistical parameters of correlations shown in our plots: Pearson correlation coefficient , Spearman’s rank correlation coefficient, , and -value — probability of the null hypothesis that the two variables have zero correlation). This is most likely explained by the modest size of our sample compared to that of Manara et al. (2016).
On the other hand, we do find an appreciable correlation between the disk size and mass, as Figure 1b demonstrates — more extended disks typically have larger dust masses. The best fit bisector regression (Isobe et al., 1990) describing this correlation is (with and measured in AU and yr), but there is significant scatter around it. This relation may seem to suggest that the values of the disk surface density at the outer edge should be roughly the same. This could raise a worry that , interpreted by Ansdell et al. (2016) as the outer 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 for the observed sample spans almost two orders of magnitude, thanks to the large scatter in Figure 1b.
In Figure 2 we show the histograms for the values of computed through equation (4) for different thermodynamic assumptions, as shown on the panel. One can see that different methods of calculating the outer disk temperature 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 given by equations (6) and (7)) or approximately uniform (for K).
The most important feature of these distributions is their broad range. Irrespective of the prescription, we find that in our sample of 26 disks the values of span more than two orders of magnitude — from to . This spread is hardly compatible with the simple idea of a single angular momentum transport mechanism setting the value of , as one would then expect a narrowly peaked distribution of values. Nor could it be several physical mechanisms operating in different systems (e.g. different instabilities driving the transport), as then one would expect to see more substructure in the distribution of .
As we show in §4.1, this spread is unlikely to be forced by the intrinsic scatter or observational errors in our sample. Thus, we are left to hypothesize that there may be some real physical reasons for this behavior of in different systems, and we try to identify them next.
3.1. Dependence of on global disk properties
What we are calculating via equation (4), given the observables, is the value of at , which determines the global evolution of the disk. For this reason it is natural to seek possible connection of with the global variables characterizing disk on scales .
In Figure 3 we plot computed for K versus the disk mass , its radial extent , , which characterizes the surface density at , and , which appears in equation (4) together with . The errors on were calculated quadratically from the uncertainties of the observables, as follows from equation (4). It is clear that these plots do not reveal significant correlations of with these global variables. This is also confirmed by the quantitative metrics of the possible relations between the pairs of variables shown in Table 2.
One may wonder whether this lack of correlation with the global disk parameters is forced by our simple assumption about the thermal state of the disk, represented by the equation (5). To address this issue, in Figure 4 we show the analog of Figure 3a, i.e. vs. , but calculated for given by equations (6) and (7). One can see that, again, no correlation is present in the data, implying that this result is robust with regard to our assumptions about the disk temperature structure. To summarize, we do not find any clear dependence of on the most obvious global characteristics of the disk.
3.2. Dependence of on
Effective viscosity computed via equation (4) depends not only on the global disk characteristics, but also on the central mass accretion rate . In Figure 5 we plot the effective viscosity for K vs. . This Figure clearly reveals a strong correlation (Pearson coefficient ) between and . Simple linear bisector regression (Isobe et al., 1990) results in a best fit line
(with measured in yr) which is consistent with a linear dependence. This relation links the broad distribution of seen in Figure 2 with the spread of in our sample.
The - correlation is robust with respect to our assumptions about , as further demonstrated in Figure 6. There we again observe that and the effective viscosity parameter computed using equations (6) and (7) are strongly correlated, despite the different assumptions about disk thermodynamics.
Correlation persists even if use the CO-based gas masses , available for 18 objects in our sample, instead of the dust-based masses when computing . This is illustrated in Figure 7. The spread in the values of measured this way is considerably larger than in Figures 5-6, and the best fit line is significantly offset from the relation (8), illustrating the problem with the CO-based disk masses (Ansdell et al., 2016; Miotello et al., 2016a).
If we were to take the CO-based masses at face value, we would conclude from Figure 7 that, in the framework of the viscous model for the disk evolution based on equations (1)-(3), some systems require . As such values of are unlikely, this could, again, demonstrate the problem with the disk mass determinations based on the CO line emission.
The existence of a tight correlation between and is a non-trivial and rather unexpected result. Indeed, if the angular momentum transport in the disk were effected by a mechanism characterized by a unique value of , then Figure 5 would look very differently, with clustering around a well-defined value regardless of , and the slope of relation being close to zero, as illustrated by the cyan line in this Figure. The corresponding to this line is for illustrative purposes, although this value has been suggested by the past studies (Hartmann et al., 1998). In that case the variation of would have been exactly compensated by the variation of . Figure 3d shows that the latter variable (proportional to ) exhibits essentially no correlation with , unlike entering the expression (4) for in an identical fashion. Thus, the very fact that a strong - correlation exists tells us something interesting.
We also explored the behavior of the dimensional kinematic viscosity , which plays the role of a diffusion coefficient for viscous spreading of the disk, see Figure 8. One possible advantage of using instead of is that its determination does not involve assumptions about the thermodynamic properties of the disk. One can see that is also strongly correlated with (although the spread around the best fit line is larger than in Figures 5-6), while at the same time being independent of either or . This again suggests that there is a certain causal relation between and the inferred disk viscosity.
3.3. Dependence of on stellar properties
Having found correlation of with , which is 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.
This lack of correlation is not surprising from the physical point of view, as one may expect only relatively weak effect of (e.g. through local shear, proportional to ) or (on which the disk temperature might depend) on the global disk properties.
Having established a close relation between the mass accretion rate onto the central star and the inferred value of on the global scale of the disk, we now seek to understand the implications of this finding. When doing this, it is also important to keep in mind the lack of any significant correlations of with other obvious characteristics of the system, be it global (like or ) or local, stellar (e.g. , or ).
There are different ways, in which such a correlation could emerge. First, it may result from various systematic effects related to the measurement of the observables (§4.1). Second, there may be a physical reason for the correlation. This would be the case if, for example, some processes related to accretion of gas onto the stellar surface are able to influence the value of globally, on scales (§4.2). Alternatively, the value of -parameter may depend on some yet unidentified property of the protoplanetary disk, resulting in observed spread, and giving rise to a variation of (§4.3). Third, the - relation (8) may simply reflect the way, in which enters the determination of in equation (4), with being, in fact, largely unrelated to the global disk characteristics. This would be the case if the central were decoupled from the global accretion rate set by the disk properties, e.g. as a result of some instability operating in the inner disk, or mass accumulation in a dead zone (§4.4). Decoupling would also be natural if the angular momentum transport in protoplanetary disks does not have a diffusive character (§4.5) and is not characterised by equations (1)-(3).
We now examine each of these possibilities in detail.
4.1. Observational biases
Our calculation of involves several observables — , , — and we need to make sure that the origin of the - correlation is not related to the possible systematic biases in their measurement. We do this next for each of these variables.
4.1.1 Uncertainty in
of stellar is a challenging task, which was accomplished in Alcalá et al. (2014, 2016) by measuring the UV excess above the stellar photospheric emission. A variety of factors, including the differences between the stellar evolution tracks computed by different groups, contribute to the uncertainty in the subsequent derivation of , which we conservatively adopted to be about 0.4 dex (Alcalá et al., 2016). However, it is not easy to see how they could enforce a systematic (and not random) correlation (8).
One way to do this might involve the unobserved portion of the accretion luminosity, which could skew the determination in a systematic way. Indeed, it may be the case that in many systems most of the accretional energy is re-emitted in the (high energy) spectral region inaccessible to ground-based instruments. In that case the accretion luminosity measured from the ground would account for only a small fraction of the bolometric accretion flux. If gas accretion is mediated by stellar magnetosphere, which truncates the disk, then one may expect (Calvet & Gullbring, 1998) the discrepancy in determination to correlate with the virial temperature of the gas striking the stellar surface in free fall — the higher would shift emission from accretion shock to shorter wavelengths and result in a more severe underestimate of . This deviation of from its true value would then lead to an underestimate of inferred through equation (4), perfectly correlated with the biased estimate of .
To address this issue, in Figure 10 we plot vs. the virial temperature calculated using stellar parameters from Table 1. One can see no correlation of a kind suggested above, with systems having higher not showing systematically lower values of . The two variables appear completely uncorrelated in our sample. This suggests that the determination of does not suffer from the bias related to the unobserved accretion luminosity.
4.1.2 Uncertainty in
In this work we also implicitly assumed that , obtained in Ansdell et al. (2016) by fitting a Gaussian to the observed intensity pattern, is the true outer disk radius, which encloses its full mass. One may worry that, in fact, this radius corresponds to the sensitivity limit of ALMA and in reality the disk extends beyond , so that both and underestimate their true values. However, Figure 3c demonstrates that this is not the case: the values of proportional to the surface brightness of the outer disk do not cluster around a single value (which could be interpreted as the sensitivity limit of observations) but rather extend over almost two orders of magnitude.
A potentially more serious issue with may arise in systems with different apparent sizes of the gas and dust disks. Evidence for this discrepancy has been found recently in TW Hya (Andrews et al., 2012), IM Lup (Cleeves et al., 2016), HD 97048 (Walsh et al., 2016), with the dust continuum emission being radially more centrally concentrated by a factor of 2-3 than the gaseous disk emitting in CO lines. This has been interpreted as the evidence for the radial inward drift of solids in these disks (Birnstiel & Andrews, 2014), which decouples radial distributions of the gas and dust surface densities. If this interpretation is correct, then the dust masses would still properly reflect the full disk mass, but the size of the main mass reservoir (gas disk) would be underestimated by a factor of several. Although this issue should be further explored observationally for our Lupus sample, we believe that it is unlikely to affect our main conclusions for the following reasons.
First, the inferred depends on rather weakly, e.g. as if K, see equation (4). Thus, a possible underestimate of by a factor of would not explain the broad distribution of the inferred values of . Second, it is not clear that the gas disk sizes based on CO measurements represent the radii where most of gas mass is concentrated (which is what the actual should correspond to). Because of the optical thickness of the CO lines, it is generally believed that the CO isotopologues are better tracers of the gas mass distribution than the CO molecule. And the sizes of regions emitting in CO and CO lines tend to be less discrepant with the dust continuum-based radii than the ones based on CO emission (Schwarz et al., 2016; Cleeves et al., 2016). This statement seems to hold in our sample too (Ansdell et al., 2016), based on the disk images obtained using different tracers.
4.1.3 Uncertainty in
Finally, we discus the effect of the uncertainty in the disk mass measurement. Miotello et al. (2016a) derived more accurate dust continuum-based masses of the disks from the sample of Ansdell et al. (2016) using detailed radiative transfer calculations of the thermal structure of the disk (instead of assuming a single K as in Ansdell et al. 2016). They found that Ansdell et al. (2016) systematically overestimate by about a factor of 2 for . However, this bias would simply uniformly shift our - relation, without affecting its scatter or slope. A similar effect would be caused by the possibility of an inward drift of solids (see §4.1.2), which tends to decrease the gas-to-dust ratio in the disk region probed by the sub-mm continuum measurements. However, such bias would just shift down (roughly uniformly) the disk mass enclosed within the dust disk radius, without breaking the correlation of increasing the spread of .
Moreover, equation (4) remains valid even if and do not characterize the full disk: as long as accounts for the full disk mass enclosed within some radius , their values can be used instead of and for the determination of via equation (4).
Another possibility for the origin of - correlation is that the central has a direct physical effect on . It is difficult to see how this connection can be realized in practice, since is set by the disk physics on global scales, while is a local property, characterizing the innermost region of the disk.
One possiblity for establishing this connection is if the accretion energy release at the stellar surface has a direct impact on the value of on global scales. This may be the case if the value of depends on the degree of ionization (as may be expected for the non-ideal MRI), and the accretional luminosity plays a major role in determining the ionization balance in the outer disk. If that were the case, one would expect to see a correlation between the global and the accretion energy flux at .
Figure 11 demonstrates that such correlation does indeed exist. However, it shows larger scatter around the best fit line than the correlation in Figure 5. This would not be expected if it were rather than alone being the real culprit behind the - correlation. Moreover, it is also unlikely that the spectral range used for inferring (longward of 310 nm, Alcalá et al. 2016) dominates the ionization balance of the disk. Nor is it clear that the accretion energy release provides major contribution to the flux of ionizing photons impinging on the disk (Glassgold et al., 2000). Furthermore, it is not obvious why this physical mechanism should give rise to an dependence with a slope so close to unity.
One final argument against a physical effect of the central on the global value of is that the disk with dependence, suggested by the best fit in Figure 11, should have a rather unusual structure. Indeed, inside the disk should converge to a constant structure, meaning that , where we used equation (2) for and took . Since the disk temperature does not increase with , this would mean that should be an increasing function of . This conclusion is hardly compatible with our understanding of the protoplanetary disk structure.
For all these reasons we do not find the direct physical effect of stellar on to be a plausible explanation of the - correlation.
Physical connection between and may also emerge in the direction opposite to that considered in §4.2, with directly affecting . Such connection is rather natural in light of the equation (4). However, in the conventional picture has a unique value, which is incompatible with the distribution shown in Figure 2. Thus, in the observed sample the broad range of has to be caused by some additional environmental parameter, which controls angular momentum transport and allows to vary over almost three orders of magnitude. And then - correlation would naturally emerge from equation (4), with the distribution of directly translating into the broad range of the values.
The hidden parameter controlling cannot be one of the global disk variables — , , global surface density — as Figure 3 shows no correlation of with them. Cooling time of the disk, which directly depends on these global disk characteristics, also cannot be the controlling parameter. This likely excludes the vertical shear instability (Urpin & Brandenburg, 1998), which sensitively depends on the local cooling time (Stoll & Kley, 2014; Lin & Youdin, 2015), from being a candidate for driving the viscous evolution of the protoplanetary disks.
At the same time, there is a variety of possible controlling parameters if transport in the disk is effected by the non-ideal MRI. They include (but are not limited to) ionization fraction on scales (Jin, 1996; Fleming et al., 2000; Bai & Stone, 2011), strength of the magnetic field in the disk (Bai & Stone, 2011), or its geometry (Simon et al., 2013b, a). All of these physical characteristics are difficult to determine observationally at the moment. Nevertheless, we believe that this way of producing - correlation is more plausible than the one outlined in §4.2.
4.4. Decoupling of stellar from the global mass accretion rate
A correlation between the inferred and would also naturally emerge if the measured through stellar accretional luminosity is, in fact, unrelated to the global mass accretion rate .
In the standard viscous disk theory the two should be equal, as demonstrated by the equation (3). However, if stellar is somehow decoupled from the global accretion rate, then equation (4) would naturally result in a strong linear correlation between the inferred (unrelated to the real global ) and . This would be true even if the real set by the physics of the angular momentum transport on scales takes on a unique value. Errors in measuring could lead to this situation, but they would need to be very dramatic (potentially exceeding two orders of magnitude), which is unlikely, as we showed in §4.1.
A decoupling between and (by more than two orders of magnitude, to explain the range of inferred ) requires a modification of the simple picture of viscous disk accretion. It can arise, for example, if some instability operates in the inner regions of protoplanetary disks, dramatically modulating local compared to its global value set on scales . The characteristic timescale for such variability should be substantial for it to have escaped detection until now. One may suspect FUor and EXor outbursts (Audard et al., 2014) to be the known realizations of such an instability. However, one would then expect the distribution of to be bimodal, with most disks being in quiescence and having low , and a small population of disks undergoing an outburst and having high inferred (Audard et al., 2014).
Another way of decoupling stellar from the global accretion rate is if the viscous mass flow towards the star accumulates in a substantial mass reservoir at some intermediate radii, e.g. in a dead zone (Gammie, 1996). This reservoir should be able to accumulate large amounts of mass, comparable to the total disk mass at the start of its evolution. This may be difficult to realize on timescales comparable to the disk lifetime (Myrs), necessitating periodic deposition of mass from the reservoir onto the star, and making this scenario similar to the aforementioned instability in the inner disk.
Alternatively, gas reaching the inner regions of the disk may be lost in a wind (photoevaporative or MHD-driven) or consumed by vigorously accreting planets. If that were the case, then in many systems the mass loss rate should be matching the global accretion rate set on large scales, with only a small amount of mass reaching the star. The wind is likely to also affect the angular momentum budget of the disk in a non-trivial manner, a possibility that we consider next.
4.5. Non-viscous evolution of the protoplanetary disks
One final, very intriguing possibility, is that the angular momentum and mass transport in the protoplanetary disks has a non-viscous (non-diffusive) character. In this case equations (1)-(3) do not hold, and - correlation emerges simply as a consequence of calculating through equation (2), with no real physical meaning for . Also, stellar may have little to do with with the global disk parameters, although the work of Manara et al. (2016) does show evidence for a correlation between and (which is not obvious in our sample of resolved disks).
Such non-viscous transport may be effected in protoplanetary disks by the magnetically-controlled winds (Blandford & Payne, 1982; Konigl, 1989). Outflows from the disks of YSOs are a well studied observational phenomenon (Frank et al., 2014). Recently self-consistent launching of the magnetocentrifugal winds has been observed in simulations of magnetized accretion disks (Suzuki & Inutsuka, 2009; Bai & Stone, 2013), adding support to this possibility.
Another potential driver of the non-diffusive evolution of the protoplanetary disks could be the density waves excited by massive perturbers, e.g. planets or stellar companions (Goodman & Rafikov, 2001; Rafikov, 2002, 2016; Dong et al., 2016). Global spiral waves have been observed recently in several protoplanetary disks in scattered light, e.g. in SAO 206462 (Garufi et al., 2013), MWC 758 (Benisty et al., 2015), HD 100453 (Wagner et al., 2015), etc.
We believe that in light of the perceived difficulty of the known local turbulent transport mechanisms to drive the protoplanetary disk evolution on Myr timescales (Turner et al., 2014), the non-diffusive mechanisms for driving disk evolution should be considered very seriously. Our work may thus provide strong indirect evidence in favor of this possibility.
4.6. Comparison with previous studies
There has been a handful of studies trying to understand viscous evolution of the protoplanetary disks based on observational data. Using the mean properties of a large sample of the protoplanetary disks Hartmann et al. (1998) have concluded that their effective viscosity should be narrowly clustered around . In this work we determine for individual objects, and find a much broader distribution of , extending over more than two orders of magnitude (Figure 2). This difference suggests that care should be taken when making inferences based on the averaged properties of the sample.
Some studies of viscous evolution of the protoplanetary disks have tried to verify the equation (3) with set equal to the age of the system (Hartmann et al., 1998; Jones et al., 2012; Manara et al., 2016). Identification of with the age of the central star is a valid procedure as long as exceeds the viscous time of the disk at its initial radius . This assumption is similar to our implicit assumption that the current disk size exceeds its initial size, , so that viscous evolution enters the self-similar regime and the memory of initial conditions gets erased. However, the problem with using this metric of disk evolution is that the determination of ages of the young stars is notoriously difficult (Soderblom et al., 2014).
Despite this drawback Hartmann et al. (1998) and Jones et al. (2012) have tried to verify that (up to a constant factor of order unity) using observational data. Both studies found significant deviations (up to two orders of magnitude) from this simple relation. In Figure 12 we show the characteristic accretion time for the objects in our sample, plotted against and . Dotted lines show the range of ages for the Lupus objects in our sample (Alcalá et al., 2016) (we do not attempt to use uncertain ages of the individual objects). It is clear that accretion times of many objects fall outside this range, by more than an order of magnitude in some cases. This agrees with the conclusions of Hartmann et al. (1998) and Jones et al. (2012).
In Figure 12a one can also see a strong anti-correlation between and , with a much weaker statistical connection for the (see Figure 12b). This demonstrates the key role of for the accretion time, just as we found for in §3.2.
Jones et al. (2012) were not able to account for the discrepancy between and even using sophisticated disk models including the effects of dead zones, photoevaporation, planet formation, etc. Instead, they concluded that it follows from the systematic errors in the determination of . However, we find that the systematic underestimate (or overestimate) of would result in a uniform overestimate (underestimate) of , but would not explain the emergence of the - correlation. Our results suggest that the discrepancy between and is more likely to be caused by the decoupling of the central from the global mass accretion rate computed based on the standard theory of viscous disk evolution, as described in §4.4-4.5.
We also note that the scenario, in which is controlled by a yet unidentified variable (§4.3), is not expected to produce significant deviations from relation (which is insensitive to in the self-similar regime). This may argue against this scenario for the protoplanetary disks, although more work is certainly needed to resolve this question.
Our approach bypasses the issue of uncertain stellar ages by simply ignoring them altogether. Instead, we use spatial information to gain insight on the physical mechanisms responsible for the angular momentum transport in the protoplanetary disks by measuring . Past efforts (Hartmann et al., 1998; Jones et al., 2012) did not have the ability to do that because they lacked the accurate information on the sizes of disks in individual objects. Thus, our work represents an independent way of testing the theory of viscous evolution of the protoplanetary disks.
In this work we explored viscous evolution of the protoplanetary disks. Using observational sample of 26 disks resolved with ALMA with measured masses (based on sub-mm continuum) and central accretion rates we derived the values of the dimensionless viscosity parameter , with the goal of constraining the mechanism of the angular momentum transport in the disk. Our findings can be summarized as follows.
The distribution of inferred values of extends over more than two orders of magnitude, from to , with no obvious preferred value inside this interval.
We found no correlation of with either the global disk parameters — mass, size, surface density — or stellar parameters — luminosity, mass, radius.
The main finding of this work is the discovery of a strong linear correlation between and central mass accretion rate , which is robust with regard to the thermodynamic assumptions about the disk. This correlation persists even if we use the CO-based gas masses for computing , and holds not only for but also for the dimensional kinematic viscosity on global scales.
These results suggest that a simple picture, in which viscous evolution of the protoplanetary disks is driven by a physical process (e.g. MRI) with a single, well-defined value of , is too simplistic and must be modified. We find that observational errors and biases cannot account for the observed - correlation, and seek other explanations. We find it unlikely that gas accretion onto the stellar surface can have a direct effect on (e.g. through the accretional energy release) on scales of order the disk size (tens to hundreds of AU). We propose three other possibilities for explaining - correlation, which effectively assume that either or are decoupled from the global characteristics (mass, size) of the disk. In that case equation (2) naturally leads to a linear relation between and . These possibilities are as follows.
The value of in every disk is controlled by some yet unobserved variable, variation of which is responsible for the broad range of . This, in turn, is the main cause of the variation of . In the case of accretion driven by the (non-ideal) MRI the role of such control parameter may be played by the disk ionization, as well as the strength or geometry of the magnetic field in the disk.
Stellar may be decoupled from the global mass accretion rate by some instability operating in the inner disk, or mass accumulation in a dead zone, or a wind with high mass loss rate. In this case the inferred values of do not characterize the global disk evolution.
Finally, disk evolution may have a non-diffusive (non-viscous) character, in which case the inferred has no physical meaning. This may be the case if mass accretion in protoplanetary disks is driven by e.g. magnetocentrifugal winds or spiral density waves.
Future work aimed at expanding the sample of resolved protoplanetary disks with well measured masses and accretion rates will help us to identify the physical reason behind the observed - correlation.
- Alcalá et al. (2014) Alcalá, J. M., Natta, A., Manara, C. F., et al. 2014, A&A, 561, A2
- Alcalá et al. (2016) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2016, ArXiv e-prints, arXiv:1612.07054
- Andrews et al. (2012) Andrews, S. M., Wilner, D. J., Hughes, A. M., et al. 2012, ApJ, 744, 162
- Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46
- Audard et al. (2014) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, Protostars and Planets VI, 387
- Bai & Stone (2011) Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144
- Bai & Stone (2013) —. 2013, ApJ, 769, 76
- Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
- Benisty et al. (2015) Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6
- Birnstiel & Andrews (2014) Birnstiel, T., & Andrews, S. M. 2014, ApJ, 780, 153
- Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
- Calvet & Gullbring (1998) Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802
- Calvet et al. (2000) Calvet, N., Hartmann, L., & Strom, S. E. 2000, Protostars and Planets IV, 377
- Chandrasekhar (1960) Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
- Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
- Cleeves et al. (2016) Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110
- Dong et al. (2016) Dong, R., Zhu, Z., Fung, J., et al. 2016, ApJ, 816, L12
- Fleming et al. (2000) Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464
- Frank et al. (2014) Frank, A., Ray, T. P., Cabrit, S., et al. 2014, Protostars and Planets VI, 451
- Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
- Gammie (2001) —. 2001, ApJ, 553, 174
- Garufi et al. (2013) Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105
- Glassgold et al. (2000) Glassgold, A. E., Feigelson, E. D., & Montmerle, T. 2000, Protostars and Planets IV, 429
- Goodman & Rafikov (2001) Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793
- Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385
- Hillenbrand (2005) Hillenbrand, L. A. 2005, ArXiv Astrophysics e-prints, astro-ph/0511083
- Isobe et al. (1990) Isobe, T., Feigelson, E. D., Akritas, M. G., & Babu, G. J. 1990, ApJ, 364, 104
- Jin (1996) Jin, L. 1996, ApJ, 457, 798
- Jones et al. (2012) Jones, M. G., Pringle, J. E., & Alexander, R. D. 2012, MNRAS, 419, 925
- Konigl (1989) Konigl, A. 1989, ApJ, 342, 208
- Latter (2016) Latter, H. N. 2016, MNRAS, 455, 2608
- Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17
- Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
- Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
- Manara et al. (2016) Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3
- Miotello et al. (2016a) Miotello, A., van Dishoeck, E. F., Kama, M., & Bruderer, S. 2016a, A&A, 594, A85
- Miotello et al. (2016b) Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2016b, ArXiv e-prints, arXiv:1612.01538
- Rafikov (2002) Rafikov, R. R. 2002, ApJ, 569, 997
- Rafikov (2015) —. 2015, ApJ, 804, 62
- Rafikov (2016) —. 2016, ApJ, 831, 122
- Schwarz et al. (2016) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91
- Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
- Simon et al. (2013a) Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ, 775, 73
- Simon et al. (2013b) Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ, 764, 66
- Soderblom et al. (2014) Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, Protostars and Planets VI, 219
- Stoll & Kley (2014) Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77
- Suzuki & Inutsuka (2009) Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJ, 691, L49
- Turner et al. (2014) Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
- Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399
- Velikhov (1959) Velikhov, E. P. 1959, JETP, 36, 1398
- Wagner et al. (2015) Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJ, 813, L2
- Walsh et al. (2016) Walsh, C., Juhász, A., Meeus, G., et al. 2016, ApJ, 831, 200
- Wardle & Koenigl (1993) Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218
|Name||,||,||,||SMA, arcsec||,||,||, pc|
|Sz65||0.64 [0.2, 1.5]||150|
|Sz68||0.68 [0.2, 1.5]||150|
|Sz69||0.034 [0.018, 0.07]||150|
|Sz71||0.096 [0.07, 0.3]||150|
|Sz83||1.5 [0.48, 4.0]||150|
|Sz84 (td)||0.11 [0.06, 0.22]||150|
|Sz90||0.056 [0.035, 0.1]||200|
|Sz98||0.066 [0.04, 0.1]||200|
|Sz108B||0.65 [0.2, 1.5]||200|
|Sz114||0.096 [0.065, 0.28]||200|
|Sz129||0.046 [0.03, 0.09]||150|
|Sz130||0.036 [0.011, 0.05]||150|
|MYLup (td)||0.083 [0.05, 0.21]||150|
|SSTc2d-J154508.9-341734||0.77 [0.25, 2.0]||150|
|SSTc2d-J160002.4-422216||0.14 [0.11, 0.7]||150|
|2MASSJ-16085324-3914401||0.034 [0.02, 0.07]||200|
|SSTc2d-J161029.6-392215||0.16 [0.1, 0.56]||200|
|2MASS-J16081497-3857145||0.022 [0.01, 0.045]||200|
Notes For each object (Name) we list the luminosity , radius , mass , and accretion rate of its parent star, the outer semi-major axis (SMA) and dust and gas masses ( and , the latter with upper and lower limits) of the disk, as well as the distance to the object . Data come from Alcalá et al. (2014, 2016), Ansdell et al. (2016), Miotello et al. (2016a). (td) near the object name indicates a transitional disk.
|-||9a||0.439 (0.354)||0.414 (0.341)||0.035 (0.095)|
|-||9b||0.448 (0.326)||0.373 (0.295)||0.061 (0.153)|
|-||9c||0.45 (0.312)||0.335 (0.252)||0.094 (0.223)|
Notes. Variables: combination of physical parameters for which the correlation is assessed; Figure: number of the Figure in which this correlation is illustrated; : Pearson correlation; : Spearman’s rank correlation coefficient; -value: probability of a null hypothesis that the two variables are completely uncorrelated. Values in parentheses correspond to the sample with the near-brown dwarf object 2MASS-J16081497-3857145 excluded.