# The emergence of the $M_D - \dot{M}_*$ correlation in the magnetohydrodynamic wind scenario

Luigi Zallio<sup>1,2,*</sup>①, Giovanni Rosotti<sup>1,3</sup>①, Benoît Tabone<sup>4</sup>①, Leonardo Testi<sup>5</sup>①, Giuseppe Lodato<sup>1</sup>①, and Alice Somigliana<sup>2,6</sup>①

Received 8 May 2024 / Accepted 30 October 2024

## ABSTRACT

*Context.* There is still much uncertainty around the mechanism that rules the accretion of proto-planetary disks. In recent years, magnetohydrodynamic (MHD) wind-driven accretion has been proposed as a valid alternative to the more conventional viscous accretion. In particular, winds have been shown to reproduce the observed correlation between the mass of the disk $M_D$ and the mass accretion rate onto the central star $\dot{M}_*$, but this has been done only for specific conditions. It is not clear whether this implies fine tuning or if it is a general result.

*Aims.* We investigated under which conditions the observed correlation between the mass of the disk $M_D$ and the mass accretion rate onto the central star $\dot{M}_*$ can be obtained.

*Methods.* We present mainly analytical calculations, supported by Monte Carlo simulations. We also perform a comparison with the observed data to test our predictions.

*Results.* In the absence of a correlation between the initial mass $M_0$ and the initial accretion timescale $t_{acc,0}$, we find that the slope of the $M_D-\dot{M}_*$ correlation depends on the value of the spread of the initial conditions of masses and lifetimes of disks. Then we clarify the conditions under which a disk population can be fitted with a single power law. Moreover, we derive an analytical expression for the spread of $\log(M_D/\dot{M}_*)$ valid when the spread of $t_{acc}$ is taken to be constant. In the presence of a correlation between $M_0$ and $t_{acc,0}$, we derive an analytical expression for the slope of the $M_D-\dot{M}_*$ correlation in the initial conditions of disks and at late times. In this new scenario, we clarify under which conditions the disk population can be fitted by a single power law, and we provide empirical constraints on the parameters ruling the evolution of disks in our models.

*Conclusions.* We conclude that MHD winds can predict the observed values of the slope and the spread of the $M_D-\dot{M}_*$ correlation under a broad range of initial conditions. This is a fundamental expansion of previous works on the MHD paradigm, exploring the establishment of this fundamental correlation beyond specific initial conditions.

**Key words.** protoplanetary disks

## 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.,

<footnote n="1">¹ Dipartimento di Fisica 'Aldo Pontremoli', Università degli Studi di Milano, via G. Celoria 16, 20133 Milano, Italy</footnote>

<footnote n="2">² European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="3">³ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="5">⁵ Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40122 Bologna, Italy</footnote>

<footnote n="6">⁶ Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinersts. 1, 81679 München, Germany</footnote>

<footnote>* Corresponding author; luigi.zallio@unimi.it</footnote>
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-\dot{M}_*$ correlation). In the viscous paradigm, the $M_D-\dot{M}_*$ 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-\dot{M}_*$ correlation. In Sect. 2, we show the condition under which it is possible to observe the $M_D-\dot{M}_*$ correlation under the assumption that the initial accretion timescale $t_{acc,0}$<fnref n="1" /> is not correlated with $M_0$. We then analytically derive the slope of the $M_D-\dot{M}_*$ 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_0$ and $t_{acc,0}$: we introduce a correlation in the form $M_0 \propto t_{acc,0}^{\phi}$, and we analytically derive the slope of the $M_D-\dot{M}_*$ correlation in the initial lifetime conditions of the disks and at late times. Then, by adding a spread in $M_0$, 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/\dot{M}_*)$. Finally, we present and investigate the $\phi-\sigma_{M_0}$ degeneracy. In Sect. 5, we summarize our conclusions.

## 2. The $M_D-\dot{M}_*$ correlation in the absence of a $M_0 - 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

$$ \frac{\partial \Sigma}{\partial t} = \frac{3}{r} \frac{\partial}{\partial r} \left\{ \frac{1}{r \Omega} \frac{\partial}{\partial r} (r^2 \alpha_{SS} \Sigma c_s^2) \right\} + \frac{3}{2r} \frac{\partial}{\partial r} \left\{ \frac{\alpha_{DW} \Sigma c_s^2}{\Omega} \right\} - \frac{3 \alpha_{DW} \Sigma c_s^2}{4(\lambda - 1) r^2 \Omega}, \quad (1) $$

where $\alpha_{SS}$ is the Shakura-Sunyaev $\alpha$-parameter, while $\alpha_{DW}$ is its analog for a MHD wind, and $\lambda$ is the magnetic lever arm. In particular, $\alpha_{DW}$ describes the amount of angular momentum extracted by the wind, while $\lambda$ 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 $\alpha_{DW}$ and $\lambda$ to be constant across the disk, and $\lambda$ to be constant in time, it is possible to solve analytically the master equation, obtaining the surface density profile

$$ \Sigma(r, t) = \Sigma_c(t) (r/r_c)^{-1+\xi} e^{-r/r_c}, \quad (2) $$

where $\xi = 1/[2(\lambda - 1)]$, and $r_c$ is the disk characteristic radius. Tabone et al. (2022a) present solutions where $\alpha_{DW}$ increases as the disk mass decreases, namely as

$$ \alpha_{DW}(t) = \alpha_{DW}(0) (\Sigma_c(t)/\Sigma_c(0))^{-\omega}, \quad (3) $$

so that $\alpha_{DW}$ accelerates the evolution of the disk. Here $\omega$ is a phenomenological parameter between 0 and 1 that describes the unknown dissipation of the magnetic field. In particular, $\omega$ describes the tendency of the magnetization to increase as the disk mass decreases. For example, the case $\omega = 1$ mimics a disk where the magnetic field strength is constant over time, while the case $\omega = 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_{acc,0}}\right)^{1/\omega}, \quad (4) $$

and the mass accretion rate $\dot{M}_*$ evolves as

$$ \dot{M}_*(t) = \frac{M_0}{2t_{acc,0}(1 + f_{M,0})} \left(1 - \frac{\omega t}{2t_{acc,0}}\right)^{1/\omega - 1}, \quad (5) $$

where $f_{M,0}$ is the initial mass ejection-to-accretion ratio

$$ f_{M,0} = (r_c/r_{in})^{\xi} - 1. \quad (6) $$

In the end, $M_D(t)$ and $\dot{M}_*(t)$ are controlled by the four independent parameters $M_0$, $f_{M,0}$, $t_{acc,0}$ and $\omega$. In this work, if not otherwise specified we set $f_{M,0} = 2.2$ and $\omega = 0.8$.

In this paper, we test the emergence of the $M_D-\dot{M}_*$ correlation. In principle, the accretion onto the central star $\dot{M}_*$ 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.

<footnote n="1">¹ From Tabone et al. (2022a): the initial accretion timescale is a generalization of the viscous timescale and corresponds to the time that would be required to accrete a fluid particle located initially at $r_c(t = 0)/2$ to the inner region of the disk with an accretion velocity equal to its initial value.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1. Isochrones of a population of disks with $t = 5$ Myr for three values of $M_0$ in the $M_D-\dot{M}_*$ plane. The colored dots are Monte Carlo simulations of $10^5$ proto-planetary disks evolving following Eqs. (4), (5) with $\omega = 0.8$, where $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 = 0.1, 1.0, 10 M_\odot$.</figcaption></figure>

## 2.1. The emergence of the $M_D-\dot{M}_*$ 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-\dot{M}_*$ plane of the sample of disks that have the same initial mass $M_0$, 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 $\dot{M}_*$ are evaluated using Eqs. (4), (5): $t_{acc,0}$ follows a lognormal distribution centered<fnref n="2" /> on the natural logarithm of 1 Myr with a spread of 0.26 dex, while $M_0 = 0.1, 1.0, 10 M_\odot$. 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-\dot{M}_*$ 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-\dot{M}_*$ correlation we set $M_0$ 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-\dot{M}_*$ 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_{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-\dot{M}_*$ correlation (Manara et al. 2016) is a power law whose slope is close to unity. In order to systematically determine the values of $\sigma_{M_0}$ and $\sigma_{t_{acc,0}}$ that lead to a

correlation between $M_D$ and $\dot{M}_*$, we use the coefficient of determination $R^2$. Combining a visual inspection of the $M_D-\dot{M}_*$ plane with the evaluation of $R^2$ from the best power law fit for each disk population, we derived a conservative “rule of thumb”: $R^2 \ge 0.5$. If $R^2 \ge 0.5$, the disk population can be safely fitted in the $M_D-\dot{M}_*$ 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-\dot{M}_*$ plane with a single power law, and that we investigate under which conditions we can derive a linear correlation between $M_D$ and $\dot{M}_*$ (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 \pm 0.12)$, which changes to $r = (0.7 \pm 0.1)$ when the upper limits are considered in the analysis.

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

Tabone et al. (2022b) show that there is no requirement to fine-tune $\omega$ to reproduce the observed accretion rates and disk masses, although explaining the rapidity of disk dispersal differs values of $\omega \sim 0$; in particular, if $0.5 \le \omega < 1$ these latter experience a steeper drop before being dispersed and in this way reproduce the observed distribution. Hence, we chose $\omega = 0.8$, which is bracketed by the $\omega = 0.5$ and $\omega = 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 $\omega = 0.2$ and $\omega = 1$, highlighting how the results of this work also hold for different values of $\omega$.

From Fig. 3, we note that our criterion $R^2 \ge 0.5$ is very roughly described by

$$ \frac{\sigma_{M_0}}{\sigma_{t_{acc,0}}} > 1. \qquad (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-\dot{M}_*$ correlation is $\gamma = 1$ (for the analytical calculations, see Appendix B).

In Fig. 4, we show the $M_D-\dot{M}_*$ correlation in the regime described by Eq. (7). Here, $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered on the natural logarithm of 0.1 $M_\odot$ 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 $\gamma = 1$. It is important to note that 4 shows that if $\sigma_{t_{acc,0}} > \sigma_{M_0}$, we would not be observing a linear correlation between $M_D$ and $\dot{M}_*$; we would instead see the boomerang isochrone shape shown in Fig. 1, which likewise depends on the value of $\omega$.

We note that Eq. (7) is only an approximate criterion. The exact criterion under which the $M_D-\dot{M}_*$ correlation is a power law with a slope of unity is more complex and depends on the age and on the value of $\omega$. In particular, we find that for low values of $\omega$ (see Appendix A) the correlation is harder to get,

<footnote n="2">² We recall that, for a Lognormal distribution, the mean value is different from the central value of the distribution. Therefore, in this work we always indicate the median values of our lognormal distributions and the spreads of their underlying normal distributions.</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 2.** $M_D-\dot{M}_*$ plane for different values of $\sigma_{M_0}$. The colored dots in every plot are Monte Carlo simulations of $10^5$ proto-planetary disks with $\omega = 0.8$ and $t = 5$ Myr, evolving following Eqs. (4), (5), where $M_0$ follows a lognormal distribution centered on the natural logarithm of $0.1 M_\odot$ 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, $t_{acc,0}$ follows a lognormal distribution centered on the natural logarithm of 1 Myr, with a spread of 0.26 dex.

<figure><img src="image_3.png" /></figure>

**Fig. 3.** $\sigma_{M_0} - \sigma_{t_{acc,0}}$ plane with the values of $R^2$ for different values of age $t$. The colored dots in every plot are Monte Carlo simulations of $10^5$ proto-planetary disks with $\omega = 0.8$, evolving following Eqs. (4), (5), where $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered on the natural logarithm of $0.1 M_\odot$ and 1 Myr, respectively. The hatched contours represent the regions where $R^2 = 0$, while the dashed line represents the points where $\sigma_{M_0} = \sigma_{t_{acc,0}}$.
<figure><img src="image_4.png" /><figcaption>Fig. 4. $M_D-\dot{M}_*$ correlation with $\omega = 0.8$ and $t = 2$ Myr. The colored dots are Monte Carlo simulations of $10^5$ proto-planetary disks evolving following Eqs. (4), (5) where $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered, respectively, on 0.1 $M_\odot$ 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 $\gamma = 1$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. $\sigma_{t_{acc,0}}-\sigma_{M_0}$ plane with the values of the slope $\gamma$ and the spread of the $M_D-\dot{M}_*$ correlation $\sigma$. Shown in orange are the values $\gamma \in [0.5, 0.9]$, while in light blue the values $\sigma \in [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 $\gamma$ and $\sigma$ for Monte Carlo disk populations obtained with the corresponding specific values of $\sigma_{t_{acc,0}}$ and $\sigma_{M_0}$. Every Monte Carlo consists of $10^4$ disks ruled by Eqs. (4), (5), where $\omega = 0.8$, $t = 2$ Myr, while $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered on the natural logarithm of 0.1 $M_\odot$ and 2 Myr, respectively. The cross shows the case investigated by Tabone et al. (2022b), while the dashed line highlights the value $\sigma_{t_{acc,0}} = 0.4$.</figcaption></figure>

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_0$, a boomerang-shaped isochrone appears in the $M_D-\dot{M}_*$ plane;

2. In the limit $\sigma_{M_0} \gg \sigma_{t_{acc,0}}$, the slope of the $M_D-\dot{M}_*$ 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-\dot{M}_*$ 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-\dot{M}_*$ correlation and its spread.

## 2.2. The slope and the spread of the $M_D-\dot{M}_*$ correlation as a function of $\sigma_{M_0}$ and $\sigma_{t_{acc,0}}$

In order to constrain the values of $\sigma_{M_0}$ and $\sigma_{t_{acc,0}}$ that can reproduce the observed slope and spread of the $M_D-\dot{M}_*$ correlation, we performed a linear regression of the disks in the $M_D-\dot{M}_*$ plane for different combinations of $\sigma_{M_0}$ and $\sigma_{t_{acc,0}}$. In Fig. 5, we show the results of our analysis. Here, the colored stripes in the $\sigma_{t_{acc,0}} - \sigma_{M_0}$ plane represent the regions of the parameter space that predict values of the slope $\gamma$ and of the spread $\sigma$ of the $M_D-\dot{M}_*$ 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 $\gamma \in [0.5, 0.9]$, while in light blue the region $\sigma \in [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/\dot{M}_*$. In this plot, $t = 2$ Myr, $\omega = 0.8$, and $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered on the natural logarithm of 0.1 $M_\odot$ 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-\dot{M}_*$ correlation for $\sigma_{M_0} \in [0.50, 1.50]$ and $\sigma_{t_{acc,0}} \in [0.40, 0.75]$.

We note that the results obtained here are valid for a fixed age $t$, a value of $\omega = 0.8$, and the fixed median values $M_0 = 0.1 M_\odot$ and $t_{acc,0} = 2$ Myr. As we show in Appendix A, changing $\omega$ will change the derived values of $\sigma_{M_0}$ and $\sigma_{t_{acc,0}}$, but it remains true that it is possible to retrieve the observed slope and the spread of the $M_D-\dot{M}_*$ 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_{acc,0}} \approx 0.4$ dex needs to be adopted to fit the observed decline of disk fraction with cluster age.

## 2.3. The spread of $\log(M_D-\dot{M}_*)$ 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/\dot{M}_*$ is a fundamental observable for disk population studies. Considering Eqs. (4), (5), we can evaluate the dimensionless spread of the accretion timescale when $\dot{M}_* \propto M_D$. Using the propagation of errors (see more details in Appendix C) on $\log(M_D/\dot{M}_*)$, the dimensionless spread is

$$ \sigma(\log(M_D/\dot{M}_*)) \approx \frac{1}{2\langle t_{acc,0} \rangle - \omega\langle t \rangle} \quad (8) \\ \sqrt{(2\langle t_{acc,0} \rangle \sigma_{t_{acc,0}})^2 + (\omega\langle t \rangle \sigma_t)^2}, $$
<figure><img src="image_6.png" /><figcaption>Fig. 6. Spread of the $M_D-\dot{M}_*$ correlation as a function of $\sigma_{t_{acc}}$, with $\sigma_t = 0.43$ dex. The values of $\sigma(\log(M_D/\dot{M}_*))$ are evaluated as the standard deviation from Monte Carlo simulations of $10^5$ protoplanetary disks with $M_0$, $t_{acc}$, and $t$ following a lognormal distribution centered respectively in 0.1 $M_\odot$, 5 Myr, and 1 Myr.</figcaption></figure>

where we use the notation $\langle t_{acc,0} \rangle$ and $\langle t \rangle$ to remark this is the median value for the population, and we also allow for the possibility of an age spread $\sigma_t$ in the population. This result is equivalent to the spread of the $M_D-\dot{M}_*$ correlation for a slope mildly different than 1, but we note that if the $M_D-\dot{M}_*$ correlation differs substantially from linearity, then Eq. (8) is not valid anymore.

The fundamental assumption underlying Eq. (8) is that $\sigma_{t_{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_{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/\dot{M}_*)$ is equal to that of the $M_D-\dot{M}_*$ 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-\dot{M}_*$ correlation, as further discussed in Sect. 4.1.

In Fig. 6, we show the spread of $\log(M_D/\dot{M}_*)$ when $t$ follows a lognormal distribution centered on the natural logarithm of 1 Myr, and $\sigma_t = 0.43$ dex. The black line is computed as the standard deviation from the Monte Carlo simulations of $10^5$ disks that follow a lognormal distribution in $t$, $t_{acc}$ and $M_0$. In particular, the value of $\sigma(\log(M_D/\dot{M}_*))$ is evaluated for $10^2$ values of $\sigma_{t_{acc}}$ that span [0, 1.1] dex, while $\sigma_t = 0$ dex. Here, $M_0$ 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_\odot$ for $M_0$. The spread of $\log(M_D/\dot{M}_*)$ is independent of the spread in $M_0$ under the conditions for which Eq. (8) holds, which is that the slope of the $M_D-\dot{M}_*$ correlation is close to unity. In Fig. 5 we see, as we get close to $\gamma = 1$, that the cyan band representing the spread of the $M_D-\dot{M}_*$ correlation becomes vertical, while as we move away from those values the dependence on $\sigma_{t_{acc}}$ 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_0$ 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-\dot{M}_*$ 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-\dot{M}_*$ correlation. Finally, we provided an approximated analytical expression for the spread on the accretion timescale $\log(M_D/\dot{M}_*)$. In the following, we study the emergence of the $M_D-\dot{M}_*$ correlation when $M_0$ is correlated to $t_{acc,0}$.

### 3. The $M_D-\dot{M}_*$ correlation in the presence of a $M_0-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_0$. However, this is not necessarily the case. In the following, we assume a correlation between $M_0$ and $t_{acc,0}$ in the form

$$M_0 \propto t_{acc,0}^{\phi} \qquad (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_{acc,0}}{\omega} \qquad (10)$$

as the disk dispersal time.

#### 3.1. The $M_D-\dot{M}_*$ correlation at initial time

When taking a population of disks without spread in $M_0$, in the initial conditions (i.e., $t/t_{disp} \ll 1$) we find that the slope is $\gamma = 1 - 1/\phi$ (see Appendix D.1 for more details).

In Fig. 7, we show the $M_D-\dot{M}_*$ correlation for a disk population (dots) of $10^5$ disks ruled from Eqs. (4), (5), where $M_0$ is given by Eq. (9). Here, $\phi = 0.8$, $\omega = 0.8$, $t = 0$ yr, and $t_{acc,0}$ follows a lognormal distribution centered on the natural logarithm of 2 Myr with a spread of 0.13 dex. The black line shows the analytical $M_D-\dot{M}_*$ 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_0$, close to dispersal (i.e., $t/t_{disp} \sim 1$) we find that $\gamma = 1 - \omega$ (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_0$ 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_0$ (as they are linked to each other). Thus, when $t/t_{disp} \sim 1$, the correlation between $M_0$ 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
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Slope of the $M_D-\dot{M}_*$ correlation with $\phi = 0.8$, $\omega = 0.8$, and $t = 0$ yr. The dots are Monte Carlo simulations of $10^5$ proto-planetary disks ruled from Eqs. (4) and (5) where $M_0$ follows Eq. (9) and $t_{acc,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).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Slope of the $M_D-\dot{M}_*$ correlation with $\phi = 0.8$, $\omega = 0.8$, and $t = 3$ Myr. The colored dots are Monte Carlo simulations of $10^5$ proto-planetary disks ruled from Eqs. (4) and (5), where $M_0$ follows Eq. (9) and $t_{acc,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).</figcaption></figure>

any possibility of detecting the initial correlation between $M_0$ and $t_{acc,0}$.

In Fig. 8, we show the $M_D-\dot{M}_*$ correlation for $\phi = 0.8$, $\omega = 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-\dot{M}_*$ 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 \lesssim 10^{-2} M_\odot$ (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-\dot{M}_*$ correlation in the above-cited region of the $M_D-\dot{M}_*$ plane, showing that when $t \sim t_{disp}$, the slope approaches the value $\gamma = 1 - \omega$. The right part of Fig. 8 instead shows the disks that are still evolving, and hence still near the condition $t/t_{disp} \ll 1$; we see that the slope here inverts

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** $M_D-\dot{M}_*$ correlation with $\phi = 1$, $\omega = 0.8$, and $t = 3.5$ Myr. The dots are Monte Carlo simulations of $10^5$ proto-planetary disks evolving following Eqs. (4) and (5), where $M_0$ follows Eq. (9) and $t_{acc,0}$ follows a lognormal distribution centered on the natural logarithm of 5 Myr with a spread of 0.13 dex. In addition, $M_0$ is then generated following a lognormal distribution with the values of $t_{acc,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}$.</figcaption></figure>

its trend, in agreement with the value $1 - 1/\phi$ that we derived in Sect. 3.1.

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

### 3.3. The role of $\sigma_{M_0}$ in the $M_D-\dot{M}_*$ correlation

If $t_{acc,0}$ is generated following a lognormal distribution with a spread and $M_0 \propto t_{acc,0}^\phi$ is satisfied, it is reasonable to think that $M_0$ should also have a spread of its own (i.e., the spread in $M_0$ 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_0$, we initially evaluated $M_0$ with Eq. (9); then $M_0$ was extracted with a lognormal distribution, with fixed spread, centered on the previously determined values of $M_0$. In Fig. 9, once the spread in $M_0$ is added, we see the appearance of the linear $M_D-\dot{M}_*$ relation. Here we show a population of disks identical to the one in Fig. 8, but we add a spread in $M_0$ 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-\dot{M}_*$ correlation with a value different than unity after the introduction of the correlation $M_0 \propto t_{acc,0}^\phi$. It is also possible to tune the value of the $M_D-\dot{M}_*$ slope, changing the values of $\phi$ 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_0$.

### 3.4. Describing the $M_D-\dot{M}_*$ correlation by a single power law

It is relevant to emphasize that, after the introduction of a spread in $M_0$, the bent shape of the $M_D-\dot{M}_*$ correlation shown in Fig. 8 can still be recovered, but only for very low values of $\sigma_{M_0}$.<fnref n="3" /> Fig- ure 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{cases} \dot{M}_* = k_1 + a M_D^{-b}, & M_D < x_0 \\ \dot{M}_* = k_2 + c (M_D - x_0)^d, & M_D \ge x_0 \end{cases} \qquad (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_0$. 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 \qquad (12) $$

to the $\dot{M}_* - M_D$ plane for different values of $\sigma_{M_0}$, namely $\sigma_{M_0} \in [0, 0.022]$ 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_0, k_1, a, b, c, d$), while Eq. (12) on only two ($a, \zeta$), 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 [0, 0.01]$ dex, hence showing that for very low values of $\sigma_{M_0}$, the $M_D-\dot{M}_*$ correlation is better fitted with a double power law. However, as little spread in $M_0$ 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<fnref n="4" />.

From this we can conclude that when a little spread in $M_0$ is added, the $M_D-\dot{M}_*$ 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-\dot{M}_*$ correlation. We show the same plot of Fig. 10 in Appendix F for $\sigma_{M_0} \in [0, 0.43]$ 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-\dot{M}_*)$

To assess whether MHD winds can reproduce both disk lifetimes and the observed spread in $\log(M_D/\dot{M}_*)$, we conducted the following analysis.

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** BIC as function of $\sigma_{M_0}$ in a range $\sigma_{M_0} \in [0, 0.022]$ dex. The BIC is obtained for Eq. (11) (in red) and Eq. (12) (in black) from the fit to the $\dot{M}_*-M_D$ plane obtained for different values of $\sigma_{M_0}$, while the simulated disk populations have the same set of parameters described in Fig. 8.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**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 $t_{disp}$. From the fit, we obtain $t_{disp} = 2.37$ Myr and $\sigma_{t_{disp}} = 0.37$ dex.</figcaption></figure>

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_{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. \qquad (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 $t_{disp}$ follows a lognormal distribution with a median value

<footnote n="3">³ Some tests were conducted on this, but we found it redundant to show the resulting plots. In particular, the bent shape shown in Fig. 8 can be recovered for $\sigma_{M_0} < 0.01$ dex. In the following, we show that as a minimum amount of $\sigma_{M_0}$ is added, the $M_D-\dot{M}_*$ correlation is well described by a single power law.</footnote>

<footnote n="4">⁴ A double power law can always be re-conducted to a single power law: the fact that the curves perfectly superpose means that the double power law became a single power law.</footnote>
$\mu_{t_{disp}} = 2.37$ Myr and a spread $\sigma_{t_{disp}} = 0.37$ dex. After some simple algebra deriving from Eq. (10), we find that the median value of the accretion timescale is $\mu_{t_{acc,0}} = 0.95$ Myr and the spread is $\sigma_{t_{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 $\sigma_t = 0^\circ$ dex, obtaining $\sigma(\log(M_D/\dot{M}_*)) = 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_\odot$, in agreement with Testi et al. (2022). For L1688 and Lupus, we obtained values of $\sigma(\log(M_D/\dot{M}_*))$ that are closely comparable with that reported above, more specifically $\sigma(\log(M_D/\dot{M}_*))_{L1688} = 0.86$ dex, $\sigma(\log(M_D/\dot{M}_*))_{Lupus} = 1.08$ dex. As we moved to the older region of ChaI, we obtained the higher spread $\sigma(\log(M_D/\dot{M}_*))_{ChaI} = 1.62$ dex, a value that increases even more for USco, the oldest region of the sample, for which we obtained $\sigma(\log(M_D/\dot{M}_*))_{USco} = 2.02$ dex. The reason behind these higher spreads is that to evaluate $\sigma(\log(M_D/\dot{M}_*))$, we assumed that the distribution of $\log(M_D/\dot{M}_*)$ 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/\dot{M}_*)$ 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 $\sigma_t = 0$ dex, and we neglected accretion variability and mass uncertainty, so the value of $\sigma(\log(M_D/\dot{M}_*))$ that we predicted must be considered as a lower limit.

## 4.2. Empirical constraints on $\sigma_{M_0}$ and $\phi$

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

In Fig. 12, we show the $\phi-\sigma_{M_0}$ plane with the corresponding values of $\gamma$. The coloured regions in the $\phi-\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 $\phi$ and $\sigma_{M_0}$. In a similar fashion to the analysis performed in Lodato et al. (2017), we show in orange the region $\gamma \in [0.5, 0.9]$, while in light blue the region $\sigma \in [0.55, 0.71]$ dex. In this case $t = 2$ Myr, $\omega = 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 $\gamma$ tends toward unity as $\sigma_{M_0}$ increases, regardless of the value of $\phi$. However, there is a

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** $\phi-\sigma_{M_0}$ plane with the values of the slope $\gamma$ and the spread of the $M_D-\dot{M}_*$ correlation $\sigma$. Shown in orange are the values $\gamma \in [0.5, 0.9]$, while in light blue the values $\sigma \in [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 $\gamma$ and $\sigma$ for Monte Carlo disk populations obtained with the corresponding specific values of $\phi$ and $\sigma_{M_0}$. Every Monte Carlo consists of $10^4$ disks evolving according to Eqs. (4), (5), where $M_0 \propto t_{acc,0}^\phi$, $\omega = 0.8$, $t = 2$ Myr, and $t_{acc,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).</figcaption></figure>

whole region of the $\sigma_{M_0} - \phi$ plane that allows values of $\gamma \neq 1$; this confirms the existence of a degeneracy between $\phi$ 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 $\gamma$ and $\sigma$; we see that only values of $\phi \in [-2, 1]$ and $\sigma_{M_0} \in [0, 1.75]$ dex are allowed. Figure 12 leads us to understand that for higher values of $\phi$ we need higher values of $\sigma_{M_0}$ in order to retrieve the observed $M_D-\dot{M}_*$ correlation. This can be explained as follows: the stronger the correlation between $M_0$ and $t_{acc,0}$ is, the bigger the spread in $M_0$ 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 $\sigma_{M_0}$ 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_0$ 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_0$ and $t_{acc,0}$.

From Fig. 12, we understand that we can obtain values of the slope of the $M_D-\dot{M}_*$ 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-\dot{M}_*$ correlation in the MHD wind-driven scenario under a broad range of conditions, a result that shows that the $M_D-\dot{M}_*$ correlation can also be easily obtained for this evolutionary paradigm.

# 5. Conclusions

In this paper, we investigated under which conditions the $M_D-\dot{M}_*$ correlation emerges in the MHD wind-driven evolutionary scenario.

<footnote n="5">⁵ For the purposes of the current work, we assume that there is no spread in stellar ages.</footnote>
We derived and analyzed the parameters that rule the appear-
ing and tuning of the $M_D-\dot{M}_*$ 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_{acc,0}}$, we obtain the correlation between the mass of the disk $M_D$ and the accretion onto the central star $\dot{M}_*$. This holds as long as $\sigma_{t_{acc,0}} < \sigma_{M_0}$, a conservative criterion that is compatible with the general rule of thumb $R^2 \ge 0.5$. In the opposite case, the $M_D-\dot{M}_*$ 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-\dot{M}_*$ correlation comparable to that found in Manara et al. (2016) varying $\sigma_{M_0}$ and $\sigma_{t_{acc,0}}$. Hence, it is possible to derive the observed slope and spread of the $M_D-\dot{M}_*$ 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/\dot{M}_*)$. The predictions arising from this formula are in agreement with the observed order of magnitude of $\sigma(\log(M_D/\dot{M}_*))$ evaluated for young populations of disks (Testi et al. 2022);

(iv) We introduced a correlation $M_0 \propto t_{acc,0}^{\phi}$ to understand if we can still derive the observed slope and spread of the $M_D-\dot{M}_*$ correlation. We analytically derived a slope $\gamma = 1 - 1/\phi$ in the initial conditions, and a slope $\gamma = 1 - \omega$ close to dispersal.

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

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

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-\dot{M}_*$ 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-\dot{M}_*$ 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-\dot{M}_*$ 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 Skłodowska-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.

References

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20

Andrews, S. M. 2020, ARA&A, 58, 483

Ansdell, M. 2020, in Five Years After HL Tau: A New Era in Planet Formation (HLTAU2020), 29

Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46

Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, VizieR Online Data Catalog: ALMA survey of protoplanetary disks in sigma Ori (Ansdell+, 2017), VizieR On-line Data Catalog: J/AJ/153/240

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Bai, X.-N. 2016, ApJ, 821, 80

Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142

Barenfeld, S. A., Carpenter, J. M., Sargent, A. I., Isella, A., & Ricci, L. 2017, ApJ, 851, 85

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Carr, J. S., Tokunaga, A. T., & Najita, J. 2004, ApJ, 603, 213

Cazzoletti, P., Manara, C. F., Liu, H. B., et al. 2019, A&A, 626, A11

Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698

Cox, E. G., Harris, R. J., Looney, L. W., et al. 2017, ApJ, 851, 83

Eisner, J. A., Arce, H. G., Ballering, N. P., et al. 2018, ApJ, 860, 77

Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72

Ferreira, J. 1997, A&A, 319, 340

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, 3rd edn.

Hartmann, L., Hinkle, K., & Calvet, N. 2004, ApJ, 609, 906

Ilee, J. D., Fairlamb, J., Oudmaijer, R. D., et al. 2014, MNRAS, 445, 3723

Lesur, G. 2021, J. Plasma Phys., 87, 205870101

Lodato, G., Scardoni, C. E., Manara, C. F., & Testi, L. 2017, MNRAS, 472, 4700

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Manara, C. F., Testi, L., Natta, A., & Alcalá, J. M. 2015, A&A, 579, A66

Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3

Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127

Manara, C. F., Natta, A., Rosotti, G. P., et al. 2020, A&A, 639, A58

Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in Astronomical Society of the Pacific Conference Series, 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 539

Mann, R. K., Di Francesco, J., Johnstone, D., et al. 2014, ApJ, 784, 82

Morbidelli, A., & Raymond, S. N. 2016, J. Geophys. Res. (Planets), 121, 1962

Najita, J., Carr, J. S., Glassgold, A. E., Shu, F. H., & Tokunaga, A. T. 1996, ApJ, 462, 919

Najita, J. R., Doppmann, G. W., Carr, J. S., Graham, J. R., & Eisner, J. A. 2009, ApJ, 691, 738

Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125

Pringle, J. E. 1981, ARA&A, 19, 137

Rosotti, G. P., Clarke, C. J., Manara, C. F., & Facchini, S. 2017, MNRAS, 468, 1631

Sanchis, E., Testi, L., Natta, A., et al. 2021, A&A, 649, A19

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Somigliana, A., Testi, L., Rosotti, G., et al. 2023, ApJ, 954, L13

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74

Tabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022a, MNRAS, 512, 2290

Tabone, B., Rosotti, G. P., Lodato, G., et al. 2022b, MNRAS, 512, L74

Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98
## Appendix A: The $\sigma_{M_0}-\sigma_{t_{acc,0}} - R^2$ plane for different values of $\omega$

In Figs. A.1, A.2 we show the $\sigma_{M_0} - \sigma_{t_{acc,0}}$ plane for $\omega = 0.2$ and $\omega = 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 $\omega$. Moreover, we note that, in Fig. A.1 for evolved populations our general criterion $R^2 \ge 0.5$ does not correspond anymore necessarily to the region $\sigma_{M_0} > \sigma_{t_{acc,0}}$. The reason for this is shown in Fig. A.3: for evolved populations, when $\omega \to 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-\dot{M}_*$ plane.

## Appendix B: The slope of the $M_D-\dot{M}_*$ 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_{acc,0}(1 + f_{M,0})} \left(1 - \frac{\omega t}{2t_{acc,0}}\right)^{-1} \quad (B.1) $$

Taking the hint from the work of Manara et al. 2016, it is possible to investigate the analytical behavior of disks evolving according to Eq. (B.1), assuming that

$$ \dot{M}_* \propto M_D^{\gamma} \quad (B.2) $$

Performing the logarithmic derivative on Eq. (B.2), the slope $\gamma$ results in the constant

$$ \gamma = \frac{M_D}{\dot{M}_*} \frac{d\dot{M}_*}{dM_D} = 1. \quad (B.3) $$

## Appendix C: The spread of the $M_D-\dot{M}_*$ correlation as function of the spread in $t$ and $t_{acc}$: analytics

Starting from Eqs. (4), (5), the ratio $M_D/\dot{M}_*$ reads

$$ \frac{M_D}{\dot{M}_*} = (1 + f_{M,0}) (2t_{acc,0} - \omega t). \quad (C.1) $$

Neglecting the covariance term between $t$ and $t_{acc,0}$, the spread of $\log(M_D/\dot{M}_*)$ correlation is obtained with a simple propagation of errors as

$$ \sigma(\log(M_D/\dot{M}_*)) = \left[ \sigma_{t_{acc,0}}^2 \left( \frac{\partial(\log(M_D/\dot{M}_*))}{\partial \log(t_{acc,0})} \right)^2 + \sigma_t^2 \left( \frac{\partial(\log(M_D/\dot{M}_*))}{\partial \log(t)} \right)^2 \right]^{1/2}, \quad (C.2) $$

which gives Eq. (8) after some simple algebra. We note that the quantity $\sigma_{t_{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-\dot{M}_*$ correlation when $M_0 \propto t_{acc,0}^\phi$: analytics

### D.1. The slope for $t/t_{disp} \ll 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} \ll 1$, $t_{acc,0} \propto M_D^{1/\phi}$. (D.1)

Including this result in Eq. (5), after the trivial substitution of Eq. (9) in Eq. (5), we obtain

$$ \dot{M}_* \propto M_D^{1-1/\phi}. \quad (D.2) $$

The logarithmic derivative gives the slope in the form

$$ \gamma = \frac{M_D}{\dot{M}_*} \frac{d\dot{M}_*}{dM_D} = 1 - \frac{1}{\phi}. \quad (D.3) $$

### D.2. The slope for $t \sim t_{disp}$

We derive the slope in the final life-time conditions of a disk starting by choosing $\epsilon$ to be arbitrarily small. Then, inserting

$$ \frac{t}{t_{disp}} = 1 - \epsilon \quad (D.4) $$

in Eqs. (4), (5) after the previous substitution of Eq. (9), $\epsilon$ is given by

$$ \epsilon \propto \left( \frac{M_D}{t_{acc,0}^\phi} \right)^\omega. \quad (D.5) $$

After some algebra, it is obtained that

$$ \dot{M}_* \propto t_{acc,0}^{\phi\omega-1} M_D^{1-\omega}, \quad (D.6) $$

which leads to a slope of the form

$$ \gamma = \frac{M_D}{\dot{M}_*} \frac{d\dot{M}_*}{dM_D} = 1 - \omega. \quad (D.7) $$

## Appendix E: The spread of the $M_D-\dot{M}_*$ correlation when $M_0 \propto t_{acc,0}^\phi$

### E.1. The spread when $M_0 \propto t_{acc,0}^\phi$ and $t \ll t_{disp}$

Following Sect. 2.3, the calculations for the spread when $M_0 \propto t_{acc,0}^\phi$ lead to

$$ \sigma(\log(M_D/\dot{M}_*)) = \sigma_{t_{acc}}. \quad (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_{acc,0}^\phi$ and $t \sim t_{disp}$

A similar result is obtained when $t \sim t_{disp}$, obtaining

$$ \sigma(\log(M_D/\dot{M}_*)) = \epsilon \sigma_{t_{acc}}. \quad (E.2) $$

It is clear that the more $t$ approaches the value of $t_{disp}$, the latter is the value of $\epsilon$. In the limit for $\epsilon \to 0$, Eq. (E.2) gives $\sigma(\log(M_D/\dot{M}_*)) \to 0$: a valid result that reflect the fact that when the disk is about to die, the slope of the $M_D-\dot{M}_*$ correlation inevitably tends toward the value $1-\omega$.
<figure><img src="image_13.png" /></figure>

**Fig. A.1.** $\sigma_{M_0} - \sigma_{t_{acc,0}}$ plane with the values of $R^2$ for different values of age $t$. The colored dots in every plot are Monte Carlo simulations of $10^5$ proto-planetary disks with $\omega = 0.2$, evolving following Eqs. (4), (5), where $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered on the natural logarithm of 0.1 $M_\odot$ and 1 Myr, respectively. The hatched contours represent the regions where there is no correlation between $M_D$ and $\dot{M}_*$, while the dashed line represents the points where $\sigma_{M_0} = \sigma_{t_{acc,0}}$.

**Fig. A.2.** $\sigma_{M_0} - \sigma_{t_{acc,0}}$ plane with the values of $R^2$ for different values of age $t$. The colored dots in every plot are Monte Carlo simulations of $10^5$ proto-planetary disks with $\omega = 1.0$, evolving following Eqs. (4), (5), where $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered on the natural logarithm of 0.1 $M_\odot$ and 1 Myr, respectively. The hatched contours represent the regions where there is no correlation between $M_D$ and $\dot{M}_*$, while the dashed line represents the points where $\sigma_{M_0} = \sigma_{t_{acc,0}}$.
<figure><img src="image_14.png" /></figure>

**Fig. A.3.** $M_D-\dot{M}_*$ correlation with $\omega = 0.2$ and $t = 2.5$ Myr. The colored dots are Monte Carlo simulations of $10^5$ proto-planetary disks evolving following Eqs. (4), (5) where $M_0$ and $t_{acc,0}$ follow a lognormal distribution centered, respectively, on the natural logarithm of $0.1 M_\odot$ 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 $\gamma = 0.72$.

<figure><img src="image_15.png" /></figure>

**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 $\dot{M}_* - M_D$ 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 $\dot{M}_* - M_D$ loses the double power law behavior in favor of a single power law written in the form of Eq. (12).