© The Authors 2024

1 Introduction

Proto-planetary disks are the birth environments of planets. In recent decades, significant progress has been made in the field of submillimeter radio astronomy, in particular thanks to the development of the Atacama Large Millimeter/Sub-millimeter Array (ALMA) observatory. Thanks to a transformational improvement in the capabilities of the new generation of telescopes, we now know that disks are substructured, showing gaps, spirals, kinks, and many more features (e.g., Andrews 2020).

In addition to high resolution observations, the study of disk demographics has also advanced significantly. In recent years, it has become possible to collect global disk properties such as submillimeter fluxes (e.g., Mann et al. 2014; Ansdell et al. 2016, 2017; Barenfeld et al. 2016; Pascucci et al. 2016; Cox et al. 2017; Eisner et al. 2018; Cazzoletti et al. 2019; Cieza et al. 2019; Ansdell 2020), disk radii (e.g., Barenfeld et al. 2017; Ansdell et al. 2018; Sanchis et al. 2021), and rates of mass accretion onto the central stars (e.g., Manara et al. 2015; Manara et al. 2017; Alcalá et al. 2017; Manara et al. 2020). This astonishing amount of new data allowed us to start investigating the statistical distributions of disk properties, which make it possible to constrain evolutionary models (e.g., Manara et al. 2023 for a review). However, despite all the progress that has been made, we are still far from a paradigm that can completely describe the formation and evolution of disks (e.g., Morbidelli & Raymond 2016).

Understanding how protoplanetary disks accrete material onto the central star and evolve is of cardinal importance for developing a standard model of planet formation and evolution. Two main physical scenarios have been proposed to explain accretion onto the central star: viscous accretion (e.g., Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), which redistributes the angular momentum within the disk (e.g., Pringle 1981; Frank et al. 2002), and magnetohydrodynamic (MHD) wind-driven accretion, where a vertical magnetic field launches a wind that extracts angular momentum from the disk (e.g., Blandford & Payne 1982; Ferreira 1997; Lesur 2021). Despite all the work put in developing these paradigms, there is still significant uncertainty regarding which scenario best describes the evolution of proto-planetary disks, as their key predictions (e.g., the different time evolution of disk radii) are difficult to observe on a statistically significant sample of observations.

In this context, observations have allowed us to find a correlation between the mass of the proto-planetary disk and the mass accretion rate onto the central star (e.g., Manara et al. 2016; hereafter called the M_D − Ṁ_* correlation). In the viscous paradigm, the M_D − Ṁ_* correlation naturally arises with a slope close to unity, in agreement with observations (e.g., Lodato et al. 2017; Rosotti et al. 2017); on the other hand, it is still not clear which conditions can produce it in the MHD scenario. In particular, it is uncertain whether it would be a general feature of wind-driven populations or if specific initial conditions would be needed.

Recently, 1D disk evolution models have been proposed to describe the effect of MHD disk-winds on the long-term evolution of disks (Suzuki et al. 2016; Bai 2016; Tabone et al. 2022a). In particular, Tabone et al. (2022a) developed analytical solutions of global disk quantities (i.e., the surface density, the characteristic disk radius, the mass of the disk, and the mass accretion rate) for the MHD wind-driven scenario, and offered some insights about retrieving the correlation. In Tabone et al. (2022b), they show that for a particular set of initial parameters, the MHD model can reproduce the observed correlation in the Lupus region, the spread around this trend, and the decline of the disk fraction. In particular, they qualitatively show that a linear relationship can be obtained assuming no correlation between disk mass and accretion timescale t_acc,0. Starting from the solutions of Tabone et al. (2022a), in this work we conduct an extensive analysis to discuss the universality of this result and whether it implies fine tuning.

The scope of the present work is to go beyond the work of Tabone et al. (2022a) to quantify under which initial conditions it is possible to reproduce the observed M_D–Ṁ_* correlation. In Sect. 2, we show the condition under which it is possible to observe the M_D−Ṁ_* correlation under the assumption that the initial accretion timescale t_acc,0 ¹ is not correlated with M₀. We then analytically derive the slope of the M_D−Ṁ_*, correlation and its spread as a function of the spread of the disk age t and t_acc,0. In Sect. 3, we relax the assumption of no correlation between M₀ and t_acc,0: we introduce a correlation in the form $M_0\propto t_{\text{acc},0}^{\varphi }$, and we analytically derive the slope of the M_D−Ṁ_* correlation in the initial lifetime conditions of the disks and at late times. Then, by adding a spread in M₀ , we show that we can retrieve the observed value of the slope under a broad range of conditions. In Sect. 4, we show that the spread derived in Sect. 2 agrees with the observed spread of log(M_D /Ṁ_*,). Finally, we present and investigate the $\varphi -{\sigma }_{M_0}$ degeneracy. In Sect. 5, we summarize our conclusions.

2 The M_D-Ṁ_*, correlation in the absence of a M₀ − t_acc,0 correlation

Tabone et al. (2022a) find that a proto-planetary disk evolving under the effect of MHD wind torque is governed by the master equation $$\begin{array}{cc}\frac{\partial \text{Σ}}{\partial t}=& \frac{3}{r}\frac{\partial }{\partial r}\left\{\frac{1}{r\text{Ω}}\frac{\partial }{\partial r}\left(r^2{\alpha }_{SS}\text{Σ}{c}_s^2\right)\right\}+\frac{3}{2r}\frac{\partial }{\partial r}\left\{\frac{{\alpha }_{\text{DW}}\text{Σ}{c}_s^2}{\text{Ω}}\right\}\\ & -\frac{3{\alpha }_{\text{DW}}\text{Σ}{c}_s^2}{4(\lambda -1)r^2\text{Ω}},\end{array}$$(1)

where α_SS is the Shakura-Sunyaev α -parameter, while α_DW is its analog for a MHD wind, and λ is the magnetic lever arm. In particular, α_DW describes the amount of angular momentum extracted by the wind, while λ determines the mass-loss rate of the wind.

When the viscous torque is neglected (i.e., when MHD winds alone describe the disk evolution and dispersal), assuming α_DW and λ to be constant across the disk, and λ to be constant in time, it is possible to solve analytically the master equation, obtaining the surface density profile $$\text{Σ}(r,t)={\text{Σc}}_{}(t){\left(r/r_{\text{c}}\right)}^{-1+\xi }{\text{e}}^{-r/r_{\text{c}}},$$(2)

where ξ = 1/[2(λ − 1)], and r_c is the disk characteristic radius. Tabone et al. (2022a) present solutions where α_DW increases as the disk mass decreases, namely as $${\alpha }_{DW}(t)={\alpha }_{DW}(0){\left({\text{Σ}}_c(t)/{\text{Σ}}_c(0)\right)}^{-\omega },$$(3)

so that α_DW accelerates the evolution of the disk. Here ω is a phenomenological parameter between 0 and 1 that describes the unknown dissipation of the magnetic field. In particular, ω describes the tendency of the magnetization to increase as the disk mass decreases. For example, the case ω = 1 mimics a disk where the magnetic field strength is constant over time, while the case ω = 0 represents a disk with constant magnetization, so that it is dispersed on the same timescale as the mass dissipation. In this set of solutions the disk mass evolves as $$M_D(t)=M_0{\left(1-\frac{\omega t}{2t_{\text{acc,}0}}\right)}^{1/\omega },$$(4)

and the mass accretion rate Ṁ_*, evolves as $${\dot{M}}_{\ast }(t)=\frac{M_0}{2t_{\text{acc,}0}\left(1+f_{M,0}\right)}{\left(1-\frac{\omega t}{2t_{\text{acc,}0}}\right)}^{1/\omega -1},$$(5)

where f_M,0 is the initial mass ejection-to-accretion ratio $$f_{M,0}={\left(r_{\text{c}}/r_{in}\right)}^{\xi }-1.$$(6)

In the end, M_D (t) and Ṁ_*(t) are controlled by the four independent parameters M₀ , f_M,0, t_acc,0 and ω. In this work, if not otherwise specified we set f_M,0 = 2.2 and ω = 0.8.

In this paper, we test the emergence of the M_D−Ṁ_* correlation. In principle, the accretion onto the central star Ṁ_* happens in the innermost region of the disk, where this latter is expected to be turbulent (e.g., Najita et al. 1996, 2009; Carr et al. 2004; Hartmann et al. 2004; Ilee et al. 2014), while in the outer region we assume that the disk is non-turbulent. However, from the modeling point of view, a small turbulent inner disk will always simply re-adjust to the mass accretion rate to which it is fed by the outer wind-driven disk, since locally the viscous timescale is smaller than the disk age. In particular, the accretion rate onto the central star will always quickly adjust to the accretion rate at the inner edge of the non-turbulent region. Hence, we do not explicitly model an inner turbulent region.

Fig. 1 Isochrones of a population of disks with t = 5 Myr for three values of M0 in the MD − Ṁ* plane. The colored dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4), (5) with ω = 0.8, where tacc,0 follows a lognormal distribution centered on the natural logarithm of 1 Myr with a spread of 0.26 dex, while M0 = 0.1, 1.0, 10 M⊙.

2.1 The emergence of the M_D − Ṁ_* correlation

When studying populations of disks of similar age, it is of great interest to consider the concept of isochrones (introduced for the first time in Lodato et al. 2017), defined as the location in the M_D − Ṁ_* plane of the sample of disks that have the same initial mass M₀ , same age t, and different accretion timescales t_acc,0 .

In Fig. 1, we show the shape of the isochrones of a population of disks. Every dot is a proto-planetary disk, for which the value of M_D and Ṁ_* are evaluated using Eqs. (4), (5): t_acc,0 follows a lognormal distribution centered² on the natural logarithm of 1 Myr with a spread of 0.26 dex, while M₀ = 0.1, 1.0, 10 M_⊙. We recover the analytical shape of the isochrone provided in Tabone et al. (2022a). As shown, these curves cannot intrinsically represent a disk population that follows the observed M_D − Ṁ_* correlation: the boomerang-shaped isochrones cannot be represented by a power law in the form ${\dot{M}}_{*}\propto M_D^{1-\omega }$. These peculiar curves show that disks with high mass and low accretion rates are the ones with the longest t_acc,0 , while disks with the lowest mass and low accretion rates are the ones with shortest t_acc,0.

To visualize the M_D − Ṁ_* correlation we set M₀ to follow a lognormal distribution rather than taking it as constant. We show in Fig. 2, the results for different values of its spread. We find that the boomerang shape of isochrones, as defined above, is visible in the M_D − Ṁ_* plane only if the spread of the initial accretion timescale is large enough. Roughly, we can say it tends to be visible for ${\sigma }_{t_{\text{acc},0}}>{\sigma }_{M_0}$, where t_acc,0 is the initial accretion timescale of the disk, while in the opposite case the points follow a power law behavior.

We note that the observed M_D − Ṁ_* correlation (Manara et al. 2016) is a power law whose slope is close to unity. In order to sis- tematically determine the values of ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$ that lead to a correlation between M_D and Ṁ_*, we use the coefficient of determination R². Combining a visual inspection of the M_D − Ṁ_* plane with the evaluation of R² from the best power law fit for each disk population, we derived a conservative “rule of thumb”: R² ≳ 0.5. If R² ≳ 0.5, the disk population can be safely fitted in the M_D − Ṁ_* with a single power law. We note that this general rule is only informative of the goodness of fit for a disk population in the M_D − Ṁ_* plane with a single power law, and that we investigate under which conditions we can derive a linear correlation between M_D and Ṁ_* (see below). We also note that our criterion is in agreement with the correlation coefficient found in Manara et al. (2016) of r = (0.56 ± 0.12), which changes to r = (0.7 ± 0.1) when the upper limits are considered in the analysis.

In Fig. 3 we show the value of R² (the coefficient of determination) measured on a large grid of models. In these plots we set M₀ and t_acc,0 to follow a lognormal distribution centered on the natural logarithm of 0.1 M_⊙ and 1 Myr, respectively, and we chose ω = 0.8. The hatched contours represent the regions where R² = 0, and the dashed line represents the points where ${\sigma }_{M_0}={\sigma }_{t_{\text{acc},0}}$. Figure 3 shows the regions of the ${\sigma }_{M_0}={\sigma }_{t_{\text{acc},0}}$ plane that allow alow or high value of R². Starting from these plots, we can extract different combinations of values for ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$, and we can visualize which of them, at different ages and for different values of ω, yield a power law correlation in the M_D − Ṁ_* plane. We chose to span from very young (t = 10⁵ yr) to evolved disk populations (t = 2.5 × 10⁶ yr) to fully explore the parameter space.

Tabone et al. (2022b) show that there is no requirement to fine-tune ω to reproduce the observed accretion rates and disk masses, although explaining the rapidity of disk dispersal disfavors values of ω ~ 0; in particular, if 0.5 ≲ ω < 1 these latter experience a steeper drop before being dispersed and in this way reproduce the observed distribution. Hence, we chose ω = 0.8, which is bracketed by the ω = 0.5 and ω = 1 solutions shown by Tabone et al. (2022b). For completeness, in Appendix A, we show the same plots as in Fig. 3, but for the extreme cases of ω = 0.2 and ω = 1, highlighting how the results of this work also hold for different values of ω.

From Fig. 3, we note that our criterion R² ≳ 0.5 is very roughly described by $$\frac{{\sigma }_{M_0}}{{\sigma }_{t_{\text{acc,0}}}}>1.$$(7)

In general, if Eq. (7) is satisfied, the population of disks that we consider can be described, in first approximation, as having constant t_acc,0. For such a population, the slope of the M_D − Ṁ_* correlation is γ = 1 (for the analytical calculations, see Appendix B).

In Fig. 4, we show the M_D − Ṁ_* correlation in the regime described by Eq. (7). Here, M₀ and t_acc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M_⊙ with a spread of 1.5 dex and on the natural logarithm of 2 Myr with a spread of 0.15 dex, respectively. The gray line is the fit of a single power law, which results in γ = 1. It is important to note that 4 shows that if ${\sigma }_{t_{\text{acc},0}}>{\sigma }_{M_0}$ , we would not be observing a linear correlation between M_D and Ṁ_* ; we would instead see the boomerang isochrone shape shown in Fig. 1, which likewise depends on the value of ω.

We note that Eq. (7) is only an approximate criterion. The exact criterion under which the M_D − Ṁ_* correlation is a power law with a slope of unity is more complex and depends on the age and on the value of ω. In particular, we find that for low values of ω (see Appendix A) the correlation is harder to get, in particular at late time. For detailed studies we thus recommend using Figs. 3, A.1, and A.2 to identify the regions of the parameter space that give a correlation instead of Eq. (7).

We can summarize the results obtained as follows:

1. In the case in which there is no spread in M₀ , a boomerangshaped isochrone appears in the M_D − Ṁ_* plane;

2. In the limit ${\sigma }_{M_0}>>{\sigma }_{t_{\text{acc},0}}$, the slope of the M_D − Ṁ_* relation is exactly 1, a result that can also be obtained analytically;

3. In the general case, we find that the correlation is recovered when the spread in initial disk mass is typically larger than the spread in initial accretion timescale. This rule of thumb can be refined using Figs. 3, A.1, and A.2.

However, observations provide us with more than a correlation; they give the spread around this correlation and the power law index. In particular, observational results (e.g., Manara et al. 2016; Testi et al. 2022) consistently find values of the M_D − Ṁ_* correlation that are <1 along with a large spread. Therefore, in the next section, we study under which conditions we can retrieve the observed value of the slope of the M_D − Ṁ_* correlation and its spread.

Fig. 2 MD − Ṁ* plane for different values of ${\sigma }_{M_0}$. The colored dots in every plot are Monte Carlo simulations of 105 proto- planetary disks with ω = 0.8 and t = 5 Myr, evolving following Eqs. (4), (5), where M0 follows a lognormal distribution centered on the natural logarithm of 0.1 M⊙ with a spread of (a) 0.04 dex, (b) 0.22 dex, (c) 0.30 dex, (d) 0.65 dex. The black dots show the “boomerang” shape of the isochrones for a disk population with the same parameters as in the other disk populations, but with ${\sigma }_{M_0}=0$. For the four plots, tacc,0 follows a lognormal distribution centered on the natural logarithm of 1 Myr, with a spread of 0.26 dex.

Fig. 3 ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$ plane with the values of R2 for different values of age t. The colored dots in every plot are Monte Carlo simulations of 105 proto-planetary disks with ω = 0.8, evolving following Eqs. (4), (5), where M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M⊙ and 1 Myr, respectively. The hatched contours represent the regions where R2 = 0, while the dashed line represents the points where ${\sigma }_{M_0}={\sigma }_{t_{\text{acc},0}}$.

Fig. 4 MD − Ṁ* correlation with ω = 0.8 and t = 2 Myr. The colored dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4), (5) where M0 and tacc,0 follow a lognormal distribution centered, respectively, on 0.1 M⊙ with a spread of 1.5 dex and on 2 Myr with a spread of 0.15 dex. The gray line represents the best fit of a single power law, for which γ = 1.

2.2 The slope and the spread of the M_D − Ṁ_* correlation as a function of ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$

In order to constrain the values of ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$ that can repro- duce the observed slope and spread of the M_D − Ṁ_* correlation, we performed a linear regression of the disks in the M_D −Ṁ_* plane for different combinations of ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$. In Fig. 5, we show the results of our analysis. Here, the colored stripes in the ${\sigma }_{t_{\text{acc},0}}-{\sigma }_{M_0}$ plane represent the regions of the parameter space that predict values of the slope y and of the spread σ of the M_D − Ṁ_* correlation in agreement with the observed values. In a similar fashion to the analysis performed in Lodato et al. (2017), we show in orange the region γ ∈ [0.5,0.9], while in light blue the region σ ∈ [0.55,0.71] dex. In contrast with Lodato et al. 2017 and Tabone et al. 2022b, we neglected the effect of accretion variability and the uncertainties on the mass estimates. This means that the populations tend to underestimate the spread in M_D/Ṁ_*. In this plot, t = 2 Myr, ω = 0.8, and M₀ and t_acc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M_⊙ and 2 Myr, respectively. The cross represents the work of Tabone et al. (2022b). Fig. 5 shows that it is possible to retrieve the observed values of the slope and the spread of the M_D− Ṁ_* correlation for ${\sigma }_{M_0}\in [0.50,1.50]$ and ${\sigma }_{t_{\text{acc}.0}}\in [0.40,0.75]$.

We note that the results obtained here are valid for a fixed age t, a value of ω = 0.8, and the fixed median values M₀ = 0.1 M_⊙ and t_acc,0 = 2 Myr. As we show in Appendix A, changing ω will change the derived values of ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$ , but it remains true that it is possible to retrieve the observed slope and the spread of the M_D− Ṁ_* correlation for a wide range of spreads in disk mass and accretion timescales in the initial conditions. We also note that in Fig. 5 not all the populations fit the disk fraction. As discussed below and in Tabone et al. (2022b), a typical value of ${\sigma }_{t_{\text{acc},0}}\simeq 0.4$ dex needs to be adopted to fit the observed decline of disk fraction with cluster age.

Fig. 5 ${\sigma }_{t_{\text{acc},0}}-{\sigma }_{M_0}$ plane with the values of the slope γ and the spread of the MD − Ṁ* correlation σ. Shown in orange are the values γ ∈ [0.5, 0.9], while in light blue the values σ ∈ [0.55, 0.71] dex, as they confine the region of values compatible with the values found in Manara et al. (2016). The colored regions show the value of γ and σ for Monte Carlo disk populations obtained with the corresponding specific values of ${\sigma }_{t_{\text{acc},0}}$ and ${\sigma }_{M_0}$. Every Monte Carlo consists of 104 disks ruled by Eqs. (4), (5), where ω = 0.8, t = 2 Myr, while M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M and 2 Myr, respectively. The cross shows the case investigated by Tabone et al. (2022b), while the dashed line highlights the value ${\sigma }_{t_{\text{acc},0}}=0.4$.

2.3 The spread of log(M_D− Ṁ_*) as a function of the spread in t and t_acc

Testi et al. (2022) and Somigliana et al. (2023) show that the spread of the accretion timescale M_D/Ṁ_* is a fundamental observable for disk population studies. Considering Eqs. (4), (5), we can evaluate the dimensionless spread of the accretion timescale when Ṁ_* ∝ M_D. Using the propagation of errors (see more details in Appendix C) on log(M_D/Ṁ_*), the dimensionless spread is $$\begin{array}{l}\sigma \left(\log \left(M_D/{\dot{M}}_{*}\right)\right)\simeq \frac{1}{2\langle t_{\text{acc},0}\rangle -\omega \langle t\rangle }\\ \sqrt{{\left(2\langle t_{\text{acc},0}\rangle {\sigma }_{t_{\text{acc},},0}\right)}^2+{\left(\omega \langle t\rangle {\sigma }_t\right)}^2,}\end{array}$$(8)

where we use the notation 〈t_acc,0〉 and 〈t〉 to remark this is the median value for the population, and we also allow for the possibility of an age spread σ_t in the population. This result is equivalent to the spread of the M_D − Ṁ_* correlation for a slope mildly different than 1, but we note that if the M_D − Ṁ_* correlation differs substantially from linearity, then Eq. (8) is not valid anymore.

The fundamental assumption underlying Eq. (8) is that ${\sigma }_{t_{\text{acc}}}$ is taken to be constant in time. In reality, this value should be computed only over the population of surviving disks. It is not necessarily true that it stays constant as, while evolving with time, a disk population loses the fastest evolving disks; as a consequence, ${\sigma }_{t_{\text{acc}}}$ should decrease as time passes. However, as noted in the work of Somigliana et al. (2023), the spread on the accretion timescale decreases slowly in time in the MHD wind-driven scenario (a factor of 0.1 dex within 2 Myr). Therefore, Eq. (8) can be considered valid for young star-forming regions. Moreover, we note that the spread of log(M_D /Ṁ_*) is equal to that of the M_D − Ṁ_* correlation when the slope of the correlation is exactly 1. Therefore, in first approximation, the expression reported in Eq. (8) gives us an estimate of the order of magnitude of the spread of the M_D − Ṁ_* correlation, as further discussed in Sect. 4.1.

In Fig. 6, we show the spread of log(M_D /Ṁ_*) when t follows a lognormal distribution centered on the natural logarithm of 1 Myr, and σ_t = 0.43 dex. The black line is computed as the standard deviation from the Monte Carlo simulations of 10⁵ disks that follow a lognormal distribution in t, t_acc and M₀. In particular, the value of σ(log(M_D /Ṁ_*)) is evaluated for 10² values of ${\sigma }_{t_{\text{acc}}}$ that span [0, 1.1] dex, while σ_t = 0 dex. Here, M₀ and t_acc follow a lognormal distribution centered, namely, on the natural logarithm of 5 Myr for t_acc and on the natural logarithm of 0.1 M_⊙ for M₀. The spread of log(M_D/Ṁ_*) is independent of the spread in M₀ under the conditions for which Eq. (8) holds, which is that the slope of the M_D − Ṁ_* correlation is close to unity. In Fig. 5 we see, as we get close to γ = 1, that the cyan band representing the spread of the M_D − Ṁ_* correlation becomes vertical, while as we move away from those values the dependence on ${\sigma }_{M_0}$ appears.

As shown in Testi et al. (2022) and Somigliana et al. (2023), the spread on the accretion timescale is an important observable that can be easily obtained from the data. We show a comparison with observations in Sect. 4.1.

Until now in the current paper, we considered disk populations for which M₀ and t_acc,0 are not correlated. First, we derived a general criterion for which we can retrieve a population of disks in the M_D − Ṁ_* plane that can be fitted with a single power law. Second, we studied under which initial conditions it is possible to reproduce the observed slope and spread of the M_D − Ṁ_* correlation. Finally, we provided an approximated analytical expression for the spread on the accretion timescale log(M_D /Ṁ_*). In the following, we study the emergence of the M_D − Ṁ_* correlation when M₀ is correlated to t_acc,0.

Fig. 6 Spread of the MD − Ṁ* correlation as a function of ${\sigma }_{t_{\text{acc}}}$, with σt = 0.43 dex. The values of σ(log(MD/Ṁ*)) are evaluated as the standard deviation from Monte Carlo simulations of 105 protoplanetary disks with M0 , tacc , and t following a lognormal distribution centered respectively in 0.1 M⊙, 5 Myr, and 1 Myr.

3 The M_D − Ṁ_* correlation in the presence of a M₀−t_acc,0 correlation

So far we have considered the initial accretion timescale t_acc,0 to be uncorrelated with the initial mass of the disk M₀. However, this is not necessarily the case. In the following, we assume a correlation between M₀ and t_acc,0 in the form $$M_0\propto t_{\text{acc},0}^{\varphi }.$$(9)

Inserting Eq. (9) in Eqs. (4), (5), it is not possible to find an analytical derivation of the slope in the form of Eq. (B.3). However, it is possible to study the analytical behavior of the slope in the initial and final stages of the evolution of the disk, starting from the definition (following Tabone et al. 2022a) of $$t_{disp}=\frac{2t_{\text{acc},0}}{\omega }$$(10)

as the disk dispersal time.

3.1 The M_D − Ṁ_* correlation at initial time

When taking a population of disks without spread in M₀ , in the initial conditions (i.e., t/t_disp ≪1) we find that the slope is γ = 1 −1/ф (see Appendix D.1 for more details).

In Fig. 7, we show the M_D − Ṁ_* correlation for a disk population (dots) of 10⁵ disks ruled from Eqs. (4), (5), where M₀ is given by Eq. (9). Here, ф = 0.8, ω = 0.8, t = 0 yr, and t_acc,0 follows a lognormal distribution centered on the natural logarithm of2 Myr with a spread of 0.13 dex. The black line shows the analytical M_D − Ṁ_* correlation (Eq. (D.2)), which agrees remarkably well with the Monte Carlo simulations.

3.2 The slope close to disk dispersal

When considering a population of disks without spread in M₀ , close to dispersal (i.e., t/t_disp ~ 1) we find that γ = 1 − ω (see Appendix D.2 for more details). This is an unexpected result, as the slope appears to reach a specific value that is independent from the initial correlation between M₀ and t_acc,0. This occurs when a disk population is close to dispersal because all the disks also share a similar t_acc,0 and, as a consequence, a similar M₀ (as they are linked to each other). Thus, when t/t_disp ~ 1, the correlation between M₀ and t_acc,0 is lost and the memory of the initial conditions is forgotten. Tabone et al. (2022a) find a similar result, even without explicitly introducing a correlation. This supports the previous statement that, taking a population of disks that is about to be dispersed, it is only possible to observe the final evolutionary track of the disk (for which ${\dot{M}}_{*}\propto M_D^{1-\omega }$), thus losing any possibility of detecting the initial correlation between M₀ and t_acc,0.

In Fig. 8, we show the M_D−Ṁ_* correlation for ϕ = 0.8, ω = 0.8, and t = 3 Myr, where the other parameters are the same as in Fig. 7. We see in this case that the M_D−Ṁ_* correlation has a bent shape that could be better described by a double power law. We are interested in the left part of the plot, where M_D ≲ 10⁻² M_⊙ (i.e., when the disk lifetime is similar to the disk dispersal time). The black line shows Eq. (D.6): the slope of the theoretical line (black) coincides with the M_D−Ṁ_* correlation in the above-cited region of the M_D−Ṁ_* plane, showing that when t ~ t_{dis p}, the slope approaches the value γ = 1 − ω. The right part of Fig. 8 instead shows the disks that are still evolving, and hence still near the condition t/t_disp ≪ 1; we see that the slope here inverts its trend, in agreement with the value 1 − 1 /ϕ that we derived in Sect. 3.1.

In Appendix E, we show the calculation of the spread of the M_D−Ṁ_* correlation when $M_0\propto t_{\text{acc},0}^{\varphi }$. The (trivial) results are in continuity with what we describe in Sect. 2.3.

Fig. 7 Slope of the MD−Ṁ* correlation with ϕ = 0.8, ω = 0.8, and t = 0 yr. The dots are Monte Carlo simulations of 105 proto-planetary disks ruled from Eqs. (4) and (5) where M0 follows Eq. (9) and tacc,0 follows a lognormal distribution centered on the natural logarithm of 2 Myr with a spread of 0.13 dex. The black line represents Eq. (D.2).

Fig. 8 Slope of the MD−Ṁ* correlation with ϕ = 0.8, ω = 0.8, and t = 3 Myr. The colored dots are Monte Carlo simulations of 105 protoplanetary disks ruled from Eqs. (4) and (5), where M0 follows Eq. (9) and tacc,0 follows a lognormal distribution centered on the natural logarithm of 2 Myr with a spread of 0.13 dex. The black line represents Eq. (D.6).

Fig. 9 MD−Ṁ* correlation with ϕ = 1, ω = 0.8, and t = 3.5 Myr. The dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4) and (5), where M0 follows Eq. (9) and tacc,0 follows a lognormal distribution centered on the natural logarithm of 5 Myr with a spread of 0.13 dex. In addition, M0 is then generated following a lognormal distribution with the values of tacc,0 following Eq. (9) with a spread of 0.43 dex. The black line represents Eq. (D.2), the gray line is the best-fitting power law for the population, resulting in ${\dot{M}}_{*}\propto M_D^{0.83}$, while the red line shows Eq. (B.2). For illustrative purposes, we also plot in blue the best-fitting power law for a population for which ${\sigma }_{M_0}=0.26$ dex, which results in ${\dot{M}}_{*}\propto M_D^{0.65}$.

3.3 The role of ${\sigma }_{M_0}$ in the M_D−Ṁ_* correlation

If t_acc,0 is generated following a lognormal distribution with a spread and $M_0\propto t_{\text{acc},0}^{\varphi }$ is satisfied, it is reasonable to think that M₀ should also have a spread of its own (i.e., the spread in M₀ reflects the fact that different stars living in the same stellar population have different masses, not only different initial timescales). To ensure that our populations have a spread in M₀, we initially evaluated M₀ with Eq. (9); then M0 was extracted with a lognormal distribution, with fixed spread, centered on the previously determined values of M₀. In Fig. 9, once the spread in M₀ is added, we see the appearance of the linear M_D−Ṁ_* relation. Here we show a population of disks identical to the one in Fig. 8, but we add a spread in M₀ of 0.43 dex. The addition of ${\sigma }_{M_0}$ hides the track shown in Fig. 8. Moreover, the slope becomes steeper as ${\sigma }_{M_0}$ increases, going toward the linear behavior, which is plotted in red. As a reference, we plotted in gray the best-fitting power law for the population, resulting in a slope of 0.83; in blue the best-fitting power law for a population with ${\sigma }_{M_0}=0.26$ dex, resulting in a slope of 0.65; while in black we report the slope obtained in Eq. (D.3).

We draw two main conclusions from this. First, it is still possible to derive the slope of the M_D−Ṁ_* correlation with a value different than unity after the introduction of the correlation $M_0\propto t_{\text{acc},0}^{\varphi }$. It is also possible to tune the value of the M_D−Ṁ_* slope, changing the values of ϕ and ${\sigma }_{M_0}$: we discuss this in Sect. 4.2. Second, the addition of ${\sigma }_{M_0}$ deletes the canonical signature of the slope shown in Figs. 7 and 8. This leads to a dramatic observational conclusion: it is not possible to directly observe the expected analytical signature of the slope because of the natural existence of spread in M₀.

3.4 Describing the M_D−Ṁ_* correlation by a single power law

It is relevant to emphasize that, after the introduction of a spread in M₀, the bent shape of the M_D−Ṁ_* correlation shown in Fig. 8 can still be recovered, but only for very low values of ${\sigma }_{M_0}$.³ Figure 8 clearly shows that a single power law cannot properly fit the evolutionary track of the disks: a double power law seems much more appropriate. Hence, we conducted a test in order to quantitatively understand when a single power law starts to fit a disk population.

First, we introduced a broken double power law in the form $$\{\begin{array}{ll}{\dot{M}}_{*}=k_1+a{M}_D^{-b},\hfill & M_D<x_0\hfill \\ {\dot{M}}_{*}=k_2+c{\left(M_D-x_0\right)}^d,\hfill & \hfill \\ \hfill & M_D\ge x_0\hfill \end{array}$$(11)

where $k_2=k_1+a{x}_0^{-b}$ is not a free parameter as the piecewise function is defined to be continuous in M_D = x₀. Then, after setting t = 3 Myr, we fitted both Eq. (11) and the equation of a single power law in the form $${\dot{M}}_{*}=a{M}_D^{\zeta }$$(12)

to the Ṁ_*−M_D plane for different values of ${\sigma }_{M_0}$, namely ${\sigma }_{M_0}\in \left[0,0.022\right]$ dex. The other parameters for the disk population are the same as in Fig. 9. Since Eq. (11) is defined based on six free-parameters (x₀, k₁, a, b, c, d), while Eq. (12) on only two (a, ζ), we performed the fit using the Bayesian Information Criterion (BIC), with the Python package RegscorePy.bic.

In Fig. 10, we show that the BIC is minimized from the double power law in the regime ${\sigma }_{M_0}\in \left[0,0.01\right]$ dex, hence showing that for very low values of ${\sigma }_{M_0}$, the M_D−Ṁ_* correlation is better fitted with a double power law. However, as little spread in M₀ is added, the red and the black curve become indistinguishable, highlighting the fact that those disk populations can be fitted both with a single or a double power law⁴.

From this we can conclude that when a little spread in Mo is added, the M_D−Ṁ_* correlation can be described by a single instead of a double power law. Thus, it is reasonable to use a single power law to fit the M_D−Ṁ_* correlation. We show the same plot of Fig. 10 in Appendix F for ${\sigma }_{M_0}\in \left[0,0.43\right]$ dex.

4 Comparison with the data. Empirical constraints on the spreads

4.1 MHD wind-driven accretion can reproduce the observed spread of log(M_D−Ṁ_*)

To assess whether MHD winds can reproduce both disk lifetimes and the observed spread in log(M_D−Ṁ_*), we conducted the following analysis.

So far, the spread in t_acc,0 has been considered as a free parameter. As shown by Tabone et al. (2022b), though with a different ansatz for the distribution of t_acc,0, the observed disk fraction sets both ${\sigma }_{t_{\text{acc},0}}$ and the median values of t_acc,0. Following the work of Tabone et al. (2022b), we fitted the disk fraction (represented as the fraction of sources that show IR excess) reported in Fedele et al. (2010) and references therein. We first defined our disk fraction as the ratio of disks that, at a certain time t, satisfies the condition $$t_{disp}>t.$$(13)

Then we performed a fit using the complementary cumulative distribution function (CCDF) assuming a lognormal distribution of t_disp. In this way, we find the parameters that describe the probability for a disk population to satisfy Eq. (13).

In Fig. 11, we show the results of the fit. We find that tdisp follows a lognormal distribution with a median value ${\mu }_{t_{disp}}=2.37$ Myr and a spread ${\mu }_{t_{disp}}=0.37$ dex. After some sim-pie algebra deriving from Eq. (10), we find that the median value of the accretion timescale is ${\mu }_{t_{\text{acc},0}}=0.95$ Myr and the spread is ${\sigma }_{t_{\text{acc},0}}=0.37$ dex. With these values, we evaluated the spread of the accretion timescale using Eq. (8) with t = 1.5 Myr and σ_t = 0⁵ dex, obtaining σ(log(M_D/Ṁ_*)) = 1.01 dex.

We compared this value with the spread of the accretion timescale evaluated from the data found in Testi et al. (2022) for the four populations of L1688, Lupus, ChaI, and USco. We removed the upper limits from the sample and considered only the systems for which M_* > 0.15 M_⊙, in agreement with Testi et al. (2022). For L1688 and Lupus, we obtained values of σ(log(M_D/Ṁ_*)) that are closely comparable with that reported above, more specifically σ(log(M_D/Ṁ_*))_L1688 = 0.86 dex, σ(log(M_D/Ṁ_*))_Lupus = 1.08 dex. As we moved to the older region of ChaI, we obtained the higher spread σ(log(M_D/Ṁ_*))_ChaI = 1.62 dex, a value that increases even more for USco, the oldest region of the sample, for which we obtained σ(log(M_D/Ṁ_*)_{US co} = 2.02 dex. The reason behind these higher spreads is that to evaluate σ(log(M_D/Ṁ_*)), we assumed that the distribution of σ(log(M_D/Ṁ_*) is normal. This assumption is valid only for young populations of disks (e.g., L1688, Lupus, and at most ChaI); as the population evolves in time, we lose disks, and the distribution of σ(log(M_D/Ṁ_*) significantly deviates from the starting distribution.

We hence conclude that the spread predicted using Eq. (8) can reproduce well the spread of the accretion timescale for young populations of disks. This result is in agreement with the analysis of Tabone et al. (2022a), but in this work we assumed a lognormal distribution to fit for t_disp, as already performed by Somigliana et al. (2023), and we applied the results to different stellar regions. We furthermore recall that for our calculation we set σ_t = 0 dex, and we neglected accretion variability and mass uncertainty, so the value of σ(log(M_D/Ṁ_*)) that we predicted must be considered as a lower limit.

Fig. 10 BIC as function of ${\sigma }_{M_0}$ in a range ${\sigma }_{M_0}\in \left[0,0.022\right]\text{\hspace{0.17em}dex}$. The BIC is obtained for Eq. (11) (in red) and Eq. (12) (in black) from the fit to the Ṁ*−MD plane obtained for different values of ${\sigma }_{M_0}$, while the simulated disk populations have the same set of parameters described in Fig. 8.

Fig. 11 Disk fraction as a function of the cluster age (data from Fedele et al. 2010 and references therein), with the fit of the CCDF of tdisp. From the fit, we obtain tdisp = 2.37 Myr and ${\sigma }_{t_{\text{disp }}}=0.37$ dex.

4.2 Empirical constraints on ${\sigma }_{M_0}$ and ϕ

We conclude this work by showing the existence of a degeneracy between ϕ and ${\sigma }_{M_0}$. As pointed out in Sect. 3.3, the slope of the M_D−Ṁ_* correlation is affected by both ϕ and ${\sigma }_{M_0}$. From Monte Carlo simulations it is possible to investigate the $\varphi -{\sigma }_{M_0}$ space in order to constrain for which values of ϕ and ${\sigma }_{M_0}$ it is possible to retrieve a slope γ ∊ [0.5 ,0.9], with a spread on the M_D−Ṁ_* correlation σ ∊ [0.55 ,0.71] dex (Manara et al. 2016). To do this, we simply evaluated γ and the spread of the M_D−Ṁ_* correlation with a linear regression of the disks residing in the M_D−Ṁ_* plane. Each Monte Carlo simulation consists of 10⁴ disks evolving according to Eqs. (4), (5), where $M_0\propto t_{\text{acc},0}^{\varphi }$.

In Fig. 12, we show the $\varphi -{\sigma }_{M_0}$ plane with the corresponding values of γ. The coloured regions in the $\varphi -{\sigma }_{M_0}$ plane show the value of the slope and of its spread obtained via linear regression of a whole population of disks generated with specific values of ϕ and ${\sigma }_{M_0}$. In a similar fashion to the analysis performed in Lodato et al. (2017), we show in orange the region γ ∊ [0.5 , 0.9], while in light blue the region σ ∊ [0.55 , 0.71] dex. In this case t = 2 Myr, ω = 0.8, and t_acc,0 follows a lognormal distribution centered on the values derived in Sect. 4.1. For comparison, we marked with a cross the region covered in Tabone et al. (2022b).

Following Sect. 3.3, we found that γ tends toward unity as ${\sigma }_{M_0}$ increases, regardless of the value of ϕ. However, there is a whole region of the ${\sigma }_{M_0}-\varphi$ plane that allows values of γ ≠ 1; this confirms the existence of a degeneracy between ϕ and ${\sigma }_{M_0}$. Moreover, following Fig. 12, we can constrain the ensemble of values of ${\sigma }_{M_0}$ that allow us to have the observed values of γ and σ; we see that only values of ϕ ∊ [−2,1] and ${\sigma }_{M_0}\in \left[0,1.75\right]$ dex are allowed. Figure 12 leads us to understand that for higher values of ϕ we need higher values of ${\sigma }_{M_0}$ in order to retrieve the observed M_D−Ṁ_* correlation. This can be explained as follows: the stronger the correlation between M₀ and t_acc,0 is, the bigger the spread in M₀ must be in order to break the bent shape shown in Fig. 8.

To conclude, we report that the range of values we found for σM₀ is in agreement with the results found in Tabone et al. (2022a) for modeling the Lupus population, although there they explored only the case with no correlation between M₀ and t_acc,0. Moreover, we note that in this section we show a procedure to constrain the parameters that rule the evolution of MHD wind-driven disks in the presence of a correlation between M₀ and t_acc,0.

From Fig. 12, we understand that we can obtain values of the slope of the M_D−Ṁ_* correlation which are compatible with Manara et al. (2016), regardless of the initial conditions of the disk populations. We hence conclude that it is possible to derive the slope of the M_D−Ṁ_* correlation in the MHD wind-driven scenario under a broad range of conditions, a result that shows that the M_D−Ṁ_* correlation can also be easily obtained for this evolutionary paradigm.

Fig. 12 $\varphi -{\sigma }_{M_0}$ plane with the values of the slope γ and the spread of the MD−Ṁ* correlation σ. Shown in orange are the values γ ∊ [0.5,0.9], while in light blue the values σ ∊ [0.55,0.71] dex as they confine the region of values compatible with the values found in Manara et al. (2016). The colored regions show the value of γ and σ for Monte Carlo disk populations obtained with the corresponding specific values of ϕ and ${\sigma }_{M_0}$. Every Monte Carlo consists of 104 disks evolving according to Eqs. (4), (5), where $M_0\propto t_{\text{acc},0}^{\varphi }$, ω = 0.8, t = 2 Myr, and tacc,0 follows a lognormal distribution centered on the values retrieved from the fit in Sect. 4.1. The cross shows the case investigated by Tabone et al. (2022b).

5 Conclusions

In this paper, we investigated under which conditions the M_D−Ṁ_* correlation emerges in the MHD wind-driven evolutionary scenario.

We derived and analyzed the parameters that rule the appearing and tuning of the M_D−Ṁ_* correlation starting from the analytical solutions given by Tabone et al. (2022b). From our analysis we reach the following conclusions:

• (i)

  Introducing a spread ${\sigma }_{t_{\text{acc},0}}$, we obtain the correlation between the mass of the disk M_D and the accretion onto the central star Ṁ_*. This holds as long as ${\sigma }_{t_{\text{acc}.0}}<{\sigma }_{M_0}$, a conservative criterion that is compatible with the general rule of thumb R² ≳ 0.5. In the opposite case, the M_D−Ṁ_* correlation tends to the boomerang-shaped isochrone introduced by Lodato et al. (2017) for the viscous case and studied by Tabone et al. (2022a), and Tabone et al. (2022b) for the wind case, which look significantly different from the observed correlation;

• (ii)

  It is possible to derive a slope and a spread of the M_D−Ṁ_* correlation comparable to that found in Manara et al. (2016) varying ${\sigma }_{M_0}$ and ${\sigma }_{t_{\text{acc},0}}$. Hence, it is possible to derive the observed slope and spread of the M_D−Ṁ_* correlation under a range of initial conditions;

• (iii)

  Neglecting the time evolution of the accretion timescale, we derived an analytical expression for the spread of the accretion timescale log(M_D−Ṁ_*). The predictions arising from this formula are in agreement with the observed order of magnitude of σ(log(M_D−Ṁ_*)) evaluated for young populations of disks (Testi et al. 2022);

• (iv)

  We introduced a correlation $M_0\propto t_{\text{acc},0}^{\varphi }$ to understand if we can still derive the observed slope and spread of the M_D−Ṁ_* correlation. We analytically derived a slope γ = 1 − 1/ϕ in the initial conditions, and a slope γ = 1 − ω close to dispersal.

  Furthermore, we found that the addition of a minimal spread in the initial disk mass distribution $\left({\sigma }_{M_0}\right)$ affects the slope of the M_D−Ṁ_* correlation. In particular, we showed that we can fit the M_D−Ṁ_* correlation with a single power law because of the existence of ${\sigma }_{M_0}$;

• (v)

  We studied the degeneracy between ϕ and ${\sigma }_{M_0}$, and we confirmed that we can obtain values of the slope and of the spread of the M_D−Ṁ_* correlation in agreement with Manara et al. (2016) for a broad range of values on both ${\sigma }_{M_0}$ and ϕ.

This work is the natural follow-up of the work of Tabone et al. (2022a) and Tabone et al. (2022b). We investigated in greater detail the emergence of the M_D−Ṁ_* correlation and demonstrated that the slope of the correlation is controlled by the correlation among the parameters describing the initial conditions (physically speaking, the star formation process). With this, we reached the fundamental conclusion that MHD wind-driven disks can reproduce the observed slope of the M_D−Ṁ_* correlation under a broad range of initial conditions. Therefore, no fine tuning of the model parameters is needed in the MHD scenario, but reproducing the correlation does exclude some regions of the parameter space, as shown by Figs. 5, and 12. While in this paper we focused only on the MHD wind scenario, our work implies that the observed M_D−Ṁ_* correlation cannot be used to discriminate between MHD winds and the viscous model. However, the correlation still constrains the model parameters for each scenario, as shown, for example, by our work and Lodato et al. (2017).

Data availability

All data used in this work are publicly available and can be found at the corresponding cited works. The scripts that have been used for this work will be shared under reasonable request to the corresponding author.

Acknowledgements

We are thankful to the anonymous referee who helped improve the clarity of the paper. The authors acknowledge support from the European Union (ERC Starting Grant DiscEvol, project number 101039651), from Fondazione Cariplo, grant No. 2022-1217, from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 823823 (Dustbusters RISE project), and from the ERC Synergy Grant “ECOGAL” (project ID 855130). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. We thank Carlo Manara for the very useful comments he provided for ensuring the correctness of this work. Moreover, we thank Giacomo Lucertini for the useful discussions in statistics, Edoardo Merli for IT support, Davide Giovagnoli for interesting dialogues in mathematics, Enrico Ragusa and Chiara Scardoni for the many chats in physics.

Appendix A The ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc}.0}}-R^2$ plane for different values of ω

In Figs. A.1, A.2 we show the ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc}.0}}$ plane for ω = 0.2 and ω = 1.0, respectively. All the other parameters are the same as the ones presented in Fig. 3. We highlight that the “rule of thumb” derived in Sect. 2 holds also for these “extreme” values of ω. Moreover, we note that, in Fig. A.1 for evolved populations our general criterion R² ≳ 0.5 does not correspond anymore necessarily to the region ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc}.0}}$. The reason for this is shown in Fig. A.3: for evolved populations, when ω → 0, the disks collect around the evolved region of the “boomerang-shaped” isochrone, thus being effectively fitted with a single power law. However, we note that our criterion is very general, therefore any simulated population based on Figs. A.1, A.2, 3, must be checked directly in the M_D−M_* plane.

Fig. A.1 ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$ plane with the values of R2 for different values of age t. The colored dots in every plot are Monte Carlo simulations of 105 proto-planetary disks with ω = 0.2, evolving following Eqs. (4), (5), where M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M⊙ and 1 Myr, respectively. The hatched contours represent the regions where there is no correlation between MD and Ṁ*, while the dashed line represents the points where ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$.

Fig. A.2 ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$ plane with the values of R2 for different values of age t. The colored dots in every plot are Monte Carlo simulations of 105 proto-planetary disks with ω = 1.0, evolving following Eqs. (4), (5), where M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of o.1 M⊙ and 1 Myr, respectively. The hatched contours represent the regions where there is no correlation between MD and Ṁ*, while the dashed line represents the points where ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$.

Fig. A.3 MD−Ṁ* correlation with ω = 0.2 and t = 2.5 Myr. The colored dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4), (5) where M0 and tacc,0 follow a lognormal distribution centered, respectively, on the natural logarithm of 0.1 M⊙ with a spread of 0.8 dex and on the natural logarithm of 1 Myr with a spread of 0.4 dex. The gray line represents the best fit to the population, resulting in a slope γ = 0.72.

Appendix B The slope of the M_D−Ṁ_* correlation when t_acc,0 is constant: analytics

Considering t_acc,0 as a constant value, Eq. (5) can be seen as $${\dot{M}}_{*}=\frac{M_D}{2t_{\text{acc},0}\left(1+f_{M,0}\right)}{\left(1-\frac{\omega t}{2t_{\text{acc},0}}\right)}^{-1}.$$(B.1)

Taking the hint from the work of Manara et al. 2016, it is possibile to investigate the analytical behavior of disks evolving according to Eq. (B.1), assuming that $${\dot{M}}_{*}\propto M_D^{\gamma }.$$(B.2)

Performing the logarithmic derivative on Eq. (B.2), the slope γ results in the constant $$\gamma =\frac{M_D}{{\dot{M}}_{*}}\frac{ \text{d}{\dot{M}}_{*}}{ \text{d}{M}_D}=1.$$(B.3)

Appendix C The spread of the M_D−Ṁ_* correlation as function of the spread in t and t_acc: analytics

Starting from Eqs. (4), (5), the ratio M_D−Ṁ_* reads $$\frac{M_D}{{\dot{M}}_{*}}=\left(1+f_{M,0}\right)\left(2t_{\text{acc},0}-\omega t\right).$$(C.1)

Neglecting the covariance term between t and t_acc,0, the spread of log(M_D−Ṁ_*) correlation is obtained with a simple propagation of errors as $$\begin{array}{l}\sigma \left(\log \left(M_D/{\dot{M}}_{*}\right)\right)=[{\sigma }_{t_{\text{acc },0}}^2{\left(\frac{\partial \left(\log \left(M_D/{\dot{M}}_{*}\right)\right)}{\partial \log \left(t_{\text{acc},0}\right)}\right)}^2+\\ {{\sigma }_t^2{\left(\frac{\partial \left(\log \left(M_D/{\dot{M}}_{*}\right)\right)}{\partial \log (t)}\right)}^2]}^{1/2},\end{array}$$(C.2)

which gives Eq. (8) after some simple algebra. We note that the quantity ${\sigma }_{t_{\text{acc}.0}}$ is the spread of the accretion timescale. This quantity, due to dispersal of the fastest evolving disks, decreases as a function of time. However, for our calculations, we keep it fixed to its initial value.

Appendix D The slope of the M_D−Ṁ_* correlation when $M_0\propto t_{\text{acc},0}^{\varphi }$: analytics

D.1 The slope for t/t_disp << 1

Injecting Eq. (9) into Eq. (4), we find that, in the initial conditions of the evolution of the disk, namely when t/t_disp << 1, $$t_{\text{acc,0}}\propto M_D^{1/\varphi }.$$(D.1)

Including this result in Eq. (5), after the trivial substitution of Eq. (9) in Eq. (5), we obtain $${\dot{M}}_{\ast }\propto M_D^{1-1/\varphi }.$$(D.2)

The logarithmic derivative gives the slope in the form $$\gamma =\frac{M_D}{{\dot{M}}_{\ast }}\frac{\text{d}{\dot{M}}_{\ast }}{\text{d}{M}_D}=1-\frac{1}{\varphi }.$$(D.3)

D.2 The slope for t ~ t_disp

We derive the slope in the final life-time conditions of a disk starting by choosing ϵ to be arbitrarily small. Then, inserting $$\frac{t}{t_{disp}}=1-ϵ$$(D.4)

in Eqs. (4), (5) after the previous substitution of Eq. (9), ϵ is given by $$ϵ\propto {\left(\frac{M_D}{t_{\text{acc},0}^{\varphi }}\right)}^{\omega }.$$(D.5)

After some algebra, it is obtained that $${\dot{M}}_{\ast }\propto t_{\text{acc},0}^{\varphi \omega -1}{M}_D^{1-\omega },$$(D.6)

which leads to a slope of the form $$\gamma =\frac{M_D}{{\dot{M}}_{\ast }}\frac{\text{d}{\dot{M}}_{\ast }}{\text{d}{M}_D}=1-\omega .$$(D.7)

Appendix E The spread of the M_D−Ṁ_* correlation when $M_0\propto t_{\text{acc},0}^{\varphi }$

E.1 The spread when $M_0\propto t_{\text{acc},0}^{\varphi }$ and t << t_disp

Following Sect. 2.3, the calculations for the spread when $M_0\propto t_{\text{acc},0}^{\varphi }$ lead to $$\sigma \left(\log \left(M_D/{\dot{M}}_{\ast }\right)\right)={\sigma }_{t_{\text{acc}}}.$$(E.1)

This shows the trivial fact that, in the initial conditions, the spread of the correlation depends uniquely on the spread of t_acc,0, a result already shown in Fig. 6.

E.2 The spread when $M_0\propto t_{\text{acc},0}^{\varphi }$ and t ~ t_disp

A similar result is obtained when t ~ t_disp, obtaining $$\sigma \left(\log \left(M_D/{\dot{M}}_{\ast }\right)\right)=ϵ{\sigma }_{t_{\text{acc}}}.$$(E.2)

It is clear that the more t approaches the value of t_disp, the littler is the value of ϵ. In the limit for ϵ → 0, Eq. (E.2) gives σ(log(M_d/Ṁ_*) → 0: a valid result that reflect the fact that when the disk is about to die, the slope of the M_d−Ṁ_* correlation inevitably tends toward the value 1 − ω.

Fig. E.1 BIC as function of ${\sigma }_{M_0}$ in a range ${\sigma }_{M_0}\in [0,0.43]$ dex. The BIC is obtained for Eq. (11) (in red) and Eq. (12) (in black) from the fit to the Ṁ* − MD plane obtained for different values of ${\sigma }_{M_0}$, while the simulated disk populations have the same set of parameters described in Fig. 8.

Appendix F The likelihood and the BIC in a wider ${\sigma }_{M_0}$ range

In Fig. E.1 we show the same results of Fig. 10, in a range ${\sigma }_{M_0}\in [0,1]$. As can be seen, for significant values of ${\sigma }_{M_0}$, the Ṁ_* − M_D loses the double power law behavior in favor of a single power law written in the form of Eq. (12).

References

All Figures

	Fig. 1 Isochrones of a population of disks with t = 5 Myr for three values of M0 in the MD − Ṁ* plane. The colored dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4), (5) with ω = 0.8, where tacc,0 follows a lognormal distribution centered on the natural logarithm of 1 Myr with a spread of 0.26 dex, while M0 = 0.1, 1.0, 10 M⊙.
In the text

Fig. 2 MD − Ṁ* plane for different values of ${\sigma }_{M_0}$. The colored dots in every plot are Monte Carlo simulations of 105 proto- planetary disks with ω = 0.8 and t = 5 Myr, evolving following Eqs. (4), (5), where M0 follows a lognormal distribution centered on the natural logarithm of 0.1 M⊙ with a spread of (a) 0.04 dex, (b) 0.22 dex, (c) 0.30 dex, (d) 0.65 dex. The black dots show the “boomerang” shape of the isochrones for a disk population with the same parameters as in the other disk populations, but with ${\sigma }_{M_0}=0$. For the four plots, tacc,0 follows a lognormal distribution centered on the natural logarithm of 1 Myr, with a spread of 0.26 dex.

In the text

Fig. 3 ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$ plane with the values of R2 for different values of age t. The colored dots in every plot are Monte Carlo simulations of 105 proto-planetary disks with ω = 0.8, evolving following Eqs. (4), (5), where M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M⊙ and 1 Myr, respectively. The hatched contours represent the regions where R2 = 0, while the dashed line represents the points where ${\sigma }_{M_0}={\sigma }_{t_{\text{acc},0}}$.

In the text

	Fig. 4 MD − Ṁ* correlation with ω = 0.8 and t = 2 Myr. The colored dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4), (5) where M0 and tacc,0 follow a lognormal distribution centered, respectively, on 0.1 M⊙ with a spread of 1.5 dex and on 2 Myr with a spread of 0.15 dex. The gray line represents the best fit of a single power law, for which γ = 1.
In the text

Fig. 5 ${\sigma }_{t_{\text{acc},0}}-{\sigma }_{M_0}$ plane with the values of the slope γ and the spread of the MD − Ṁ* correlation σ. Shown in orange are the values γ ∈ [0.5, 0.9], while in light blue the values σ ∈ [0.55, 0.71] dex, as they confine the region of values compatible with the values found in Manara et al. (2016). The colored regions show the value of γ and σ for Monte Carlo disk populations obtained with the corresponding specific values of ${\sigma }_{t_{\text{acc},0}}$ and ${\sigma }_{M_0}$. Every Monte Carlo consists of 104 disks ruled by Eqs. (4), (5), where ω = 0.8, t = 2 Myr, while M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M and 2 Myr, respectively. The cross shows the case investigated by Tabone et al. (2022b), while the dashed line highlights the value ${\sigma }_{t_{\text{acc},0}}=0.4$.

In the text

	Fig. 6 Spread of the MD − Ṁ* correlation as a function of ${\sigma }_{t_{\text{acc}}}$, with σt = 0.43 dex. The values of σ(log(MD/Ṁ*)) are evaluated as the standard deviation from Monte Carlo simulations of 105 protoplanetary disks with M0 , tacc , and t following a lognormal distribution centered respectively in 0.1 M⊙, 5 Myr, and 1 Myr.
In the text

	Fig. 7 Slope of the MD−Ṁ* correlation with ϕ = 0.8, ω = 0.8, and t = 0 yr. The dots are Monte Carlo simulations of 105 proto-planetary disks ruled from Eqs. (4) and (5) where M0 follows Eq. (9) and tacc,0 follows a lognormal distribution centered on the natural logarithm of 2 Myr with a spread of 0.13 dex. The black line represents Eq. (D.2).
In the text

	Fig. 8 Slope of the MD−Ṁ* correlation with ϕ = 0.8, ω = 0.8, and t = 3 Myr. The colored dots are Monte Carlo simulations of 105 protoplanetary disks ruled from Eqs. (4) and (5), where M0 follows Eq. (9) and tacc,0 follows a lognormal distribution centered on the natural logarithm of 2 Myr with a spread of 0.13 dex. The black line represents Eq. (D.6).
In the text

Fig. 9 MD−Ṁ* correlation with ϕ = 1, ω = 0.8, and t = 3.5 Myr. The dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4) and (5), where M0 follows Eq. (9) and tacc,0 follows a lognormal distribution centered on the natural logarithm of 5 Myr with a spread of 0.13 dex. In addition, M0 is then generated following a lognormal distribution with the values of tacc,0 following Eq. (9) with a spread of 0.43 dex. The black line represents Eq. (D.2), the gray line is the best-fitting power law for the population, resulting in ${\dot{M}}_{*}\propto M_D^{0.83}$, while the red line shows Eq. (B.2). For illustrative purposes, we also plot in blue the best-fitting power law for a population for which ${\sigma }_{M_0}=0.26$ dex, which results in ${\dot{M}}_{*}\propto M_D^{0.65}$.

In the text

	Fig. 10 BIC as function of ${\sigma }_{M_0}$ in a range ${\sigma }_{M_0}\in \left[0,0.022\right]\text{\hspace{0.17em}dex}$. The BIC is obtained for Eq. (11) (in red) and Eq. (12) (in black) from the fit to the Ṁ*−MD plane obtained for different values of ${\sigma }_{M_0}$, while the simulated disk populations have the same set of parameters described in Fig. 8.
In the text

	Fig. 11 Disk fraction as a function of the cluster age (data from Fedele et al. 2010 and references therein), with the fit of the CCDF of tdisp. From the fit, we obtain tdisp = 2.37 Myr and ${\sigma }_{t_{\text{disp }}}=0.37$ dex.
In the text

Fig. 12 $\varphi -{\sigma }_{M_0}$ plane with the values of the slope γ and the spread of the MD−Ṁ* correlation σ. Shown in orange are the values γ ∊ [0.5,0.9], while in light blue the values σ ∊ [0.55,0.71] dex as they confine the region of values compatible with the values found in Manara et al. (2016). The colored regions show the value of γ and σ for Monte Carlo disk populations obtained with the corresponding specific values of ϕ and ${\sigma }_{M_0}$. Every Monte Carlo consists of 104 disks evolving according to Eqs. (4), (5), where $M_0\propto t_{\text{acc},0}^{\varphi }$, ω = 0.8, t = 2 Myr, and tacc,0 follows a lognormal distribution centered on the values retrieved from the fit in Sect. 4.1. The cross shows the case investigated by Tabone et al. (2022b).

In the text

Fig. A.1 ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$ plane with the values of R2 for different values of age t. The colored dots in every plot are Monte Carlo simulations of 105 proto-planetary disks with ω = 0.2, evolving following Eqs. (4), (5), where M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of 0.1 M⊙ and 1 Myr, respectively. The hatched contours represent the regions where there is no correlation between MD and Ṁ*, while the dashed line represents the points where ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$.

In the text

Fig. A.2 ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$ plane with the values of R2 for different values of age t. The colored dots in every plot are Monte Carlo simulations of 105 proto-planetary disks with ω = 1.0, evolving following Eqs. (4), (5), where M0 and tacc,0 follow a lognormal distribution centered on the natural logarithm of o.1 M⊙ and 1 Myr, respectively. The hatched contours represent the regions where there is no correlation between MD and Ṁ*, while the dashed line represents the points where ${\sigma }_{M_0}-{\sigma }_{t_{\text{acc},0}}$.

In the text

Fig. A.3 MD−Ṁ* correlation with ω = 0.2 and t = 2.5 Myr. The colored dots are Monte Carlo simulations of 105 proto-planetary disks evolving following Eqs. (4), (5) where M0 and tacc,0 follow a lognormal distribution centered, respectively, on the natural logarithm of 0.1 M⊙ with a spread of 0.8 dex and on the natural logarithm of 1 Myr with a spread of 0.4 dex. The gray line represents the best fit to the population, resulting in a slope γ = 0.72.

In the text

	Fig. E.1 BIC as function of ${\sigma }_{M_0}$ in a range ${\sigma }_{M_0}\in [0,0.43]$ dex. The BIC is obtained for Eq. (11) (in red) and Eq. (12) (in black) from the fit to the Ṁ* − MD plane obtained for different values of ${\sigma }_{M_0}$, while the simulated disk populations have the same set of parameters described in Fig. 8.
In the text