Protoplanetary disc ‘isochrones’ and the evolution of discs in the M˙-Md plane

Abstract

In this paper, we compare simple viscous diffusion models for the disc evolution with the results of recent surveys of the properties of young protoplanetary discs. We introduce the useful concept of ‘disc isochrones’ in the accretion rate–disc mass plane and explore a set of Monte Carlo realization of disc initial conditions. We find that such simple viscous models can provide a remarkable agreement with the available data in the Lupus star forming region, with the key requirement that the average viscous evolutionary time-scale of the discs is comparable to the cluster age. Our models produce naturally a correlation between mass accretion rate and disc mass that is shallower than linear, contrary to previous results and in agreement with observations. We also predict that a linear correlation, with a tighter scatter, should be found for more evolved disc populations. Finally, we find that such viscous models can reproduce the observations in the Lupus region only in the assumption that the efficiency of angular momentum transport is a growing function of radius, thus putting interesting constraints on the nature of the microscopic processes that lead to disc accretion.

planets and satellites: formation, planet-disc interaction – protoplanetary discs

1 INTRODUCTION

Understanding the evolution of protoplanetary discs is essential in order to study the process of star and planet formation. Since the pioneering studies of Lynden-Bell & Pringle (1974), viscous evolution has been considered as playing a fundamental role in this context, although several uncertainties in the mechanisms leading to disc evolution still remain. First, the physical origin of disc’viscosity’ is still debated, with the magneto-rotational instability (MRI) (Balbus 2003) as the main candidate in the inner, hotter parts of the disc, and the gravitational instability (Kratter & Lodato 2016) as the main candidate during the early phases of disc evolution, when the disc is still cold and massive. Recently, non-ideal magneto-hydrodynamical simulations have found that a large portion of the disc is actually stable to the MRI, and new models based on the effects of magnetic disc winds have been proposed as a way to remove, rather than redistribute, angular momentum in the disc (e.g. Bai 2016). Secondly, it is clear that other, non-viscous mechanisms play a role in disc clearing at later times, such as photoevaporation (Clarke, Gendrin & Sotomayor 2001) and planet-formation. Finally, it is well known that protoplanetary accretion proceeds ‘in bursts’, such that relatively long periods of quiescence, with moderate/low accretion rates (of the order of ∼10^−8–10 M_⊙ yr⁻¹), are followed by sudden bursts of accretion, during which the accretion rate can rise by a factor of a few (e.g. Costigan et al. 2014) or even by several orders of magnitude (Audard et al. 2014).

In the context of viscous disc evolution, the analytical self-similar solutions of Lynden-Bell & Pringle (1974) have often been used to describe to a first approximation the evolution of surface density, mass and accretion rate in a protoplanetary disc. While clearly being a simplification, and certainly not accounting for the above mentioned accretion bursts, such models are able to describe the average behaviour of protoplanetary discs in simple terms. One key parameter in these models is represented by the disc viscosity ν, whose origin, as mentioned above, is still unknown, and it is often parametrized in terms of the Shakura & Sunyaev (1973) prescription whereby

\begin{eqnarray} \nu =\alpha c_{\rm s}H, \end{eqnarray}

(1)

where c_s is the sound speed of the gas in the disc, H = c_s/Ω is the disc thickness, Ω is the disc angular velocity (assumed to be Keplerian) and α is a dimensionless scale parameter, whose value is variously estimated to be in the range 0.001–0.1 (Hartmann 1998; King, Pringle & Livio 2007).

Viscous accretion discs have been used to reproduce the general observed relationship between mass accretion rates and stellar mass. Muzerolle et al. (2005) and Natta, Testi & Randich (2006) completed the first spectroscopic surveys of mass accretion rates in young stellar objects with discs in the Taurus and ρ-Ophiuchus star forming regions. Dullemond, Natta & Testi (2006) analysed those samples in the framework of viscous accretion discs models and showed that populations of discs derived from a set of initial conditions consistent with observed prestellar cores properties could explain the observed relationship between mass accretion rates and stellar masses and between disc masses and stellar mass. Dullemond et al. (2006) also showed that their models predict a very tight correlation between the disc mass and the mass accretion rate, which was not observed in the limited set of measurements available then.

Recently, more extensive spectroscopic surveys to measure mass and accretion rates in discs in large and complete samples with an homogeneous method, such as in the Lupus or Chamaeleon I star forming regions (Alcalá et al. 2014, 2017; Manara et al. 2016a, 2017), offer a way to test such evolutionary models in a statistically coherent way, in particular when combined with millimetre-interferometry surveys of the same regions (Ansdell et al. 2016; Pascucci et al. 2016). A pioneering analysis of this kind was attempted by Hartmann (1998) in order to give estimates in α, combining the measured mass accretion rates with the age of the targets to derive how the former decreased with time. Jones, Pringle & Alexander (2012) have performed an analysis of how well do the simple self-similar models of Lynden-Bell & Pringle (1974) reproduce the evolution of accretion discs, arguing that in order to reproduce the observed mass, accretion rates and age of a sample of accretion discs, significant deviations from the self-similar models were required. The analysis of Jones et al. (2012) was limited by several factors. Their sample was collected from a number of sources, and thus resulted to be highly inhomogeneous, and the uncertainties in the determination of almost all disc observables was very high.In particular, the age of young stars is notoriously challenging to measure, as it depends on pre-main-sequence evolutionary tracks, which are very uncertain at young ages (Baraffe et al. 2002; Baraffe & Chabrier 2010; Soderblom et al. 2014). Recently, Rosotti et al. (2017) expanded the analysis by Jones et al. (2012) to study how several effects, such as the presence of dead zones in the disc, planet formation and external photoevaporation affect the evolution of discs. Also, Rafikov (2017) has used similarity solutions to estimate the viscosity based on observed disc mass, accretion rates and disc radii in Lupus.

In this paper, we perform a similar analysis but with two important differences. First, from the observational point of view, we rely on the homogeneous sample of mass accretion rates and disc masses in the Lupus star forming region obtained by Alcalá et al. (2014, 2017) and Ansdell et al. (2016) and already combined and analysed by Manara et al. (2016b). Secondly, we introduce the concept of protoplanetary disc ‘isochrones’, which is the locus of points in the |$\dot{M}\text{--}M_{\rm d}$| plane occupied by a number of sources assumed to have the same age. Analogously to the case of stellar isochrones, this analysis can test viscous evolutionary disc models and return an estimate of the age of a given cluster of protostars, thus alleviating one of the uncertainties of previous analyses.

The paper is organized as follows. In Section 2, we introduce the concept of disc isochrones, we recall the self-similar solutions and provide an analytical expression to evaluate them. In Section 3, we perform some Monte Carlo simulations, where we produce a synthetic population of protoplanetary discs and place them in the |$\dot{M}\text{--}M_{\rm d}$| plane, to check the usefulness of the isochrones. In Section 4, we apply our modelling to the data obtained by Manara et al. (2016b) in the Lupus region. In Section 6, we draw our conclusions.

2 DISC ‘ISOCHRONES’

2.1 The self-similar solutions

The (viscous) evolution of a protoplanetary disc is determined by the function Σ(R, t) that describes the surface density of the disc as a function of time t and radius R. The simplest case that is often used for protoplanetary discs is represented by the so-called ‘self-similar’ solutions, obtained by Lynden-Bell & Pringle (1974), under the assumption that disc viscosity has a simple power-law dependence on radius as follows:

\begin{eqnarray} \nu =\nu _0(R/R_0)^\gamma , \end{eqnarray}

(2)

where R₀ is a scaleradius, ν₀ is the value of viscosity at R₀ and γ is a free index. The self-similar solution is

\begin{eqnarray} \Sigma (R,t) = \frac{M_0}{2\pi R_0^2}(2-\gamma )\left(\frac{R}{R_0}\right)^{-\gamma } T^{-\eta }\exp \left(-\frac{(R/R_0)^{(2-\gamma )}}{T}\right),\nonumber\\ \end{eqnarray}

(3)

where M₀ is the disc mass at t = 0, η = (5/2 − γ)/(2 − γ), T = 1 + t/t_ν and the viscous time is defined as

\begin{eqnarray} t_{\nu }=\frac{R_0^2}{3(2-\gamma )^2\nu _0}. \end{eqnarray}

(4)

From equation (3), one thus sees that the scaleradius R₀ represents the initial truncation radius of the disc. The disc truncation radius, according to equation (3), evolves with time as R_out = R₀T^{1/(2 − γ)}. The disc mass is obtained by simply integrating equation (3) over radius:

\begin{eqnarray} M_{\rm d}(t)=M_0T^{(1-\eta )}. \end{eqnarray}

(5)

Note that in order to have a disc mass that decreases with time, we need to require γ < 2. In this scenario, the only mass that accretes on to the star is the mass lost by viscous spreading by the disc and it is thus given by the opposite of the time derivative of M_d(t):

\begin{eqnarray} \dot{M}=-\frac{\mbox{d}M_{\rm d}}{\mbox{d}t}=(\eta -1)\frac{M_0}{t_{\nu }}T^{-\eta }. \end{eqnarray}

(6)

A frequently used parameter when describing disc evolution is the ‘disc lifetime’, defined as the ratio of disc mass to accretion rate (Jones et al. 2012; Rosotti et al. 2017):

\begin{eqnarray} t_{\rm disc}=\frac{M_{\rm d}(t)}{\dot{M}_{\rm acc}(t)}=2(2-\gamma )(t+t_{\nu }). \end{eqnarray}

(7)

It is worth noting that, apart from a constant factor, the disc lifetime t_disc is a measure of the age of the disc, if t ≫ t_ν, and a measure of its viscous time, if t ≪ t_ν.

2.2 Isochrones for the self-similar solution

We define the isochrone of a population of protoplanetary discs as the locus in the |$\dot{M}\text{--}M_{\rm d}$| plane of all the discs of the same age t and same initial disc mass M₀, defined parameterically as a function of the viscous time t_ν. The isochrones for the self-similar solution depends thus on the free parameter γ. By combining equations (5) and (6),it is possible to write the isochrone, for a given initial mass M₀, explicitly as

\begin{eqnarray} \dot{M}=\frac{M_{\rm d}}{2(2-\gamma )t}\left[1-\left(\frac{M_{\rm d}}{M_0}\right)^{2(2-\gamma )}\right]. \end{eqnarray}

(8)

Some example of isochrones, for the case γ = 3/2, are shown in Fig. 1. The solid lines refer to an age of t = 10⁶ yr, for different log (M₀/ M_⊙) = −4, −3.5, −3, −2.5 and −2. The dashed lines refer to t = 10⁵ yr for log (M₀/ M_⊙) = −4 and −2, respectively.

Example of protoplanetary disc isochrones, for γ = 3/2. The solid lines show a series of isochrones at an age t = 106 yr, with different initial disc masses log (M0/ M⊙) = −4, −3.5, −3, −2.5 and −2, respectively. The dashed lines refer to t = 105 yr for log (M0/ M⊙) = −4 and −2, respectively. The red dotted lines show instead a few examples of evolutionary tracks, referring to M0 = 10−2 M⊙ and log (tν/yr) = 6, 5 and 4 (from right- to left-hand side in the plot).

Figure 1.

Example of protoplanetary disc isochrones, for γ = 3/2. The solid lines show a series of isochrones at an age t = 10⁶ yr, with different initial disc masses log (M₀/ M_⊙) = −4, −3.5, −3, −2.5 and −2, respectively. The dashed lines refer to t = 10⁵ yr for log (M₀/ M_⊙) = −4 and −2, respectively. The red dotted lines show instead a few examples of evolutionary tracks, referring to M₀ = 10⁻² M_⊙ and log (t_ν/yr) = 6, 5 and 4 (from right- to left-hand side in the plot).

Open in new tab Download slide

The isochrone can be divided into two regions: for low masses and accretion rates (on the left-hand side in Fig. 1), there is a linear relation between |$\dot{M}$| and M_d, corresponding to relatively evolved discs, with t ≫ t_ν, such that both quantities are simple power-laws with respect to time. On the right-hand side of Fig. 1, we find instead those systems that have not yet had time to evolve and for which the mass is still close to the initial mass M₀.

The disc isochrone should not be confused with the evolutionary track of a single disc with a given t_ν, which can be easily obtained by eliminating the age t between equations (5) and (6):

\begin{eqnarray} \dot{M}=\frac{1}{2(2-\gamma )}\frac{M_0}{t_{\nu }}\left(\frac{M_{\rm d}}{M_0}\right)^{5-2\gamma }, \end{eqnarray}

(9)

which defines a steeper relation in the |$\dot{M}\text{--}M_{\rm d}$| plane, which (unlike the isochrone in the self-similar part), in turn, depends on the parameter γ. A few examples of evolutionary tracks are shown in Fig. 1 as red dotted lines, referring to M₀ = 10⁻² M_⊙ and log (t_ν/yr) = 6, 5 and 4 (from the right- to left-hand side in the plot).

Note that the various isochrones evaluated at the same age, but with different initial masses, all tend to the same linear asymptotic relation, whose normalization depends on the isochrone age. Therefore, if we consider an ensemble of discs with different viscous times and initial masses, we would expect to find a tight correlation between mass and accretion rate for low disc masses. On the contrary, for higher disc masses, we would probe those discs that have evolved the less and thus keep ‘memory’ of their initial condition. In this case, the distribution of discs in the plane should maintain the original scatter and thus show a looser correlation. In order to check this behaviour, we have run some Monte Carlo realization of initial conditions, which we describe in the next section.

3 ‘Disc POPULATION SYNTHESIS’ MODELS

We have simulated initial conditions for a sample of 100 protoplanetary discs. In the following, unless otherwise stated, we will always consider the case where γ = 3/2. We have assumed a lognormal distribution of initial masses and viscous times, assuming the following parameters: 〈log (t_ν/yr)〉 = 7, |$\sigma _{t_\nu }=1$|⁠, 〈log (M₀/ M_⊙)〉 = −2.5, |$\sigma _{M_0}=0.5$|⁠, where |$\sigma _{t_\nu }$| and |$\sigma _{M_0}$| are the Gaussian widths for log t_ν and log M₀, respectively. We then evolved the population up to a fixed age t and plotted the resulting values of disc mass and accretion rates in Fig. 2, for three different ages: t = 0.1, 1 and 10 Myr, shown with the blue, green and cyan diamonds. In Fig. 2, we also plot for comparison the corresponding isochrones, assuming log (M₀/ M_⊙) = −1.5. This plot confirms the predictions that we made in the previous section: at t = 10 Myr most discs have evolved to the self-similar part of the viscous diffusion evolution and indeed display a tight correlation in the |$\dot{M}\text{--}M_{\rm d}$| plot. At earlier ages, instead, there is progressively more scatter, especially for high disc masses.

$Results of a Monte Carlo ‘disc population synthesis’ model. A Monte Carlo realization of 100 initial conditions for disc evolution, assuming a lognormal distribution in M0 and tν has been produced (see text for details). The plot shows the location of the disc population in the $\dot{M}\text{--}M_{\rm d}$ plane for t = 10 Myr (cyan diamonds), t = 1 Myr (green diamonds) and t = 0.1 Myr (blue diamonds). The three lines show the corresponding isochrones (with the same colour scheme), assuming log (M0/ M⊙) = −1.5.$

Figure 2.

Results of a Monte Carlo ‘disc population synthesis’ model. A Monte Carlo realization of 100 initial conditions for disc evolution, assuming a lognormal distribution in M₀ and t_ν has been produced (see text for details). The plot shows the location of the disc population in the |$\dot{M}\text{--}M_{\rm d}$| plane for t = 10 Myr (cyan diamonds), t = 1 Myr (green diamonds) and t = 0.1 Myr (blue diamonds). The three lines show the corresponding isochrones (with the same colour scheme), assuming log (M₀/ M_⊙) = −1.5.

Open in new tab Download slide

An important corollary of the fact that the scatter is more pronounced for high disc masses with respect to the lower disc masses (opening in a ‘fan’-like way below the corresponding isochrone) is that if one wants to obtain a correlation from the disc population, the resulting correlation will be less steep than linear, a feature that is observed, for example, in the Lupus sample (Manara et al. 2016b), where |$\dot{M}\propto M_{\rm d}^a$|⁠, with a ≈ 0.7, and even when combining the disc masses and mass accretion rates in the Lupus and in the Chamaeleon I regions (Mulders et al. 2017). However, due to the limited sensitivity of current surveys, the disc mass range spanned by observations is too little to indicate any mass dependent spread in the correlation.

Additionally, we would expect that, for an evolved sample, where most of the discs have evolved to t ≫ t_ν, the scatter in the disc mass–accretion rate plane should be small, while it should progressively increase for younger systems. We are thus left with an interesting possibility of disentangling the degeneracy between age and viscous time inherent from equation (7): while a measure of the average disc lifetime for a sample of protoplanetary discs gives a measure of the sum t + t_ν, a measure of the scatter in the correlation gives us a measure of t/t_ν. We will perform such an analysis on the Lupus data in the following Section.

4 AN APPLICATION TO THE LUPUS STAR FORMING REGION

In Fig. 3, we plot the same data as in Fig. 2, but with the observed data from the Lupus survey (Manara et al. 2016b) overlaid with red circles. The mass accretion rates shown here have been derived measuring the excess emission in the Balmer continuum region from broad-band flux-calibrated VLT/X-Shooter spectra (Alcalá et al. 2014, 2017), while the disc masses are obtained from the 0.87-mm continuum flux (Ansdell et al. 2016) and converted into total disc mass assuming a gas-to-dust ratio of 100. Note that the viscous disc models give a prediction for the gas disc masses. We have chosen to use dust masses converted into gas masses using a standard gas-to-dust ratio because accurate measurements of gas masses based on CO emission are debated (Miotello et al. 2016). Here, we include in the analysis all the discs detected at mm-wavelengths. We considered the measured values for the mass accretion rate for all these targets, including those where the measured accretion rates are compatible with the chromospheric emission (see discussion by Alcalá et al. 2017; Manara et al. 2017). The discs that are undetected at mm-wavelengths are excluded, but this does not alter the overall scatter of the data along the best fit (see fig. 1 of Manara et al. 2016b). It is apparent that the Lupus data show a significant scatter, much larger than the one expected from the models at t ≈ 1 Myr, which is of the order of the average age of the Lupus stars. In the following, we assume the age of this region to be t_Lupus = 1.6 Myr, with a spread of ≈0.3 dex, based on the estimates by Comerón (2008) or Alcalá et al. (2017), ranging from 1 to 3 Myr. However, it should be kept in mind that Fig. 3 has been produced by evolving our simulated disc population to a fixed given age, while in reality, as mentioned above, we know that there is a significant spread in ages as well in the observed sample.

Figure 3.

Same as Fig. 2, but with the data from the Lupus survey (Manara et al. 2016b) overlaid with red circles.

Open in new tab Download slide

We have thus re-run our disc population synthesis model (using 300 simulated discs) by allowing not only a spread in initial conditions (as above), but also a spread in ages. Additionally, we have also simulated an observational error on the measured masses and accretion rates by adding an additional stochastic spread of 0.1 dex for the final disc masses and 0.45 dex for the final accretion rates from the model.

We have then run a large number of models, varying the average viscous time 〈log t_ν〉 and the average age 〈log t〉, and keeping a fixed Gaussian width |$\sigma _{t_\nu }$| equal to either 1 or 1.2 dex (see below) and a given width σ_t, taken to be equal to the known 0.3 dex dispersion in Lupus ages.

In order to compare these models to the data, we can define the average disc lifetime for a given population as t_disc = 10^τ yr, where |$\tau =\langle (\log M-\log {\dot{M}})\rangle$|⁠. The scatter of the population is then defined as |$s=\langle (\log M-\log {\dot{M}}-\tau )^2\rangle$|⁠. We then compute t_disc and s for each population (either observed or obtained from our population synthesis models) and retain only those models that match with a given precision the observational parameters. We have computed t_disc,obs and s_obs of the data, obtaining t_disc,obs = 2.5 Myr and s_obs = 0.32. We have then computed t_disc and s of each model we have run, and thanks to the comparison between these values and the observed ones, we found the models that match the data.

The results of this modelling are shown in Fig. 4, where we show in purple the combinations of 〈log t_ν〉 and 〈log t〉 that match within 10 per cent (which is the typical uncertainty in the derived parameters) the observed average t_disc in Lupus, which is 〈log (t_disc,obs/yr)〉 = 6.4, and in cyan the combinations that match the observed scatter. The left-hand panel refers to the case where |$\sigma _{t_\nu }=1$|⁠, while the right-hand panel refers to |$\sigma _{t_\nu }=1.2$|⁠. The horizontal dashed line indicates the known average age in Lupus (≈1.6 Myr), while the yellow vertical line indicates the expected viscous time for a disc whose initial truncation radius was R₀ = 100 au, assuming α = 0.1 and H/R = 0.1, evaluated at R = R₀. It can be seen that the observed Lupus age falls within our confidence region, marginally for |$\sigma _{t_\nu }=1$|⁠, and more comfortably for |$\sigma _{t_\nu }=1.2$|⁠.

$Combination of disc population synthesis parameters that match within 10 per cent the observed value of the average disc lifetime $\langle \log (t_{\rm \text{disc},obs}/\mbox{yr})\rangle = 6.4$ (purple points) and the observed value of the scatter in the $\dot{M}\text{--}M_{\rm d}$ plane (cyan points). The left-hand panel refers to $\sigma _{t_\nu }=1$, while the right-hand panel refers to $\sigma _{t_\nu }=1.2$. The horizontal dashed line indicates the known average age in Lupus, while the yellow vertical line indicates the expected viscous time for a disc whose initial truncation radius was R0 = 100 au, assuming α = 0.1 and H/R = 0.1, evaluated at R = R0.$

Figure 4.

Combination of disc population synthesis parameters that match within 10 per cent the observed value of the average disc lifetime |$\langle \log (t_{\rm \text{disc},obs}/\mbox{yr})\rangle = 6.4$| (purple points) and the observed value of the scatter in the |$\dot{M}\text{--}M_{\rm d}$| plane (cyan points). The left-hand panel refers to |$\sigma _{t_\nu }=1$|⁠, while the right-hand panel refers to |$\sigma _{t_\nu }=1.2$|⁠. The horizontal dashed line indicates the known average age in Lupus, while the yellow vertical line indicates the expected viscous time for a disc whose initial truncation radius was R₀ = 100 au, assuming α = 0.1 and H/R = 0.1, evaluated at R = R₀.

Open in new tab Download slide

The interpretation of the purple region in Fig. 4 is straightforward. Models on the left horizontal purple branch correspond to t ≫ t_ν and thus the average age must reproduce the observed average disc lifetime (see equation 7). However, these models provide a scatter that is too small compared to the observed one. On the contrary, the vertical purple line in the right part of the plot corresponds to t ≪ t_ν, for which it is the viscous time that reproduces the observed disc lifetime. Such models provide a scatter that is much larger than the observed scatter. In order to provide a match to both the average disc lifetime and the observed scatter, we have to sit in an intermediate class of models, for which the age and viscous times are comparable. In particular, our confidence regions for the age and viscous time in Lupus, corresponding to the intervals that provide a match to the observed data with the required accuracy, are

\begin{eqnarray} \langle \log (t/\mbox{yr})\rangle & = & 5.9\pm 0.3, \end{eqnarray}

(10)

\begin{eqnarray} \langle \log (t_\nu /\mbox{yr})\rangle & = & 5.8\pm 0.4, \end{eqnarray}

(11)

for the |$\sigma _{t_\nu }=1$| case, and

\begin{eqnarray} \langle \log (t/\mbox{yr})\rangle & = & 6.05\pm 0.25, \end{eqnarray}

(12)

\begin{eqnarray} \langle \log (t_\nu /\mbox{yr})\rangle & = & 5.5\pm 0.5, \end{eqnarray}

(13)

for the |$\sigma _{t_\nu }=1.2$| case. Note that up to now we did not have to make any assumption on the distribution of initial truncation radii (which are built in our model through the relation between viscous time and truncation radius, equation 4). At the same time, we did not have to define any specific model for the disc viscosity, for example, in terms of α. If we want to translate the average viscous times into average initial truncation radii R₀ for the discs, we have to assume some relationship between viscosity and radius. If we assume an α viscosity, with α = 0.1 and H/R = 0.1, evaluated at R₀, we obtain 〈log (R₀/au)〉 = 2.3 ± 0.2, for |$\sigma _{t_\nu }=1$|⁠, and log (R₀/au)〉 = 2.1 ± 0.2, for |$\sigma _{t_\nu }=1.2$|⁠. Clearly these estimates of initial disc radii reflect our choice of α and are expected to be smaller for smaller α, since, for a given t_ν, we have that R₀ ∝ α^1/2.

In principle, we could have obtained solutions that bracket the observed Lupus age better by allowing an even larger |$\sigma _{t_\nu }$|⁠, and thus moving the cyan points more to the upper left-hand corner of the plots in Fig. 4, but this would have implied extremely small values for the initial radii of our simulated disc population, with 1σ values for the radii distribution corresponding to R₀ ≈ 15 au, for the case α = 0.1. While there is some observational evidence that Class I protostars are surrounded by very compact discs (Miotello et al. 2014), which might be related to magnetic braking during the main infall phase (Hennebelle et al. 2016), in order to test whether this is a reasonable value for very young discs, systematic measurements of disc sizes in very young protostars are needed.

In Fig. 5, we show our ‘best-fit’ model for |$\sigma _{t_\nu }=1$|⁠, with 〈log (t/yr)〉 = 5.9 and 〈log (t_ν/yr)〉 = 5.8, and further assuming 〈log (M₀/ M_⊙)〉 = −2.2 and |$\sigma _{M_0}=0.2$|⁠. The Lupus data, with error bars, are shown in red, while the diamonds show our disc population synthesis. It can be seen that our modelling does reproduce the observed population well.

$Distribution in the $\dot{M}\text{--}M_{\rm d}$ plane of our best-fitting model, for $\sigma _{t_\nu }=1$, with 〈log (t/yr)〉 = 5.9 and 〈log (tν/yr)〉 = 5.8, and further assuming 〈log (M0/ M⊙)〉 = −2.2 and $\sigma _{M_0}=0.2$. The Lupus data, with error bars, are shown in red, while the diamonds show our disc population synthesis.$

Figure 5.

Distribution in the |$\dot{M}\text{--}M_{\rm d}$| plane of our best-fitting model, for |$\sigma _{t_\nu }=1$|⁠, with 〈log (t/yr)〉 = 5.9 and 〈log (t_ν/yr)〉 = 5.8, and further assuming 〈log (M₀/ M_⊙)〉 = −2.2 and |$\sigma _{M_0}=0.2$|⁠. The Lupus data, with error bars, are shown in red, while the diamonds show our disc population synthesis.

Open in new tab Download slide

5 DISCUSSION

5.1 Changing the value of γ

In the sections above, we have always assumed that γ = 3/2. This has the nice property that for this value of γ:

\begin{eqnarray} t_{\rm disc}=t+t_\nu . \end{eqnarray}

(14)

Here, we discuss what happens if we allow γ to have different values. Before showing the results of the disc population synthesis calculations, we can make some analytical predictions. In the general case

\begin{eqnarray} t_{\rm disc}=2(2-\gamma )(t+t_\nu )>2(2-\gamma )t, \end{eqnarray}

(15)

or, put in other terms, the ratio of the observed average age of a given population to the average disc lifetime must satisfy

\begin{eqnarray} \frac{\langle t\rangle }{\langle t_{\rm disc}\rangle }<\frac{1}{2(2-\gamma )}. \end{eqnarray}

(16)

This condition has been first pointed out by Rosotti et al. (2017). The above equation can be turned around to give a requirement on γ for a given observed age and disc lifetime as follows:

\begin{eqnarray} \gamma >2-\langle t_{\rm disc}\rangle /2\langle t\rangle . \end{eqnarray}

(17)

If we evaluate this for Lupus, for which 〈t_disc〉 ≈ 2.5 Myr and 〈t〉 ≈ 1.6 Myr, we obtain γ ≳ 1.2. We should thus expect that models with smaller values of γ cannot be able to match the observed data.

In Fig. 6, we show the result of our modelling, assuming γ = 1 (left-hand panel) and γ = 0.5 (right-hand panel). With these values of γ, the purple line lies always below the dashed line indicating the age of Lupus, and thus our models will always require an age that is much smaller than the Lupus age, even assuming the extreme case that t ≫ t_ν.

Figure 6.

Same as Fig. 4, but for γ = 1 (left-hand panel) and γ = 0.5 (right-hand panel).

Open in new tab Download slide

5.2 Systematic uncertainties

The analysis above is based on observational data that suffer from often quite strong systematic uncertainties. The most important systematics are in the estimates of disc masses, which are based on dust measurement and thus depend on the assumed dust opacity and dust-to-gas ratio, and in the estimates of stellar ages, for which pre-main-sequence tracks are very uncertain at young ages. Disc mass uncertainties, in turn, affect the measurement of the disc lifetime t_disc ∝ M_d. This implies that our conclusions based on equation (17) above concerning the dependence of viscosity on radius can certainly be affected by such uncertainties.

However, we consider it unlikely that the disc masses are very strongly underestimated in our sample. The disc masses in our sample span a large range but are already relatively close to marginal gravitational stability. This can be seen, for example, in fig. 10 of Pascucci et al. (2016), which shows the disc-to-star mass ratio of a sample of discs obtained in the same way as in our sample, but for the Chamaleon region (that shares very similar properties to Lupus in terms of disc mass and accretion rate distributions). Thus, systematic uncertainties in disc masses would probably have the effect of decreasing the value of t_disc. At the same time, Soderblom et al. (2014) discuss the uncertainties in stellar ages and conclude that these would go in the direction of increasing the stellar age. Both effects would thus make the constraint implied by equation (17) even more stringent than γ ≳ 1.2 implied by the data as they stand.

5.3 Implications on disc evolution

In order to reproduce the observed Lupus data, as mentioned above, our models, which are based on very standard viscous accretion models, have to satisfy two less standard requirements: (i) that the average viscous time is a factor of a few larger than previously thought, and (ii) that the viscosity (and thus the surface density profile) is a steeper function of radius than previously thought.

Disc dispersal is thought to occur within a few Myr, based on the fraction of stars with infrared excess in star formation regions of different ages (Haisch, Lada & Lada 2001; Fedele et al. 2010). Disc dispersal occurs through a combination of various effects, such as planet formation, but most likely due to the interplay of viscous evolution and photoevaporation (Alexander et al. 2014), which is not included in our models. Pure viscous models usually result in dispersal time-scales that are much longer than observed, so that a reduction in the viscous time (as predicted by us) would go in the direction of alleviating the discrepancy. Clearly, a more complete model including photoevaporation would be needed in order to assess the implications of our conclusions on the fraction of stars with infrared excess at any given time.

For what concerns the steepness of the viscosity law and of the surface density profile, observational data are not conclusive. The best observations of surface density profile in the gas component are provided by the nearby TW Hya disc (Zhang et al. 2017). These authors obtain a value of γ = 0.9 ± 0.4 for TW Hya. Now, while their best-fitting value is lower than our prediction (that, we remind is a prediction for the average population and not for any individual system), our required value of γ ≳ 1.2 lies comfortably within their uncertainties. Another benchmark can be done instead by comparing our value of γ with the slope of the Minimum-Mass Solar Nebula, which is γ = 1.5, perfectly consistent with our estimates.

6 CONCLUSIONS

In this paper, we have introduced two tools that can be useful when interpreting the data obtained from current and future surveys of protoplanetary disc properties in selected star formation regions: the concept of ‘protoplanetary disc isochrone’, for example, in the |$\dot{M}\text{--}M_{\rm d}$| plane and the ‘disc population synthesis’ models.

We have specified our isochrones and disc synthesis models for the simplest class of disc evolutionary models: the self-similar solutions of Lynden-Bell & Pringle (1974). By analysing these models we reach the following conclusions:

We expect that a population of discs with similar age should show a tight correlation in the |$\dot{M}\text{--}M_{\rm d}$| plane only if the age of the system is significantly larger than the average viscous time, and that the scatter should be more prominent for larger disc masses, which are less evolved.
As a consequence of the point above, we conclude that viscous disc evolution should result in a slope in the |$\dot{M}\text{--}M_{\rm d}$| relation that is, in general, shallower than 1, as observed in the Lupus star forming region (Manara et al. 2016b), as well as combining the data from this region with those from a similar data set in the similarly young Chamaleon I region (Mulders et al. 2017), unless the discs have ages much larger than their viscous time.
Our model, despite being the simplest conceivable viscous evolution model, provides a remarkable fit to the data coming from the Lupus survey.
The slope and scatter in the |$\dot{M}\text{--}M_{\rm d}$| relation are a function of disc age, approaching a slope of 1 and a small scatter for older systems. This prediction of our modelling can be easily tested by looking at the disc populations at different ages. The recent survey of the Chamaleon region (Pascucci et al. 2016) does not reveal any difference with respect to Lupus, given that their ages are very similar. Recently, Ansdell et al. (2017) have performed an ALMA survey of the older σ −Ori region, but we still lack sufficiently accurate measurements of accretion rates for this region to properly test our predictions. Additionally, in σ-Orionis, external photoevaporation may alter the disc evolution at late times (Rigliaco et al. 2009).
It is easy to demonstrate analytically that, given the observed ages and disc lifetimes in Lupus, a fit with viscous model can only be done if the viscosity is a relatively steep function of radius, with ν ∝ R^γ, and γ ≳ 1.2. This has important implications for what concerns the mechanism of angular momentum redistribution in discs. Indeed, adopting the standard α prescription for viscosity, we have
\begin{eqnarray} \nu \propto \alpha TR^{3/2}\propto \alpha R, \end{eqnarray}
(18)
where in the last equality we have made the assumption that T ∝ R^−1/2, as usually observed. A slope larger than 1 (as required to reproduce the Lupus data with our modelling) requires that α is a growing function of radius, thus implying an increasing efficiency of angular momentum transport with increasing radius.

Mulders et al. (2017) have recently performed a similar analysis, using models with constant α but with a time dependent viscosity (Chambers 2009). Their analysis shows, consistent with our findings, that short viscous time-scales do not reproduce the observed spread in the mass accretion–disc mass relation.

Our modelling is obviously simplified and can be further expanded in several directions. From the theoretical point of view, we have just assumed the simplest viscous evolution for the discs, neglecting several important effects such as planet formation, episodic and variable accretion (Audard et al. 2014), and photoevaporation. Rosotti et al. (2017) have considered the effects of both internal and external photoevaporation, of dead zones and of planet formation, and have shown that internal photoevaporation, planet formation or dead zones would lead to a steepening of the isochrone in the low-disc-mass regime, while external photoevaporation would result in the opposite effect. These effects, when considered, would generate a spread along the relation predicted by viscous evolution. From the observational point of view, the obvious step to do is to test our models on older star-forming regions, for which we have provided clear theoretical predictions.

Acknowledgements

We thank an anonymous referee for a thorough reading of our manuscript. We thank the Munich Institute for Astro- and Particle Physics for hospitality during the ‘Protoplanetary Disks and Planet Formation and Evolution’ programme. We also thank Ilaria Pascucci and Cathie Clarke for interesting discussions. This work was partly supported by the Italian Ministero dell'Istruzione, Università e Ricerca, through the grant Progetti Premiali 2012 iALMA (CUP C52I13000140001). CFM acknowledges an European Space Agency (ESA) Research Fellowship.

REFERENCES

Alcalá

J. M.

et al. ,

2014

A&A

561

Month:	Total Views:
September 2017	16
October 2017	18
November 2017	28
December 2017	20
January 2018	32
February 2018	8
March 2018	8
April 2018	7
May 2018	2
June 2018	1
July 2018	8
August 2018	8
September 2018	1
October 2018	3
November 2018	2
December 2018	3
January 2019	3
February 2019	6
March 2019	9
April 2019	4
May 2019	2
June 2019	6
July 2019	3
August 2019	1
September 2019	2
October 2019	15
November 2019	15
December 2019	6
January 2020	6
February 2020	1
March 2020	3
April 2020	6
May 2020	3
June 2020	8
July 2020	6
September 2020	5
October 2020	3
November 2020	8
December 2020	11
January 2021	10

Article Contents

Abstract

1 INTRODUCTION

2 DISC ‘ISOCHRONES’

2.1 The self-similar solutions

2.2 Isochrones for the self-similar solution

3 ‘Disc POPULATION SYNTHESIS’ MODELS

4 AN APPLICATION TO THE LUPUS STAR FORMING REGION

5 DISCUSSION

5.1 Changing the value of γ

5.2 Systematic uncertainties

5.3 Implications on disc evolution

6 CONCLUSIONS

Acknowledgements

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Protoplanetary disc ‘isochrones’ and the evolution of discs in the M˙-M_d plane