On the time evolution of the M_d−M_⋆ and Ṁ–M_⋆ correlations for protoplanetary discs: the viscous time-scale increases with stellar mass

ABSTRACT

Large surveys of star-forming regions have unveiled power-law correlations between the stellar mass and the disc parameters, such as the disc mass |$M_{\mathrm{d}} \!-\! {M_{\star }}$| and the accretion rate |$\dot{M} \!-\! {M_{\star }}$|⁠. The observed slopes appear to be increasing with time, but the reason behind the establishment of these correlations and their subsequent evolution is still uncertain. We conduct a theoretical analysis of the impact of viscous evolution on power-law initial conditions for a population of protoplanetary discs. We find that, for evolved populations, viscous evolution enforces the two correlations to have the same slope, λ_m = λ_acc, and that this limit is uniquely determined by the initial slopes λ_{m, 0} and λ_{acc, 0}. We recover the increasing trend claimed from the observations when the difference in the initial values, δ₀ = λ_{m, 0}−λ_{acc, 0}, is larger than 1/2; moreover, we find that this increasing trend is a consequence of a positive correlation between the viscous time-scale and the stellar mass. We also present the results of disc population synthesis numerical simulations, that allow us to introduce a spread and analyse the effect of sampling, which show a good agreement with our analytical predictions. Finally, we perform a preliminary comparison of our numerical results with observational data, which allows us to constrain the parameter space of the initial conditions to λ_{m, 0} ∈ [1.2, 2.1], λ_{acc, 0} ∈ [0.7, 1.5].

accretion, accretion discs, planets and satellites: formation, protoplanetary discs

1 INTRODUCTION

Protoplanetary discs are the cradle of planets. Their evolution and dispersal strongly impact the outcome of planet formation; they especially affect the extent and availability of planetesimals, the building blocks of planets (Morbidelli et al. 2012; Mordasini et al. 2015).

Discs also serve as a mass reservoir for the central accreting protostar. For accretion to take place, material stored in the disc needs to lose most of its angular momentum. The trigger to this process is conventionally identified as a macroscopic viscosity; the pioneering work of Lynden-Bell & Pringle (1974), based on the α prescription of Shakura & Sunyaev (1973), set the ground for numerous following studies treating accretion as a redistribution of angular momentum within the disc. Despite being by far the most widely used, viscosity is not the only accretion theory; several studies have suggested MHD winds as promising candidates to explain protoplanetary disc accretion (Lesur, Kunz & Fromang 2014; Bai 2017; Béthune, Lesur & Ferreira 2017; Lesur 2020; Tabone et al. 2022a,b). In this scenario, angular momentum is removed instead of being redistributed, leading to significant differences in the evolutionary predictions.

Thanks to the great technological development of the last decades, and in particular to the advent of facilities like the Atacama Large Millimeter Array (ALMA), observational data allow to test evolutionary models. Extended data sets collecting information on a large number of young stellar objects provide the ideal ground to test theoretical predictions; the observational focus is therefore on surveys of entire star-forming regions. A number of such surveys have been already carried out (for example; Ansdell et al. 2016, 2017; Barenfeld et al. 2016; Pascucci et al. 2016; Testi et al. 2016, 2022; Alcalá et al. 2017; Manara et al. 2017, 2020; Cieza et al. 2019; Williams et al. 2019; Sanchis et al. 2020), unveiling interesting features and patterns, such as power-law correlations between the properties of discs and their host stars.

Two crucial steps are required to test evolutionary models: first, performing numerical simulations of the different prescriptions; and secondly, a comparison of these numerical results with data. Identifying key predictions for each model allows to distinguish between different scenarios. The viscous case shows a characteristic behaviour, known as viscous spreading: as part of the disc mass loses angular momentum and drifts inwards, accreting the protostar, another portion of the disc gains the same amount of angular momentum instead, expanding towards larger radii. Therefore, the radial extent of the disc increases with time, despite the ongoing accretion. This prediction does not apply to the MHD winds scenario (Trapman et al. 2022); its removal of angular momentum causes disc radii to decrease as evolution proceeds. Whether the available data in the Lupus star-forming region agree with the viscous spreading predictions has been debated in the literature (Sanchis et al. 2020; Trapman et al. 2020; Toci et al. 2021). Unfortunately, detecting viscous spreading using dust radius as a tracer is non-trivial: Rosotti et al. (2019a) showed that this method requires significantly deep observations, targeting a fraction as high as 95|${{\,\rm per\,cent}}$| of the total flux.

The disc mass–accretion rate correlation provides another diagnostic criterion. In the purely viscous case, where the self-similar solution (Lynden-Bell & Pringle 1974) holds, such a correlation naturally stems from the analytical prescriptions for M_d and |$\dot{M}$|⁠; in particular, it is expected to have slope of ∼1 (Hartmann et al. 1998; Lodato et al. 2017; Mulders et al. 2017; Rosotti et al. 2017) and a spread decreasing in time (Lodato et al. 2017), as the age of the population reaches and then outgrows the viscous time-scale. The analysis of the first data from Lupus and Chameleon (Manara et al. 2016b) showed a possible agreement with this prediction. However, Tabone et al. (2022b) showed that an MHD disc winds model could be tuned to reproduce the |$M_\mathrm{d} \!-\! \dot{M}$| correlation equally well, both in slope and spread – making it more challenging to distinguish between these two models solely using this relation. Moreover, there are some behaviours that cannot be fully explained by any of the two scenarios alone: an example is the Upper Scorpius star-forming region (Manara et al. 2020), where the observational data show a large scatter in the |$M_{\mathrm{d}} \!-\! \dot{M}$| relation (see also Testi et al. 2022). Additional mechanisms have been invoked to explain these inconsistencies, such as internal and external photoevaporation (Rosotti et al. 2017; Somigliana et al. 2020; Sellek, Booth & Clarke 2020b) and dust evolution (Sellek, Booth & Clarke 2020a); however, understanding the right combination of processes to retrieve the observed spread is non-trivial.

These similar behaviours, and the difficulties in observing the key prediction of viscous spreading, may sound discouraging. However, surveys also provide additional information – namely, the correlation between stellar and disc parameters. Many independent works (Muzerolle et al. 2003; Mohanty, Jayawardhana & Basri 2005; Natta, Testi & Randich 2006; Herczeg & Hillenbrand 2008; Alcalá et al. 2014; Kalari et al. 2015; Manara et al. 2017) have found a correlation between the disc accretion rate |$\dot{M}$| and the stellar mass M_⋆, with a slope of ∼1.8 ± 0.2; on the other hand, the disc mass versus stellar mass correlation appears to steepen with time (Ansdell et al. 2017). Some attempts to explain both the |$\dot{M} \!-\! M_{\star }$| correlation (Alexander & Armitage 2006; Clarke & Pringle 2006; Dullemond, Natta & Testi 2006; Ercolano et al. 2014) and the M_d−M_⋆ correlation (Pascucci et al. 2016; Pinilla, Pascucci & Marino 2020) have been made; none the less, it is still not clear whether these correlations are determined by the initial conditions, or rather established later on in the disc lifetime as a consequence of the evolutionary processes. Whether the time evolution of these correlations can be understood in the context of viscously evolving disc populations has not been investigated so far.

The numerical counterpart of star-forming regions surveys is population syntheses, i.e. generating and evolving synthetic populations of discs through numerical methods. Performing population syntheses is particularly useful to test evolutionary models, as well as the impact of different physical effects on the diagnostic quantities of interest (such as disc masses and radii). Population syntheses have been carried out already, investigating different aspects of evolution and dispersal of protoplanetary discs (Lodato et al. 2017; Somigliana et al. 2020; Sellek et al. 2020a,b); however, none of them included both a proper Monte Carlo drawing of the involved parameters and the correlations between disc and stellar properties. The usual assumptions were a fixed stellar mass and a linear span of the parameter space, which make a good first approximation but lack the spread and statistics that really make population syntheses a powerful tool. In this paper, we employ a new and soon to be released python code, Diskpop, which performs a coherent population synthesis of protoplanetary discs and sets the basis for future developments, taking into account more and more physical effects acting on discs.

In this work, we aim to study the dependency of disc properties on the stellar mass M_⋆, discussing its implications from the evolutionary point of view. These observed correlations are most likely linked to both evolution and initial conditions, and our goal is to disentangle between the two. We investigate the case where the correlations are already present as initial conditions. With this first paper we, focus on setting up the framework for population syntheses: we limit our case study to purely viscous discs, but the natural progression is to include more evolutionary predictions to compare (and hopefully distinguish) between each other. The structure of the paper is as follows: In Section 2, we report and discuss state of the art of the observational evidence on disc masses, accretion rates, and radii; in Section 3, we present our analytical considerations on the effects of viscous evolution on power-law initial correlations between the disc parameters and the stellar mass; in Section 4, we discuss our population synthesis model, including both the implementation in Diskpop and the numerical results; we compare our results with observational data and discuss their implications in Section 5, and finally, we draw the conclusions of this work in Section 6.

2 SUMMARY OF OBSERVATIONAL EVIDENCE

2.1 Disc mass

Numerous surveys of star-forming regions, where disc masses were determined by observing the sub-mm continuum emission of the dust component (Ansdell et al. 2016, 2017; Barenfeld et al. 2016; Pascucci et al. 2016; Testi et al. 2016, 2022; Sanchis et al. 2020), have highlighted a power-law correlation between the disc mass M_dust (of the dusty component) and the stellar mass M_⋆. This relationship can be parametrized as linear in the logarithmic plane,

$$\begin{eqnarray} \log _{10} \left(\frac{M_{\mathrm{dust}}}{M_{\oplus }} \right) = \lambda _{\mathrm{m}, \mathrm{obs}} \log _{10} \left(\frac{M_*}{\mathrm{M}_{\odot }} \right) + q_{\mathrm{obs}} + \varepsilon _{\mathrm{obs}}, \end{eqnarray}$$

(1)

with slope λ_{m, obs} and intercept q_obs.¹ε_obs is a Gaussian random variable with mean 0, representing the scatter of the correlation.

The values of the parameters λ_{m, obs} and q_obs can be determined from observations by fitting the correlation between M_dust and M_⋆. Ansdell et al. (2017) and Testi et al. (2022) (hereafter A17 and T22, respectively) performed this fit for different star-forming regions; the results are shown in Table 1 (same as table 4 in A17 and table H.1 in T22), where Δ_obs represents the intrinsic dispersion (namely, the standard deviation of the distribution of ε_obs).

Table 1.

Fitted values of λ_m,obs, q_obs and Δ_obs (see Section 2.1 for details) for different star-forming regions as performed by A17 (top table) and T22 (bottom table).

Region	Age [Myr]	λm, obs	qobs	Δobs
Taurus	1–2	1.7 ± 0.2	1.2 ± 0.1	0.7 ± 0.1
Lupus	1–3	1.8 ± 0.4	1.2 ± 0.2	0.9 ± 0.1
Cha I	2–3	1.8 ± 0.3	1.0 ± 0.1	0.8 ± 0.1
σ Orionis	3–5	2.0 ± 0.4	1.0 ± 0.2	0.6 ± 0.1
Upper Sco	5–11	2.4 ± 0.4	0.8 ± 0.2	0.7 ± 0.1

Region	Age [Myr]	λm, obs	qobs	Δobs
Corona	0.6	1.3 ± 0.5	0.4 ± 0.4	1.1 ± 0.7
Taurus	0.9	1.5 ± 0.2	1.1 ± 0.1	0.8 ± 0.3
L1668	1	1.5 ± 0.2	1.0 ± 0.1	0.8 ± 0.3
Lupus	2	1.7 ± 0.3	1.4 ± 0.2	0.7 ± 0.3
Cha I	2.8	1.6 ± 0.3	1.1 ± 0.2	0.7 ± 0.4
Upper Sco	4.3	2.2 ± 0.3	0.8 ± 0.2	0.7 ± 0.3

Table 1 shows that the mean value of q_obs tends to get lower and lower with time. This behaviour, which is more visible in the top panel, is an intrinsic characteristic of the standard viscous scenario: as discs get older, part of their mass is lost due to the ongoing accretion of the central protostar, which eventually depletes the disc. However, it is worth pointing out that some star-forming regions appear not to follow this trend (Cazzoletti et al. 2019; Williams et al. 2019), and the reason behind that is still unclear.

On the other hand, the mean value of λ_{m, obs} appears to be increasing with time, implying a steepening of the correlation between the (dust) disc mass and the stellar mass. Investigating the mathematical origin, physical meaning, and expected evolution of such a trend is one of the main goals of this paper, and will be addressed in Section 3.

2.2 Disc accretion rate

The main signatures of young stars accreting material from their surrounding disc can be found in their spectra. Gas falling on to the stellar surface along the magnetic field lines (Calvet & Gullbring 1998) causes an excess emission, particularly visible in the UV area of the spectrum (and especially in the Balmer continuum, see Gullbring et al. 1998). Characteristic emission line profiles are also typical indicators of accretion. Modelling the Balmer continuum excess in the spectra of young stars, and fitting emission line profiles, provide effective ways of measuring accretion rates.

Numerous surveys focusing on different star-forming regions have targeted accretion rates (Muzerolle et al. 2003; Natta et al. 2004; Mohanty et al. 2005; Dullemond et al. 2006; Herczeg & Hillenbrand 2008; Rigliaco et al. 2011; Manara et al. 2012, 2016a, 2017, 2020; Alcalá et al. 2014, 2017; Kalari et al. 2015; Testi et al. 2022). Many of them have found a power-law correlation between the accretion rate and the stellar mass, |$\dot{M} \propto {M_{\star }}^{\lambda _{\mathrm{acc}, \mathrm{obs}}}$|⁠. The best-fitting value of λ_{acc, obs} ≈ 1.8 ± 0.2 seems to be roughly constant throughout different regions, suggesting that it could be independent on age (unlike the M_d−M_⋆ correlation, see Section 2.1). On the contrary, Manara et al. (2012) (hereafter M12) do see an increasing trend of λ_{acc, obs} with the age of the population, a trend similar to the one claimed by A17 with respect to λ_{m, obs}. There is, however, a significant difference between this latter work and the others mentioned: M12 analysed a sample of ∼700 stars in the single star-forming region of the Orion Nebula Cluster, determining the isochronal age of each object in the region independently. On the other hand, A17 and T22 considered different star-forming regions and assumed all of the objects in each of them to be coeval – determining therefore the mean age of each sample. The ages of young stellar objects are usually determined by comparing their position on the Hertzsprung–Russell (H–R) diagram with theoretical isochrones (see Soderblom et al. 2014 for a review): however, this method comes with a series of caveats (Preibisch 2012). In particular, translating a spread in luminosity into a spread in ages is not straightforward due to a number of factors that can impact the shape of the H–R diagram, such as measurement uncertainties and variations of the accretion processes. Determining a mean age for a whole star-forming region absorbs part of this uncertainty, which is the reason why the approach of M12 is less used. Nonetheless, their results intriguingly show a behaviour of λ_{acc, obs} similar to that of λ_{m, obs}, and therefore represent a case study worth considering. In this work, when compared to observational data (see Section 5.1), we will consider the fits for λ_{acc, obs} obtained from M12 as well as T22. Table 2 summarizes the fitted values from both works.

Table 2.

Fitted values of the |$\dot{M} - M_{\star }$| slope, λ_{acc, obs} (see Section 2.2 for details), by M12 (top table) and T22 (bottom table).

Age [Myr]	λacc, obs
0.8	1.15 ± 2.00
1	1.26 ± 2.02
2	1.61 ± 2.06
5	2.08 ± 2.13
8	2.32 ± 2.16
10	2.43 ± 2.18

Region	Age [Myr]	λacc, obs
L1668	1.0	1.8 ± 0.5
Lupus	2.0	1.6 ± 0.3
Cha I	2.8	2.3 ± 0.3
Upper Sco	4.3	1.5 ± 0.8

As discussed for λ_{m, obs}, there is in principle no theoretical reason for the correlation of the accretion rate with the stellar mass to steepen or flatten in time. However, if we assume the viscous framework to hold, we do expect the accretion rate to be a decreasing function of time; in particular, Hartmann et al. (1998) showed that the viscous evolution implies |$\dot{M} \propto t^{- \eta }$|⁠, where η ∼ 1.5.

2.3 Disc radius

Evolutionary models predict that disc sizes to vary with time. In particular, radial drift is expected to influence the size of the dust discs (Weidenschilling 1977): large dust grains drift inwards, eventually disappearing, while small grains are left behind and follow the motion of the gaseous component. In the viscous framework, gas is subject to the so-called viscous spreading: as a consequence of the conservation of angular momentum, the accretion on to the central star leads to a radial expansion of the disc. Once the large grains are removed, discs are expected to be wide and faint (Rosotti et al. 2019b) and may be challenging to observe. These caveats must be taken into account when discussing disc radii, and, in practice, may severely limit the ability to observe viscous spreading, even when it does take place (Toci et al. 2021).

The radial size of the dust component in discs is measured by analysing the extent of the millimetric thermal continuum emission. In the Ophiucus (Cox et al. 2017; Cieza et al. 2019), Lupus (Ansdell et al. 2016; Tazzari et al. 2017; Andrews et al. 2018; Hendler et al. 2020) and Taurus (Long et al. 2019; Kurtovic et al. 2021) star-forming regions, this method has been widely employed. Andrews et al. (2018) found evidence of a correlation between the dust disc radius and the stellar mass, namely R_d ∝ M_⋆^0.6; recent works (Andrews et al. 2018; Sanchis et al. 2021) have also found a correlation between the disc radius and the disc dust mass, which could be used to derive additional correlations with the stellar mass. However, as measuring radii requires to spatially resolve the discs, the sample of objects with measured R_d is smaller than that with measured M_d; moreover, both those measurements carry significant uncertainties due to optically thick emission. This makes it not convenient to prefer this relation to R_d − M_⋆.

Using gas tracers, such as the rotational line emission of the ¹²CO molecule, one can also measure the gas disc size (R_CO). As these observations are very time consuming, less data are available for the gaseous component (Barenfeld et al. 2016; Ansdell et al. 2018; Sanchis et al. 2021). In their work, Ansdell et al. (2018) did not see any correlation between the gas disc size and the stellar masses in Lupus, but their sample is biased towards the highest mass discs around the highest mass stars. At the present time, there is no strong evidence of a correlation between the disc gas size and the stellar mass; however, if a correlation exists, it is probably positive (see Appendix  A).

3 ANALYTICAL CONSIDERATIONS

In the previous section, we have analysed the observational evidence for correlations between the stellar mass and the three major disc properties – mass, accretion rate, and radius. We have presented the possibility, discussed in previous works, to describe these correlations (with the possible exception of R_d − M_⋆) with power-laws. In this section, we conduct a theoretical analysis of these correlations and their evolution in the purely viscous framework. We aim at determining the initial conditions needed to recover the observational findings.

We start with two assumptions:

• Self-similar discs: we assume that all discs in the population only evolve under the effect of the viscosity ν. ν can be modelled as a power-law of the disc radius, |$\nu = \nu _{\rm c} \left(\frac{R}{R_{\rm c}} \right)^{\gamma }$|⁠, where ν_c = ν(R = R_c) and R_c is the exponential cutoff radius. The value of ν is prescribed following Shakura & Sunyaev (1973) as ν = αc_sH (where α is a dimensionless parameter, c_s is the sound speed, and H is the height of the disc). This implies that the self-similar solution by Lynden-Bell & Pringle (1974) describes the evolution of the gas component in protoplanetary discs,

  $$\begin{eqnarray} \Sigma (R, t) &=& \frac{M_{\mathrm{d}}(0)}{2 \pi {R_{\rm c}}^2}(2-\gamma) \left(\frac{R}{R_{\rm c}} \right)^{-\gamma } T^{-\eta } \\&&\times \,\exp \left(- \frac{(R/R_c)^{2-\gamma }}{T} \right), \end{eqnarray}$$

  (2)

  where η = (5/2 − γ)/(2 − γ), T = 1 + t/t_ν and t_ν is the viscous time-scale, defined as t_ν = R_c²/(3ν_c), on which viscous processes leading to evolution of discs take place.
• Power-law initial conditions: we assume power-law correlations with the stellar mass as initial conditions for M_d, |$\dot{M}$|⁠, and R_c. The value of the slopes can vary as evolution takes place, and we will refer to the evolved values as λ_m, λ_acc, and ζ. The initial correlations are set as follows:

  $$\begin{eqnarray} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}M_\mathrm{d} (0) \propto {M_{\star }}^{\lambda _{\mathrm{m}, 0}} \\ \dot{M} (0) \propto {M_{\star }}^{\lambda _{\mathrm{acc}, 0}} \\ R_c (0) \propto {M_{\star }}^{\zeta _0}. \end{array}\right. \end{eqnarray}$$

  (3)

The self-similar solution provides analytical expressions for both the disc mass and accretion rate

$$\begin{eqnarray} M_\mathrm{d} (t) &=& M_{\mathrm{d}, 0} \left(1 + \frac{t}{t_{\nu }} \right)^{1 - \eta },\\ \dot{M} (t) &=& (\eta - 1) \frac{M_{\mathrm{d}, 0}}{t_{\nu }} \left(1 + \frac{t}{t_{\nu }} \right)^{- \eta }. \end{eqnarray}$$

(4),(5)

Note that as the ratio |$M_{\mathrm{d}} / \dot{M}$| has the dimension of a time, by choosing the initial M_d and |$\dot{M}$|⁠, we completely specify their time evolution. Depending on the age of the disc with respect to its viscous time-scale, set by the radius of the disc itself, discs can be considered young (t ≪ t_ν) or evolved (t ≫ t_ν). In these two limits, it is possible to simplify equations (4) and (5). We divide the following discussion based on these two evolutionary scenarios, as they lead to different results and theoretical expectations (Sections 3.1 and 3.2).

Equations (4) and (5) make clear that the time evolution of the slopes depends on how t_ν scales with stellar mass. Because at t = 0 the initial accretion rate can be written as |$\dot{M}_0 \propto M_{\mathrm{d}, 0}/t_{\nu }$|⁠, a scaling of t_ν is already implicitly assumed given a choice of λ_{m, 0} and λ_{acc, 0}. We can see from the definition of the viscous time-scale that the dependence on stellar mass is contained in any of the three parameters α, c_s, or H. In this paper, we assume α to be a constant across all stellar masses; moreover, we also consider it as constant in radius and time, and fix its value to 10⁻³. However, note that the assumption that α does not depend on stellar mass does not affect our general results on the evolution of the slopes (see the end of Section 3.2), since this depends only on the scaling of viscous time-scale with stellar mass. On the other hand, this assumption does affect how aspect ratio H/R and disc radius scale with stellar mass. Regarding the aspect ratio, this quantity can be defined as

$$\begin{eqnarray} \frac{H}{R} = \frac{c_s}{v_k},\\ \end{eqnarray}$$

where v_k is the keplerian velocity. Assuming a radial temperature profile T ∝ R^−1/2, then c_s ∝ R^−1/4 so that H/R will scale as R^1/4. Note that this implies γ = 1, which will be the case from now on. We can parametrize the dependency on M_⋆ as a power-law with exponent β

$$\begin{eqnarray} \frac{H}{R} = \frac{H}{R}\big |_{R=R_0} \left(\frac{R}{R_0} \right)^{1/4} \left(\frac{M_{\star }}{\mathrm{M}_{\odot }} \right)^{\beta }. \end{eqnarray}$$

(6)

In this work, we have considered three possible values for β, which correspond to different physical situations:

• H/R does not depend on the stellar mass, which implies β = 0;

• c_s does not depend on the stellar mass, implying β = −1/2;

• T ∝ M_⋆^0.15 (derived by radiative transfer models, see Sinclair et al. 2020), which implies β = −0.425. Physically, this means that, while radiative transfer predicts that discs around lower mass stars should be colder due to their lower luminosity, this is a sub-dominant effect in setting the scale-height: the weaker gravity of these stars drives most of the change in scale-height.

Of these three possibilities, the third one is the more realistic; none the less, the resulting value of β is very close to −1/2, meaning that we can approximate it to the second case, which makes mathematical calculations more straightforward. On the other hand, despite not being realistic, the first scenario is usually assumed for the sake of simplicity, and for this reason, we decided to include it in this work. In conclusion, we considered β to be either equal to 0 or 1/2.

3.1 Young populations: t ≪ t_ν

If the age of the population is much smaller than its viscous time-scale, no evolution has taken place yet. This means that we are still observing the initial condition, and that the parameters λ_m, λ_acc, and ζ coincide with λ_m,0, λ_acc,0, and ζ₀. However, in the viscous scenario these parameters are always linked to one another; for young discs, Equation (5) can be written as |$\dot{M}(0) \approx M_{\mathrm{d}}(0)/t_{\nu }$|⁠. In our notation, M_d represents the total disc mass, which is given by |$99 {{\,\rm per\,cent}}$| gas mass and |$1 {{\,\rm per\,cent}}$| dust mass; since the observed power-law correlation (equation 1) refers to the dust mass, we need to translate M_dust to M_d dividing by the dust-to-gas ratio, ε, which is assumed to be constant in time and M_⋆. The initial, theoretical disc mass can therefore be written as

$$\begin{eqnarray} M_{\mathrm{d}}(0) = K_1 \left(\frac{M_{\star }}{\mathrm{M}_{\odot }} \right)^{\lambda _{\mathrm{m}, 0}}, \end{eqnarray}$$

(7)

where the normalization constant K₁ is given by 10^qM_⊕/ε. Taking into account the dependence of the aspect ratio of the disc H/R on the stellar mass through equation (6), and using a generic value for β, we can write the initial accretion rate in the self-similar scenario as

$$\begin{eqnarray} \dot{M}(0) = K_2 \left(\frac{M_{\star }}{\mathrm{M}_{\odot }} \right)^{\lambda _{\mathrm{m}, 0} + 2 \beta + 1/2 - \zeta _0}, \end{eqnarray}$$

(8)

defining the normalization constant K₂ as

$$\begin{eqnarray} K_2 = \frac{10^q M_{\oplus }}{\varepsilon } \sqrt{\frac{G \mathrm{M}_{\odot }}{R_{c, \odot }}} \frac{3}{2} \alpha \left(\frac{H}{R}\big |_{R=R_{c, \odot }} \right)^2, \end{eqnarray}$$

(9)

where R_{c, ⊙} is the value of the cut-off radius R_c for a solar-type star. The only free parameter in equation (8) is the stellar mass M_⋆: all of the other parameters (λ_{m, 0}, ζ₀, R_{c, ⊙}, β) can be determined by comparing this expression to the observational data. In particular, since |$\dot{M} \propto {M_{\star }}^{\lambda _{\mathrm{acc}, 0}}$|⁠, assuming that we can determine the values of λ_{m, 0} and λ_{acc, 0} from observations and that α is known, we could constrain the value of ζ₀,

$$\begin{eqnarray} \zeta _0 \equiv \lambda _{\mathrm{m}, 0} + 2 \beta + \frac{1}{2} - \lambda _{\mathrm{acc}, 0}. \end{eqnarray}$$

(10)

On the other hand, R_{c, ⊙} can be set imposing the normalization constant K₂ to equal the typical value for the accretion rate of a solar-type star. Based on the observation of Manara et al. (2017), |$\log _{10}(\dot{M} / \mathrm{M}_{\odot } \text{yr}^{-1}) = -8.44$|⁠, which leads to R_{c, ⊙} = 1.6 × 10¹⁴ cm ≈ 10 au (assuming α = 10⁻³, as stated above). This is a reasonable initial disc radius; various recent papers have also assumed radii of the same order of magnitude (Trapman et al. 2019; Rosotti et al. 2019a,b; Toci et al. 2021).

3.2 Evolved populations: t ≫ t_ν

If the age of the population is larger than its viscous time-scale, discs can be considered evolved. In principle, at this stage, the initial conditions are not observed anymore – meaning that we expect to see an evolution of both the slopes and spreads of the correlations. In particular, if t ≫ t_ν, the self-similar accretion rate reduces to

$$\begin{eqnarray} \dot{M} (t) = \frac{1}{2} \frac{M_\mathrm{d}(t)}{t}; \end{eqnarray}$$

(11)

as t does not depend on the stellar mass, we can see that equation (11) implies that the dependency on M_⋆ in evolved discs is the same for |$\dot{M}$| and M_d. Moreover, the distributions of the two quantities, M_d and |$\dot{M}$|⁠, are expected to resemble each other more and more as evolution takes place and the condition given in equation (11) is reached: this means that, eventually, we expect both the slope and the spread of the two correlations to reach the same value.

3.2.1 Slopes

As we did in the young populations scenario, expressing M_d(t) through the first assumption in equation (3) and using the definition of t_ν leads to

$$\begin{eqnarray} M_\mathrm{d} (t) = K_3 \left(\frac{M_{\star }}{\mathrm{M}_{\odot }} \right)^{\lambda _{\mathrm{m}, 0}+\zeta _0/2-\beta -1/4} t^{-1/2}, \end{eqnarray}$$

(12)

where K₃ is a normalization constant that does not depend on the stellar mass or the age, defined by

$$\begin{eqnarray} K_3 = M_{\mathrm{d}, \odot } \left(\frac{{R_{c, \odot }}^3}{G \mathrm{M}_{\odot }} \right)^{1/4} \frac{1}{\sqrt{3 \alpha }} \left(\frac{H}{R}\big |_{R=R_{c, \odot }} \right)^{-1}; \end{eqnarray}$$

M_{d, ⊙} represents the initial disc mass for a solar-type star. Denoting the evolved slopes, at t ≫ t_ν, as λ_m and λ_acc, this implies that

$$\begin{eqnarray} \lambda _{\mathrm{m}} \equiv \lambda _{\mathrm{acc}} = \lambda _{\mathrm{m}, 0} + \frac{\zeta _0}{2} - \beta - \frac{1}{4}. \end{eqnarray}$$

(13)

There are two ways in which equation (13) can be satisfied:

• λ_{m, 0} = λ_{acc, 0}: if the two slopes start from the same value, they do not evolve with time.

• λ_{m, 0} ≠ λ_{acc, 0}: if the initial values of the two slopes are different, an evolution of their values must take place due to accretion processes.

In the first scenario, as we do not expect any evolution with time of the two slopes, we can use the initial condition (see equation 10) to determine the value of ζ₀, finding ζ₀ = 2β + 1/2. As expected, substituting this value of ζ₀ in equation (13) leads to λ_m = λ_acc = λ_{m, 0}. This implies that observing the same slope for the two correlations M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$| is not enough to claim anything on the age of the population. It could either be evolved, in which case we are observing a direct consequence of viscous evolution, or it could be young, in which case we would be observing the initial condition. This endorses what we suggested earlier, that observed correlations can be due either to the initial conditions or the evolution or, most likely, a combination of the two. It would be possible to disentangle between the two possibilities if additional information was provided, for example

• observations at earlier ages (to try and see a change in the slopes with time);

• measurements of disc radii R_i, which is in principle observable (although with some caveats, see Toci et al. 2021) to determine the viscous time-scale and give an estimate of the evolutionary stage of the population. Note that there is a degeneracy in the determination of the viscous time-scale, as it depends on α, R_c, and |$H/R \big |_{R = R_c}$|⁠;

• the |$\dot{M} \!-\! M_{\mathrm{d}}$| correlation: a large spread implies that the population cannot be considered evolved yet (Lodato et al. 2017).

The other possibility is that the two initial values of the slopes are different. If that is the case, we can define the difference between the initial values as δ₀ = λ_{m, 0} − λ_{acc, 0}. This means that the initial condition for ζ₀ (equation 10) is given by ζ₀ = 1/2 + 2β + δ₀, which can be used in equation (13) to give

$$\begin{eqnarray} \lambda _{\mathrm{m}} = \lambda _{\mathrm{acc}} = \frac{3 \lambda _{\mathrm{m}, 0} - \lambda _{\mathrm{acc}, 0}}{2}. \end{eqnarray}$$

(14)

Note that equation (14) is not symmetric in λ_{m, 0} and λ_{acc, 0}, meaning that the impact of the initial condition for the M_d−M_⋆ correlation is more long-lived than that for |$\dot{M} \!-\! M_{\star }$|⁠. Whether λ_{m, 0} and λ_{acc, 0} grow steeper or shallower with time, towards the limit value determined by equation (14), depends on the sign of δ₀:

• A steepening of the correlation, λ_m > λ_{m, 0}, is achieved if λ_{m, 0} > λ_{acc, 0} (δ₀ > 0); this condition also implies λ_acc > λ_{acc, 0}.

• On the other hand, a flattening of the correlation, λ_m < λ_{m, 0}, is obtained if λ_{m, 0} < λ_{acc, 0} (δ₀ < 0); as in the previous case, this also implies λ_acc < λ_{acc, 0}.

In conclusion, there are two possible scenarios: either both slopes increase towards the common value for evolved populations, if λ_{m, 0} > λ_{acc, 0}, or they both decrease towards it, it λ_{m, 0} < λ_{acc, 0}. In both cases, the difference between the two slopes δ = λ_m − λ_acc will tend towards zero. This argument stems from purely mathematical considerations, but can be easily understood from the physical point of view. As the initial accretion time-scale of the disc can be written as the ratio |$M_{\mathrm{d}} / \dot{M}$|⁠, it will itself show a power-law correlation with the stellar mass. In particular, its slope will be the difference of those of M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$|⁠; in our notation, its initial value is δ₀. A positive correlation (i.e. a positive δ₀) means that discs around more massive stars will evolve slower than discs around less massive stars, causing a steepening of the correlations between disc and stellar parameters. The contrary argument applies for a negative δ₀, which leads to flattening correlations. As the evolution proceeds, δ = λ_m − λ_acc → 0, eventually reaching a homologous depletion of discs around more and less massive stars. Fig. 1 illustrates this concept in form of sketch.

Figure 1.

Illustrative sketch to show how steepening slopes are a consequence of a positive correlation between viscous time-scales and stellar masses. The dots represent artificial data sets with their linear fits displayed as solid lines. The left-hand panel shows an early evolutionary stage, where |$t_{\nu } \propto {M_{\star }}^{\delta _0}$| with δ₀ > 0: as discs around more massive stars have a longer viscous time-scale, in the time interval Δ_t = t₂ − t₁ their mass is less depleted than that of discs around less massive stars, which in turn have a short viscous time-scale. This implies that, after Δ_t, the correlation between M_d and M_⋆ is steeper. On the other hand, the right-hand panel represents a later evolutionary stage, where δ has reached the limit value of 0. In that case, the viscous time-scale does not depend on the stellar mass anymore, and this implies a homologous depletion of discs around stars of all masses. The same argument applies for |$\dot{M}$|⁠.

These conclusions on the evolutionary behaviour of λ_m and λ_acc only depend on the scaling of the viscous time-scale with the stellar mass. This is the reason why our assumption that α does not depend on M_⋆ does not influence this result. It should be noted, however, that if α did depend on the stellar mass, we would find a different relation between λ_{m, 0}, λ_{acc, 0}, and ζ₀, as the viscous time-scale is set by α and R_c.

3.2.2 Spread

As we have mentioned before, not only the mean slopes but also the spreads of the M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$| correlations are expected to reach the same limit value in the viscous scenario. Lodato et al. (2017) have shown that the |$\dot{M} \!-\! M_{\mathrm{d}}$| correlation is expected to tighten in time, leading to a decreasing spread; on the other hand, as the stellar masses hardly change, we expect the M_d−M_⋆ correlation to follow the same behaviour as that of the distribution of disc masses, which we discuss below.

If we assume the initial disc masses and radii to be distributed log-normally, we can show that viscous evolution preserves the lognormal shape of the distributions. Moreover, the spread of the M_d distribution at t ≫ t_ν writes (see Appendix  B for derivation)

$$\begin{eqnarray} \sigma ^2_M(t) = {\sigma ^2_M}(0) + (1 - \eta)^2 {\sigma ^2_{t_{\nu }}}(0); \end{eqnarray}$$

(15)

|$\sigma _{t_{\nu }}$| is the spread of the distribution of viscous time-scales, which in our case, given the direct proportionality between the viscous time-scale and radius, is very close to that of radii σ_R. Equation (15) shows that the dispersion in the distribution of disc masses increases with disc evolution.

Conversely, we expect the dispersion of the distribution of accretion rates to decrease with time. As the initial accretion rate can be written as |$\dot{M} \propto M_{\mathrm{d}, 0}/t_{\nu }$|⁠, the initial dispersion |$\sigma _{\dot{M}}(0)$| is given by |${\sigma ^2_{\dot{M}}}(0) = {\sigma ^2_{M}}(0) + {\sigma ^2_{t_{\nu }}}(0)$|⁠. At evolved times, instead, we have already noticed that the two distributions will coincide: this means that also |$\sigma _{\dot{M}}(t)$| will be given by equation (15). For γ = 1, (1 − η)² = 1/4, meaning that |$\sigma _{\dot{M}}(t) \gt \sigma _{\dot{M}}(0)$|⁠.

3.3 Summary: implications of the different scenarios

If we assume a power-law dependence of the disc mass M_d and accretion rate |$\dot{M}$| on the stellar mass M_⋆ as initial condition for a population of discs, viscous theory predicts four different evolutionary scenarios for the observed correlations at later times:

• young populations (t ≪ t_ν): evolution has not taken place yet, therefore the observed slopes λ_m and λ_acc are the same as the initial values λ_{m, 0} and λ_{acc, 0};

• old populations (t ≫ t_ν): evolution has taken place and both slopes have had the time to reach the evolved value, λ_m = λ_acc = (3λ_{m, 0} − λ_{acc, 0})/2. This can take place either via

  • both slopes already starting at the final value λ, if λ_{m, 0} = λ_{acc, 0} (implying δ₀ = λ_{m, 0} − λ_{acc, 0} = 0);

  • a steepening of both slopes towards λ_m = λ_acc, if λ_{m, 0} > λ_{acc, 0} (implying δ₀ > 0);

  • a flattening of both slopes towards λ_m = λ_acc, if λ_{m, 0} < λ_{acc, 0} (implying δ₀ < 0).

Table 3 summarizes all of the possible theoretical scenarios, with each line representing a set of parameters. The steps to determine a set of values are as follows:

• First, the initial values of the slopes of the correlations M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$|⁠, λ_{m, 0}, and λ_{acc, 0}, are chosen. This determines δ₀, defined as the difference between the two.

• β, the slope of the correlation H/R − M_⋆ (equation 6), is chosen.

• Using equation (10), ζ₀ – the slope of the correlation R_d − M_⋆ – is dletermined.

Table 3.

Summary of the different possible theoretical scenarios divided by the relative values of λ_{m, 0} and λ_{acc, 0} (see text for details, especially equation 10 and Section 3.3). The definition of each parameter appearing in the table is reminded in the top grey box. Each line represents a set of parameters. As discussed in the text, we dismiss the cases where β = 0 (which would imply an aspect ratio independent on M_⋆) and ζ₀ < 0 (implying disc radii decreasing with M_⋆): for visual purposes, cells corresponding to these two conditions have a grey background. The only line which does not include any grey cells is the one that is not influenced by any of these restraints, i.e. the most physically reasonable scenario, and is framed in green.

However, as we discussed above, not all parameters can assume any possible theoretical value while still maintaining physical relevance. In particular, β (the slope of the correlation of H/R with M_⋆) is unlikely to be zero, and ζ₀ (the initial slope of the correlation between R_c and M_⋆) is unlikely to be negative. These conditions are visualized in Table 3 as cells with a grey background. Note that the assumption of a constant α does influence this argument: the increasing, decreasing and constant behaviour of λ_m and λ_acc based on the initial conditions is not affected, but the physical likelihood of the different scenarios is. Some lines in Table 3 show the same value of β but different intervals for δ₀; this is because, based on the value of δ₀, different possibilities for ζ₀ may arise, in particular, it is meaningful to split the different possibilities if they include a change in the sign of ζ₀. Given that each line in the table represents a set of parameters, lines that contain one or more grey cells turn out to be discarded.

The only line that is not influenced by these restrictions, framed in green, represents the physically meaningful scenario. It is interesting to note that the requirement for both β and ζ₀ to assume reasonable values determines the difference between the initial slopes λ_{m, 0} and λ_{acc, 0}. Specifically, it prescribes δ₀ to be greater than 1/2: this implies λ_{m, 0} > λ_{acc, 0}, making the steepening slopes scenario the most reasonable one. Moreover, as we discussed in the previous section, δ₀ > 0 also implies a positive initial correlation between the accretion time-scale and the stellar mass. The observational claim by A17 is intriguingly matching our theoretical consideration; a preliminary comparison with data and a discussion of its implications is performed in Section 5.1.

4 POPULATION SYNTHESIS

In this section, we discuss the numerical counterpart of the theoretical analysis presented in Section 3. In particular, we describe the population synthesis model that we have implemented in Diskpop: we briefly present the physical framework, and we show the numerical results as well as discuss their agreement with the theoretical predictions.

4.1 Numerical methods – Diskpop

Developing a population synthesis model requires to implementing a numerical code to generate and evolve a synthetic population of protoplanetary discs. In this paper, we have used the python code Diskpop that we developed and that we will release soon. In this section, we briefly describe its basic functioning; for a more detailed presentation of it, we refer to the upcoming paper.

Diskpop generates an ensemble (which we term population) of N discs using the following scheme:

• Determines N stellar masses: this is achieved performing a random draw from an input probability distribution. In Diskpop, we use the initial mass function proposed by Kroupa (2001);

• Determines the N mean values of initial disc masses and radii: this is where we account for correlations between stellar and disc parameters. Following the initial correlations listed in equation (3), we choose the values of λ_{m, 0} and λ_{acc, 0} and evaluate ζ₀ using equation (10); the initial mean disc mass and radius are then computed, for each of the N stellar masses from step (i), using the prescriptions (7) and (8). Note that, by doing so, the correlations are intrinsically set to hold only for the mean value of the relevant parameters;

• Draws the N disc masses and radii: for this, the user needs to set two distributions (usually normal or lognormal) and their initial spread in dex (usually between 0.5 and 1.5 dex). After these parameters are determined, the draw is performed using the mean values computed in step (ii).

The other relevant parameters besides M_⋆, M_d, and R are fixed in our model. In particular, the aspect ratio of the disc at 1 au and the α parameter for the Shakura & Sunyaev (1973) prescription are set to be |$\frac{H}{R} \big |_{R = 1 \text{ au}} = 0.03$| for a 1 M_⊙ star and α = 10⁻³. As for the value of α, it is still very much debated how high, or low, it should be or even if it should be constant in time or across the discs in a population at all (Rafikov 2017); some works have also suggested a dependence of α on the radial position in the disc, namely a lower value in the inner disc with respect to the outer disc (Liu et al. 2018). None the less, assuming a constant value of around 10⁻³ leads to reproducing the observed evolution. A proper study of the α value is out of the scope of this paper.

Once the population is initialized, we want to evolve it using the chosen prescription. In the purely viscous case, which is the focus of this first paper, we can either use the analytical self-similar solution (2) or solve the viscous evolution equation

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} = \frac{3}{R} \frac{\partial }{\partial R} \left(R^{1/2} \frac{\partial }{\partial R} (\nu \Sigma R^{1/2}) \right). \end{eqnarray}$$

(16)

To do so, Diskpop employs the python code presented in Booth et al. (2017). In this work, we have used the numerical approach, but as long as an analytical solution is available, there is no difference in the two methods.

The raw numerical solution is an array of values for the gas surface densities Σ at the chosen ages, from which we compute the other quantities of interest: the disc mass is defined as the integral of the surface density from the inner (R_in) to the outer (R_out) radius,

$$\begin{eqnarray} M_{\mathrm{d}} = \int _{R_{in}}^{R_{out}} 2 \pi r \Sigma \, {\rm d}r, \end{eqnarray}$$

(17)

and R as the radius enclosing 68 |${{\,\rm per\,cent}}$| of the total disc mass.

4.2 Numerical simulations and comparison with theoretical expectations

In this section, we show the results obtained from simulating populations of N = 100 protoplanetary discs using Diskpop, evolved via numerical solution of equation (16). Following the scheme of Section 3, we divide the following discussion on the basis of the evolutionary regime considered. Our aim is to compare population synthesis results with the analytical prescription that we derived in the previous section.

4.2.1 Young populations: t ≪ t_ν

The discussion in Section 3.1 shows that, if the considered population is much younger than its viscous time-scale, we do not expect any significant evolution of the initial conditions. It is worth making a couple of considerations on the order of magnitude of the viscous time-scale. For γ = 1, t_ν scales with the disc radius (as well as α⁻¹, which is fixed in this work though), and given that our initial disc sizes are of the order of 10 au, the typical viscous time-scale will be shorter than 1 Myr; this is the case for more than |$98 {{\,\rm per\,cent}}$| of the discs in the synthetic population. This constraint stems from the need to reproduce the mean observed mass and accretion rate, as discussed in Section 3.1. Such a short viscous time-scale means that young populations would be very challenging to observe. The earliest data available for populations of discs usually correspond to ages around 1 Myr (e.g. Lupus and ρ Ophiucus); if the distribution of t_ν reproduced the viscous time-scales of these regions, already the youngest populations would be too old for t ≪ t_ν to hold. For this reason, in the following, we only show numerical results for the evolved populations scenario. As mentioned earlier, all of this discussion holds for the particular viscous time-scales that derive from the chosen parameters in our simulations: there is in principle the possibility of obtaining longer viscous time-scales (lower α, larger R_c), which would reverse this argument. However, from the observational point of view, this would require larger disc masses or lower accretion rates.

4.2.2 Evolved populations: t ≫ t_ν

In the case of evolved populations, which are observed at ages longer than their viscous time-scales, we present the numerical results for all of the three scenarios discussed in Section 3.3. Fig. 2 shows the mean fitted slopes of the M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$| correlations (λ_m and λ_acc, respectively) versus the age of the population. As summarized in Section 3.3, from the theoretical point of view, we expect both λ_m and λ_acc to reach the same evolved value, either through a steepening, a flattening, or a constant evolution. Each of these three possibilities corresponds to a choice of the initial parameters, and is determined by the sign of their difference: δ₀ = λ_{m, 0} − λ_{acc, 0}. In Fig. 2, we present our results, where the cyan and pink lines refer to λ_m and λ_acc, respectively, obtained in each of these scenarios:

• the left-hand panel corresponds to λ_{m, 0} = 1.3, λ_{acc, 0} = 1.9⇒δ₀ < 0, which is expected to produce a flattening of the slopes towards the evolved value;

• the middle panel shows the evolution if λ_{m, 0} = λ_{acc, 0} = 1.7⇒δ₀ = 0, which corresponds to the theoretical expectation of constant slopes;

• finally, in the right-hand panel we set λ_{m, 0} = 2.1, λ_{acc, 0} = 1.5⇒δ₀ > 0, which should lead to an increasing trend of the slopes.

Figure 2.

Time evolution of the mean values of the slopes of the M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$| correlations, λ_m and λ_acc (cyan and pink, respectively), as obtained running Diskpop and fitting the results using linmix. The black dashed line represents the limit value of both slopes for t → +∞ (see equation 14). Each panel corresponds to a different choice in the initial values λ_{m, 0} and λ_{acc, 0}. Left-hand panel: λ_{m, 0} = 1.3, λ_{acc, 0} = 1.9, δ₀ = λ_{m, 0} − λ_{acc, 0} < 0; middle panel: λ_{m, 0} = λ_{acc, 0} = 1.7, δ₀ = 0; right-hand panel: λ_{m, 0} = 2.1, λ_{acc, 0} = 1.5, δ₀ > 0. The numerical behaviour follows the theoretical prediction derived in Section 3: λ_{m, 0} < λ_{acc, 0} (left-hand panel) leads to a flattening trend, λ_{m, 0} = λ_{acc, 0} (middle panel) to a constant slope, and λ_{m, 0} > λ_{acc, 0} to a steepening trend. Moreover, the theoretical limit value is recovered for every choice of parameters. In every simulation, H/R|_{R = 1au} = 0.33 for M_⋆ = 1M_⊙, α = 10⁻³, β = −0.5, and both initial conditions M_d(0) and R_c(0) follow a lognormal distribution.

All simulations have the same values of α = 10⁻³, H/R|_{R = 1au} = 0.03 for M_⋆ = 1M_⊙, β = −1/2. The black dashed line shows the theoretical evolved value, computed based on the initial values of the slopes as per equation (14). We can see that all of the three simulations do match our analytical evolutionary predictions: the limit value for evolved populations is recovered, and the flattening, constant, and increasing trends are recovered. The middle panel shows a small evolution at very early ages, but it is most likely due to numerical effects and can be neglected. Note that the synthetic populations shown in the plots are evolved up to 20 Myr; while we do not expect to observe disc-bearing protostars at those very late ages (due to the disc removal processes at play that this work does not include – see Mamajek 2009), we evolved our synthetic populations up to such a long age to make sure no transient affected our results. Moreover, the evolved value is reached at very late ages (∼10 Myr), when we expect processes other than viscous evolution to have taken place.

4.3 Numerical simulations introducing a spread

The real power of disc population syntheses is the possibility of introducing a spread in the initial conditions of the population: in this section, we explore the influence of such a spread on our results. Fig. 3 shows the output of Diskpop: the left and right-hand panels display the correlation M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$|⁠, respectively, at four subsequent timesteps as per the legend. Each dot represents a disc in the population; the fit was performed using the python package linmix.

Figure 3.

M_d−M_⋆ (left-hand panel) and |$\dot{M} \!-\! M_{\star }$| (right-hand panel) correlations obtained with Diskpop at four consequent time steps from top to bottom, as shown in the legend. The dots represent discs in the population, and the numerical fit is overplotted. These simulations were performed with the same parameters as the right-hand panel of Fig. 2 (δ₀ > 0), with the addition of a spread in the initial conditions of σ_{M, 0} = 0.65 dex, and σ_{R, 0} = 0.52 dex.

We also analysed the time evolution in the spreads of both correlations. The spread is determined as the standard deviation of the vertical distances of every point from the fitted relation. Fig. 4 shows the results with input values σ_M(0) = 0.65 dex, σ_R(0) = 0.52 dex (chosen to be consistent with the observed values in the data sets from A17 and T22, referring to lognormal distributions). The light blue and pink lines show the evolution of the dispersion of the M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$| correlations respectively. As we pointed out in Section 3.2.2, we expect the spread of the two correlations to reach the same evolved value, represented in the figure by the black dashed line. The numerical limit of 0.71 dex is very close to the theoretical one of 0.7; the slight difference is most likely due to the fact that the numerical dispersion in t_ν also accounts for the dispersion in stellar masses, which is not considered in the analytical calculation.

Figure 4.

Time evolution of the spread of the M_d−M_⋆ (light blue) and |$\dot{M} \!-\! M_{\star }$| (pink) correlations. The black dashed line shows the analytical spread of the M_d distribution for t ≫ t_ν as per equation (15). Both spreads reach the same value after some million years of evolution, in agreement with M_d and |$\dot{M}$| reaching a unity correlation in the logarithmic plane; this final spread slightly differs from the analytical estimation for evolved populations due to the addition of a spread in the stellar masses.

5 DISCUSSION

In Section 3, we have performed a theoretical analysis of the viscous evolution of power-law initial conditions. Through it, we have determined three possible physical scenarios which differ by the relative values of the initial slopes λ_{m, 0} and λ_{acc, 0}. We then analysed the output of our population synthesis code Diskpop and compared the evolutionary behaviour with our theoretical expectations (Section 4). In this section, we perform a preliminary comparison with some of the available observational data (Section 5.1) and discuss the implications of our findings (Section 5.2).

5.1 Preliminary comparison to observational data

After we tested Diskpop, recovering the analytical prediction for the viscous evolution of λ_m and λ_acc, we performed a preliminary comparison using three sets of observational data (M12; A17; T22)

There are some relevant caveats to these comparisons. First of all, we did not include any observational effects, biases, or sensitivity limitations in our simulations. However, the goal of this comparison is not to fit any data set, nor to perfectly reproduce the observations: our aim is just to have a qualitative idea of the parameter space, and to see whether the general trend shown by observations can (at least in principle) be reproduced from the theoretical point of view.

Figs 5–7 show the data comparison for the M_d−M_⋆ (Fig. 5) and the |$\dot{M} \!-\! M_{\star }$| (Figs 6 and 7) correlation. The top panel of Fig. 5 refers to data from A17, the bottom panel of Figs 5 and 6 from T22, and Fig. 7 from M12. Each panel in these figures contains two sets of simulations, represented by the solid and dashed lines respectively, corresponding to two different sets of parameters (λ_{m, 0} and λ_{acc, 0}, see legend), with a spread of σ_M(0) = 0.65 dex and σ_R(0) = 0.52 dex. The darker shaded areas around the solid and dashed lines represent the 25th and 75th percentiles for the fitted slopes, obtained from 100 different statistical realisations of the same simulations. The initial conditions on λ_m and λ_acc were chosen to qualitatively match the observed slopes at later and younger ages respectively; the other parameters in the simulations are fixed; in particular, H/R|_{R = 1au} = 0.33 for M_⋆ = 1 M_⊙, α = 10⁻³, and β = −0.5. The mean observational values are represented by diamonds, and the vertical and horizontal lines show their uncertainties. As we already discussed in Section 2, there is a key difference between the data by M12 and all of the other observations: M12 focused on a single star-forming region, dividing the protostars based on their isochronal age group, while data from both A17 and T22 instead refer to different star-forming regions, each with its own mean age. For this reason, the horizontal bars in Fig. 7 do not refer to the errors in age, but rather to the extent of the age intervals considered when performing the fits. Because of these considerable differences from the observational point of view, we divided Figs 6 and 7.

Figure 5.

Comparison of the numerical results for λ_m obtained from Diskpop with the observational data from A17 (top panel) and T22 (bottom panel). The solid and dashed lines in each panel represent the results of simulations performed with two different sets of parameters, as shown in the legend. The diamonds correspond to the observational mean values, while the vertical and horizontal lines represent the 1σ error bars. The two sets of parameters displayed are able to match (from a qualitative point of view) the youngest and oldest populations, respectively; they represent the limits of the parameter space. Both are chosen to fall in the most physically meaningful scenario, with δ > 1/2. The darker shaded regions around the solid and dashed lines represent the 25th and 75th percentile for the fitted slopes, obtained from 100 different statistical realizations of the same simulation.

Figure 6.

Same as Fig. 5, but referring to λ_acc, showing both the numerical results from Diskpop and the observational data from T22. See the caption of Fig. 5 for detail on annotations.

Figure 7.

Same as Fig. 6, but compared to the observational data from M12. See the caption of Fig. 5 for detail on annotations. Note that the horizontal bars in this figure do not represent the errors on the ages, but rather the extent of the age interval that was considered in performing the fit (see text for details).

It is worth recalling that, out of the two diagnostics M_d and |$\dot{M}$|⁠, the second is more readily translated in the observational space. On the other hand, disc mass measurements always refer to dust content; given that our simulations only include the gaseous component, all of the observed masses shown here have been multiplied by the dust-to-gas ratio (set to the usual value of 100, see Bohlin, Savage & Drake 1978). If Diskpop included dust evolution, the numerical disc mass M_d would be affected and so would the M_d − M_⋆ correlation and its time evolution; on the other hand, we do not expect a significant influence of the presence of dust on the accretion rate. At this stage, we also assumed spread-less initial conditions for the disc masses and radii. For these reasons, as well as not including any observational effects or biases, the comparison shown in Figs 5–7 is qualitative.

Nonetheless, our comparison does suggest some interesting ideas. A17 claimed to have found an increasing trend of λ_m with the age of the sample, which can also apply to T22 – even if it is not as evident. Furthermore, in T22 the mean value for λ_acc seems to be increasing as well: the visible exception of the last point (Upper Scorpius) must be treated carefully, as it represents a highly incomplete sample – hence the large error bars (see the original paper for more detail). Due to the method adopted to determine the age of the subgroups in the Orion Nebula Cluster, M12 found significant errors on the fitted slopes; despite not being statistically significant, the mean measured values are increasing, and these data show a general agreement with the steepening trend. The interesting message of Figs 5–7 is that the observed increasing trend can be reproduced from the theoretical point of view, as the simulations show. Despite including a spread of around 0.5 dex in both the disc mass and radius, the analytical expectations are well recovered: given the choices of parameters (in particular, δ₀ > 1/2), our simulations show the steepening evolution of both λ_m and λ_acc. The main result is that, if this is the case, such observed correlations would be naturally explained in the purely viscous framework by

• assuming power-law correlations as initial conditions;

• imposing a positive correlation between the viscous time-scale and the stellar mass.

Figs 5–7 also show the effect of a spread in the initial conditions on the evolution of the slopes. At t = 0, the median slopes (solid and dashed lines for the two initial conditions) and the input values λ_{m, 0} and λ_{acc, 0} coincide with a precision of ∼10⁻³ and ∼10⁻², respectively; moreover, the evolutionary path of both λ_m and λ_acc is not changed. Therefore, including a spread in the initial condition only influences the starting point of the evolution and does not affect the fundamental results that we have discussed.

As we mentioned above, we chose the parameters in the two different runs of Diskpop to qualitatively match the youngest and most evolved values for both λ_{m, obs} and λ_{acc, obs}: this identifies a region of the parameter space for λ_{m, 0} and λ_{acc, 0}, which lead to an evolution that roughly matches the observations (within their 1σ error bars). These two sets of parameters are given by λ_{m, 0} = 2.1, λ_{acc, 0} = 1.5 and λ_{m, 0} = 1.2, λ_{acc, 0} = 0.7; the interval in the parameter space is

$$\begin{eqnarray} \lambda _{\mathrm{m}, 0} \in [1.2, 2.1] \quad \lambda _{\mathrm{acc}, 0} \in [0.7, 1.5]. \end{eqnarray}$$

(18)

Figs 5–7 unveil another feature worth pointing out. Different values of λ_{m, 0} and λ_{acc, 0} can change the starting point of the curves obtained, but cannot change their shape, i.e. the rate of steepening or flattening of the slopes is not influenced by the initial parameters in the purely viscous case. We will discuss this further in the next section.

5.2 Evolution of the slopes

As we discussed above, the key promising aspect worth pointing out is that there is indeed a theoretical scenario which is able to reproduce the observed steepening of the correlations. Furthermore, this case also corresponds to the most likely physical regime, which leads to a positive correlation between the disc radius and the stellar mass. If these analytical requirements are matched, our analysis finds that power-law initial conditions for the correlations λ_m and λ_acc can only be steepened by viscous evolution. This implies that the observed correlations could be easily explained in terms of initial conditions and viscous evolution, without invoking any additional physical process.

On the other hand, we also found that changing the initial values λ_{m, 0} and λ_{acc, 0} raises or lowers the evolution curve, adding only a small modification to the rate at which the exponents of the correlation steepen. Indeed, even if the difference δ₀ = λ_{m, 0} − λ_{acc, 0} does impact the viscous time-scales, with these choices of parameters they still remain of the order of 1 Myr; as the bulk of the evolution in our models happens on comparable time-scales, the effect of these changes is minimal.

Initial conditions within the parameters space specified in equation (18), and with δ₀ > 1/2, are able to qualitatively match the observed slopes within 1σ (although note that, as the error bars are large, this is not strongly constraining on any of the three models). However, the youngest and oldest populations seem to be better described by different simulations. This can be explained in two different ways:

• Initial conditions are not the same for every star-forming region. This would mean that different curves must be considered for different regions to match their evolution;

• Viscous evolution alone is not enough to explain the whole extent of star-forming regions.

The second possibility is very much reasonable; a number of physical effects are now thought to affect disc evolution, sooner or later in their lifetime. We can mention in particular MHD winds, which are expected to play a role alongside viscosity (Tabone et al. 2022a), internal and external photoevaporation (Clarke, Gendrin & Sotomayor 2001; Alexander, Clarke & Pringle 2004; Owen, Clarke & Ercolano 2012), or dust evolution (Sellek et al. 2020a); these effects are complicated to include in the mathematical analysis, but may be implemented in numerical codes to test whether they do lead to any modifications to the steepening rates of λ_m and λ_acc. Dust evolution in particular would likely affect the computed disc mass, although it should not have a considerable influence on the accretion rate. This expected behaviour is encouraging, as Fig. 5 shows that the shape of the time evolution of λ_acc (right-hand panel) qualitatively matches the observed data way better than that of λ_m (left and central panels). An ideal candidate to show the signature of more complex evolutionary patterns would be Upper Scorpius, which is notorious for being difficult to explain through classic viscous models (Manara et al. 2020; Trapman et al. 2020) unless adding additional physics, such as dust evolution (Sellek et al. 2020a).

6 CONCLUSIONS

In this paper, we investigated the correlation between the stellar mass and disc properties, M_d−M_⋆ and |$\dot{M} \!-\! M_{\star }$|⁠. Assuming power-law correlations as initial conditions (⁠|$M_{\mathrm{d}}(0) \propto {M_{\star }}^{\lambda _{\mathrm{m}, 0}}$|⁠, |$\dot{M}(0) \propto {M_{\star }}^{\lambda _{\mathrm{acc}, 0}}$|⁠), as well as viscous evolution, we obtained analytical equations to describe how the exponents of these correlations should evolve with time. Our main findings are the following:

• In the viscous picture, the two correlations follow the same trend: they both either grow steeper or shallower with time. Given enough time, they tend to the same value, which is determined by the initial slopes or, in other terms, by the scaling of the viscous time-scale with the stellar mass. If we further make the plausible assumption that H/R ∝ M_⋆^β, with β = −1/2 (as supported by radiative transfer models), α does not depend on M_⋆ and δ₀ = λ_{m, 0} − λ_{acc, 0} > 1/2, both correlations steepen with time, as tentatively suggested by observations. This is a consequence of viscous time increasing with stellar mass. Moreover, the spread of the two correlations is also bound to evolve to reach the same value, determined by the initial spreads in disc masses and viscous time-scales.

• Motivated by this early success, we attempted a comparison of our predictions with the available data. To this end, we built a numerical tool, Diskpop, to evolve synthetic disc populations and include effects that cannot be tackled in the analytical approach (such as a spread in the initial conditions and the effect of random sampling).

• With the two sets of initial conditions λ_{m, 0} = 1.2, λ_{acc, 0} = 0.7 and λ_{m, 0} = 2.1, λ_{acc, 0} = 1.5, we obtain two limit curves which match either the youngest or the oldest region in the samples (within their error bars). This means that the range λ_{m, 0} ∈ [1.2, 2.1], λ_{acc, 0} ∈ [0.7, 1.5] as exponents in the initial conditions allows to span the observed parameter space.

• Changing the initial condition can either raise or lower the curve, without significantly modifying its shape: This implies that our model cannot exactly match both the youngest and the oldest populations. This behaviour can imply two consequences: either the initial condition is not the same for every star-forming region, or viscosity alone is not enough to explain the observed long-term evolution [as already discussed in previous works, e.g. Mulders et al. (2017), Manara et al. (2020), Testi et al. (2022)]. Of course, these two possibilities can also be valid at the same time. In particular, we expect dust evolution and planet formation, followed by disc-planet interaction, to play a crucial role: affecting the disc mass M_d, it will likely change the time evolution of λ_m. However, it would probably not affect λ_acc, which intriguingly already resembles the observations better.

In our future work, we aim at further developing Diskpop, increasing its complexity adding relevant physical effects such as MHD winds, photoevaporation and dust evolution; this will allow us to test the influence of these phenomena on the evolution of populations, and to obtain numerical results more suitable for a comparison with observational data.

ACKNOWLEDGEMENTS

We thank an anonymous referee for their valuable comments that helped us improve the clarity of this paper. We thank Richard Alexander, Barbara Ercolano, Til Birnstiel, and Kees Dullemond for the 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), by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Ref no. 325594231 FOR 2634/1 TE 1024/1-1, by the DFG cluster of excellence Origins (www.origins-cluster.de), and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). GR acknowledges support from an STFC Ernest Rutherford Fellowship (grant number ST/T003855/1). MT has been supported by the UK Science and Technology Research Council (STFC) via the consolidated grant ST/S000623/1. GL, GR, LT, CFM, and CT have been supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 823823 (RISE DUSTBUSTERS project).

DATA AVAILABILITY

All the data sets generated and analysed during this study are available from the corresponding author on reasonable request.

Footnotes

1

Our notation slightly differs from the ones previously used by Pascucci et al. (2016) and Ansdell et al. (2017): In the first paper, the slope (here λ_{m, obs}) is α and the intercept (here q_obs) is β, while the second one uses the opposite convention. Moreover, we named the intrinsic dispersion Δ_obs instead of δ to avoid confusion with another parameter that we define later in our paper.

2

Using the standard convention for linear relations, y = α + βx.

APPENDIX A: R_d − M_⋆ CORRELATION

Fig. A1 shows a plot of the gas (CO) disc radius as a function of the stellar mass using data by Sanchis et al. (2021), fitted using linmix. The best-fitting parameters² are α = 2.14 ± 0.08, β = 0.32 ± 0.13, with a spread σ = 0.32 ± 0.16 dex. As discussed in Section 2.3, the correlation is positive but the fit is not particularly strong; the correlation coefficient is ρ_XY = 0.32, which implies weak correlation.

Figure A1.

Radius obtained from CO observations [data by Sanchis et al. (2021)] as a function of the stellar mass in the Lupus star-forming region.

APPENDIX B: EVOLUTION OF THE DISC MASS DISTRIBUTION IN THE SELF-SIMILAR SCENARIO

Here, we demonstrate that an initially lognormal distribution of disc masses keeps being lognormal at long times under the assumption of self-similar viscous evolution and estimate the relevant distribution parameters.

We start by defining some useful quantities: let m₀ = log M_d,0 and τ_ν = log t_ν be the log of the initial disc mass and viscous time of a specific disc in a population. Now, the initial distributions of disc masses and viscous times are assumed to be lognormal, so that:

$$\begin{eqnarray} \frac{\partial N}{\partial m_0}=N_1\exp [-(m_0-\bar{m}_0)^2/2\sigma _M^2], \end{eqnarray}$$

(B1)

$$\begin{eqnarray} \frac{\partial N}{\partial \tau _\nu }=N_2\exp [-(\tau _\nu -\bar{\tau }_\nu)^2/2\sigma _{t_\nu }^2], \end{eqnarray}$$

(B2)

where |$\bar{m}_0$| and |$\bar{\tau }_\nu$| are the average values of the initial disc mass and viscous time (in logarithmic sense), σ_M and |$\sigma _{t_\nu }$| are their dispersions, and N₁ and N₂ are normalization constants. Note that the two lognormal probability distributions are independent, in the m₀ − τ_ν space, and their means and variances may depend on the stellar mass. For a self-similar disc in the limit t ≫ t_ν, the disc mass at time t is M_d(t) = M_{d, 0}(t/t_ν)^a, where a = 1 − η (see equation 4). In logarithmic form, this is:

$$\begin{eqnarray} m=m_0+a\tau -a\tau _\nu , \end{eqnarray}$$

(B3)

where m = log M_d(t) and τ = log t. We can invert equation (B3), obtaining

$$\begin{eqnarray} \tau _{\nu ,m}=\frac{m_0}{a}-\frac{m}{a}+\tau , \end{eqnarray}$$

(B4)

so that τ_{ν, m} is the viscous time for which a disc with initial mass m₀ evolves into m after a time τ.

The distribution of disc masses at time t can then be obtained by integrating over the distribution of initial disc masses and viscous times, under the condition that equation (B4) is satisfied

$$\begin{eqnarray} \frac{\partial N}{\partial m}=\int \int \frac{\partial N}{\partial m_0}\frac{\partial N}{\partial \tau _\nu }\delta (\tau _\nu -\tau _{\nu ,m}){\rm d}\tau _\nu {\rm d} m_0, \end{eqnarray}$$

(B5)

where δ is the Dirac function. The integral over viscous times can be readily be done, and, after inserting the initial distributions of disc masses and viscous times and using equation (B4), we obtain

$$\begin{eqnarray} \frac{\partial N}{\partial m} &=& N_1 N_2 \int \exp \Bigg\lbrace -\frac{(m_0-\bar{m}_0)^2}{2\sigma _M^2}\\&&-\,\left[ \frac{(m_0-m)}{a}+(\tau -\bar{\tau }_\nu)\right]^2\frac{1}{2\sigma _{t_\nu }^2} \Bigg\rbrace {\rm d} m_0. \end{eqnarray}$$

(B6)

Let us define |$\bar{m}=\bar{m}_0+a\tau -a\bar{\tau }_\nu$|⁠, which we will see is the average disc mass at time τ. We also define |$\tilde{m}=m-\bar{m}$| and |$\tilde{m}_0=m_0-\bar{m}_0$|⁠. The exponent in equation (B6) can be then rewritten as (neglecting the factor −1/2)

$$\begin{eqnarray} \mathcal {M}=\frac{\tilde{m}_0^2}{\sigma _M^2}+\left[\frac{\tilde{m}_0}{a}-\frac{\tilde{m}}{a}\right]^2\frac{1}{\sigma _{t_\nu }^2}. \end{eqnarray}$$

(B7)

We can now easily rewrite the exponent |$\mathcal {M}$| in such a way to isolate the dependence on m₀ (that we integrate upon). After some straightforward algebra, we get

$$\begin{eqnarray} \mathcal {M}=\left(A\tilde{m}_0-\frac{B^2}{A}\tilde{m}\right)^2 + \frac{A^2-B^2}{A^2/B^2}\tilde{m}^2, \end{eqnarray}$$

(B8)

where

$$\begin{eqnarray} B^2=\frac{1}{a^2\sigma ^2_{t_\nu }} \end{eqnarray}$$

(B9)

and

$$\begin{eqnarray} A^2=B^2+\frac{1}{\sigma ^2_M}. \end{eqnarray}$$

(B10)

Now, integrating equation (B6) over |$\tilde{m}_0$|⁠, the whole first term on the RHS in equation (B8) translates into a normalization constant, leaving us with just a Gaussian distribution for |$\tilde{m}$|⁠. This means that the disc masses at time t are distributed log-normally, with mean |$\bar{m}=\bar{m}_0+a\tau -a\bar{\tau }_\nu$| and with dispersion

$$\begin{eqnarray} \sigma _M(t)=\frac{A^2/B^2}{A^2-B^2}=\sigma _M^2+a^2\sigma _{t_\nu }^2. \end{eqnarray}$$

(B11)