The time evolution of $M_{\mathrm{d}}/\dot M$ in protoplanetary discs as a way to disentangle between viscosity and MHD winds

As the classic viscous paradigm for protoplanetary disk accretion is challenged by the observational evidence of low turbulence, the alternative scenario of MHD disk winds is being explored as potentially able to reproduce the same observed features traditionally explained with viscosity. Although the two models lead to different disk properties, none of them has been ruled out by observations - mainly due to instrumental limitations. In this work, we present a viable method to distinguish between the viscous and MHD framework based on the different evolution of the distribution in the disk mass ($M_{\mathrm{d}}$) - accretion rate ($\dot M$) plane of a disk population. With a synergy of analytical calculations and 1D numerical simulations, performed with the population synthesis code \texttt{Diskpop}, we find that both mechanisms predict the spread of the observed ratio $M_{\mathrm{d}}/\dot M$ in a disk population to decrease over time; however, this effect is much less pronounced in MHD-dominated populations as compared to purely viscous populations. Furthermore, we demonstrate that this difference is detectable with the current observational facilities: we show that convolving the intrinsic spread with the observational uncertainties does not affect our result, as the observed spread in the MHD case remains significantly larger than in the viscous scenario. While the most recent data available show a better agreement with the wind model, ongoing and future efforts to obtain direct gas mass measurements with ALMA and ngVLA will cause a reassessment of this comparison in the near future.


INTRODUCTION
The gaseous component of protoplanetary disks has traditionally been described as undergoing viscous accretion (Lynden-Bell & Pringle 1974, Pringle 1981. In recent years, however, a growing observational evidence is challenging this picture, as the low levels of turbulence detected in protoplanetary disks appear incompatible with the observed evolution (Pinte et al. 2016, Flaherty et al. 2018, Rosotti 2023. The best alternative to the classic viscous scenario is currently provided by MHD disk winds, originally proposed by Blandford & Payne (1982). This model has gained increasing popularity in the recent years, as several studies (see Lesur 2021 for a review) have shown it to reproduce the key evolutionary features of protoplanetary disks; moreover, Tabone et al. (2022a) have developed a simple analytical parametrization, making it a valid alternative to the viscous theory.
A compelling question is which of these mechanism, or which combination of the two, drives angular momentum transport in protoplanetary disks (Manara et al. 2023). Answering this question has proven to be a surprisingly difficult task: even though the two models are in principle well distinguishable through their characteristic theoretical predictions, the observational counterpart is lagging behind (e.g., Rosotti et al. 2019b, Ilee et al. 2022. A good example of this problem is viscous spreading, a fundamental feature of viscous evolution that causes the gaseous component of disks to expand in radius as they evolve. As MHD evolution does not arXiv:2309.04496v1 [astro-ph.EP] 7 Sep 2023 show a similar behavior (Zagaria et al. 2022b), it would in principle be a good candidate for disentangling between the two predictions: however, the high sensitivity required to detect it has until now represented a limit. While Class 0 objects are widely accepted to be born small (< 60 au: Maury et al. 2019, also supported by the numerical experiments of, e.g., Lebreuilly et al. 2021) and grow wider in the first 1-2 Myr of evolution (Najita & Bergin 2018), whether the radius of Class II disks increases or decreases with time is widely debated. Dust continuum radii are observed to be shrinking with time (Hendler et al. 2020, Zagaria et al. 2022b, as an effect of radial drift, while gas observations (Ansdell et al. 2018, Sanchis et al. 2021, Toci et al. 2021, Long et al. 2022 have covered too small of a sample at too low sensitivities to draw firm conclusions. The advent of ALMA Band 1 (Carpenter et al. 2020) and the next-generation VLA (ngVLA, Tobin et al. 2018) in the near future will allow to perform surveys of protoplanetary disks at unprecedentedly long wavelengths, which will play a crucial role in determining the leading evolutionary mechanism. At the same time, finding novel approaches to tackle this problem is crucial to obtain significant results.
In this Letter, we suggest a new method to distinguish between the two models from the population perspective: through a joint theoretical and population synthesis approach, we investigate the time evolution of disks in the disk mass -accretion rate plane, proving it to be a good approach for our goal. This work is structured as follows: in Section 2 we describe the evolutionary prescriptions that we adopt and we discuss their numerical implementation. In Section 3 we present our results and we compare them with the observations. Finally, in Section 4 we discuss the implications of these results and draw our conclusions.

Secular evolution
The simulations presented in this work have been carried out using the 1D Python population synthesis code Diskpop. For a detailed description of the code, as well as its public release, we refer to our upcoming paper (Somigliana et al. in prep; earlier implementations of the code, its basic assumptions and features have been described in Rosotti et al. 2019a, Rosotti et al. 2019b, Toci et al. 2021, Somigliana et al. 2022. The viscous and MHD evolution are implemented following Lynden-Bell & Pringle (1974) and Tabone et al. (2022a) respectively. In this section we briefly present both models, referring to the original papers for a deeper discussion.
In the viscous case, we solve the classic evolution equation following the prescription by Shakura & Sunyaev (1973), the viscosity ν is modeled as α SS c s H, where α SS is a dimensionless parameter, c s is the sound speed, and H is the height of the disk. Furthermore, assuming the viscosity to be a power-law of the disk radius for ease of solving the equation, and R c is a scale radius), the analytical solution by Lynden-Bell & Pringle (1974) holds.
In the MHD case instead (Tabone et al. 2022a), the evolution equation is given by where Ω is the keplerian orbital frequency, λ is the magnetic lever arm parameter, and α DW is a magnetic equivalent of α SS . Equation (2) is a generalization of Equation (1) if the gas surface density evolves not only because of the viscous torque (first term on the RHS) but also because of the effects of MHD disk winds, which extract angular momentum and induce a mass loss (second and third term on the RHS respectively). Assuming that both λ and α DW are constant across the disk, and that α DW ∝ Σ c −ω (where Σ c = Σ(R = R c )), Equation (2) can be solved analytically (see Tabone et al. 2022a).

Isochrones
Isochrones are defined as the curves described by a population of objects of the same age in a given plane. In the case of protoplanetary disks, isochrones in the M d −Ṁ plane have been the focus of recent studies (Lodato et al. 2017, Somigliana et al. 2020. For viscously evolving disks (Lodato et al. 2017), the isochrone readsṀ the only free parameter in Equation (3) is the initial disk mass M 0 , which only sets the starting point of the isochrone. Nonetheless, at late stages (when M d ≪ M 0 ) all disks in a population are bound to reach the same locus on the M d −Ṁ plane: while this happens at different ages for each disk, depending on its viscous timescale t ν = R c 2 /(3(2 − γ) 2 ν c ), a fully evolved population (t → +∞) will necessarily sit on the theoretical isochrone of the corresponding age.
For MHD disks, the isochrone is defined as (Tabone et al. 2022a) Equation (4) depends not only on M 0 , but also on the equivalent of t ν in the MHD winds case, the initial accretion timescale t acc,0 , through f M,0 (determined by the disk radius -see Tabone et al. 2022a for details).
The interpretation of the isochrones in the two models is therefore different: while the viscous curves for all disks in a population lie on top of each other (except at the early stages, when M d ∼ M 0 ), MHD evolution never loses memory of the initial conditions. This is because, depending on whether we fix M 0 or t acc,0 , we can define two types of isochrones for an MHD population. As a result, disks with a different M 0 will occupy an area of the M d −Ṁ plane, rather than sitting on a single curve, and this will be the case even for evolved populationswhich means that it is not possible to use the isochrones to obtain age estimates for disk populations. Based on this argument, we investigate whether the evolution of a population of disks in the M d −Ṁ plane could carry tangible signatures of the evolutionary model.

Population synthesis
In this work we adopt a population synthesis approach, which consists of generating and evolving a synthetic population of protoplanetary disks via numerical methods. We employ the Python tool Diskpop, which we expanded from our previous work (Somigliana et al. 2022) to include MHD disk wind evolution. In this section, we present a brief outline of the workflow, referring to the upcoming code release for a detailed description of the methods and the implementation.
First, we generate N ∼ 100 stars, whose masses M ⋆ follow the Kroupa (2001) initial mass function. We then assemble a Young Stellar Object (YSO) by assigning a disk to each star: to determine the initial mass and radius of said disk, we assume that the initial disk mass and accretion rate scale as power-laws of the stellar mass (M d ∝ M ⋆ λm andṀ ∝ M ⋆ λacc ). In our previous work (Somigliana et al. 2022) we have demonstrated how λ m,0 ∈ [0.7, 1.5] and λ acc,0 ∈ [1.2, 2.1] can reproduce the slopes of observed correlations of disk properties with stellar mass at later ages; we refer to that paper for a detailed discussion. We determine M d andṀ for each disk drawing from a log-normal distribution, centered in the mean value computed via the power-law correla-tions and with a width (σ) of choice; R d is then derived from considerations onṀ (see Somigliana et al. 2022 for details). The other relevant quantities besides M ⋆ , M d and R d are fixed in our model: Table 1 shows the parameters that we used in the simulations presented in this work, based on the disc evolution studies of Lodato et al. (2017) and Tabone et al. (2022a) for viscosity and MHD winds respectively. While a detailed study of the parameters space is outside of the scope of this work, we have tested two more combinations of parameters (shown by Tabone et al. (2022b) to reproduce the Lupus star-forming region) and we found that our results are independent on the particular combination chosen. Once the population of YSOs is generated, it is evolved following the viscous or MHD prescription via a 1D implementation of the models described in Section 2.1. Although Diskpop allows to numerically solve the evolution equations, in this work we have used the analytical solutions to Equation (1) and (2); it is therefore important to note that our results depend on the assumptions needed to obtain such solutions (e.g., the power-law scaling of viscosity with the disk radius).
It is crucial to point out that disk dispersal is an intrinsic feature of MHD winds, but not of viscous evolution. Our code includes an observational effect by considering as dispersed disks with masses lower than 10 −6 M ⊙ ; this simulates a dispersal effect even in the viscous scenario, which would otherwise generate disks with infinite lifetime, that do not match the observed disk fraction (see Appendix C). This problem is usually solved in the literature by adding other physical effects to the purely viscous model, such as internal photoevaporation (see e.g. Hollenbach et al. 1994, Clarke et al. 2001, Owen et al. 2011, Picogna et al. 2019, Emsenhuber et al. 2023. In order to account for the statistical effect of reducing our sample throughout the evolution caused by disk dispersal, we performed 100 simulations for both setups described in Table 1 and then considered not only the median evolution of the interesting quantities, but also the interval between the 25th and 75th percentile (see Section 3).

RESULTS
In this Section, we show the results of the evolution of viscous and MHD populations of protoplanetary disks in the M d −Ṁ plane: in particular, we consider the ratio of the two quantities (hereafter t lt , disk lifetime -see Jones et al. 2012). We first discuss the expected evolution of the distribution of disk lifetimes from an analytical point of view (paragraph 3.1), and then we confirm our theoretical results through Diskpop simu-  (2)  lations (paragraph 3.2); finally, we compare our results with the observations (paragraph 3.3).

Disk lifetimes distribution
In the traditional viscous picture (Dullemond et al. 2006, Lodato et al. 2017, disks lie on the theoretical isochrone (Equation 3) at a given age t if their initial viscous timescale t ν,0 is much shorter than t; as evolution proceeds, more and more disks reach this stage and therefore the population converges around the cor-responding isochrone. As a consequence, the spread around the isochrones decreases with time: eventually, once the population is fully self-similar (i.e., its age is larger than all of the viscous timescales), the spread will be vanishingly small and the correlation between M d andṀ will be perfectly linear. This trend is illustrated in the top panel of Figure 1: the solid lines show three theoretical isochrones at different ages, while the dots represent a synthetic population of 100 disks obtained with Diskpop evolving in time with the same color coding. The aforementioned convergence to the theoretical isochrone starts as early as 1 Myr, while at 10 Myr the population is almost fully evolved and closely resembles the theoretical curve. From this argument, we can expect the moments of the distribution of t lt to evolve in the viscous case as follows: (i) the mean value of t lt will converge towards the actual age of the region, (ii) the spread will decrease until t ν < t for every disk in the population, (iii) the skewness will increase. For a more detailed discussion on the expected and observed evolution of the skewness, we refer to Appendix A.
The bottom panel of Figure 1 shows a synthetic population of disks evolved via MHD winds in the M d −Ṁ plane. As discussed in paragraph 2.2, the evolved population does not converge to the same isochrone: the large spread at all ages is such that making a prediction on the time evolution of the distribution of t lt is not as straightforward as for a viscous population. Tabone et al. (2022b) have shown that, assuming an exponential distribution of t acc,0 (which is determined fitting the observed disk fraction), the distribution of t lt does not depend on time; however, this result is specific of the exponential distribution. If we consider a different distribution of t acc,0 , that of t lt for an evolved population may depend on time: this is the case for our choice of a log-normal distribution of t acc,0 , which can still reproduce both the disk and accretion fraction (see Appendix B).  Figure 2 shows the time evolution of the mean (top) and width (bottom) of the distributions of t lt for the viscous (blue) and MHD (yellow) models. The lighter shades of both models include an additional observational uncertainty, σ obs , that we implemented by adding an extra spread on the disk mass and the accretion rate, of 0.1 dex and 0.45 dex respectively (as an estimate of the observational uncertainty, see Manara et al. 2023, Testi et al. 2022. As stated in Section 2, we performed 100 runs for each simulation: the solid line represents the median, while the shaded areas around it show the 25th-75th percentile intervals. As the MHD model removes disks more effectively, the sample size decreases more than in the viscous case, making the statistical fluctuations between different simulations larger: this leads the yellow lines to have broader shaded areas. Considering the mean values of the distributions, adding σ obs only slightly shifts the curves for both the viscous and MHD case, resulting in a negligible difference. The two evolutionary models differ at early stages (< 1 Myr), but soon reach a common behavior that makes them indistinguishable within the 25th-75th percentile intervals. On the other hand, the widths of the distribution (bottom panel) significantly differ from one case to the other. The viscous case without additional uncertainty (darker blue) steeply decreases, as expected from viscous theory (Lodato et al. 2017) and discussed in paragraph 3.1. This is not the case for the MHD prescription (orange): while the general trend is still decreasing, it is not as steep as the viscous, and ultimately does not tend to zero but rather to an evolved value determined by the initial conditions.
The convolution with observational uncertainty in the viscous case (light blue) significantly shifts the curve up, as well as modifying its shape. The total width of the distribution is the root sum squared of the intrinsic spread (σ int ) and the observational uncertainty (σ obs ), σ tot = σ int 2 + σ 2 obs . The intrinsic spread σ int , given by the initial conditions, tends to zero as discussed above: therefore, we expect the final width to tend to σ obs , which is exactly what we recover. This causes the evolved population to have a significantly larger spread than that predicted by theory. On the contrary, despite still being shifted at larger values as an effect of the additional uncertainty, in the MHD case (yellow) the shape of the curve is not dramatically modified. This is because σ int is comparable to σ obs at all times, which makes this argument strongly dependent on the initial condition: as the total spread is given by √ σ int 2 + σ obs 2 , the behavior of the MHD case will only be significantly different from the viscous case if σ int is non negligible with respect to σ obs . In our previous work (Somigliana et al. 2022) we have shown how initial spreads of 0.65 dex and 0.52 dex for M d and R d respectively are able to reproduce the observed spreads around the correlations with the stellar masses; therefore, we set these values for the MHD simulation, while we choose a bigger spread of 1 dex for the viscous case, as it can better reproduce the observed values (see 3.3).
As mentioned in Section 2.3, the purely viscous model does not account for disk dispersal. Without exploring the whole parameter space, which is beyond the scope of this Letter, we have run a test model with photoevaporation, assuming the standard model of Owen et al. (2010), with a mass-loss rate of 10 −10 M ⊙ yr −1 following the latest constraints (Alexander et al. 2023). The mean and the width of the distribution of t lt increase with respect to the purely viscous case, but the difference is minimal and becomes negligible including the and width (bottom) for the viscous (blue) and MHD (yellow) models, including observational uncertainties, with the observations (gray diamonds). The shaded areas are as in Figure 2, while the gray bars represent the interval between the 16th and 84th percentiles (top) and the uncertainty on the width (bottom). While both models overestimate the mean values (see text for details), the evolution of the width of the distribution suggests a better match with the MHD model. observational uncertainty; therefore, our conclusions are not affected.

Comparison with the observations
In paragraph 3.2 we have shown the viscous and MHD predictions for the time evolution of the mean and width of the distribution of t lt ; in this paragraph, we compare our results with observations of different star-forming regions. We used the table 1 compiled by Manara et al. (2023) for Taurus, Lupus, Chameleon I and Upper Sco, and the data by Testi et al. (2022) for L1688 (to limit the contamination from sub-populations with different ages in the Ophiuchus complex).
Before commenting on the comparison itself, it is important to note that our simulations do not include dust evolution, making our definition of disk mass solely based on the gas content of disks; on the other hand, the observed disk masses rely on sub-mm fluxes, tracing the dust content instead. As the bulk of disk masses is in the 1 The table is available at http://ppvii.org/chapter/15/. gaseous phase, inferring the total mass from dust observations requires to i) constrain the dust-to-gas ratio in disks and ii) assume optically thin emission; however, as the accuracy of these assumptions is debated, the community is striving towards obtaining more reliable disk mass estimates (see Bergin et al. 2013, McClure et al. 2016Veronesi et al. 2021 for dynamical measurements; Anderson et al. 2022, Trapman et al. 2022 for a combination of gaseous tracers). The results of the ALMA Large Programs AGE-PRO and DECO will further contribute to this goal; moreover, the advent of the ALMA Band 1 and ngVLA will allow to move to longer wavelengths, where dust emission is less optically thick (Tazzari et al. 2021). In light of these forthcoming developments, our work can be considered a prediction that will be interpreted to its full potential with the results of this observational effort. The data comparison presented in the following is therefore intended as a state-of-the-art, which we anticipate to revise in the near future. Figure 3 shows the result of our comparison: the mean and width of the distribution are shown in the top and bottom panel respectively, and both include the viscous and MHD (blue and yellow line, as in Figure 2 and 5) numerical evolution. The gray diamonds represent the observed star-forming regions. None of the two evolutionary mechanisms reproduces the observed mean values, which are systematically lower. A potential reason for this mismatch could be an underestimation of disk masses; a difference of a factor as little as 3 in the observed masses would be sufficient to explain the discrepancy with the models -confirming the need to repeat this comparison with more accurate disk masses estimates. Moreover, Zagaria et al. (2022a) have shown how taking stellar multiplicity into account can explain the high accretors in Upper Sco; we expect this effect to shift the theoretical prediction to lower values of t lt for evolved populations. Dust growth and evolution prescriptions, which were not included in this work, are also likely to play a role as they can better explain the observed disk mass -accretion rate correlation (Sellek et al. 2020). The width of the distribution, on the other hand, provides more interesting results. The viscous prediction manages to marginally recover observed values at the earliest evolutionary stages, but as such values increase in time, the discrepancy with the viscous expectation grows larger and larger. This result was already anticipated by Manara et al. (2020) (see also Manara et al. 2023). It should be kept in mind that our viscous simulations have a σ int of 1 dex for both the disk mass and radius (see Table 1); as large as the intrinsic spread can be, the steeply decreasing viscous trend will always evolve the width of the distribution to σ obs . The MHD simulation instead falls within the error bars of the earliest observed star-forming region, up until ages on ∼ 2.5 Myr. There is an increasing discrepancy for more evolved populations, up until around 20% for Upper Sco; however, the oldest populations also represent the less complete samples, and therefore they carry a significant bias that should be kept in mind when comparing with simulations. Moreover, there are caveats to our own simulations, as in the viscous case we neglect disk dispersal mechanisms (such as internal or external photoevaporation, e.g. Malanga et al. in prep.) and only consider a detection threshold in disk masses.

DISCUSSION AND CONCLUSIONS
In this work, we have investigated how the time evolution of the distribution of a population of disks in the M d −Ṁ plane is impacted by the evolutionary model, considering the viscous and MHD prescriptions respectively. We have presented a combination of analytical considerations and numerical simulations, performed through the 1D population synthesis code Diskpop, in the case of a log-normal distribution of initial accretion timescales (which reproduces both the disk and accretion fraction). We find that, while the mean of the distribution of t lt = M d /Ṁ is not significantly impacted by the chosen model, the expected behavior of the width shows considerable differences depending on the evolutionary prescription; when including the observational biases in the form of additional uncertainty, this distinctive behavior is maintained.
Our predictions will be exploited to their full potential through a comparison with the results of the current observational effort to obtain direct estimates of disk gas masses; for the time being, we compare our evolutionary trends with the latest available observational data (based on dust observations) in different star-forming regions. We find that the purely viscous case only manages to marginally reproduce the observations at the earliest ages, while the MHD curve resembles them better. Based on these results, we suggest the analysis of these distributions as a viable method to disentangle between the viscous and MHD evolutionary models; our data comparison hints at a better agreement with the MHD model.

ACKNOWLEDGMENTS
We thank an anonymous referee for their comments that helped us improving the clarity of the manuscript. This The skewness of a distribution, defined as the third standardized moment, measures the asymmetry of the distribution about its mean. As we mentioned in paragraph 3.1, alongside the mean value and the width, in the viscous case we expect also the skewness of the distribution of t lt to evolve in time; in this Appendix we discuss this theoretical expectation and show the results of our numerical simulations.  Figure 1. Full dots represent disks whose initial viscous timescale is shorter than the age of the population, and that can therefore be considered evolved.
The left panel of Figure 4 shows a population of viscously evolving disks (dots) at three subsequent ages, as well as the corresponding theoretical isochrones (solid lines). Full dots represent disks whose initial viscous timescale is shorter than the age of the population, which as a whole can therefore be considered evolved: from viscous theory, such disks are expected to have reached the self-similar condition and lie on the analytical isochrone, i.e., to show a linear correlation between the disk mass and the accretion rate. On the other hand, empty dots represent not-yet-evolved disks, which lie below the theoretical isochrone. As the population evolves, more disks satisfy the t ν < t condition, as can be visualized by the increasing number of full dots in Figure 4; this implies that more disks lie on the theoretical isochrone, bringing the population on the M d −Ṁ plane closer to a line. While this causes the width of the distribution of t tl to decrease with time, the skewness on the other hand increases -as we show in the right panel of Figure 4, which represents the corresponding histograms at all ages. This skewing effect is due to the fact that younger disks, which do not lie on the isochrone yet, have a t lt longer than the actual age of the region, and therefore contribute to positively skew the distribution -while evolved disks, which make up the bulk of the population, cluster close to the mean value. Figure 5 shows the evolution of the skewness of a population of disks generated and evolved with Diskpop with the same color coding and shaded areas as Figure 1; the left panel represents the case with no observational uncertainty, where the viscous distribution (blue) gets more and more skewed as expected, growing by a factor of 2 between 0.1 and 10 Myr. On the other hand, the MHD distribution (orange) remains symmetrical within the 25th-75th percentile for the whole evolution, resulting in a factor 3 difference from the viscous model for evolved populations. As significant as this theoretical difference is, including the observational biases (right panel) completely smooths it out: the two expected observed behaviors are indistinguishable once convoluted with the additional observational uncertainties.
In conclusion, while the evolution of the skewness makes an interesting theoretical argument stemming from the different interpretation of isochrones in the two models, it does not provide a reliable method to compare viscosity and MHD from the observational point of view.

B. TIME EVOLUTION OF THE DISTRIBUTION OF t lt
As t lt depends on t acc,0 as t lt = (1 + f M )(2t acc,0 − ωt), the evolved distribution of t lt is determined by the choice of initial distribution of t acc,0 : Tabone et al. (2022b) have shown that, choosing an exponential distribution for t acc,0 , the corresponding distribution of t lt reads where f M is defined in Tabone et al. (2022a) and τ = 2.5 Myr to fit the disk fraction, f D (t) = exp (−t/τ ). As f D is only a normalization factor, (B1) still have an exponential shape; moreover, it does not depend on time, as well as its mean value. On the other hand, if we pick a log-normal distribution for t acc,0 , we can still reproduce both the disk and the accretion fraction (see Appendix C) but in that case the evolved distribution of t lt becomes where µ and σ are the mean value and width of the initial log-normal distribution. Notice that Equation (B2) is not a log-normal in t lt ; moreover, it does depend on time, and so does its mean value and spread.

C. IMPACT OF INTERNAL PHOTOEVAPORATION
As mentioned in the main paper, disk dispersal in an intrinsic feature of MHD winds. These models manage to reproduce both the disk and accretion fraction, defined as the fraction of young stars with infrared excess (Hernández et al. 2007) and accreting (i.e., withṀ > 10 −11 M ⊙ yr −1 following Fedele et al. 2010) objects respectively, as shown by the orange lines in Figure 6. On the other hand, purely viscous models do not account for disk dispersal. This leads to a mismatch between the predicted and observed disk and accretion fraction, represented by the blue lines in Figure 6: the disk fraction is almost constant to 1, the little decrease being due to the observational threshold that we introduced in our simulations (considering dispersed disks with masses lower than 10 −6 M ⊙ , see Section 2.3), while the accretion fraction does decrease, but not enough to match the observed values. This problem is usually overcome in the literature by including internal photoevaporation, a two-timescale process that introduces a disk dispersal mechanism, allowing to reproduce the observations as shown by the purple lines in Figure 6. We ran the test simulation presented in this Appendix using the standard photoevaporative model of Owen et al. (2012), with a mass-loss rate of 10 −1− M ⊙ yr −1 , consistent with the latest constraints (Alexander et al. 2023).  Figure 2. The dashed blue lines show the exponential fits to the data. Following the original paper, we define the accretion fraction as the fraction of sources with accretion rate higher than 10 −11 M⊙/yr. Our choice of a log-normal distribution of initial accretion timescales for the MHD model reproduces both the disk and accretion fraction, as does the exponential distribution chosen by Tabone et al. (2022b). The viscous model does not reproduce any of the fractions due to the lack of a disk dispersal mechanism, while including internal photoevaporation allows to recovered the observed behavior.
Once internal photoevaporation kicks in, it lowers the accretion rates for a given disc mass, introducing therefore a spread in the M d −Ṁ plane (Somigliana et al. 2020); therefore, it could in principle affect the conclusions of this work. However, we have tested that the mean and width of the t lt distribution in the presence of photoevaporation do not significantly deviate from the purely viscous prediction; without observational spread the photoevaporative case lies between the viscous and MHD models, and becomes indistinguishable from the viscosity when the observational spread is included.