# The evolution of the $M_d-M_\star$ and $\dot{M}-M_\star$ correlations traces protoplanetary disc dispersal

Alice Somigliana<sup>1,2,*</sup>①, Leonardo Testi<sup>3,4</sup>, Giovanni Rosotti<sup>5</sup>①, Claudia Toci<sup>1</sup>, Giuseppe Lodato<sup>5</sup>, Rossella Anania<sup>5</sup>①, Benoît Tabone<sup>6</sup>①, Marco Tazzari<sup>4</sup>①, Ralf Klessen<sup>7,8</sup>①, Ugo Lebreuilly<sup>9</sup>①, Patrick Hennebelle<sup>9</sup>, and Sergio Molinari<sup>10</sup>

Received 16 May 2024 / Accepted 1 July 2024

## ABSTRACT

Observational surveys of entire star-forming regions have provided evidence of power-law correlations between the disc-integrated properties and the stellar mass, especially the disc mass ($M_d \propto M_\star^{\lambda_m}$) and the accretion rate ($\dot{M} \propto M_\star^{\lambda_{acc}}$). Whether the secular disc evolution affects said correlations is still a matter of debate: while the purely viscous scenario has been investigated, other evolutionary mechanisms could have a different impact. In this paper, we study the time evolution of the slopes $\lambda_m$ and $\lambda_{acc}$ in the wind-driven and viscous-wind hybrid case and compare it to the purely viscous prediction. We use a combination of analytical calculations, where possible, and numerical simulations performed with the 1D population synthesis code *Diskpop*, which we also present and release to the community. Assuming $M_d(0) \propto M_\star^{\lambda_{m,0}}$ and $\dot{M}(0) \propto M_\star^{\lambda_{acc,0}}$ as initial conditions, we find that viscous and hybrid accretion preserve the power-law shape of the correlations, while evolving their slope; on the other hand, magneto-hydrodynamic winds change the shape of the correlations, bending them in the higher or lower end of the stellar mass spectrum depending on the scaling of the accretion timescale with the stellar mass. However, we show how a spread in the initial conditions conceals this behaviour, leading to power-law correlations with evolving slopes as in the viscous and hybrid case. We analyse the impact of disc dispersal, intrinsic in the wind model and due to internal photoevaporation in the viscous case: we find that the currently available sample sizes (~30 discs at 5 Myr) introduce stochastic oscillations in the slopes' evolution, which dominate over the physical signatures. We show that we could mitigate this issue by increasing the sample size: with ~140 discs at 5 Myr, corresponding to the complete Upper Sco sample, we would obtain small enough error bars to use the evolution of the slopes as a proxy for the driving mechanism of disc evolution. Finally, from our theoretical arguments, we discuss how the observational claim of steepening slopes necessarily leads to an initially steeper $M_d-M_\star$ correlation with respect to $\dot{M}-M_\star$.

**Key words.** accretion, accretion disks – planets and satellites: formation – protoplanetary disks – stars: pre-main sequence

## 1. Introduction

The secular evolution of protoplanetary discs is deeply inter-
twined with both the planet formation process (Morbidelli et al.
2012) and the accretion onto the central protostar (Hartmann
et al. 1998). Planetesimals, the building blocks of planets, form
and evolve within the disc following the dynamics of either the
gaseous or solid component, depending on their relative size and
their coupling (or lack thereof) with the gas particles; on the
other hand, the protostar is fed by the disc itself, through the
accretion of material that loses angular momentum and drifts
inwards. The ideal ground to explore the connection between
protoplanetary discs and their host stars is provided by large
surveys of entire star-forming regions, targeting the properties of

both discs and protostars; the last decade has seen a significant
observational effort in the direction of these population-level
studies, also thanks to the advent of facilities such as the
Atacama Large Millimeter Array (ALMA) (see the PPVII
reviews by Manara et al. 2023; Miotello et al. 2023).

Disc masses and accretion rates are arguably the most stud-
ied integrated disc properties. Accretion rates are inferred from
the spectra of the central stars, which show an excess emis-
sion (especially prominent in the UV) when accretion is taking
place; surveys performed across different star-forming regions
(Muzerolle et al. 2003; Natta et al. 2004; Mohanty et al. 2005;
Dullemond et al. 2006; Herczeg & Hillenbrand 2008; Rigliaco
et al. 2011; Manara et al. 2012; Alcalá et al. 2014; Manara
et al. 2016; Alcalá et al. 2017; Manara et al. 2017; Venuti et al.
2019; Manara et al. 2020) agree on the presence of a power-law

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

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

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

<footnote n="4">⁴ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="5">⁵ Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy</footnote>

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

<footnote n="7">⁷ Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany</footnote>

<footnote n="8">⁸ Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany</footnote>

<footnote n="9">⁹ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="10">¹⁰ Istituto Nazionale di Astrofisica-IAPS, Via Fosso del Cavaliere 100, 00133 Roma, Italy</footnote>

<footnote>* Corresponding author; alice.somigliana@eso.org</footnote>
correlation between the accretion rate and the stellar mass, $\dot{M} \propto M_{\star}^{\lambda_{\text{acc,obs}}}$ (following the notation of Somigliana et al. 2022). On the other hand, disc masses have traditionally been determined from observations of the sub-millimetre continuum emission of the solid component of discs; due to the large number of assumptions involved in converting fluxes into total disc masses (see Miotello et al. 2023), one of the current main goals of the protoplanetary disc community is the accurate determination of total disc masses – both from dynamical constraints (Veronesi et al. 2021; Lodato et al. 2023) and direct measurements of the total gas content (e.g. Bergin et al. 2013; Anderson et al. 2022; Trapman et al. 2022). Despite the systematic uncertainties involved in their determination, dust-based disc masses also seem to show a power-law correlation with the stellar mass, $M_{\text{d}} \propto M_{\star}^{\lambda_{\text{m,obs}}}$, across different star-forming regions (Ansdell et al. 2016, 2017; Barenfeld et al. 2016; Pascucci et al. 2016; Testi et al. 2016, 2022; Sanchis et al. 2020).

The existence of the disc mass-stellar mass and accretion rate-stellar mass correlations is now generally accepted; however, there is no consensus on the physical reason behind their establishment and their evolution with time. While the $\dot{M}-M_{\star}$ correlation appears to have a roughly constant slope<fnref n="1" /> of $\lambda_{\text{acc,obs}} \approx 1.8 \pm 0.2$ (as first suggested by Muzerolle et al. 2003 and supported by many of the following works mentioned above), the $M_{\text{d}}-M_{\star}$ correlations is claimed to be steepening with time (Ansdell et al. 2017), from the lowest $\lambda_{\text{m,obs}}(t \sim 1 \text{ Myr}) = 1.7 \pm 0.2$ (Taurus) to the highest $\lambda_{\text{m,obs}}(t \sim 5 \text{ Myr}) = 2.4 \pm 0.4$ (Upper Sco). Whether these correlations reflect the initial conditions of disc populations, or are rather a product of the secular evolution, is still under debate. Both possibilities have been discussed for the $\dot{M}-M_{\star}$ correlation: Alexander & Armitage (2006) assumed it to hold as an initial condition, favouring the correlation to be present in young populations, whereas Dullemond et al. (2006) derived it from a simple model of disc formation from a rotating collapsing core, which provided an explanation for evolved disc populations. At the same time, the claimed increase in the slope of $M_{\text{d}}-M_{\star}$ does suggest an evolutionary trend; Somigliana et al. (2022) have found that, assuming power-law correlations between both $M_{\text{d}}$ and $\dot{M}$ and the stellar mass as initial conditions, secular evolution can indeed alter the slopes of the correlations themselves (see Section 3 for details). However, their analysis was limited to the standard viscous evolution paradigm, whereas the driving mechanism of accretion is far from being constrained (see Manara et al. 2023 for a review).

The traditional viscous accretion model prescribes a macroscopic viscosity as the cause of redistribution of angular momentum within the disc (Lynden-Bell & Pringle 1974; Pringle 1981). In this scenario, while part of the material loses angular momentum and moves radially closer to the star, other material gains the same amount of angular momentum and moves further away, increasing the disc size. The viscous paradigm can explain many key features of disc evolution, but it cannot account for disc dispersal – as determined from the observational evidence of an exponentially decreasing fraction of both disc-bearing (Hernández et al. 2007) and accreting (Fedele et al. 2010) sources in star-forming regions with time; furthermore, the low levels of turbulence detected in discs (Pinte et al. 2016; Flaherty et al. 2018; Rosotti 2023) appear incompatible with the observed evolution. While the discrepancy in the disc and accretion fraction can be mended considering mechanisms such as internal or external photoevaporation

(Alexander et al. 2014; Winter et al. 2018), which effectively clear discs on timescales comparable with the observed decline, the tension between the expected and observed amount of turbulence does not appear to have been solved yet. On the other hand, the magneto-hydrodynamic (MHD) disc winds’ scenario offers a promising alternative. Pioneering work (Blandford & Payne 1982; Ferreira 1997) supported by recent numerical simulations (e.g. Béthune et al. 2017) demonstrated that MHD winds launched from the disc surface have the net effect of removing angular momentum as a consequence of the extraction of material; MHD wind-driven accretion can even lead to disc dispersal (Armitage et al. 2013; Tabone et al. 2022b). Following disc evolution at a population level in numerical simulations remains out of reach; however, 3D core-collapse simulations have shown how non-ideal magnetohydrodynamics and ambipolar diffusion play a fundamental role in shaping the resulting population of early-type young stellar objects (Lebreuilly et al. 2021, 2024). While some 3D studies of isolated disc formation have attempted to bridge the gap between Class 0/I and Class II stages (Machida & Hosokawa 2013; Hennebelle et al. 2020; Xu & Kunz 2021a,b; Machida & Basu 2024; Mauxion et al. 2024), the high numerical cost of the simulations for 3D population synthesis does not allow one to follow the evolution of the discs up to very evolved stages where they can be considered isolated from the surrounding environment. MHD wind-driven disc populations can however be modelled in 1D using simple prescriptions as proposed by Suzuki et al. (2016) or Tabone et al. (2022a). Detecting characteristic signatures of either of the two evolutionary prescriptions is a compelling issue (Long et al. 2022; Alexander et al. 2023; Somigliana et al. 2023; Trapman et al. 2023; Coleman et al. 2024).

In the context of the evolution of the correlations between the disc properties and the stellar mass, while the purely viscous scenario has been extensively studied by Somigliana et al. (2022), the wind-driven paradigm remains unexplored; with this paper, we address this deficiency and investigate the impact of MHD wind-driven evolution on the $M_{\text{d}}-M_{\star}$ and $\dot{M}-M_{\star}$ correlations, with a particular focus on their time evolution. We also extend the work of Somigliana et al. (2022) by including internal photoevaporation to the viscous framework. We employ numerical simulations of populations of protoplanetary discs, performed with the population synthesis code Diskpop, which we also introduce and release to the community.

The paper is structured as follows: in Section 2, we present Diskpop and describe its main features, set up and solution algorithm; in Section 3, we discuss the time evolution of the $M_{\text{d}}-M_{\star}$ and $\dot{M}-M_{\star}$ correlations in three evolutionary scenarios from the theoretical perspective; in Section 4, we show the impact of a spread in the initial conditions and dispersal mechanisms on the evolution of the slopes, and we present the numerical results obtained from realistic disc population synthesis; in Section 5 we interpret the implications of our findings in the context of the observational determination of the slopes, and finally in Section 6 we draw the conclusions of this work.

## 2. Numerical methods: Diskpop

In this section, we present the 1D population synthesis code Diskpop<fnref n="2" />. We describe the master equation for the secular

<footnote n="2">² Diskpop and the output analysis library popcorn can be installed via the Python Package Index, pip install diskpop and pip install popcorn_diskpop. The full documentation and tutorials are available at https://alicesomigliana.github.io/diskpop-docs/index.html. If you use Diskpop in your work, please cite this paper (Somigliana et al. 2024).</footnote>

<footnote n="1">¹ Throughout this work, we use ‘slope’ as a synonym for the power-law index, referring to the correlations in the logarithmic plane.</footnote>
evolution of discs (Section 2.1), the initial conditions to generate a synthetic population (Section 2.2), the solution algorithm (Section 2.3), and the user interface and output (Section 2.4). For a more detailed description, we refer to the code documentation; for a validation of the code, see Appendix C.

## 2.1. Master equation

The master equation of protoplanetary disc evolution,

$$ \frac{\partial \Sigma}{\partial t} = \frac{3}{r} \frac{\partial}{\partial r} \left[ \frac{1}{\Omega r} \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} - \dot{\Sigma}_{\text{photo}}, \quad (1) $$

describes the time evolution of the gas surface density in the most general framework, where $\Sigma$ is the gas surface density, $\Omega$ the Keplerian orbital frequency, $\alpha_{SS}$ the Shakura & Sunyaev (1973) $\alpha$ parameter, $\alpha_{DW}$ the MHD equivalent of $\alpha_{SS}$ (Tabone et al. 2022a), $c_s$ the sound speed, and $\lambda$ the magnetic lever arm parameter, which quantifies the ratio of extracted to initial specific angular momentum. The four terms on the right hand side (RHS) refer to (i) the viscous torque, whose strength is parameterised by $\alpha_{SS}$, (ii) the wind-driven accretion, which corresponds to an advection term, parameterised by $\alpha_{DW}$, (iii) mass loss due to MHD disc winds, parameterised by $\lambda$ and (iv) mass loss due to other physical phenomena (in our case, we consider internal and external photoevaporation). Depending on the values of the specific parameters, Equation (1) can describe a purely viscous ($\alpha_{DW} = 0$), purely MHD wind-driven ($\alpha_{SS} = 0$) or hybrid ($\alpha_{SS}, \alpha_{DW} \neq 0$) evolution, with ($\dot{\Sigma}_{\text{photo}} \neq 0$) or without ($\dot{\Sigma}_{\text{photo}} = 0$) the influence of photoevaporation. In the following, we briefly describe the various evolutionary scenarios and the available analytical solutions.

**Viscously evolving discs.** In the case of purely viscous evolution, the MHD winds parameter $\alpha_{DW}$ asset to zero. If we also neglect the influence of photoevaporation, Equation (1) reduces to the first term on the RHS and its solution depends on the functional form of the effective viscosity, parameterised as $\nu = \alpha_{SS}c_sH$ (where $H$ is the vertical height of the disc). A popular analytical solution for viscous discs is the Lynden-Bell & Pringle (1974) self-similar solution, which assumes viscosity to scale as a power-law of the radius ($\nu \propto R^\gamma$).

**MHD winds-driven evolution.** There are two classes of analytical solutions to Equation (1) in the MHD wind-driven scenario, associated with a specific prescription of $\alpha_{DW}$ (Tabone et al. 2022a). We briefly describe their key features, and refer to the original paper for their derivation and an in-depth discussion.

1. The simplest class of solutions (so-called 'hybrid solutions'), which highlight the main features of wind-driven accretion in comparison to the viscous model, assume a constant $\alpha_{DW}$ with time; these solutions depend on the value of $\psi \equiv \alpha_{DW}/\alpha_{SS}$, which quantifies the relative strength of the radial and vertical torque.

2. Another class of solutions, which describe the unknown evolution of the magnetic field strength, assume a varying $\alpha_{DW}$ with time. To obtain these, Tabone et al. (2022a) parameterised $\alpha_{DW}(t) \propto \Sigma_c(t)^{-\omega}$, with $\Sigma_c = M_d(t)/2\pi R_c^2(t)$ (where $R_c$ is a characteristic radius) and $\omega$ as a free parameter, and neglect the radial transport of angular momentum ($\alpha_{SS} = 0$).

**Photoevaporation.** The generic $\dot{\Sigma}_{\text{photo}}$ term in Equation (1) allows to account for photoevaporative processes, both internal

and external. The exact form of $\dot{\Sigma}_{\text{photo}}$ depends on the specific model considered; therefore, the availability (or lack thereof) of analytical solutions needs to be considered case by case.

Diskpop allows populations of discs to be evolved analytically. In particular, as of this release, it includes implementations of the Lynden-Bell & Pringle (1974) self-similar solution and all the analytical solution proposed by Tabone et al. (2022a). In the cases where Equation (1) cannot be solved analytically, the code relies on the solution algorithm described in Section 2.3.

## 2.2. Initial conditions and parameters

Every Diskpop simulation begins with the generation of a synthetic population of Young Stellar Objects (YSOs). Each YSO constitutes a star and a disc, whose key initial parameters (stellar mass, disc mass, accretion rate, disc radius, evolutionary parameters $\alpha_{SS}, \alpha_{DW}, \lambda, \omega...$) can be set by the user. In the following, we describe the standard case where we consider the stellar masses to be distributed according to an Initial Mass Function (IMF) and correlating with the disc mass and radius, and briefly mention the other possible choices; for a deeper discussion, we refer to the Diskpop documentation.

Diskpop assembles YSOs by determining their parameters as follows:

1. Stellar mass $M_*$: determined following the Kroupa (2001) IMF. Other possible choices are a constant mass for all the stars in the population, or a set of custom stellar masses.

2. Initial disc mass $M_d$, accretion rate $\dot{M}$: determined from log-normal distributions of given width and mean value. In the standard case, Diskpop considers an initial power-law correlation between the initial $M_d$ and $\dot{M}$ and the stellar mass (see Section 3.1 for a detailed discussion), where the normalisation at $1 M_*$, the slope and the scatter around the power-laws are free parameters. If the correlations with the stellar mass are neglected, the user sets the mean value and spread of the distributions.

3. Accretion parameters ($\alpha_{SS}, \alpha_{DW}, \lambda, \omega$): global properties of the whole population, given as input from the user. By setting the parameters controlling accretion, Diskpop determines the disc radius $R_d$ and accretion timescale $t_{acc}$ – which are instead disc-specific and linked to the disc mass and accretion rate.

4. Internal photoevaporation parameters ($\dot{M}_{\text{wind}}, L_X$): the total photoevaporative mass-loss rate, $\dot{M}_{\text{wind}}$, can either be set by the user or computed from the stellar X luminosity $L_X$ as (Owen et al. 2012)

$$ \dot{M}_{\text{wind}} = 6.25 \times 10^{-9} \times \left(\frac{M_{\star}}{M_{\odot}}\right)^{-0.068} \left(\frac{L_X}{10^{30} \text{ergs}^{-1}}\right)^{1.14} M_{\odot} \text{yr}^{-1}. $$

The surface mass-loss profile $\dot{\Sigma}_{\text{photo}}$ (Equation (B2) in Owen et al. 2012) is then scaled so that $\int 2\pi R \dot{\Sigma}_{\text{photo}} dR$ is equal to $\dot{M}_{\text{wind}}$. As for $M_d$ and $\dot{M}$, $\dot{M}_{\text{wind}}$ (or equivalently $L_X$) is extracted from a log-normal distribution whose mean is determined assuming power-law correlations with the stellar mass, while the normalisation at $1 M_*$, the slope and the width of the distribution are free parameters.

5. External photoevaporation parameters (FUV): FUV flux experienced by each disc, in units of $G_0^3$. This parameter can be set to any value accessible in the FRIEDv2 grid

<footnote n="3">³ $G_0$ stands for the Habing unit (Habing 1968), the flux integral over the range of wavelengths [912–2400] Å weighted by the average value in the solar neighbourhood ($1.6 \times 10^{-3} \text{ erg s}^{-1} \text{ cm}^{-2}$).</footnote>
of mass loss rate (Haworth et al. 2023), spanning from 1 to $10^5$ G$_0$.

## 2.3. Solution algorithm

After generating the initial population of YSOs as described above, Diskpop proceeds to evolve it by integrating the master equation (1). Our solution algorithm employs an operator splitting method: the original equation is separated into different parts over a time step, and the solution to each part is computed separately. Then, all the solutions are combined together to form a solution to the original equation. We split Equation (1) into five different pieces, related to viscosity, wind-driven accretion onto the central star, wind-driven mass loss, internal and external photoevaporation respectively. Furthermore, Diskpop includes the possibility to trace the dust evolution in the disc, which is split in radial drift and dust diffusion. In the following, we describe the solution algorithm for each process.

1. Viscous accretion: the standard viscous solver is based on the freely available code by Booth et al. (2017). We assume a radial temperature profile $T \propto R^{-1/2}$, which results in $c_s \propto R^{-1/4}$ and $H/R \propto R^{1/4}$. We note that this implies $v \propto R$ (i.e. $\gamma = 1$), which will be the case from now on. We assume $H/R = 1/30$ at 1 AU and a mean molecular weight of 2.4. We refer to the original paper for details on the algorithm.

2. Wind-driven accretion: the second term in Equation (1) is effectively an advection term. The general form of the advection equation for a quantity $q$ with velocity $v$ is $\partial_t q(x, t) + v \partial_x q(x, t) = 0$; in the case of wind-driven accretion, the advected quantity is $R\Sigma$, while the advection (inwards) velocity is given by $v_{DW} = (3\alpha_{DW}Hc_s)/2R$. We solve the advection equation with an explicit upwind algorithm (used also for dust radial drift).

3. Wind-driven mass loss: the mass loss term (third in Equation (1)) does not involve any partial derivative, and therefore is simply integrated in time multiplying by the time step.

4. Internal photoevaporation: effectively, internal photoevaporation (implemented through the model of Owen et al. 2012) is another mass loss term – therefore, as above, its contribution is computed with a simple multiplication by the time step. Once the accretion rate of the disc drops below the photoevaporative mass loss rate, a gap opens in the disc at the radius of influence of photoevaporation: in the model of Owen et al. (2012), the prescription changes depending on the radial location in the disc, with respect to the gap itself. Later, the gap continues to widen; when it eventually becomes larger than the disc, we stop the evolution and consider the disc as dispersed.

5. External photoevaporation: for a given stellar mass and FUV flux experienced by the disc, the mass loss rate arising from external photoevaporation is obtained, at each radial position, from a bi-linear interpolation of the FRIEDv2 grid (Haworth et al. 2023) using the disc surface density at each radial cell. The outside-in depletion of material is implemented following the numerical approach of Sellek et al. (2020): we define the ‘truncation’ radius, $R_t$, as the position in the disc corresponding to the maximum photoevaporation rate (which is related to the optically thin to thick transition of the wind), and we remove material from each grid cell at $R > R_t$ weighting on the total mass outside this radius. The mass loss attributed to the cell $i$ can be written as:

$$ \dot{M}_{\text{ext},i} = \dot{M}_{\text{tot}} \frac{M_i}{M(R > R_t)}, \quad (2) $$

where $M_i$ is the mass contained in the cell $i$, and $\dot{M}_{\text{tot}}$ is the total mass loss rate outside the truncation radius.

6. Dust evolution<fnref n="4" />: based on the two populations model by Birnstiel et al. (2012) and the implementation of Booth et al. (2017). We consider the dust grain distribution to be described by two representative sizes, a constant monomer size and a time-dependent larger size, which can grow up to the limit imposed by the fragmentation and radial drift barriers. We evolve the dust fraction of both sizes following Laibe & Price (2014), and also include a diffusive term: the diffusion comes from the coupling with the turbulent gas, which has the effect of mixing the dust grains, counteracting gradients in concentration (Birnstiel et al. 2010). The dust-gas relative velocities are computed following Tanaka et al. (2005) and include feedback on the gas component. We refer to Booth et al. (2017) for details on the numerical implementation. Dust evolution is included in the release of Diskpop, however the scientific results presented in this work are based on gas simulations only.

The separate pieces of Equation (1) must be solved over the same time step to be joined in a coherent solution. We calculate the time step for each process imposing the Courant-Friedrichs-Lewy (CFL) condition. The CFL condition reads $\Delta t = C \min(\Delta x/v)$ and ensures that, within one time step $\Delta t$, the material moving at velocity $v$ does not flow further than one grid spacing $\Delta x$. The Courant number $C$ must be positive and smaller than 1, with $C = 1$ corresponding to the maximum allowed timestep to keep the algorithm stable. In our implementation, we pick $C = 0.5$. We use zero gradients boundary conditions, setting the value of the first and last cell in our grid to that of the second and second to last. We solve the equation on a radial grid of $10^3$ points with power-law spacing and exponent $1/2$, extending from $3 \times 10^{-3}$ au to $10^4$ au. From the physical point of view, this choice corresponds to assuming boundary layer accretion (see e.g. Popham et al. 1993; Kley & Lin 1996) – however the difference from magnetic truncation accretion is negligible beyond $\sim 10^{-3}$ au.

After each process has been solved separately, all the pieces are put back together to compute the new surface density, from which the integrated disc quantities are then calculated. As each disc evolves independently of the others in the population, the solver can easily be run in parallel.

## 2.4. User interface and output

The user interface of Diskpop is a .json parameters file that includes all the user-dependent parameters. Aside from the number of objects in the population and the evolutionary mechanism, the user can set the chosen IMF (either Kroupa 2001, single stellar mass, or custom input file), the distributions to draw the disc parameters from (single value, flat, normal, log-normal), as well as the normalisation, slope and spread of the correlations, the times at which snapshots are generated, and the initial dust-to-gas ratio. Furthermore, the user can determine a limit disc mass: this is to be intended as a threshold below which the disc would not be detectable anymore, and is therefore considered dispersed in the simulation as well. When a disc is dispersed, the corresponding YSO turns into a Class III object consisting of the central star only.

<footnote n="4">⁴ As the dust evolution module was forked from Richard Booth’s repository, users of Diskpop who wish to use dust in their work ought to cite Booth et al. (2017) together with this paper.</footnote>
The output of Diskpop is a .hdf5 file containing the properties of both the disc and the star at all chosen time steps for each YSO in the population: this includes the stellar mass, luminosity, temperature, disc mass, accretion rate, accretion timescale, gas and dust surface density, disc radius, dust grain sizes. The output can be easily read and analysed with the dedicated library popcorn, released with the code. For a more in-depth description of the parameters, the user interface and the output, we refer to the Diskpop documentation.

### 3. The time evolution of the correlations between disc properties and stellar mass under different accretion drivers: Analytical considerations

The existence of power-law correlations between the main integrated disc properties – namely the disc mass and stellar accretion rate – and the stellar mass is supported by various surveys across a number of different star-forming regions (e.g. on $M_d - M_\star$: Ansdell et al. 2016; Barenfeld et al. 2016; Pascucci et al. 2016; Testi et al. 2016; on $\dot{M}-M_\star$: Muzerolle et al. 2003; Natta et al. 2004; Mohanty et al. 2005; Alcalá et al. 2014; Manara et al. 2016; Alcalá et al. 2017; Manara et al. 2017; Venuti et al. 2019; Manara et al. 2020; Testi et al. 2022). However, whether the establishment and subsequent evolution of said correlations is a product of the secular evolution of discs, or rather an imprint of the initial conditions, remains unclear. Somigliana et al. (2022) explored a combination of both possibilities, assuming the correlations to hold as initial conditions and investigating the impact of purely viscous evolution; we briefly recall their main theoretical results (Section 3.1) and extend their analysis to the hybrid (Section 3.2) and purely wind-driven (Section 3.3) models from the theoretical perspective. We note that the results presented in this work are based on gas simulations.

Following Somigliana et al. (2022), we assume power-law correlations between the disc properties and the stellar mass to hold as ‘initial conditions’. We focus on the disc mass $M_d$, the stellar accretion rate $\dot{M}$ and the disc radius $R_d$, and label the slopes of their correlations with the stellar mass $\lambda_m$, $\lambda_{acc}$, and $\zeta$ respectively. The initial correlations are set as follows:

$$ \begin{cases} M_d(0) \propto M_\star^{\lambda_{m,0}}, \\ \dot{M}(0) \propto M_\star^{\lambda_{acc,0}}, \\ R_d(0) \propto M_\star^{\zeta_0}. \end{cases} \quad (3) $$

To analyse the impact of secular evolution on this set of initial conditions, we analytically determine the evolved expressions for $M_d(t)$, $\dot{M}(t)$ and $R_d(t)$ in the three different scenarios. Table 1 summarises the parameters introduced in this Section.

#### 3.1. Purely viscous model

The full calculations for the purely viscous case can be found in Somigliana et al. (2022). Here, we briefly remind the main assumptions and results, and we refer to the original paper for a detailed discussion.

As mentioned in Section 2.1, assuming a power-law scaling of viscosity with the disc radius ($\nu \propto R^\gamma$) allows the viscous evolution equation to be solved analytically, recovering the so-called self-similar solution (Lynden-Bell & Pringle 1974). In this case, the disc mass and accretion rate read

$$ M_d(t) = M_{d,0} \left(1 + \frac{t}{t_v}\right)^{1-\eta}, \quad (4) $$

<table><caption>**Table 1.** Summary and description of the parameters used throughout the paper.</caption><thead><tr><th>Parameter</th><th>Description</th></tr></thead><tbody><tr><td>λ</td><td>Magnetic lever arm parameter</td></tr><tr><td>ψ</td><td>Wind-to-turbulent α ratio</td></tr><tr><td>ω</td><td>Power-law index of α<sub>DW</sub> with Σ<sub>c</sub></td></tr><tr><td>ξ</td><td>Mass ejection index</td></tr><tr><td>f<sub>M,0</sub></td><td>Initial mass ejection-to-accretion ratio</td></tr><tr><td>λ<sub>m</sub></td><td>Power-law index of M<sub>d</sub> with M<sub>★</sub> (Eq. (3))</td></tr><tr><td>λ<sub>acc</sub></td><td>Power-law index of M<sub>★</sub> (Eq. (3))</td></tr><tr><td>ζ</td><td>Power-law index of R<sub>d</sub> with M<sub>★</sub> (Eq. (3))</td></tr><tr><td>β</td><td>Power-law index of H/R with M<sub>★</sub></td></tr><tr><td>μ</td><td>Power-law index of t<sub>acc</sub> with M<sub>★</sub></td></tr><tr><td>δ</td><td>λ<sub>m</sub> - λ<sub>acc</sub></td></tr></tbody></table>

**Notes.** The top block refers to the MHD parameters defined in Tabone et al. (2022a), while the bottom block shows the slopes of the correlations between the disc properties and the stellar mass.

$$ \dot{M}(t) = (\eta - 1) \frac{M_{d,0}}{t_v} \left(1 + \frac{t}{t_v}\right)^{-\eta}, \quad (5) $$

where $\eta = (5/2 - \gamma)/(2 - \gamma)$ and the viscous timescale $t_v = R_c^2/[3(2 - \gamma)^2\nu(R = R_c)]$ at the characteristic radius $R_c$. Because $\dot{M}_0 \propto M_{d,0}/t_{v,0}$, a power-law scaling of $M_{d,0}$ and $\dot{M}_0$ with the stellar mass implies the viscous timescale $t_{v,0}$ to scale as a power-law with the stellar mass as well, which we define as $t_{v,0} \propto M_\star^{\mu_0}$; furthermore, this scaling corresponds to the difference between the scaling of the disc mass with the stellar mass and of the accretion rate with the stellar mass. Defining $\delta_0 = \lambda_{m,0} - \lambda_{acc,0}$, in this case $\mu_0 = \delta_0^5$; therefore, the scaling of $t_{v,0}$ with the stellar mass is determined by the relative values of $\lambda_{m,0}$ and $\lambda_{acc,0}$. The main results of Somigliana et al. (2022) are that (i) viscous evolution maintains the power-law shape of the correlations between the stellar mass and the disc parameters, however (ii) the slope of said correlations may evolve with time, depending on the initial conditions. This is because in a purely viscous framework, the $M_d - \dot{M}$ correlation is bound to reach a linear correlation with slope unity (Lodato et al. 2017; Rosotti et al. 2017), which implies the two quantities to have the same dependence on the stellar mass, as $M_d/\dot{M} \propto M_\star^{\lambda_m - \lambda_{acc}}$. Therefore, $\lambda_m$ and $\lambda_{acc}$ must eventually reach the same value, determined by the initial conditions as

$$ \lambda_{m,evo} = \lambda_{acc,evo} = \frac{3\lambda_{m,0} - \lambda_{acc,0}}{2}. \quad (6) $$

Depending on the sign of $\delta_0 = \lambda_{m,0} - \lambda_{acc,0}$, the initial slopes can either

- ‘steepen’, that is, $\lambda_{m,evo} > \lambda_{m,0}$, if $\delta_0 > 0$ (implying also $\lambda_{acc,evo} > \lambda_{acc,0}$);

- ‘flatten’, that is, $\lambda_{m,evo} < \lambda_{m,0}$, if $\delta_0 < 0$ (implying also $\lambda_{acc,evo} < \lambda_{acc,0}$);

- ‘remain constant’, that is, $\lambda_{m,evo} = \lambda_{m,0}$, if $\delta_0 = 0$ (implying also $\lambda_{acc,evo} = \lambda_{acc,0}$).

Because in the viscous case $\delta_0 = \mu_0$, where we note again that $\mu_0$ is the slope of the correlation between the viscous timescale and the stellar mass ($t_{v,0} \propto M_\star^{\mu_0}$), we can also interpret these

<footnote n="5">5 The definition of $\delta_0$ might seem redundant at this stage, but it will become important in the following discussion.</footnote>
scenarios from the viscous timescale perspective. If $\mu_0 > 0$, meaning that the viscous timescale increases with the stellar mass, discs around less massive stars will have shorter viscous timescales, which leads to a faster evolution, compared to discs around more massive stars, which will in turn have longer viscous timescales. This uneven evolution across the stellar mass spectrum leads to a steepening of the linear correlation, as is visualised by Somigliana et al. (2022) in Figure 1. The same reasoning, but with an opposite or constant trend, applies to the other two scenarios.

### 3.2. Hybrid model – $\omega = 0$

In the hybrid viscous and MHD winds model, the general analytical solution by Tabone et al. (2022a) gives

$$M_d(t) = M_0 \left( 1 + \frac{t}{(1 + \psi)t_{\text{acc},0}} \right)^{-(\psi + 2\xi + 1)/2}, \quad (7)$$

$$\dot{M}(t) = \dot{M}_0 \left( 1 + \frac{t}{(1 + \psi)t_{\text{acc},0}} \right)^{-(\psi + 4\xi + 3)/2}, \quad (8)$$

where $\dot{M}_0$ is defined as

$$\dot{M}_0 = \frac{\psi + 1 + 2\xi}{\psi + 1} \frac{M_0}{2t_{\text{acc},0}} \frac{1}{(1 + f_{M_0})}; \quad (9)$$

in this notation, $\psi = \alpha_{\text{DW}}/\alpha_{\text{SS}}$ represents the relative strength of MHD winds and viscosity,

$$\xi = \frac{1}{4}(\psi + 1) \left[ \sqrt{1 + \frac{4\psi}{(\lambda - 1)(\psi + 1)^2}} - 1 \right]$$

is the mass ejection index quantifying the local mass loss rate to the local accretion rate, and $f_{M,0} = (R_{c,0}/R_{in})^\xi - 1$ the dimensionless mass ejection-to-accretion ratio (with $R_{in}$ initial disc radius). If we neglect the MHD-driven mass loss ($\psi \ll 1$ and $\xi \ll 1$, which correspond to $f_{M,0} \ll 1$ as well), Equations (7) and (8) reduce to the viscous case; on the other hand, if mass loss is included, it depends on the radial extent of the disc through $f_{M,0} + 1$ – which has an impact on the initial accretion rate (Equation (9)). Because the accretion timescale $t_{\text{acc}}$ is a generalisation of $t_v$ in the MHD winds framework, the dependence of the two timescales on the stellar mass is exactly equivalent, and we keep the same notation as above: $t_{\text{acc}} \propto M_*^\mu$. However, as mentioned above $\dot{M}_0$ depends on the stellar mass not only through $M_0$ and $t_{\text{acc},0}$ as in the viscous case, but also through $f_{M,0} + 1$. As $f_{M,0} + 1 \propto R_{c,0}^\xi$, and $R_{c,0} \propto M_*^{\xi_0}$, the additional dependence will have a slope of $\xi_0\xi$. Therefore, in the MHD winds scenario we can link $\delta_0$ with $\mu_0$ as $\delta_0 = \mu_0 + \xi_0\xi$. The practical meaning of this difference is that, while in the viscous scenario the difference in slope between the two correlations depends only on the scaling of the viscous timescale with the stellar mass, in the hybrid scenario it depends also on the scaling between the disc radius and the stellar mass. It is important to note that $\xi$ is a small number, typically of the order of ~0.1, therefore the difference between the viscous and hybrid case is not particularly prominent. For evolved populations, the disc mass and accretion rate read (Tabone et al. 2022a)

$$M_d(t \gg t_{\text{acc}}) \sim M_0 \left( \frac{t}{t_{\text{acc},0}} \right)^{-(\psi + 2\xi + 1)/2}, \quad (10)$$

$$\dot{M}(t \gg t_{\text{acc}}) \sim \dot{M}_0 \left( \frac{t}{t_{\text{acc},0}} \right)^{-(\psi + 4\xi + 3)/2}; \quad (11)$$

this brings the evolved slopes $\lambda_{\text{m,evo}}$ and $\lambda_{\text{acc,evo}}$ to

$$\begin{align}
\lambda_{\text{m,evo}} &= \lambda_{\text{m,0}} + \frac{1}{2}\mu_0(\psi + 2\xi + 1), \tag{12} \\
\lambda_{\text{acc,evo}} &= \lambda_{\text{acc,0}} + \frac{1}{2}\mu_0(\psi + 4\xi + 3),
\end{align}$$

which reduces to Equation (6) in the viscous case ($\psi \ll 1$, $\xi \ll 1$). Similarly to viscosity, a hybrid secular evolution maintains the power-law shape of the correlation; moreover, Equation (12) provides a theoretical prediction for the evolved slopes: they can steepen, flatten or remain the same as the initial conditions, depending on the involved parameters. As the terms in parentheses in Equation (12) are sums of positive values, the sign of the evolved slopes depends on the sign of $\mu_0$ as in the viscous case. However, there is a difference from the viscous case: as in the hybrid scenario $\mu_0$ and $\delta_0$ do not coincide anymore, a constraint on the value of $\mu_0$ is translated into a constraint on $\delta_0 - \xi_0\xi$. In particular, the slopes will increase if $\delta_0 > \xi_0\xi$ (corresponding to $\mu_0 > 0$), whereas if $\delta_0 < \xi_0\xi$ (corresponding to $\mu_0 < 0$), the slopes will decrease; and finally, if $\delta_0 = \xi_0\xi$ (corresponding to $\mu_0 = 0$) the slopes will remain constant in time. Another difference from the viscous scenario is that $\lambda_{\text{m,evo}}$ and $\lambda_{\text{acc,evo}}$ are not expected to reach the same value anymore. The limit difference is given by $\delta_{\text{evo}} = \lambda_{\text{m,evo}} - \lambda_{\text{acc,evo}}$: substituting the values from Equation (12) one finds $\delta_{\text{evo}} = \delta_0 - \mu_0(\xi + 1)$, which can be further reduced to $\delta_{\text{evo}} = \xi(2\beta + \frac{1}{2})$ by using the definition of $\mu_0$ as the slope of the $t_{\text{acc},0} - M_*^\mu$ correlation. In this expression, $\beta$ is the slope of the correlation between the disc aspect ratio and the stellar mass, which comes from definition of the accretion timescale; with a standard $\beta = -1/2$ (a reasonable approximation of the value derived by radiative transfer simulation, see e.g. Sinclair et al. 2020), we obtain $\delta_{\text{evo}} = -\xi/2$. As $\xi$ is positive by definition, in the hybrid scenario $\delta_{\text{evo}}$ is always negative, meaning that in an evolved population, the correlation between the accretion rate and the stellar mass will necessarily be steeper than that between the disc mass and the stellar mass. However, we stress once more that $\xi$ is a small number and therefore the predicted difference is also small. Figure 1 shows the evolution of the slopes from Diskpop simulations (with no spread in the initial conditions) for the hybrid model (coloured) compared with the viscous case (grey), which matches the theoretical expectations discussed above. The list of parameters used in the simulations is available in Table A.1.

Summarising, both the hybrid and viscous secular evolution preserve the power-law shape of the correlations between the disc properties and the stellar mass. The main difference is that the hybrid model does not predict the slopes of the $M_d-M_*$ and $\dot{M}-M_*$ correlations to reach the same limit value (unlike the viscous case). The predicted difference in the evolved slopes is given by $\delta_{\text{evo}} = -\xi/2$ (of the order of 0.1). The current observational uncertainties on $\delta_{\text{evo}}$ range from 0.4 to 0.8 (Testi et al. 2022, see Table 3), 4 to 8 times larger than the predicted difference, making it not observable at this stage.

### 3.3. Pure wind – $\omega \neq 0$

The solution to the pure wind model (i.e. time-dependent $\alpha_{\text{DW}}$ through $\alpha_{\text{DW}}(t) \propto \Sigma_c(t)^{-\omega}$) by Tabone et al. (2022a) gives

$$M_d(t) = M_0 \left( 1 - \frac{\omega}{2t_{\text{acc},0}} t \right)^{1/\omega}, \quad (13)$$
<figure><img src="image_1.png" /><figcaption>Fig. 1. Time evolution of the slopes of the $M_d - M_\star$ and $\dot{M} - M_\star$ correlations, $\lambda_m$ and $\lambda_{acc}$ (blue and orange solid line, respectively) in the hybrid scenario with $\alpha_{SS} = \alpha_{DW} = 10^{-3}$ ($\psi = 1$), $\lambda = 3$, $\beta = 0.5$, resulting in $\xi = 0.11$. The dashed lines represent the expected evolved value of both slopes, as in Equation (12). For comparison, we include the viscous evolution as well, represented by grey solid (actual evolution) and dashed (expected evolved value) lines (see Somigliana et al. 2022 for a detailed discussion). The three panels show different values of $\mu_0$, slope of the $t_{acc,0} - M_\star$ correlation, which is directly linked to the difference between $\lambda_{m,0}$ and $\lambda_{acc,0}$ (see text for details): we expect the slopes to (i) decrease if $\mu_0 < 0$ (left panel), (ii) remain constant if $\mu_0 = 0$ (central panel), and (iii) increase if $\mu_0 > 0$ (right panel). Contrary to the viscous case, the slopes are not expected to reach the same value anymore, but rather settle to a constant difference of $-\xi/2$. This difference is always negative, meaning that the evolved $\dot{M} - M_\star$ correlation is always steeper than that of $M_d - M_\star$ (explaining the lines crossing in the right panel).</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2. Time evolution of $M_d - M_\star$ (left) and $\dot{M} - M_\star$ (right) in the pure wind model ($\omega \neq 0$). This plot is obtained with disc population synthesis modelling, without any spread in the initial conditions. Each dot represents a disc in the population at different ages as shown in the colour bar. The initial power-law correlation, shown in light blue, is lost as early as ~1 Myr (corresponding to ~2 $< t_{acc,0} >$ with these parameters) due to a downward bending corresponding to lower stellar masses. In this simulation, we have used $N = 100$ discs, $\alpha_{DW} = 10^{-3}$, $\lambda = 3$, $\omega = 0.25$, $\lambda_{m,0} = 2.1$, $\lambda_{acc,0} = 1.5$. The set of MHD parameters is based on Tabone et al. (2022a).</figcaption></figure>

$$ \dot{M}(t) = \frac{M_0}{2t_{acc,0}(1+f_{M,0})} \left(1 - \frac{\omega}{2t_{acc,0}}t\right)^{-1+1/\omega} ; \qquad (14) $$

in this case, the functional form of Equations (13) and (14) does not allow us to derive a simple analytical expression for $M_d(t \gg t_{acc})$ and $\dot{M}(t \gg t_{acc})$. Therefore, to explore the evolution of the correlations in the pure wind case, we fully rely on Diskpop simulations. In order to account for the impact of secular evolution only, we input perfect correlations between the disc properties and the stellar mass – that is, we do not include any spread in the initial conditions. Figure 2 shows the time evolution of $M_d$ (left panel) and $\dot{M}$ (right panel) as a function of the stellar mass, from younger (darker) to older (lighter) populations. The

input power-law correlation (light blue) corresponds to a line in the logarithmic plane; however, as early as ~1 Myr (corresponding to ~2 $< t_{acc,0} >$ for this simulation), the input correlation starts to bend downwards at lower stellar masses. This behaviour reveals a significantly different trend from the viscous and hybrid model: in the pure wind scenario, the initial power-law shape of the correlations is not preserved by the secular evolution, but rather broken. In Figure 2, the faster evolution of discs around lower mass stars is the consequence of a positive $\mu_0$, implying a positive correlation between the stellar mass and the accretion timescale; a negative correlation between $M_\star$ and $t_{acc,0}$ (i.e. a negative $\mu_0$) would lead to faster evolution of discs around higher mass stars, causing the correlation to bend towards the other end of the stellar mass spectrum (see Figure B.1). Another parameter
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Same as Figure 2 but with the addition of a spread in the initial correlations between the disc properties and the stellar mass ($\sigma_{M_d}(0) = 0.65$ dex, $\sigma_{R_d}(0) = 0.52$ dex). Despite the linear correlation being readily broken in theory (see Section 3.3), the scatter introduced by the spread in the initial conditions simulates the correlation also at more evolved ages. As an example, we show the fitted line at 5 Myr (the age of the oldest observed population) in both panels ($\log_{10}(M_d/M_\odot) = 3.1 \log_{10}(M_\star/M_\odot) - 3.3, \log_{10}(\dot{M}/M_\odot\text{yr}^{-1}) = 2.5 \log_{10}(M_\star/M_\odot) - 10.4$).</figcaption></figure>

that might impact the evolution of the correlation is the dispersal timescale, $t_{disp}$; in the wind-driven case, $t_{disp} \propto t_{acc}$, therefore $t_{disp}$ does not introduce any further dependence on the stellar mass.

# 4. Population synthesis

In Section 3 we have discussed analytical trends, and presented simulations with no spread to analyse the effect of secular disc evolution alone on the evolution of the slopes. In order to test our theoretical predictions against observational data, we need to account for both a spread in the initial conditions and disc dispersal mechanisms. In this Section, we discuss the impact of both factors on the $M_d-M_\star$ and $\dot{M}-M_\star$ slopes and run realistic population synthesis simulations, to determine whether the model-dependent evolutionary features described in Section 3 would be observable with the currently available data.

## 4.1. Effects of a spread in the initial conditions

The introduction of an observationally motivated spread in the initial conditions is crucial to produce realistic population synthesis model. In the purely viscous case, Somigliana et al. (2022) have shown how the spread does not significantly impact the evolution of either the $M_d-M_\star$ or $\dot{M}-M_\star$ correlations; the shape of the curves (grey in Figure 1) is unaffected, except for their starting point, and the statistical fluctuation – determined as the interval between the 25th and 75th percentile out of 100 realisations of numerical simulations with the same initial conditions – is of the order ~0.1 for both slopes. This is a factor two less than the smallest observational uncertainty, and therefore does not produce a detectable difference in the predicted results.

Following Somigliana et al. (2022), we set $\sigma_{M_d}(0) = 0.65$ dex and $\sigma_R(0) = 0.52$ dex (determined from Ansdell et al. 2017 and Testi et al. 2022) for the log-normal distributions of $M_d(0)$ and $R_d(0)$ in the hybrid scenario (with $\alpha_{SS} = \alpha_{DW} = 10^{-3}$ hence $\psi = 1, \lambda = 3, \omega = 0$). We find that, just as in the purely viscous case, a spread in the initial conditions only shifts the starting point of the curves (coloured in Figure 1) and does not have any significant effect on the shape of the evolution of the slopes (see Figures 5–7 in Somigliana et al. 2022). The statistical fluctuation for both slopes is again of the order of 0.1, and therefore below the observational error and not impacting our predictions.

On the other hand, wind-driven models with increasing $\alpha_{DW}$ in time and a spread in the initial conditions (Figure 3) behave quite differently from the theoretical expectation discussed in Section 3.3. As the spread introduces a stochastic component, the discs will have higher or lower masses and accretion rates with equal probability; the practical result for the initial correlations is that the bending towards lower stellar masses (approximately $\log_{10}(M_\star/M_\odot) < 0.5$ in Figure 2) is lost to the stochastic displacement of the discs in the $M_d-M_\star$ and $\dot{M}-M_\star$ planes. Actually, the resulting distribution of discs both in the $M_d-M_\star$ and $\dot{M}-M_\star$ log-log plane does simulate a linear correlation; this implies that, while the stellar mass and the disc properties ‘should not’ exhibit a linear correlation already after a few $t_{acc}$, the presence of a spread mimics such correlation, making the wind-driven scenario indistinguishable from the viscous and hybrid ones. The one feature that remains observable, despite the presence of a spread, is the removal of discs around more or less massive stars – depending on the value of $\mu_0$, as discussed in Section 3.3; in the simulation shown in Figure 3 we have set $\mu_0 > 0$, implying that discs around less massive stars evolve more rapidly and are therefore more readily dispersed, as is visualised by the lack of sources around lower stellar masses at evolved ages.

Summarising, an initial power-law correlation between the disc properties and the stellar mass would keep its power-law shape under wind-driven evolution, similarly to the viscous or hybrid case; however, the interpretation of the observed correlations is different depending on the theoretical framework. While the viscous and hybrid models ‘preserve’ an initially established correlation, making the characteristic evolution of the slopes $\lambda_m$ and $\lambda_{acc}$ a tracer of disc evolution itself, the apparent correlation observed in wind-driven populations is merely a signature of the initial spread, rather than the evolutionary mechanism at play.

## 4.2. Accounting for disc dispersal: Internal photoevaporation

Out of the three theoretical scenarios discussed so far, wind-driven evolution is the only one that manages to reproduce the disc and accretion fraction (as measured by Hernández et al. 2007 and Fedele et al. 2010 respectively) – and therefore, the only one whose predictions can reasonably be compared with observations. Traditionally, the problem of disc dispersal in viscous
populations is addressed by including internal photoevaporation (see e.g. Hollenbach et al. 1994; Clarke et al. 2001; Alexander et al. 2006a,b); in this Section, we discuss the impact of internal photoevaporation on the previously described expectations for the evolution of the slopes in the purely viscous scenario. As the lack of analytical solutions to the general equation (1) does not allow for analytical arguments, we base the following discussion on physical considerations.

Internal photoevaporation is a threshold process, that kicks in after the accretion rate drops below the photoevaporative mass-loss rate (Clarke et al. 2001). The moment where the effect of photoevaporation becomes non-negligible depends therefore on the initial accretion rate: assuming for simplicity a fixed photo-evaporation rate for the whole population, as the accretion rate scales positively with the stellar mass we can expect discs around lower mass stars to show the effects of photoevaporation earlier. Moreover, given that the disc mass also scales positively with the stellar mass, said sources correspond to the less massive ones in the population. From these considerations, we can expect discs with lower initial mass to be the first ones to be affected by photoevaporation, causing a steepening of the $M_d-M_\star$ and $\dot{M}-M_\star$ correlations. However, photoevaporation also disperses discs: removing sources from the population may alter the expected behaviour, hence the need to perform numerical simulations to understand the evolution of the correlations for a population of discs undergoing internal photoevaporation. Our simulations remove discs either because the photoevaporative gap becomes too large, or because the disc mass or accretion rate fall below a certain detectability threshold. As we mentioned above, the dispersal timescale $t_{disp}$ might also play a role, if it has a different scaling with the stellar mass with respect to the accretion timescale ($t_v$ in the viscous case). Our numerical implementation of internal photoevaporation follows Owen et al. (2010), and we assume a mass-loss rate of $4 \times 10^{-10} M_\odot \text{yr}^{-1}$ – which allows us to reproduce the observed disc fraction for the set of parameters of our simulation – for all discs in the population; therefore, no further $M_\star$ dependence is introduced, and the $t_v - M_\star$ scaling is the only one that matters. We stress that a stellar mass-dependent photoevaporative rate ‘is’ expected (Picogna et al. 2021): we explore the influence of such dependence on disc observables in an upcoming work (Malanga et al., in prep.). We discuss the results of our Diskpop simulations in the following Section.

### 4.3. What the slopes are tracing

Figure 4 shows the time evolution of the $M_d-M_\star$ slope for 15 realisations with the same initial conditions for the viscous plus photoevaporative (central panel) and wind-driven (right panel) models, both reproducing the observed disc and accretion fractions, compared to the purely viscous scenario (left panel). The initial number of discs in the populations was determined to recover a certain number of objects at 5 Myr (increasing from top to bottom), and varies in the different simulations, as the decline (or lack thereof) of the disc fraction is model-dependent (see Figure 6 of Somigliana et al. 2023).

The number of objects in the simulation displayed in the first row was set to obtain ~30 discs at 5 Myr, corresponding to the currently available sample size in Upper Sco, the oldest observed star-forming region. In the left panel, we see how the evolution of the slope in the purely viscous model is not significantly affected by a spread in the initial conditions: the single realisations resemble each other remarkably well, the only difference being the starting point of the curve (as found by Somigliana et al. 2022).

On the other hand, the photoevaporative and wind-driven models have a dissimilar behaviour: each realisation can deviate substantially from the others, as we can particularly notice by the location, amplitude and direction of the bumps. The key difference between these models and the purely viscous case is disc dispersal: the stochastic nature of the slope evolution suggests that it does not trace the underlying secular disc evolution, as in the viscous scenario, but rather carry the signatures of disc dispersal itself – making it impossible to use the evolution of the slopes as a proxy for accretion mechanisms. There are two main factor that play a role in this context:

1. Initial conditions and spread in the correlations. The exact evolution of the slopes will depend on the initial conditions, both of the disc mass-stellar mass correlations themselves and of the population-wide parameters. Furthermore, the removal of discs from the population would not impact the results of the fitting procedure only if there was perfect correlations between the disc parameters and the stellar mass; with a spread in the initial correlations, on the other hand, the results may differ depending on ‘which’ discs in the population are dispersed;

2. Small number statistics. Depending on the initial number of objects, disc removal can lead to small samples – so small that it might lead to low number statistics issues. This is the case for the top row of Figure 4, where the number of objects at 10 Myr is of the order of 10 or lower.

In this work, we have used one specific set of parameters (summarised in Table A.1), determined following Somigliana et al. (2022) (viscous model) and Tabone et al. (2022a) (hybrid and wind-driven), and we leave a deeper exploration of the parameters space to a future work. While the exact shape of the slope evolution, and therefore the accretion model signature, might depend on the initial conditions, the top panel of Figure 4 shows that with the current sample sizes the noise dominates over the physical evolution. However, the currently available sample of the oldest star-forming region (Upper Sco) with both disc masses (derived from the millimetre flux, Barenfeld et al. 2016) and accretion rate (Manara et al. 2020) estimates is highly incomplete; it is therefore worth investigating whether a higher level of completeness would help reducing the entity of the oscillations, allowing to disentangle between the different evolutionary models.

The central and bottom rows of Figure 4 show how a larger sample would impact the oscillations of the slope evolution. The simulations in the middle row are performed imposing a double sample size at 5 Myr with respect to the current one (~60 discs), while in the bottom row we assume to have the complete Upper Sco sample, totalling ~140 discs (Carpenter et al., in prep.). We remind that we focus on Upper Sco as the oldest observed star-forming region, which makes it the most affected by disc dispersal.

As expected, statistical significance increases with a larger sample, leading to a decreased impact of the oscillations on the global slope evolution; with the complete sample, in particular, we can reduce the spread in the evolution by a factor of ~2 compared to the current available data. This argument confirms the importance of larger sample sizes in discriminating between the viscous and wind-driven models, as already suggested by Alexander et al. (2023) in the context of the accretion rates distribution.

As we mentioned in Section 4.2, our simulations consider discs as dispersed if their masses or accretion rates fall below the imposed detectability threshold of $10^{-12} M_\odot \text{yr}^{-1}$. We have also included a threshold in accretion rates in post-processing to
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Time evolution of the slope of the $M_d-M_\star$ correlation for 15 statistical realisations in the purely viscous (left, $\alpha_{SS} = 10^{-3}$), photoevaporative (centre, $M_w = 4 \times 10^{-10} M_\odot$/yr, $\alpha_{SS} = 10^{-3}$), and wind-driven (right, $\omega = 0.25$, $\lambda = 3$, $\alpha_{DW} = 10^{-3}$) model. The dashed lines show the median evolution, while the shaded area represents the interval between the 25th and the 75th percentiles. The three rows show different sample sizes, increasing from top to bottom. The initial size of each population was chosen to obtain a certain number of discs at 5 Myr (the age of the oldest observed population, Upper Sco) with the different disc fractions. In the top row, we match the current size of the Upper Sco sample (~30 objects, Testi et al. 2022) with both accretion rate and disc mass measurements; the middle row shows double the current sample size (~60 objects), while the bottom row assumes a complete sample (~140 objects). While the viscous model produces a remarkably similar evolution for all simulations, the latter two show stochastic oscillations from one realisation to another, suggesting that disc dispersal impacts the observed slope more than the evolutionary model does – at least with the currently available sample sizes; increasing the number of sources significantly mitigates the oscillations. The slope of the $\dot{M}-M_\star$ correlation behaves the same way.</figcaption></figure>

account for non-accreting objects. From the observational point of view, this latter selection depends on both the instrumental sensitivity and the definition of disc itself: how Class II objects are defined, and in turn how Class III sources are removed from the observed samples, impacts the resulting slope. Summarising, with the current sample sizes, the evolution of the slopes is significantly more affected by disc dispersal than it is by secular evolution; therefore, at the state of the art, it cannot be used as a proxy to disentangle between the different evolutionary models. Increasing the sample size would allow the effect of low number statistics to be reduced, potentially allowing the different evolution of the slopes under the two evolutionary mechanisms to be observed; we further discuss this possibility in the following Section.

**5. The observational relevance of the slopes**

Observed star-forming regions have both a spread in the initial conditions in addition to some disc dispersal mechanism

(be it photoevaporation or MHD winds); as we discussed in Section 4.3, with the current sample sizes, the statistical significance of the observationally determined slopes is undermined and their evolution traces disc dispersal, rather than the accretion mechanism. In this Section, we perform a statistical analysis of our simulated slopes and compare them with the currently available measurements; furthermore, we show the relevance of measuring the slopes despite these limitations and discuss the conditions under which they allow us to put constraints on disc evolution.

**5.1. Comparison of different evolutionary models**

Figure 5 shows the comparison of the evolution of the slopes between the viscous + photoevaporative (solid line) and wind-driven (dashed line) models for both the $M_d - M_\star$ and $\dot{M}-M_\star$ correlations (top and bottom row, respectively), including the measured slopes in four star-forming regions from Testi et al. (2022) as grey dots. As we mentioned above, both models lead
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Comparison of the time evolution of $\lambda_m$ (top) and $\lambda_{acc}$ (bottom) between the viscous+photoevaporative and wind-driven case including measured slopes from four star-forming regions. To account for statistical fluctuations, each simulation combines 100 realisations of the same initial conditions: the lines show the median evolution, while the shaded area represent the interval between the 25th and 75th percentiles. The simulations in the three columns differ for the initial number of discs, determined to obtain a specific sample size at 5 Myr – currently available sample (left), double the currently available sample (centre) and the complete sample (right). Observed slopes from Testi et al. (2022).</figcaption></figure>

to disc dispersal consistently with the observed disc and accretion fraction (shown in Figure 6 of Somigliana et al. 2023); the three columns show simulations performed with a different initial number of discs, increasing from left to right, to obtain a different sample size at 5 Myr according to the predicted decline of observed discs. As in Figure 4, the number of objects at 5 Myr is ~30, ~60 and ~140 from left to right, increasing from the currently available measurements in Upper Sco to the virtually complete sample. To estimate the effect of statistical fluctuations, given by the spread in the initial conditions, we ran 100 simulations for each set up: the solid and dashed lines represent the median evolution, while the intervals between the 25th and 75th percentile are visualised by the shaded areas. The growing shaded area, particularly visible with smaller sample sizes, is representative of the decreasing amount of sources on which the fit is performed: with the current sample size (left column), that leads to ~30 discs at 5 Myr, we end up with a $1\sigma$ deviation from the median value of ~0.5–0.6. Larger sample sizes significantly reduce the scatter, leading to $\sigma \sim 0.4$ with a double sample and $\sigma \sim 0.2$ for the complete sample, reducing the current one by a factor 3. As mentioned in Section 4.2, with the currently available number of objects the dominant role in the evolution of the slopes is played by disc dispersal, which makes it difficult to trace the imprint of the secular evolution. The expected Upper Sco complete sample (right column) allows for a better separation between the two models – particularly for the $M_d - M_\star$ correlation: the expected slope in the two scenarios differs by ~0.5, while the typical uncertainty of the currently measured slopes is between 0.2 and 0.3. Larger sample sizes would further

decrease this uncertainty, allowing us to discriminate between the two models based on the slope evolution.

The observed slopes (from Testi et al. 2022) are only included in the left column of Figure 5 as they refer to the current sample size. The main source of uncertainty in the current measurements is given by Upper Sco, mainly due to the incomplete sample; moreover, it is worth pointing out that external photoevaporation is likely to play a significant role in this region (Anania et al., in prep.). This comparison is meant as a first glance of the parameters space of the observed slopes, and we anticipate a proper exploration of the initial conditions once the full sample will be available. In the following Section, we discuss the other constraints that we can put on disc evolution, besides identifying the driving accretion mechanism.

## 5.2. What the slopes are 'really' tracing

Despite not allowing to conclusively discriminate between different evolutionary scenarios with the current sample sizes, the slopes of the $M_d-M_\star$ and $M-M_\star$ correlations can still help with constraining other properties from the theoretical considerations presented in Section 3, which we summarise in Table 2. If we 'assume' an evolutionary model to begin with, and we can estimate (directly or indirectly) either $\mu_0$ or $\delta_0$, we can constrain the other parameter. When discussing the observational determination of what we have so far referred to as 'initial conditions', it is important to clarify the meaning of 'initial'. Diskpop deals with and evolves Class II, potentially Class III, objects; hence, the initial conditions we input refer to the beginning of the Class II
<table><caption>**Table 2.** Summary of the different possible theoretical scenarios described in Section 3 to visualise the relative signs of the parameters at play.</caption><thead><tr><th colspan="8">μ<sub>0</sub> = δ<sub>0</sub> - ζ<sub>0</sub>ξ</th></tr><tr><th>ξ</th><th>μ<sub>0</sub></th><th>t<sub>acc,0</sub> - M<sub>★</sub></th><th>ζ<sub>0</sub></th><th>R<sub>d</sub> - M<sub>★</sub></th><th>δ<sub>0</sub></th><th>M<sub>d</sub> - M<sub>★</sub> (a) and M&#x0303; - M<sub>★</sub> (b)</th><th>⇒</th></tr></thead><tbody><tr><td rowspan="3">0</td><td>&lt; 0</td><td>t<sub>acc,0</sub> ↓ M<sub>★</sub> ↑</td><td>Any</td><td>Any</td><td>&lt; 0</td><td>(a) shallower than (b)</td><td>=</td></tr><tr><td>0</td><td>t<sub>acc,0</sub> ↔ M<sub>★</sub> ↑</td><td>Any</td><td>Any</td><td>0</td><td>(a) as steep as (b)</td><td>=</td></tr><tr><td>&gt; 0</td><td>t<sub>acc,0</sub> ↑ M<sub>★</sub> ↑</td><td>Any</td><td>Any</td><td>&gt; 0</td><td>(a) steeper than (b)</td><td>=</td></tr><tr><td rowspan="6">(0, 1)</td><td rowspan="2">&lt; 0</td><td rowspan="2">t<sub>acc,0</sub> ↓ M<sub>★</sub> ↑</td><td>≤ 0</td><td>R<sub>d</sub> ↓↔ M<sub>★</sub> ↑</td><td>≤ 0</td><td>(a) shallower or as steep as (b)</td><td>=</td></tr><tr><td>&gt; 0</td><td>R<sub>d</sub> ↑ M<sub>★</sub> ↑</td><td>(-∞, ζ<sub>0</sub>ξ)</td><td>(a) shallower or steeper than (b)<sup>(a)</sup></td><td>≠</td></tr><tr><td rowspan="2">0</td><td rowspan="2">t<sub>acc,0</sub> ↔ M<sub>★</sub> ↑</td><td>≤ 0</td><td>R<sub>d</sub> ↓↔ M<sub>★</sub> ↑</td><td>≤ 0</td><td>(a) shallower or as steep as (b)</td><td>○</td></tr><tr><td>&gt; 0</td><td>R<sub>d</sub> ↑ M<sub>★</sub> ↑</td><td>≥ 0</td><td>(a) steeper or as steep as (b)</td><td>=</td></tr><tr><td rowspan="2">&gt; 0</td><td rowspan="2">t<sub>acc,0</sub> ↑ M<sub>★</sub> ↑</td><td>≤ 0</td><td>R<sub>d</sub> ↓↔ M<sub>★</sub> ↑</td><td>(ζ<sub>0</sub>ξ, +∞)</td><td>(a) shallower or steeper than (b)<sup>(b)</sup></td><td>○</td></tr><tr><td>&gt; 0</td><td>R<sub>d</sub> ↑ M<sub>★</sub> ↑</td><td>(ζ<sub>0</sub>ξ, +∞)</td><td>(a) steeper than (b)</td><td>=</td></tr></tbody><tfoot><tr><td colspan="8">δ<sub>evo</sub> = ξ(ζ<sub>0</sub> - μ<sub>0</sub>), ξ ≠ 0</td></tr><tr><td colspan="8">M<sub>d</sub> - M<sub>★</sub> (a) and M&#x0303; - M<sub>★</sub> (b)</td></tr><tr><td colspan="8">&lt; 0 (a) shallower than (b)</td></tr><tr><td colspan="8">0 (a) as steep as (b)</td></tr><tr><td colspan="8">&gt; 0 (a) steeper than (b)</td></tr></tfoot></table>

**Notes.** From left to right in the top table, the columns show (i) $\xi$, a proxy for the evolutionary model (viscous if $\xi = 0$, hybrid or wind-driven otherwise); (ii) $\mu_0$, the slope of the $t_{acc,0} - M_\star$ correlation and its implication on the correlation itself; (iii) $\zeta_0$, the slope of the $R_d - M_\star$ correlation, and its implication on the correlation itself; (iv) $\delta_0 = \lambda_{m,0} - \lambda_{acc,0}$, the difference between the initial slopes of the $M_d - M_\star$ and $M\ddot{M} - M_\star$ correlations and its implication on their relative steepness. The final column summarises whether the signs of $\mu_0$ and $\delta_0$ are 'necessarily the same' ($=$), 'necessarily opposite' ($\neq$) or 'can be either' ($\circ$). When discussing the implications on correlations, up(down)wards arrows represent an in(de)crease of the parameters, while horizontal arrows describes the lack of correlation. The different cell colours are purely meant to guide the eye. The top table links the initial conditions, while the bottom table summarises the implications of the evolved difference in the slopes. $^{(a)}$ $\zeta_0\xi$ is positive, therefore $\delta_0$ can either be negative (implying (a) shallower than (b)) or positive (implying (a) steeper than (b)). $^{(b)}$ $\zeta_0\xi$ is negative, therefore the same argument as in *a* holds.

phase, where the protostellar collapse is over and the disc is already formed. From the observational point of view, this means that we expect $\delta_0$ and $\mu_0$ to refer to young Class II objects – around, or younger than, approximately 1 Myr. Earlier phases like the Class 0 and I need a dedicated study, as the accretion of the protostellar envelope is expected to play a prominent role in those stages.

In the following, we discuss the constraints we can put in both directions and comment on their feasibility based on the currently available estimates of $\mu_0$ and $\delta_0$.

### 5.2.1. Constraining $\delta_0$ from $\mu_0$

Ansdell et al. (2017) claimed $\lambda_m$ to be increasing with time. As shown by Somigliana et al. (2022) and discussed in Section 3, increasing slopes imply that discs around less massive stars evolve faster than those around more massive stars; this can be interpreted in terms of increasing accretion timescale with stellar mass, which corresponds to $\mu_0 > 0$ (with $t_{acc,0} \propto M_\star^{\mu_0}$). In this section we discuss the implications of the increasing slope scenario on the initial conditions, $\lambda_{m,0}$ and $\lambda_{acc,0}$.

The top panel of Table 2 shows the relation between $\mu_0$ and $\delta_0$ in the different evolutionary models. As $\mu_0 = \delta_0 - \zeta_0\xi$, in the viscous case (corresponding to $\xi = 0$) we have $\mu_0 = \delta_0$ as mentioned in Section 3. This means that, to recover the suggested increasing slopes scenario, $\delta_0$ necessarily needs to be positive – regardless of the value of any other parameters: this translates to the initial $M_d - M_\star$ correlation being steeper than $M\ddot{M} - M_\star$. In the hybrid and wind-driven case, instead, the implication is less straightforward as a positive $\mu_0$ can lead to opposite signs of $\delta_0$: this is determined by the scaling of the disc radius with the stellar

mass, which is suggested to be (weakly) positive from observational evidences (e.g. Andrews et al. 2018). In principle, as $\delta_0 \in (\zeta_0\xi, +\infty)$, the sign of $\zeta_0$ determines whether negative values of $\delta_0$ are possible; however, as $\xi$ is a small number (0.1 in this work), only a limited area of the parameters space would lead to a negative $\delta_0$. Summarising, if we assume increasing slopes ($\mu_0 > 0$) we can constrain the sign of $\delta_0$ regardless of the evolutionary model assumed: in both cases $\delta_0$ needs to be positive, which leads to an initially steeper $M_d - M_\star$ than $M\ddot{M} - M_\star$ correlation.

### 5.2.2. Constraining $\mu_0$ from $\delta_{evo}$

Instead of assuming increasing slopes, we can start from the currently measured values of $\delta_{evo}$ and estimate $\delta_0$ in the different evolutionary models. In the viscous case, because $\delta_{evo} = 0$, we focus on the value of the single slopes instead: as $\lambda_{m,evo} = \lambda_{acc,evo} = \delta_0/2 + \lambda_{m,0}$, the measured final value of the slopes does not help in constraining $\delta_0$ as it also depends on $\lambda_{m,0}$. In the hybrid case, instead, we have $\delta_{evo} = \xi(\zeta_0 - \mu_0)$, meaning that if we can determine the sign of $\delta_{evo}$ we can constrain that of $\mu_0$ as well. While in principle the sign of $\zeta_0$ influences that of $\mu_0$, as we mentioned above $\zeta_0$ is likely a small number: therefore, effectively, $\delta_{evo}$ and $\mu_0$ have opposite signs for the vast majority of the parameters space.

Assuming that the observed disc populations can be considered evolved enough for the above arguments to hold, we can estimate $\delta_{evo}$ from the most recent and homogeneous measurements available of $\lambda_m$ and $\lambda_{acc}$ (Testi et al. 2022). The resulting values of $\delta_{obs}$ (which we label ‘observed’ as opposed to the theoretical expectation, ‘evolved’), summarised in Table 3, are
<table><caption>**Table 3.** Values of $\delta$ derived from the currently available measurements.</caption><thead><tr><td>Region</td><td>Median age [Myr]</td><td>$\delta_{\text{obs}}$</td></tr></thead><tbody><tr><td>L1668</td><td>1</td><td>-0.3 ± 0.5</td></tr><tr><td>Lupus</td><td>2</td><td>0.1 ± 0.4</td></tr><tr><td>Cha I</td><td>2.8</td><td>-0.7 ± 0.4</td></tr><tr><td>Upper Sco</td><td>4.3</td><td>0.7 ± 0.8</td></tr></tbody></table>

**Notes.** The values of $\delta$ are computed from measurements of $\lambda_m$ and $\lambda_{acc}$ reported by Testi et al. (2022).

oscillating: out of the four regions L1668, Lupus, Chameleon I and Upper Sco, we find two positive and two negative median values. Moreover, in three cases out of four the uncertainties are so large that $\delta_{\text{obs}}$ would be compatible with both a positive and a negative value. The difficulty in assessing the sign of $\delta_{\text{obs}}$ from the current measurements of the slopes make constraining $\mu_0$ from $\delta_{\text{evo}}$ not trivial. Larger sample sizes would give a better measurement of the slopes and reduce the uncertainty, leading to a more solid determination of the sign of $\delta_{\text{obs}}$ – which would possibly allow $\mu_0$ to be constrained.

Summarising, the (admittedly not robust) observational evidence pointing towards increasing accretion timescale with stellar mass allows us to constrain the initial correlations between the stellar mass and disc parameters; regardless of the evolutionary model considered, the initial slope of the $M_d - M_\star$ correlation needs to be larger than that of $\dot{M}-M_\star$. The other way around, constraining the slope of the accretion timescale – stellar mass correlation from the difference between $\lambda_m$ and $\lambda_{acc}$ at the present time, requires sample sizes larger by at least a factor two.

## 6. Conclusions

In this paper, we have investigated the impact of disc evolution models on the correlations between the stellar mass and the disc properties – especially the disc mass and the accretion rate. We have explored the purely viscous, wind-driven, viscous and wind hybrid, and photoevaporative models. Assuming power-law correlations to hold as initial conditions, $M_d(0) \propto M_\star^{\lambda_{acc,0}}$, $\dot{M}(0) \propto M_\star^{\lambda_{acc,0}}$, we performed analytical calculations (where possible) and population synthesis simulations for both evolutionary scenarios, and compared them with the purely viscous case discussed in Somigliana et al. (2022). Our main results are the following:

1. The viscous and hybrid models change the slope of the initial correlations as function of the evolutionary time, but preserve their shape. In the wind-driven model, instead, the correlations deviate from the original power-law shape: this is visualised in the logarithmic plane as a bending of the linear correlation (see Figure 1). The bending direction is towards the less or more massive stars depending on the scaling of the accretion timescale with the stellar mass (positive and negative correlation respectively);

2. The characteristic behaviour of the slopes in the wind-driven model is concealed by the presence of a spread in the initial conditions, which introduces a scatter in the correlations and makes it no longer possible to detect the bending (Figure 2). This leads to a considerably similar evolution of the correlations in the different accretion models;

3. Performing our simulations with evolutionary models that match the disc dispersal timescales (intrinsic in the wind-driven model and including internal photoevaporation in

the viscous case), we find that the evolution of the slopes is significantly impacted by the removal of discs from the population (Figure 3). Different realisations of the same simulation dramatically differ from one another, and show a stochastic behaviour with large variations (Figure 4). This has both a physical (presence of a spread in the initial conditions) and a statistical (low number of objects left after a few Myr of evolution) reason;

4. While a proper exploration of the parameters space, outside of the scope of this work, would be needed to assess the impact of the initial conditions, with the currently available sample sizes the noise dominates over the physical evolution;

5. Increasing the sample size can mitigate the effects of disc dispersal on the evolution of the slope by removing the stochastic effects. We find that, for our parameters choice, the complete sample of Upper Sco (~140 sources) at 5 Myr would reduce the oscillations enough to make the slopes a proxy for the evolutionary model (Figure 5);

6. While the currently available sample sizes do not yet allow to distinguish between the different evolutionary models, we can use them to put some constraints on the initial conditions. We find that in all evolutionary scenarios, the observational claim of increasing slopes leads to an initially steeper correlation between the disc mass and the stellar mass than between the accretion rate and the stellar mass. The other possible way, measuring the current slopes and inferring the correlation between the accretion timescale and the stellar mass from them, provides weaker constraints because of the high uncertainties in the current measurements;

7. We have presented and released the 1D Python disc population synthesis code Diskpop and its output analysis library popcorn.

In this work, we have shown how large enough samples of protoplanetary discs can provide a way of distinguishing between the evolutionary models (with a standard set of parameters) through the observation of the time evolution of the correlations between the disc properties and the stellar mass. We have shown how the stochastic fluctuations seen with the currently available observations could be significantly reduced if we had access to the complete Upper Sco sample, consisting of approximately 140 sources at 5 Myr. We strongly support the observational effort in the direction of obtaining larger amounts of data for evolved star-forming regions, and encourage the exploration of the parameters space beyond the standard case.

**Acknowledgements.** We thank an anonymous referee for their useful comments that helped us improve the clarity of the manuscript. This work was partly supported by the Italian Ministero dell'Istruzione, Università e Ricerca through the grant Progetti Premiali 2012-iALMA (CUP C52I13000140001), by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Ref no. 325594231 FOR 2634/2 TE 1024/2-1, by the DFG Cluster of Excellence Origins (www.origins-cluster.de). This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 823823 (DUSTBUSTERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECO-GAL (grant 855130) and ERC Starting Grant DiscEvol (grant 101039651). GR acknowledges support from Fondazione Cariplo, grant No. 2022-1217. 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 Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

## References

Alcalá, J. M., Natta, A., Manara, C. F., et al. 2014, A&A, 561, A2

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20
Alexander, R. D., & Armitage, P. J. 2006, ApJ, 639, L83

Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216

Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006b, MNRAS, 369, 229

Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475

Alexander, R., Rosotti, G., Armitage, P. J., et al. 2023, MNRAS, 524, 3948

Anderson, D. E., Cleeves, L. I., Blake, G. A., et al. 2022, ApJ, 927, 229

Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157

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, AJ, 153, 240

Armitage, P. J., Simon, J. B., & Martin, R. G. 2013, ApJ, 778, L14

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

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Birstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79

Birstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

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

Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994

Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485

Coleman, G. A. L., Mroueh, J. K., & Haworth, T. J. 2024, MNRAS, 527, 7588

Dullemond, C. P., Natta, A., & Testi, L. 2006, ApJ, 645, L69

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

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Haworth, T. J., Coleman, G. A. L., Qiao, L., Sellek, A. D., & Askari, K. 2023, MNRAS, 526, 4315

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594

Hernández, J., Hartmann, L., Megeath, T., et al. 2007, ApJ, 662, 1067

Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654

Kley, W., & Lin, D. N. C. 1996, ApJ, 461, 933

Kroupa, P. 2001, MNRAS, 322, 231

Laibe, G., & Price, D. J. 2014, MNRAS, 444, 1940

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2024, A&A, 682, A30

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

Lodato, G., Rampinelli, L., Viscardi, E., et al. 2023, MNRAS, 518, 4481

Long, F., Andrews, S. M., Rosotti, G., et al. 2022, ApJ, 931, 6

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

Machida, M. N., & Hosokawa, T. 2013, MNRAS, 431, 1719

Machida, M. N., & Basu, S. 2024, arXiv e-prints [arXiv:2405.08271]

Manara, C. F., Robberto, M., Da Rio, N., et al. 2012, ApJ, 755, 154

Manara, C. F., Fedele, D., Herczeg, G. J., & Teixeira, P. S. 2016, A&A, 585, A136

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

Mauxion, J., Lesur, G., & Maret, S. 2024, A&A, 686, A253

Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 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, 501

Mohanty, S., Jayawardhana, R., & Basri, G. 2005, ApJ, 626, 498

Morbidelli, A., Lunine, J. I., O'Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Annu. Rev. Earth Planet. Sci., 40, 251

Muzerolle, J., Hillenbrand, L., Briceño, C., Calvet, N., & Hartmann, L. 2003, in Brown Dwarfs, 211, 141

Natta, A., Testi, L., Muzerolle, J., et al. 2004, A&A, 424, 603

Owen, J. E., Ercolano, B., Clarke, C. J., & Alexand er, R. D. 2010, MNRAS, 401, 1415

Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880

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

Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Popham, R., Narayan, R., Hartmann, L., & Kenyon, S. 1993, ApJ, 415, L127

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

Rigliaco, E., Natta, A., Randich, S., Testi, L., & Biazzo, K. 2011, A&A, 525, A47

Rosotti, G. P. 2023, New A Rev., 96, 101674

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

Sanchis, E., Testi, L., Natta, A., et al. 2020, A&A, 633, A114

Sellek, A. D., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 492, 1279

Shakura, N. I., & Sunyaev, R. A. 1973, in IAU Symposium, 55, X- and Gamma-Ray Astronomy, eds. H. Bradt, & R. Giacconi, 155

Sinclair, C. A., Rosotti, G. P., Juhasz, A., & Clarke, C. J. 2020, MNRAS, 493, 3535

Somigliana, A., Toci, C., Lodato, G., Rosotti, G., & Manara, C. F. 2020, MNRAS, 492, 1120

Somigliana, A., Toci, C., Rosotti, G., et al. 2022, MNRAS, 514, 5927

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

Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414

Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111

Testi, L., Natta, A., Manara, C. F., et al. 2022, A&A, 663, A98

Trapman, L., Zhang, K., van't Hoff, M. L. R., Hogerheijde, M. R., & Bergin, E. A. 2022, ApJ, 926, L2

Trapman, L., Rosotti, G., Zhang, K., & Tabone, B. 2023, ApJ, 954, 41

Venuti, L., Stelzer, B., Alcalá, J. M., et al. 2019, A&A, 632, A46

Veronesi, B., Paneque-Carreño, T., Lodato, G., et al. 2021, ApJ, 914, L27

Winter, A. J., Clarke, C. J., Rosotti, G., et al. 2018, MNRAS, 478, 2700

Xu, W., & Kunz, M. W. 2021a, MNRAS, 502, 4911

Xu, W., & Kunz, M. W. 2021b, MNRAS, 508, 2142
# Appendix A: Parameters used throughout the paper

Table A.1 shows the parameters used throughout the paper. The general simulation parameters refer to the initial correlations (distribution, slopes and spreads), while the following three tables are divided by evolutionary model and display the main parameters for each of them.

<table><caption>**Table A.1.** Values of the parameters used throughout the paper, unless explicitly stated otherwise.</caption><tbody><tr><th colspan="2">General simulation parameters</th></tr><tr><td>λ<sub>m,0</sub></td><td>2.1</td></tr><tr><td>λ<sub>acc,0</sub></td><td>1.5</td></tr><tr><td>σ<sub>M<sub>d</sub></sub></td><td>0.65 dex</td></tr><tr><td>σ<sub>R<sub>d</sub></sub></td><td>0.52 dex</td></tr><tr><td>M<sub>d</sub>, R<sub>d</sub> distribution</td><td>lognormal</td></tr><tr><th colspan="2">Viscous model</th></tr><tr><td>α<sub>SS</sub></td><td>10<sup>-3</sup></td></tr><tr><th colspan="2">Photoevaporative model</th></tr><tr><td>α<sub>SS</sub></td><td>10<sup>-3</sup></td></tr><tr><td>M&#x0307;<sub>w</sub></td><td>4 × 10<sup>-10</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><th colspan="2">Wind-driven model</th></tr><tr><td>α<sub>SS</sub></td><td>10<sup>-3</sup></td></tr><tr><td>α<sub>DW</sub></td><td>10<sup>-3</sup></td></tr><tr><td>λ</td><td>3</td></tr><tr><td>ω</td><td>0.25</td></tr></tbody></table>

# Appendix B: MHD model with $\mu < 0$

As mentioned in Section 3.3, the breaking of the linear correlation between the disc properties and the stellar mass happens towards higher or lower stellar masses depending on the value of $\mu_0$. Figure 2 shows the a simulation $\mu_0 > 0$, while in Figure B.1 we show the opposite case. As $\mu_0$ is linked to $\delta_0$ through $\mu_0 = \delta_0 - \zeta_0\xi$, if $\delta_0 < \zeta_0\xi$ the resulting $\mu_0$ will be negative, leading to a specular bending of the correlation. Given that $\zeta_0\xi$ is a small number (~0.1 in our simulation), this generally requires a negative $\delta_0$. The simulation in Figure B.1 has $\lambda_{m,0} = 1.3$ and $\lambda_{acc,0} = 1.7$, resulting in $\delta_0 = -0.4$.

# Appendix C: Validation of Diskpop

Figure C.1 shows the evolution of the gas surface density as a function of the disc radius in the cases of evolution driven by (i) viscosity, (ii) viscosity and internal photoevaporation, (iii) MHD winds, and (iv) viscosity and external photoevaporation for a single disc simulated with Diskpop. The top left panel, corresponding to viscous evolution (i), shows the key feature of viscous spreading: while the global surface density decreases as a consequence of the accretion onto the central star, the radial extent of the disc increases. This is a consequence of the redistribution of angular momentum, that causes part of the disc material to move towards larger radii. The top right panel,

<figure><img src="image_6.png" /><figcaption>**Fig. B.1.** Same as Figure 2, but with a choice of initial slopes resulting in a negative $\mu_0$ ($\lambda_{m,0} = 1.3, \lambda_{acc,0} = 1.7$). The bending of the linear correlation happens towards larger disc masses, in agreement with the prediction.</figcaption></figure>

where the disc evolves under the combined effect of viscosity and internal photoevaporation (ii), shows the characteristic two-timescales behaviour (Clarke et al. 2001): the evolution is effectively viscous in the earliest stages, as long as the accretion rate is larger than the photoevaporative mass-loss rate; then, photoevaporation opens a cavity within the disc, which gets divided into an inner and an outer disc. The inner disc is less extended and therefore has a shorter viscous timescale, which means it evolves much faster and is quickly completely accreted onto the protostar; the outer disc instead keeps on evolving on timescales comparable to the original one, making photoevaporation a two-timescales process. The bottom left panel shows a disc evolving due to MHD winds (iii): the absolute value of the surface density drops faster than in the viscous model, because of the increase of $\alpha_{DW}$ in time. Furthermore, as angular momentum is removed from the wind (together with material), the disc does not spread but rather shrinks in time, as expected from the theoretical prediction (Tabone et al. 2022a). Finally, the bottom right panel shows the evolution of a disc undergoing external photoevaporation combined with viscosity (iv): the latter dominates at the earliest stages, producing the characteristic features like viscous spreading, while the effect of external photoevaporation is visible at later ages as a truncation of the disc that also halts its spreading. In this case, the disc truncation and the outside-in depletion of disc material is the consequence of the photo-dissociation of gas molecules due to the FUV radiation emitted by massive stars and experienced by the disc. The efficiency of this process depends primarily on the stellar mass and the FUV flux experienced: given a fixed FUV flux, a disc around a lower mass star will lose its material to external winds more easily compared to a disc around a higher mass star, because of the higher gravitational bond of the system. For the same reason,
<figure><img src="image_7.png" /><figcaption>Fig. C.1. Gas surface density as a function of the radius for a disc generated with Diskpop at different times as the colour bar shows. The four models are purely viscous (top left, $\alpha_{SS} = 10^{-3}$), viscous including internal photoevaporation (top right, $\alpha_{SS} = 10^{-3}$, $\dot{M}_w = 4 \times 10^{-10} M_{\odot}yr^{-1}$), wind-driven (bottom left, $\alpha_{DW} = 10^{-3}$, $\lambda = 3$, $\omega = 0.25$) and viscous including external photoevaporation (bottom right, $\alpha_{SS} = 10^{-3}$, FUV = 100 $G_0$) respectively.</figcaption></figure>

more extended and less massive discs are more prone to external truncation.

Figure C.2 shows the isochrones in the $\dot{M} - M_d$ plane at 0.1, 1 and 10 Myr for the three evolutionary models of viscosity, viscosity and photoevaporation, and MHD winds. Each dot represents a disc in the population, while the solid lines show the analytical prediction (when applicable): in the viscous case, the isochrones read (Lodato et al. 2017)

$$ \dot{M} = \frac{M_d}{2(2-\gamma)t} \left[ 1 - \left( \frac{M_d}{M_0} \right)^{2(2-\gamma)} \right], \qquad (C.1) $$

while in the MHD wind-driven scenario (Tabone et al. 2022a)

$$ \dot{M} = \frac{1}{\omega(1 + f_{M,0})t} M_d \left[ \left( \frac{M_d}{M_0} \right)^{-\omega} - 1 \right]. \qquad (C.2) $$

The left panel shows the viscous model, where the discs tend more and more towards the theoretical isochrone during their evolution (Lodato et al. 2017); the central panel includes internal photoevaporation, which has the effect of bending the isochrones once the accretion rate becomes comparable to the photoevaporative mass-loss rate (Somigliana et al. 2020); finally, the right panel shows an MHD wind-driven population, where the scatter in the $\dot{M} - M_d$ plane remains significant during the evolution - contrary to the viscous prediction (Somigliana et al. 2023).
<figure><img src="image_8.png" /><figcaption>**Fig. C.2.** Isochrones at 0.1, 1, and 10 Myr for disc populations undergoing viscous, viscous+internal photoevaporation, and wind-driven evolution (left, centre, and right panel, respectively), with the same parameters as Figure C.1. Each dot represents a disc, while the solid lines show the theoretical isochrones at the corresponding age, where available.</figcaption></figure>