© The Authors 2024

1 Introduction

The secular evolution of protoplanetary discs is deeply intertwined 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 studied integrated disc properties. Accretion rates are inferred from the spectra of the central stars, which show an excess emission (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 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_{\mathrm{d}}\propto M_{\star }{}^{{\lambda }_{\mathrm{m},\text{ 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¹ of λ_acc,obs ≈ 1.8 ± 0.2 (as first suggested by Muzerolle et al. 2003 and supported by many of the following works mentioned above), the M_d–M_* correlations is claimed to be steepening with time (Ansdell et al. 2017), from the lowest λ_m,obs(t − 1 Myr) = 1.7 ± 0.2 (Taurus) to the highest λ_m,obs(t ~ 5 Myr) = 2.4 ± 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_d–M_* does suggest an evolutionary trend; Somigliana et al. (2022) have found that, assuming power-law correlations between both M_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_d–M_* 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_d–M_* 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². We describe the master equation for the secular 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, $$\begin{array}{r}\frac{\mathrm{\partial }\mathrm{\Sigma }}{\mathrm{\partial }t}=\frac{3}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[\frac{1}{\mathrm{\Omega }r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{\alpha }_{\mathrm{S}\mathrm{S}}\mathrm{\Sigma }{c}_s{}^2\right)\right]+\frac{3}{2r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left[\frac{{\alpha }_{\mathrm{D}\mathrm{W}}\mathrm{\Sigma }{c}_s{}^2}{\mathrm{\Omega }}\right]\\ -\frac{3{\alpha }_{\mathrm{D}\mathrm{W}}\mathrm{\Sigma }{c}_s{}^2}{4(\lambda -1)r^2\mathrm{\Omega }}-{\dot{\mathrm{\Sigma }}}_{\text{photo }},\end{array}$$(1)

describes the time evolution of the gas surface density in the most general framework, where ∑ is the gas surface density, Ω the Keplerian orbital frequency, α_SS the Shakura & Sunyaev (1973) α parameter, α_DW the MHD equivalent of α_SS (Tabone et al. 2022a), c_s the sound speed, and λ 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 α_SS, (ii) the wind-driven accretion, which corresponds to an advection term, parameterised by α_DW, (iii) mass loss due to MHD disc winds, parameterised by λ 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 (α_DW = 0), purely MHD wind-driven (α_SS = 0) or hybrid (α_SS, α_DW ≠ 0) evolution, with $({\dot{\mathrm{\Sigma }}}_{\text{photo }}\ne 0)$ or without $({\dot{\mathrm{\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 α_DW isset 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 ν = α_SSc_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 (ν ∝ R^γ).

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 α_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 α_DW with time; these solutions depend on the value of ψ ≡ α_DW/α_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 α_DW with time. To obtain these, Tabone et al. (2022a) parameterised α_DW(t) ∝ ∑_c(t)^−ω, with ∑_c = M_d(t)/2πR_c²(t) (where R_c is a characteristic radius) and ω as a free parameter, and neglect the radial transport of angular momentum (α_SS = 0).

Photoevaporation. The generic ${\dot{\mathrm{\Sigma }}}_{\text{photo }}$ term in Equation (1) allows to account for photoevaporative processes, both internal and external. The exact form of ${\dot{\mathrm{\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 of a star and a disc, whose key initial parameters (stellar mass, disc mass, accretion rate, disc radius, evolutionary parameters α_SS, α_DW, λ, ω...) 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 logormal 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 (α_SS, α_DW, λ, ω): 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}{\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{s}}^{-1}}\right)}^{1.14}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}.$$

   The surface mass-loss profile ${\dot{\mathrm{\Sigma }}}_{\text{photo }}$ (Equation (B2) in Owen et al. 2012) is then scaled so that $\int 2\pi R{\dot{\mathrm{\Sigma }}}_{\text{photo}}\mathrm{d}\mathrm{R}$ 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₀³. This parameter can be set to any value accessible in the FRIEDv2 grid of mass loss rate (Haworth et al. 2023), spanning from 1 to 10⁵ G₀.

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 ∝ R^−1/2, which results in c_s ∝ R^−1/4 and H/R ∝ R^1/4. We note that this implies ν ∝ R (i.e. γ = 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 υ is ∂_tq(x, t) + υ∂_xq(x, t) = 0; in the case of wind-driven accretion, the advected quantity is R∑, while the advection (inwards) velocity is given by υ_DW = (3α_DWHc_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}}_{\mathrm{e}\mathrm{x}\mathrm{t},i}={\dot{M}}_{\mathrm{t}\mathrm{o}\mathrm{t}}\frac{M_i}{M(R>R_{\mathrm{t}})},$$(2)

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

6. Dust evolution⁴: 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 Δt = Cmin(Δx/υ) and ensures that, within one time step Δt, the material moving at velocity υ does not flow further than one grid spacing Δ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³ points with power-law spacing and exponent 1/2, extending from 3 × 10⁻³ au to 10⁴ 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 ~10⁻³ 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-togas 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.

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_*: 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 λ_m, λ_acc, and ζ respectively. The initial correlations are set as follows: $$\{\begin{array}{l}{M}_{\mathrm{d}}(0)\propto M_{\star }{}_{}^{{\lambda }_{\mathrm{m},0}},\\ \dot{M}(0)\propto M_{\star }{}^{{\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c},0}},\\ R_{\mathrm{d}}(0)\propto M_{\star }{}^{{\zeta }_0}.\end{array}$$(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 (ν ∝ R^γ) 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_{\mathrm{d}}(t)=M_{\mathrm{d},0}{(1+\frac{t}{t_{\nu }})}^{1-\eta },$$(4) $$\dot{M}(t)=(\eta -1)\frac{M_{\mathrm{d},0}}{t_{\nu }}{(1+\frac{t}{t_{\nu }})}^{-\eta },$$(5)

where η = (5/2 − γ)/(2 − γ) and the viscous timescale t_ν = R_c²/[3(2 − γ)²ν(R = R_c)] at the characteristic radius R_c. Because ${\dot{M}}_0\propto M_{\mathrm{d},0}/t_{\nu ,0}$, a power-law scaling of M_d,0 and ${\dot{M}}_0$ with the stellar mass implies the viscous timescale t_ν,0 to scale as a power-law with the stellar mass as well, which we define as t_ν,0 ∝ M_*^μ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 δ₀ = λ_m,0 − λ_acc,0, in this case μ₀ = δ₀⁵; therefore, the scaling of t_ν,0 with the stellar mass is determined by the relative values of λ_m,0 and λ_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_{\mathrm{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_{\mathrm{d}}/\dot{M}\propto M_{\star }{}^{{\lambda }_{\mathrm{m}}-{\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c}}}$. Therefore, λ_m and λ_acc must eventually reach the same value, determined by the initial conditions as $${\lambda }_{\mathrm{m},\mathrm{e}\mathrm{v}\mathrm{o}}={\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{v}\mathrm{o}}=\frac{3{\lambda }_{\mathrm{m},0}-{\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c},0}}{2}.$$(6)

Depending on the sign of δ₀ = λ_m,0 − λ_acc,0, the initial slopes can either

• ‘steepen’, that is, λ_m,evo > λ_m,0, if δ₀ > 0 (implying also λ_acc,evo > λ_acc,0);

• ‘flatten’, that is, λ_m,evo < λ_m,0, if δ₀ < 0 (implying also λ_acc,evo < λ_acc,0);

• ‘remain constant’, that is, λ_m,evo = λ_m,0, if δ₀ = 0 (implying also λ_acc,evo = λ_acc,0).

Because in the viscous case δ₀ = μ₀, where we note again that μ₀ is the slope of the correlation between the viscous timescale and the stellar mass (t_ν,0 ∝ M_*^μ0), we can also interpret these scenarios from the viscous timescale perspective. If μ₀ > 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 − ω = 0

In the hybrid viscous and MHD winds model, the general analytical solution by Tabone et al. (2022a) gives $$M_{\mathrm{d}}(t)=M_0{(1+\frac{t}{(1+\psi )t_{\mathrm{a}\mathrm{c}\mathrm{c},0}})}^{-(\psi +2\xi +1)/2},$$(7) $$\dot{M}(t)={\dot{M}}_0{(1+\frac{t}{(1+\psi )t_{\mathrm{a}\mathrm{c}\mathrm{c},0}})}^{-(\psi +4\xi +3)/2}\text{, }$$(8)

where ${\dot{M}}_0$ is defined as $${\dot{M}}_0=\frac{\psi +1+2\xi }{\psi +1}\frac{M_0}{2t_{\mathrm{a}\mathrm{c}\mathrm{c},0}}\frac{1}{(1+f_{{\mathrm{M}}_0})};$$(9)

in this notation, ψ = α_DW/α_SS represents the relative strength of MHD winds and viscosity, $$\xi =\frac{1}{4}(\psi +1)[\sqrt{1+\frac{4\psi }{(\lambda -1)(\psi +1{)}^2}}-1]$$

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)^ξ − 1 the dimensionless mass ejection-to-accretion ratio (with R_in initial disc radius). If we neglect the MHD-driven mass loss (ψ ≪ 1 and ξ ≪ 1, which correspond to f_M,0 ≪ 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_acc is a generalisation of t_ν 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_acc ∝ M_*^μ. However, as mentioned above ${\dot{M}}_0$ depends on the stellar mass not only through M₀ and t_acc,0 as in the viscous case, but also through f_M,0 + 1. As f_M,0 + 1 ∝ R_c,0^ξ, and R_c,0 ∝ M_*^ζ0, the additional dependence will have a slope of ζ₀ξ. Therefore, in the MHD winds scenario we can link δ₀ with μ₀ as δ₀ = μ₀ + ζ₀ξ. 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 ξ 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_{\mathrm{d}}(t\gg t_{\mathrm{a}\mathrm{c}\mathrm{c}})\sim M_0{\left(\frac{t}{t_{\mathrm{a}\mathrm{c}\mathrm{c},0}}\right)}^{-(\psi +2\xi +1)/2},$$(10) $$\dot{M}(t\gg t_{\mathrm{a}\mathrm{c}\mathrm{c}})\sim {\dot{M}}_0{\left(\frac{t}{t_{\mathrm{a}\mathrm{c}\mathrm{c},0}}\right)}^{-(\psi +4\xi +3)/2};$$(11)

this brings the evolved slopes λ_m,evo and λ_acc,evo to $$\begin{array}{rl}& {\lambda }_{\mathrm{m},\mathrm{e}\mathrm{v}\mathrm{o}}={\lambda }_{\mathrm{m},0}+\frac{1}{2}{\mu }_0(\psi +2\xi +1),\\ & {\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c},\mathrm{e}\mathrm{v}\mathrm{o}}={\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c},0}+\frac{1}{2}{\mu }_0(\psi +4\xi +3),\end{array}$$(12)

which reduces to Equation (6) in the viscous case (ψ ≪ 1, ξ ≪ 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 μ₀ as in the viscous case. However, there is a difference from the viscous case: as in the hybrid scenario μ₀ and δ₀ do not coincide anymore, a constraint on the value of μ₀ is translated into a constraint on δ₀ − ζ₀ξ. In particular, the slopes will increase if δ₀ > ζ₀ξ (corresponding to μ₀ > 0), whereas if δ₀ < ζ₀ξ (corresponding to μ₀ < 0), the slopes will decrease; and finally, if δ₀ = ζ₀ξ (corresponding to μ₀ = 0) the slopes will remain constant in time. Another difference from the viscous scenario is that λ_m,evo and λ_acc,evo are not expected to reach the same value anymore. The limit difference is given by δ_evo = λ_m,evo − λ_acc,evo: substituting the values from Equation (12) one finds δ_evo = δ₀ − μ₀(ξ + 1), which can be further reduced to ${\delta }_{\mathrm{e}\mathrm{v}\mathrm{o}}=\xi (2\beta +\frac{1}{2})$ by using the definition of μ₀ as the slope of the t_acc,0 − M_* correlation. In this expression, β 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 β = −1/2 (a reasonable approximation of the value derived by radiative transfer simulation, see e.g. Sinclair et al. 2020), we obtain δ_evo = −ξ/2. As ξ is positive by definition, in the hybrid scenario δ_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 ξ 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_{\star }$ correlations to reach the same limit value (unlike the viscous case). The predicted difference in the evolved slopes is given by δ_evo = −ξ/2 (of the order of 0.1). The current observational uncertainties on δ_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 − ω = 0

The solution to the pure wind model (i.e. time-dependent α_DW through α_DW(t) ∝ ∑_c(t)^−ω) by Tabone et al. (2022a) gives $$M_{\mathrm{d}}(t)=M_0{(1-\frac{\omega }{2t_{\mathrm{a}\mathrm{c}\mathrm{c},0}}t)}^{1/\omega },$$(13) $$\dot{M}(t)=\frac{M_0}{2t_{\mathrm{a}\mathrm{c}\mathrm{c},0}(1+f_{\mathrm{M},0})}{(1-\frac{\omega }{2t_{\mathrm{a}\mathrm{c}\mathrm{c},0}}t)}^{-1+1/\omega };$$(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 ≫ t_acc) and $\dot{M}(t\gg t_{\mathrm{a}\mathrm{c}\mathrm{c}})$. 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 μ₀, implying a positive correlation between the stellar mass and the accretion timescale; a negative correlation between M_* and t_acc,0 (i.e. a negative μ₀) 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 that might impact the evolution of the correlation is the dispersal timescale, t_disp; in the wind-driven case, t_disp ∝ t_acc, therefore t_disp does not introduce any further dependence on the stellar mass.

Fig. 1 Time evolution of the slopes of the Md − M* and $\dot{M}-M_{\star }$ correlations, λm and λacc (blue and orange solid line, respectively) in the hybrid scenario with αSS = αDW = 10−3 (ψ = 1), λ = 3, β = 0.5, resulting in ξ = 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 μ0, slope of the tacc,0 − M* correlation, which is directly linked to the difference between λm,0 and λacc,0 (see text for details): we expect the slopes to (i) decrease if μ0 < 0 (left panel), (ii) remain constant if μ0 = 0 (central panel), and (iii) increase if μ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 −ξ/2. This difference is always negative, meaning that the evolved $\dot{M}-M_{\star }$ correlation is always steeper than that of Md–M* (explaining the lines crossing in the right panel).

Fig. 2 Time evolution of Md–M* (left) and $\dot{M}-M_{\star }$ (right) in the pure wind model (ω ≠ 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 < tacc,0 > with these parameters) due to a downward bending corresponding to lower stellar masses. In this simulation, we have used N = 100 discs, αDW = 10−3, λ = 3, ω = 0.25, λm,0 = 2.1, λacc,0 = 1.5. The set of MHD parameters is based on Tabone et al. (2022a).

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_{\mathrm{d}}}(0)=0.65$ dex, ${\sigma }_{R_{\mathrm{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 (log10(Md/M⊙) = 3.1 log10(M*/M⊙) − 3.3, ${\log }_{10}\left(\dot{M}/M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}\right)=2.5$ log10(M*/M⊙) − 10.4).

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_* 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_* 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 σ_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 α_SS = α_DW = 10⁻³ hence ψ = 1, λ = 3, ω = 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 α_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_lo(M_*/M_⊙) < 0.5 in Figure 2) is lost to the stochastic displacement of the discs in the M_d–M_* and $\dot{M}-M_{\star }$ planes. Actually, the resulting distribution of discs both in the M_d–M_* 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 μ₀, as discussed in Section 3.3; in the simulation shown in Figure 3 we have set μ₀ > 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 λ_m and λ_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 photoevaporation 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_* 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_ν in the viscous case). Our numerical implementation of internal photoevaporation follows Owen et al. (2010), and we assume a mass-loss rate of 4 × 10⁻¹⁰ M_⊙ yr⁻¹ − 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_* dependence is introduced, and the t_ν − M_* 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_* 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⁻¹² M_⊙ yr⁻¹. We have also included a threshold in accretion rates in post-processing to 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.

Fig. 4 Time evolution of the slope of the Md–M* correlation for 15 statistical realisations in the purely viscous (left, αSS = 10−3), photoevaporative (centre, ${\dot{M}}_w=4\times {10}^{-10}{\mathrm{M}}_{\odot }/\mathrm{y}\mathrm{r},{\alpha }_{\mathrm{S}\mathrm{S}}={10}^{-3}$), and wind-driven (right, ω = 0.25, λ = 3, α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.

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_* 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 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σ deviation from the median value of ~0.5–0.6. Larger sample sizes significantly reduce the scatter, leading to σ ~ 0.4 with a double sample and σ ~ 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_* 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.

Fig. 5 Comparison of the time evolution of λm (top) and λ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).

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_* and $\dot{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 μ₀ or δ₀, 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 phase, where the protostellar collapse is over and the disc is already formed. From the observational point of view, this means that we expect δ₀ and μ₀ 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 μ₀ and δ₀.

5.2.1 Constraining δ₀ from μ₀

Ansdell et al. (2017) claimed λ_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 μ₀ > 0 (with t_acc,0 ∝ M_*^μ0). In this section we discuss the implications of the increasing slope scenario on the initial conditions, λ_m,0 and λ_acc,0.

The top panel of Table 2 shows the relation between μ₀ and δ₀ in the different evolutionary models. As μ₀ = δ₀ − ζ₀ξ, in the viscous case (corresponding to ξ = 0) we have μ₀ = δ₀ as mentioned in Section 3. This means that, to recover the suggested increasing slopes scenario, δ₀ necessarily needs to be positive – regardless of the value of any other parameters: this translates to the initial M_d–M_* correlation being steeper than $\dot{M}-M_{\star }$. In the hybrid and wind-driven case, instead, the implication is less straightforward as a positive μ₀ can lead to opposite signs of δ₀: 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 δ₀ ∈ (ζ₀ξ, + ∞), the sign of ζ₀ determines whether negative values of δ₀ are possible; however, as ξ is a small number (0.1 in this work), only a limited area of the parameters space would lead to a negative δ₀. Summarising, if we assume increasing slopes (μ₀ > 0) we can constrain the sign of δ₀ regardless of the evolutionary model assumed: in both cases δ₀ needs to be positive, which leads to an initially steeper M_d–M_* than $\dot{M}-M_{\star }$ correlation.

5.2.2 Constraining μ₀ from δ_evo

Instead of assuming increasing slopes, we can start from the currently measured values of δ_evo and estimate δ₀ in the different evolutionary models. In the viscous case, because δ_evo = 0, we focus on the value of the single slopes instead: as λ_m,evo = λ,_acc,evo = δ₀/2 + λ_m,0, the measured final value of the slopes does not help in constraining δ₀ as it also depends on λ_m,0. In the hybrid case, instead, we have = ξ(ζ₀ − μ₀), meaning that if we can determine the sign of δ_evo we can constrain that of μ₀ as well. While in principle the sign of ζ₀ influences that of μ₀, as we mentioned above ζ₀ is likely a small number: therefore, effectively, δ_evo and μ₀ 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 δ_evo from the most recent and homogeneous measurements available of λ_m and λ_acc (Testi et al. 2022). The resulting values of δ_obs (which we label ‘observed’ as opposed to the theoretical expectation, ‘evolved’), summarised in Table 3, are 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 λ_obs would be compatible with both a positive and a negative value. The difficulty in assessing the sign of λ_obs from the current measurements of the slopes make constraining μ₀ from δ_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 δ_obs − which would possibly allow μ₀ 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_* 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 λ_m and λ_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_{\mathrm{d}}(0)\propto M_{\star }^{{\lambda }_{\mathrm{m},0}},\dot{M}(0)\propto M_{\star }{}^{{\lambda }_{\mathrm{a}\mathrm{c}\mathrm{c},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 Sklodowska- Curie grant agreement No. 823823 (DUSTBUSTERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (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.

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.

Appendix B MHD model with μ < 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 μ₀. Figure 2 shows the a simulation μ₀ > 0, while in Figure B.1 we show the opposite case. As μ₀ is linked to δ₀ through μ₀ = δ₀ − ζ₀ξ, if δ₀ < ζ₀ξ the resulting μ₀ will be negative, leading to a specular bending of the correlation. Given that ζ₀ξ is a small number (~0.1 in our simulation), this generally requires a negative δ₀. The simulation in Figure B.1 has λ_m,0 = 1.3 and λ_acc,0 = 1.7, resulting in δ₀ = −0.4.

Fig. B.1 Same as Figure 2, but with a choice of initial slopes resulting in a negative μ0 (λm,0 = 1–3, λacc,0 = 1–7). The bending of the linear correlation happens towards larger disc masses, in agreement with the prediction.

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, 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 α_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, more extended and less massive discs are more prone to external truncation.

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, αSS = 10−3), viscous including internal photoevaporation (top right, ${\alpha }_{\mathrm{S}\mathrm{S}}={10}^{-3},{\dot{M}}_{\mathrm{w}}=4\times {10}^{-10}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$). wind-driven (bottom left, αDW = 10−3, λ = 3, ω = 0.25) and viscous including external photoevaporation (bottom right, αSS = 10−3, FUV = 100 G0) respectively.

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_{\mathrm{d}}}{2(2-\gamma )t}[1-{\left(\frac{M_{\mathrm{d}}}{M_0}\right)}^{2(2-\gamma )}],$$(C.1)

while in the MHD wind-driven scenario (Tabone et al. 2022a) $$\dot{M}=\frac{1}{\omega (1+f_{\mathrm{M},0})t}{M}_{\mathrm{d}}[{\left(\frac{M_{\mathrm{d}}}{M_0}\right)}^{-\omega }-1].$$(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).

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.

References

All Tables

All Figures

Fig. 1 Time evolution of the slopes of the Md − M* and $\dot{M}-M_{\star }$ correlations, λm and λacc (blue and orange solid line, respectively) in the hybrid scenario with αSS = αDW = 10−3 (ψ = 1), λ = 3, β = 0.5, resulting in ξ = 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 μ0, slope of the tacc,0 − M* correlation, which is directly linked to the difference between λm,0 and λacc,0 (see text for details): we expect the slopes to (i) decrease if μ0 < 0 (left panel), (ii) remain constant if μ0 = 0 (central panel), and (iii) increase if μ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 −ξ/2. This difference is always negative, meaning that the evolved $\dot{M}-M_{\star }$ correlation is always steeper than that of Md–M* (explaining the lines crossing in the right panel).

In the text

Fig. 2 Time evolution of Md–M* (left) and $\dot{M}-M_{\star }$ (right) in the pure wind model (ω ≠ 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 < tacc,0 > with these parameters) due to a downward bending corresponding to lower stellar masses. In this simulation, we have used N = 100 discs, αDW = 10−3, λ = 3, ω = 0.25, λm,0 = 2.1, λacc,0 = 1.5. The set of MHD parameters is based on Tabone et al. (2022a).

In the text

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_{\mathrm{d}}}(0)=0.65$ dex, ${\sigma }_{R_{\mathrm{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 (log10(Md/M⊙) = 3.1 log10(M*/M⊙) − 3.3, ${\log }_{10}\left(\dot{M}/M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}\right)=2.5$ log10(M*/M⊙) − 10.4).

In the text

Fig. 4 Time evolution of the slope of the Md–M* correlation for 15 statistical realisations in the purely viscous (left, αSS = 10−3), photoevaporative (centre, ${\dot{M}}_w=4\times {10}^{-10}{\mathrm{M}}_{\odot }/\mathrm{y}\mathrm{r},{\alpha }_{\mathrm{S}\mathrm{S}}={10}^{-3}$), and wind-driven (right, ω = 0.25, λ = 3, α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.

In the text

Fig. 5 Comparison of the time evolution of λm (top) and λ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).

In the text

	Fig. B.1 Same as Figure 2, but with a choice of initial slopes resulting in a negative μ0 (λm,0 = 1–3, λacc,0 = 1–7). The bending of the linear correlation happens towards larger disc masses, in agreement with the prediction.
In the text

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, αSS = 10−3), viscous including internal photoevaporation (top right, ${\alpha }_{\mathrm{S}\mathrm{S}}={10}^{-3},{\dot{M}}_{\mathrm{w}}=4\times {10}^{-10}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$). wind-driven (bottom left, αDW = 10−3, λ = 3, ω = 0.25) and viscous including external photoevaporation (bottom right, αSS = 10−3, FUV = 100 G0) respectively.

In the text

	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.
In the text