© The Authors 2023

1 Introduction

In recent years, the number of observed extrasolar planets has continued to rise exponentially. The growing amount of observational data offers the opportunity to put theoretical planet formation models to the test (e.g. Alessi & Pudritz 2018; Emsenhuber et al. 2021a,b). In working to understand the processes involved with planet formation, we recognise the key role played by the evolution of the protoplanetary disc, as it delivers the material from which the planets grow. Currently, most global planet formation models rely on the viscous α-disc model (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), where α is the scaling parameter for the effectiveness of angular momentum transport, whose physical origins are still under debate. Recently, the focus shifted towards magnetic fields as a promising new alternative to inducing angular momentum transport.

Protoplanetary discs are thought to be threaded by a large-scale poloidal magnetic field as a remnant of the parental molecular cloud core magnetic field. The presence of such a magnetic field is known to have strong influence on accretion processes inside the disc. Magnetorotational instability (hereafter, MRI; Balbus & Hawley 1991, 1997) was regarded as a promising mechanism for inducing turbulence and, hence, angular momentum transport in protoplanetary discs, whereas the presence of a magnetic field can also drive a magneto centrifugal wind, removing both mass and angular momentum (Blandford & Payne 1982; Königl & Salmeron 2010). Although low ionization rates are sufficient for the magnetic field to couple with the gas, non-ideal magnetohydrodynamics (hereafter, MHD) effects, such as ohmic dissipation and ambipolar diffusion, can render large swathes of the disc dead (MRI inactive). Gammie (1996) realised that ohmic dissipation would create a region near the midplane where MRI is suppressed (the so-called dead zone) and suggested that accretion is happening via MRI active layers above the midplane. When including ambipolar diffusion, the dead zone is extended to regions with low density, allowing only for a very thin layer to drive accretion (Perez-Becker & Chiang 2011b,a; Bai & Stone 2013). This leaves magneto-centrifugal winds as a promising mechanism to explain the angular momentum removal. Far-ultraviolet (FUV) radiation can lead to the efficient ionization of the upper layer of the disc, while strongly coupling gas to the magnetic field. Gas is being loaded onto the field lines and centrifugally accelerated, removing both mass and angular momentum.

Observing magnetic fields and related outflows in protoplanetary discs is a challenging task that has only been accomplished very recently. There is evidence for hourglass-shaped magnetic fields in NGC 2024 and NGC 1333 IRAS4A (Crutcher 2012). Furthermore, Harrison et al. (2021) recently measured upper limits on the magnetic field strength of AS 209, while Whelan et al. (2021) found evidence of an MHD disc wind emerging from the two accreting T Tauri stars: RU Lupi and AS 205 N.

Magnetically driven disc winds may not only have strong influence on angular momentum transport inside the protoplanetary disc, but they are also able to inhibit or even revert the type I migration of protoplanets (Ogihara et al. 2015a,b; Suzuki et al. 2016). The last of the mentioned works introduced a prescription for the removal of angular momentum and mass by MHD disc winds in a 1D disc model. However, the authors could not constrain all the parameters of their model because of the lack of observational detections of the process as well as the uncertainties in the underlying hydrodynamical simulations. Here, we aim to better constrain their prescription by looking at population-level results. We investigate the evolution of MRI-inactive (i.e. low-α) protoplanetary discs under the combined effects of magnetically driven disc winds as well as internal and external photoevaporation. We use different model settings that were suggested in Suzuki et al. (2016) and we determine which ones match protoplanetary disc observations in terms of lifetimes, stellar accretion rates, and masses. This will allow us to investigate the imprint of magnetically driven disc winds on planet populations in future works.

This paper is structured as follows: in Sect. 2 we give a full description of the evolution model and initial conditions. The results of disc population syntheses are shown in comparison with observational data in Sect. 3, followed by a discussion in Sect. 4. Our conclusions are listed in Sect. 5.

2 Methodology

We conducted 1D simulations of protoplanetary disc evolution, including the effects of magnetically driven disc winds and both internal and external photoevaporation, across a wide range of initial conditions (e.g. initial disc mass, inner disc truncation radius, ambient FUV field strength, etc.). The distributions of initial conditions are chosen to reflect the conditions found in young star forming regions. The disc evolution model is presented in Sect. 2.1 and initial conditions are given in Sect. 2.2.

2.1 Disc evolution model

The model at hand is an enhanced version of the gas disc model used in the Generation III Bern global model of planetary formation and evolution (Emsenhuber et al. 2021a). The new model includes magnetically driven disc winds according to the model of Suzuki et al. (2016) to account for the removal of angular momentum and mass originating from magnetic fields threading the disc. Furthermore, here we use a more physically motivated model for the external photoevaporation based on the FRIED grid from Haworth et al. (2018). In the following, we describe the model, with an emphasis on the newly added parts.

The protoplanetary disc is treated as a 1D radial axissymmetric structure in cylindrical coordinates (r, ϕ, z). The time evolution of the surface density Σ(r) = ∫ ρ(r, z)dz is performed by numerically solving the evolution equation: (1)

Here, we use standard notation where r is the distance from the host star, t is the time, and Ω corresponds to the angular velocity which is assumed to be Keplerian , with G being the gravitational constant, is the sound speed, and the molecular weight is µ = 2.24, hydrogen atom mass, m_H, with k_в being the Boltzmann constant, and ρ is the gas density. The subscript ()_mind denotes values at the mid-plane. and are parametrisations for the redistribution and removal of angular momentum.

The equation incorporates classical viscous diffusion, advection through angular momentum removal and sink terms for removal of mass. It is based on the equation derived in Appendix A of Suzuki et al. (2016). Here, corresponds to the effective viscosity used in the classical α -disc model (Shakura & Sunyaev 1973); it is mathematically connected to a turbulent viscosity via ¹. However, in the picture of magnetic fields, we ought to think of it as a parameter for angular momentum transport through MRI. The advection is parameterised by , which corresponds to the non-dimensional stress acting as a torque on the disc, removing angular momentum and driving accretion inside the disc. Furthermore, we have various sink terms for mass removal via magnetically driven disc winds (), internal photoevaporation () and external photoevaporation () This is similar to Kunitomo et al. (2020), although we concurrently include internal and external photoevaporation. A cold MHD wind and photoevaporation are both of different physical nature and, strictly speaking, we would expect a single wind of an intermediate nature (Bai et al. 2016). Here, we represent the view where both winds play important roles at different stages of the disc evolution (e.g. Pascucci et al. 2022).

The equation is solved on a logarithmically-spaced grid with 3400 grid cells between the inner disc edge, given by an initial condition (Sect. 2.2.5), and 1000 AU. We use an advection-diffusion algorithm derived in Appendix A of Birnstiel et al. (2010). We set a minimum value surface density of Σ_min = 10⁻⁴ g cm⁻² throughout the grid and impose Dirichlet boundary conditions Σ_min at the disc edges.

2.1.1 Thermal structure

The thermal structure of the disc is evaluated at each time evolution step of the disc, following the approach in Nakamoto & Nakagawa (1994): (2)

Here, σ_SB is the Stefan-Boltzmann constant. T_mid corresponds to the temperature at the midplane. F_rad is the transferred energy from viscous dissipation and liberated gravitational energy; T_ext accounts for external sources of heating. The Planck mean optical depth is given by τ_P = 2.4τ_R, while the Rosseland mean optical depth is expressed as τ_R = κ(ρ_mid, T_mid)∑, where the opacity κ is a function of midplane density and temperature. We use the maximum of the values given by the expressions from Bell & Lin (1994) and Freedman et al. (2014); the latter for a solar composition.

External heating is considered by disc surface irradiation-and direct irradiation through the midplane by the host star and heating from the surrounding molecular cloud; (3)

For the irradiation of the disc surface by the host star, we follow the calculations for flared discs (Ruden & Pollack 1991; Hueso & Guillot 2005): (4)

where we adopt ∂ ln(H)/∂ ln(r) = 9/7 from Chiang & Goldreich (1997). Since the inner edge of the disc is directly exposed to irradiation from the host star, the disc midplane gets additionally heated, as follows: (5)

while taking into account the optical depth through the midplane, given as . To account for heating of the disc surface by the surrounding molecular cloud, we add a constant value of T_cloud = 10 K. This is a strong simplification, since we neglected variations in the heating from the stellar cluster environment (e.g. Ndugu et al. 2018).

2.1.2 Magnetic disc winds

Magnetically driven disc winds are incorporated using the model from Suzuki et al. (2016). The model considers both removal of angular momentum (magnetic braking) and mass.

The magnetic braking results in the advective term in the evolution equation (Eq. (1)) and is parametrised by . Strength and time evolution of the magnetic field is still not well understood. Bai (2013) reported a positive dependence on the strength of the magnetic field , with torques between from local MHD simulations. We adopt the two parametrisations for presented in Suzuki et al. (2016).

Firstly, in the case where the magnetic field diffuses outwards with decreasing Σ, the torque will stay approximately constant due to the dependency of ρ on Σ and we therefore chose = const. Secondly, if the magnetic field stays constant, only ρ is dependent on Σ. This can be parametrised by (6)

where is the initial value of and Σ_init is the initial surface density at the given location.

The sink term for mass removal by magnetic disc winds is given by (7)

with C_w being a non-dimensional mass loss rate. This non-dimensional mass loss rate is constrained by the disc wind energetics C_w,e (Eqs. (9) and (11)) and confined by an upper limit of C_w,0 = 10⁻⁵ to avoid very high mass fluxes, (8)

From the conservation of total MHD energy, we can derive an energy constrain on the disc winds; speaking in terms of variables this means that the mass flux C_w can be expressed in terms of , , and other local quantities. We adopted two approaches, namely: ‘strong disc wind’ and ‘weak disc wind’, as discussed in Suzuki et al. (2016).

Strong disc wind

The strong disc wind scenario assumes that all liberated gravitational energy contributes to launching a wind, while viscous heating is transferred to radiation (i.e. disc heating). This leads to: (9) (10)

Weak disc wind

In the weak disc wind scenario, the liberated energy contributes to both launching a wind and heating of the disc. The percentage that contributes to disc heating is controlled by the parameter ϵ_rad ∊ [0, 1]. Here, we have (11) (12)

It can be seen that in the limit of no disc winds (ϵ_rad = 1 and ), we get back to the equation for the viscous dissipation rate of (Nakamoto & Nakagawa 1994; Emsenhuber et al. 2021a). For the scenario of a weak wind, we followed Suzuki et al. (2016) and adopted ϵ_rad = 0.9, meaning that only 10% of the liberated energy goes into launching the wind.

2.1.3 Photoevaporation

Protoplanetary discs are exposed to high-energy radiation from the central star and the surrounding cluster. The MHD wind model at hand describes a cold, magnetocentrifugal wind which is launched by the liberated accretion energy. The presence of high energy radiation would lead to a single wind of possibly intermediate nature, namely, magneto-thermal wind (e.g. Bai et al. 2016; Wang et al. 2019; Rodenkirch et al. 2020). However, the physical nature of such a wind has yet to be fully characterised. Therefore, we investigate the scenario, where internal photoevaporation is suppressed by the emerging magnetic wind in an early phase and dominates only the final stage of the disc evolution (Lesur et al. 2022; Pascucci et al. 2022). While such a separation in time would not work for external photoevaporation, it seems likely that the presence of an external UV radiation field would lead to a considerably enhanced mass loss rate at the outer disc compared to the cold wind scenario; this is a hypothesis that will have to be investigated with magneto-thermal models in the future.

Given these assumptions, we consider both internal photoevaporation by the host star through extreme-ultraviolet radiation (EUV; 13.6eV < hv) and external photoevaporation by far-ultraviolet radiation (FUV; 6eV < hv < 13.6eV) from the surrounding cluster. The model for internal photoevaporation is the same as used in Emsenhuber et al. (2021a), except for the EUV flux scaling and shielding by disc winds. For external photoevaporation we now use a new model based on precalculated values from Haworth et al. (2018). Recent work on internal photoevaporation showed that X-ray photoevaporation is expected to excel EUV photoevaporation by order of magnitudes (e.g. Jennings et al. 2018). Emsenhuber et al. (2023) compared the evolution of viscous disc populations including both internal X-ray and external FUV photoevaporation. They found that the high mass-loss rates of the two leads towards too short lived discs, which would be inconsistent with observations. This coupled with the fact that X-ray can penetrate much higher column densities (making the desired separation in time by radiation shielding ineffective) only motivate the inclusion of a simple EUV photoevaporation prescription.

Internal photoevaporation

We used a simple parametrisation for the internal photoevaporation following the Clarke et al. (2001) model, that is, one based on the weak stellar wind model from Hollenbach et al. (1994). EUV radiation from the host star heats the surface of the gas disc and creates a layer of ionized hydrogen with a temperature of T_II ≈ 10⁴ K and a mean molecular weight of µ_II = 0.68.² In the model from Hollenbach et al. (1994), this heated layer is expected to launch a thermal wind at sound speed beyond some distance from_t the host star, where the sound speed of the ionized gas, , is greater than the escape velocity, , and is not gravitationally bound any more (i.e. the gravitational radius; ). Similarly to Alexander & Pascucci (2012), we reduced the gravitational radius to a critical radius, r_crit, by a factor of β_II = 0.14 (see also Liffman 2003).

Following Clarke et al. (2001), we defined the base density of the wind at the critical radius as , where r_crit,14 = β_IIr_g,II/10¹⁴cm is the scaling radius and Φ₄₁ corresponds to the ionizing photon luminosity in units of 10⁴¹ s⁻¹. The ionizing photon luminosity is scaled with stellar mass according to Φ_EUV = 10^40.7(M_⋆/M_⊙)^1.5 (see Table 2, Gorti & Hollenbach 2009). This scaling law holds true for stellar masses ≲ 3M_⊙, which is the case in our populations. The base density outside the critical radius is assumed to follow the power law n₀(r) = n₀(r_crit,14) · (r/r_crit)^−5/2. Putting everything together leaves us with the sink term for the internal photoevaporation: (13)

The total mass loss rate scales with stellar mass as , which is below recent calculations from Komaki et al. (2021), who derived a stellar mass dependence of considering not only EUV but also FUV and X-ray irradiation by the host star.

Photoevaporative winds can only be launched when high energy photons manage to heat the disc beyond the critical radius r_crit. This becomes inherently difficult with a strong magnetic wind emerging from the inner disc < 1AU and it is therefore expected that internal photoevaporation is not effective in the early stage of disc evolution (Pascucci et al. 2022). In order to take this shielding effect into account, we calculate the column density of the emerging magnetic disc wind along the disc surface. This gives (14)

where υ_wind = 70 km s⁻¹ is the typical velocity of the emerging wind (Pascucci et al. 2022) and µ_wind = 2.23 is the mean molecular weight of the wind in terms of the hydrogen atom mass, m_H. The factor of 2 takes into account that the wind emerges from both sides of the disc. We assume that the EUV radiation can penetrate only a column density of < 10¹⁹ cm⁻² (<10²⁰ cm⁻² Ercolano et al. 2009; ≲10^19–20 cm⁻² Pascucci et al. 2022). For regions where n_wind(r) > 10¹⁹ cm⁻², the internal photoevaporation is suppressed (i.e., Σ_PEW,int(r) = 0). This represents a simple, but physically motivated approach to describe the interplay between MHD and photoevaporative winds.

External photoevaporation

Star formation usually occurs in groups of hundreds to thousands of stars. Nearby young massive stars emit ultraviolet radiation that can heat up the gas disc, which results in a thermally driven wind. We use the FUV Radiation Induced Evaporation of Discs (FRIED) grid from Haworth et al. (2018) to obtain mass loss rates from far-ultraviolet irradiation, which is considered to be the dominant driver of external photoevaporation (Adams et al. 2004). The grid consists of precalculated mass loss rates for various stellar masses (0.05–1.9M_⊙), FUV field strengths (10–10⁴ G₀)³, disc masses (3.2 × 10⁻⁵ − 0.2 M_⋆ of the stellar mass), and sizes (1 – 400 AU). Linear interpolation was used to retrieve mass loss rates for given stellar masses and FUV field strengths, as well as disc sizes and masses in between the grid points⁴. For the FUV field strength, we carried out a linear interpolation in loglog. For the stellar mass, we carried out a linear interpolation in lin-log (linear in stellar mass).

While stellar mass and FUV field strength are given relatively straightforward, this is not the case for the disc radius and the corresponding mass of the disc. As discussed in Sellek et al. (2020a) the effective outer radius of the disc is located at the transition between the optically thick and optically thin regime. Analogously to Sellek et al. (2020a), we evaluate the mass loss rate for all disc grid cells and set the disc radius to where the obtained mass loss rate Ṁ_PEW,ext is maximal.

The FRIED grid includes a floor value of 10⁻¹⁰ M_⊙ yr⁻¹ for the evaporation rate. This is a significant rate and prevents the proper investigation of weak FUV field environments. Nonetheless, in striving to study them, we subtracted the floor value of 10⁻¹⁰ M_⊙ yr⁻¹ from the mass-loss rate returned by the interpolation, while ensuring a negligible minimum mass-loss rate of 10⁻¹⁵ M_⊙ yr⁻¹. Further, we extended the range of possible FUV field strengths down to 1 G₀ by linearly interpolating the value returned from the grid at 10 G₀ and a new floor value of 10⁻¹⁵ M_⊙ yr⁻¹ at 1 G₀.

The 2D calculations from Haworth & Clarke (2019) suggest that the mass loss rates from the FRIED grid have to be regarded as lower limits. Their simulations furthermore illustrated that the mass loss rate is set entirely by the outer half of the disc and originates mostly from within the outer 10% of the disc outer edge, R_edge, which we define as where the outer surface density reaches Σ_min, our minimum value set throughout the grid. We therefore assume that the mass is removed uniformly from the outer 10% of the disc (β_ext = 0.9): (15)

This new model for external photoevaporation replaces our old method described in Emsenhuber et al. (2021a,b), where the mass loss rate Ṁ_PEW,ext was chosen such that the synthetic disc lifetimes fitted the observations. With the loss of this tuning parameter, disc lifetimes are entirely given by the initial conditions, inferred from observations.

2.2 Initial conditions

We used the approach of evolving large populations of protoplanetary discs, exploring a wide parameter space of initial conditions (e.g. Lodato et al. 2017; Mulders et al. 2017; Schib et al. 2021; Emsenhuber et al. 2021b; Tabone et al. 2022b). We conducted disc population syntheses for different combinations of disc wind scenarios and ambient FUV field strengths, summarised in Table 1. These scenarios cover the four possible combination of strong or weak disc winds in combination with either a constant torque (decreasing magnetic field) or increasing torque (constant magnetic field), as described in Sect. 2.1.2. Following the discussion above, all discs are assumed to be MRI inactive with (Suzuki et al. 2016). Besides these specified initial parameters, we also have a set of Monte Carlo variables that are distributed in order to reproduce the various conditions found in star-forming clusters, namely: (1) the dust-to-gas ratio (f_D/G, Sect. 2.2.2), (2) the initial disc mass (M_disc,init Sect. 2.2.3), (3) the disc’s characteristic radius (R_char, Sect. 2.2.4), (4) the disc’s inner edge (R_in, Sect. 2.2.5), (5) the stellar mass (M_⋆, Sect. 2.2.6), and (6) the FUV field strength (ℱ_FUV, Sect. 2.2.7).

The resulting distributions of the initial parameters are shown in Fig. 1. With the exception of the reseating of the disc mass with stellar mass, we assume the Monte Carlo variables to be independent of each other. However, it may well be that there is some relation between the stellar mass and the FUV field strength, for example.

Fig. 1 Adopted distributions of the Monte Carlo variables for our disc populations: dust-to-gas ratio (panel a), initial dust masses (panel b), resulting gas mass (panel c), inner disc radius (panel d), FUV field strength for weak-FUV and strong-FUV environments (panel e), and stellar mass (panel f).

2.2.1 Initial surface density

We assume an initial gas surface density similar to Veras & Armitage (2004), with a power-law index β and an exponential cut-off with characteristic radius, R_char·The inner disc is assumed to be truncated by the magnetic field of the host star through magnetospheric accretion; therefore, R_in is taken to be the corrotational radius. The initial surface density profile is set with (16)

Here, ∑_0,5.2AU is the initial surface density at 5.2AU, and β = 0.9 was chosen according to Andrews et al. (2010). The value of ∑_0,5.2AU is fully determined by the integral of Eq. (16), ignoring the decrease at the inner edge, which is negligible (Eq. (14) in Emsenhuber et al. 2021a).

2.2.2 Dust-to-gas ratio

Disc masses are usually obtained from dust continuum emission measurements. The dust-to-gas ratio, f_D/G, is a necessary quantity to convert these initial dust masses to initial gas masses. Following Emsenhuber et al. (2021b), we assume that stellar and disc metallicities [Fe/H] are identical and, thus, we used the relation (Murray et al. 2001): (17)

with the Sun’s dust-to-gas ratio being f_D/G,⊙ = 0.0149 (Lodders 2003). For the distribution of the metallicity [Fe/H], we used data from the Coralie RV search sample (Santos et al. 2005), where the metallicity is normal distributed with µ = −0.02 and σ = 0.22. We note that we limited the metallicity to −0.6 < [Fe/H] < 0.5 to avoid unrealistically high or low values in comparison to the solar neighbourhood.

Fig. 2 Initial relation of disc dust mass and characteristic radius for a synthetic population of 1000 discs, following Fig. 12 of Tobin et al. (2020).

2.2.3 Initial disc mass

We used a log-normal fit to dust disc masses of Class 1 protoplanetary discs in the Perseus star-forming region from Tychoniec et al. (2018) with log₁₀(µ/M_⊕) = 2.03 and σ =0.35 dex. As the host star masses of the sample are not known, we assume that it is representative for stellar masses of M_⋆ ~ 0.3 M_⊙ (see discussion in Tobin et al. 2016). We further assumed a simple linear correlation between the disc and stellar masses (M_disc ∝M_⋆, Raymond et al. 2007; Andrews et al. 2013). Somigliana et al. (2022) inferred a slightly steeper initial correlation from comparing observations and numerical models of viscous disc evolution . Furthermore, the relation has found to be steepening with age (Pascucci et al. 2016; Somigliana et al. 2022). We use the previously defined dust-to-gas ratio to convert the dust masses to gas masses. Again, we limited the distribution to disc masses of 4 × 10⁻⁴ M_⋆ < M_disc < 0.16M_⋆, where the upper limit ensures self-gravitational stability.

2.2.4 Characteristic disc radius

The characteristic disc radius, R_char, is inferred from the disc dust mass using the relation R_char = 70·[M_dust/100 M_⊕]^0.25, with an applied spread of 0.1 dex to achieve a similar distribution as in Fig. 12 of Tobin et al. (2020). The resulting correlation between dust masses and characteristic radii is shown in Fig. 2.

2.2.5 Inner disc edge

Following Emsenhuber et al. (2021b) once again, we assume that the inner disc is truncated by the stellar magnetic field (i.e. magnetospheric accretion). The inner radius can thus be inferred from rotation rates of young stellar objects by calculating the corotational radius r_co = (G·M_⋆/Ω)^1/3, with Ω = 2π/P being the angular velocity and P the rotation period of the central star. We used a log-normal distribution fit to observed rotation rates of young stars in the NGC 2264 open cluster from Venuti et al. (2017), log₁₀(µ/d) = 0.676 and σ = 0.306 dex. To avoid having inner radii smaller than the initial stellar radius, we set a lower bound to R_in = 1.65 × 10⁻² AU (as in Emsenhuber et al. 2021b).

2.2.6 Stellar mass

Both internal and external photoevaporation processes are sensitive to the mass of the host star. Observed protoplanetary discs usually surround stars with stellar masses below 1 M_⊙. We therefore employed a commonly used parametrization for the initial mass function by Chabrier (2003), with a log-normal distribution for M_⋆ ≤ 1 M_⊙ and a power-law decay above M_⋆ > 1 M_⊙. Our stellar mass distribution spans a range from 0.08 M_⊙ to 1.5 M_⊙, with a median mass of 0.21 M_⊙, which is in agreement with observations of stars hosting a protoplanetary disc (e.g. Rigliaco et al. 2011; Tobin et al. 2020).

2.2.7 FUV field strength

The local FUV flux, ℱ_FUV, which drives external photoevaporation, is determined by the location inside the stellar cluster. We used two approaches to characterise either strong or weak FUV field environments. For the strong FUV field approach, we used an ensemble distribution for FUV fluxes experienced by stars in a cluster, calculated by Adams et al. (2006, their Fig. 9). The distribution shows approximately log-normal behaviour with a tail towards lower values of G₀. We approached the weak FUV field by assuming a simple log-normal distribution with log₁₀ (µ/G₀) = 1 and σ = 0.5 dex. Values beyond 10⁴ G₀ are treated as the boundary value. This step is not expected to have any consequence in this case, since discs in a FUV field of 10⁴ G₀ or more have already evaporated after a few 10⁵ yr.

3 Results

For each of the eight combinations listed in Table 1, we calculated the temporal evolution of a population, comprising 1000 discs each. The resulting observables (e.g. stellar accretion rate, disc mass, disc lifetime, etc.) can thus be compared with observational data.

3.1 Secular evolution: Exemplary cases

Before going into the results of the population syntheses, we offer a brief discussion of some general characteristic features of our model. We selected a 0.1 M_⋆ disc around a 0.25 M_⊙ star, which gives an initial disc mass of 0.025 M_⊙, and we evolved it for different disc wind scenarios. Figure 3 shows the evolution of the disc in a 10 G₀ FUV field, while Fig. 4 shows the evolution of the same disc in a 5000 G₀ environment. For the eight resulting discs, we show the detailed surface density, radial mass flow, and outflow evolution. We compute the radial mass flow as: (18)

where the first term accounts for contributions of the diffusive term and the second term represents accretion driven by magnetic braking (Suzuki et al. 2016). The outflows contain the sum of the three different components (magnetic winds and internal and external photoevaporation).

In the weak FFUV case (Fig. 3), the evolution of the surface density is similar to the model of Kunitomo et al. (2020), with the exception of the outer disc truncation due to external photoevaporation. Over a large part of the disc (except for the outer disc), magnetic winds initially dominate the mass loss. With strong winds, the inner disc is exposed to a strong wind, resulting in a rapid decrease in the surface density in the inner region. This also means that only a fraction of the material reaches the inner disc, thus the stellar accretion rate always remains low (~ 10⁻¹⁰ M_⊙ yr⁻¹). With weak disc winds, the contrast of the radial mass flow in the inner and outer region is much lower, leading to larger stellar accretion rates. Furthermore, the torque strength has a strong influence on the radial mass flow. With the Σ-dependent torque prescription we use here (which has an initial torque strength of ), the radial mass flow in the disc is ≲1 × 10⁻⁹ M_⊙ yr⁻¹ ; whereas for a strong constant torque , it is larger initially. The accretion rate in the Σ-dependent case decreases less rapidly with time. For this set of initial conditions, only the combination of weak disc winds and a strong constant torque is able to produce a stellar accretion rate larger than 10⁻⁹ M_⊙ yr⁻¹. The Σ-dependent torques with larger initial values of could also produce larger stellar accretion rates (see Sect. 3.4 for some insights on higher initial values). In all but one case, an inner cavity opens up towards the end of the evolution, leading to a complete stop of accretion onto the star.

For a field of 5000 G₀ (Fig. 4), the discs are dispersed very rapidly ≲1 Myr outside-in; thus, there is not enough time for an inner cavity to open up. Here, internal photoevaporation does not play a role in disc dispersal, as the discs are quickly dispersed from the outside. As external photoevaporation has little effect in the inner region of the disc, we still find that only the combination of weak disc winds and strong constant torques produces a stellar accretion rate of 10⁻⁹ M_⊙ yr⁻¹ or more.

We note that the radial flow diagrams show a rapid outflow at the outer edge of the disc. This is due to the rapid fall off in surface density and the low-α viscosity imposed. However, we do not witness any viscous spreading of the disc, as external photoevaporation (even in the case of weak ℱ_FUV) compensates for this effect.

Fig. 3 Exemplary disc evolution shown for initial conditions: Mdisc = 0.1 M⋆ with M⋆ = 0.25 M⊙, Rin = 0.04 AU and ℱFUV = 10 G0. We show temporal evolution of the surface density, radial mass flow rate and rate of change in surface density due to outflows (i.e. magnetic disc winds, internal- and external photoevaporation) for all disc wind scenarios. The individual contributions can be distinguished roughly when considering the outflows . External photoevaporation corresponds to the outermost peak while internal photoevaporation acts at ~ 0.3 AU to a few AU. MHD winds remove mass from the whole disc and increase in strength towards the inner disc. Dotted lines show the evolution at 104 yr, dashed lines at 5 × 104 yr. Coloured lines are spaced by 105 yr.

Fig. 4 Exemplary disc evolution, analogous to Fig. 3 but for a strong ambient FUV field strength of 5000 G0.

Fig. 5 Stellar accretion rate Ṁacc vs. gas disc mass Mdisc shown for all disc wind scenarios (see Table 1) with a weak FUV field strength distribution. Evolution tracks of individual systems, coloured by ℱFUV, and snapshots at 2 Myr (green circles and 3 Myr (orange circles are displayed. Tracks of exemplary cases are highlighted with black dashed and dotted lines. We show the detailed evolution of those exemplary discs in Figs. 3 and 4. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant Mdisc/Ṁacc are shown for 0.1, 1 and 10 Myr (thin dotted lines).

3.2 Stellar accretion rate and disc dispersal processes

The stellar accretion rate is an important observable that is correlated with the disc mass and age of the system (Manara et al. 2016; Hartmann et al. 2016; Testi et al. 2022). Determining the stellar accretion rates, disc masses, and stellar ages is very challenging (see Hartmann et al. 2016, for a review). As the determination of the age is affected by large uncertainty and the choice of a time of ‘zero’ in our simulations is not necessarily in agreement, we only look at the stellar accretion rates and disc masses. The stellar accretion rate is usually derived from the accretion luminosity, assuming magnetospheric accretion (e.g. Gullbring et al. 1998; Alcalá et al. 2017), while the mass of the disc is usually derived from dust continuum measurements and with the assumption of f_D/G = 0.01 (e.g. Tychoniec et al. 2018; Manara et al. 2019).

In our simulations, the stellar accretion rate corresponds to the value of the radial mass flow (Eq. (18)) at the inner boundary of the grid (magnetospheric accretion, Sect. 2.2.5. In Figs. 5 and 6, we show the stellar accretion rate and the gas disc mass of the systems as a side-by-side comparison of the different disc wind scenarios and FUV field strength distributions. Observational data for comparison is taken from Rigliaco et al. (2011) and Ansdell et al. (2017) for σ Orionis and Manara et al. (2019) for Lupus and Chamaeleon I. In Tables 2 and 3, we give some characteristic numbers on time scales and mass-loss from different processes for each population.

Evolution tracks in the M_disc-Ṁ_acc plane show essentially two distinct behaviours of the evolutionary path.

Firstly, there are systems where the stellar accretion rate stays above ≳ 10⁻¹² M_⊙ yr⁻¹, while the disc mass declines to masses below ≲ 10⁻⁵ M_⊙, caused by large parts of the mass being removed from the outer disc by external photoevaporation, similar to the case shown in Fig. 4. This corresponds to quasi horizontal tracks. We refer to this process as ‘outside-in dispersion’ (e.g. Scally & Clarke 2001; Koepferl et al. 2013).

Secondly, the other systems show a rapid decrease in stellar accretion rate while still possessing a disc of a mass of >10⁻⁵ M_⊙. This corresponds to tracks which ultimately move nearly vertically downwards. This is caused by an inner cavity opening up, preventing stellar accretion, in a similar way to the case of the discs shown in Fig. 3. These discs are referred to as transitional discs, which show inside-out dispersal (e.g. Alexander et al. 2014).

As discussed in Sect. 3.1, there is a strong dependence of the evolutionary path on the FUV field strength. This is reflected in the results shown in Fig. 6, where discs exposed to strong FUV fields beyond approximately 3000 G₀ sustain strong mass loss from external photoevaporation. These discs are rapidly dispersed before an inner hole can open. In the presence of weaker FUV fields, the external photoevaporation leads to a lower mass-loss rate, giving magnetically driven disc winds and internal photoevaporation time to open up an inner cavity. The inner cavity in transitional discs is created by preventing the inner disc from being resupplied by inward flowing gas from the mass reservoir in the outer disc. The inner cavities in our simulations usually initially open up at less than ≲ 1 AU.

In Tables 2 and 3, we list the fraction of discs opening up an inner cavity. We define discs going through such a transitional phase as ones that are optically thin (1 > τ) at radii ≳ 1 AU, with part of the outer disc being optically thick (1 < τ). We get fractions of 35% to 63% of all discs going through a transitional phase in weak ambient FUV fields, whereas discs in a strong FUV environment are less likely to open up a cavity (11% to 34%. Discs exposed to a strong magnetic disc wind are more likely to go through a transitional phase.

To compare our synthetic populations with observations, we offer snapshots of the populations at fixed ages. For a weak ambient FUV field, we display snapshots at 2 Myr (orange circles) and 3 Myr (green circles). We note that in the presence of strong disc winds, discs are dispersed much faster. Thus, at 2 Myr, there are merely 10% to 32% left, compared to substantial 81% to 88% for weak disc winds. Hence, the number of synthetic data points depends on the magnetic disc wind scenario. Furthermore, we note that the location of the bulk in the M_disc - Ṁ_acc plane does not show strong dependence on time. For the strong FUV field distribution we show a snapshot at 4 Myr (orange circles). At that time, discs exposed to very strong FUV fields have already dispersed. Discs still present at that time are exposed to FUV fields of ≲1000 G₀ (see Fig. 7). We note again the low number of data points for strong disc winds.

Since the magnetic disc winds are strongest within <1 AU from the inner edge (Suzuki et al. 2016; Kunitomo et al. 2020), the evolution of the inner disc is strongly affected by the magnetically driven disc winds. This can be seen in the ensemble evolution of the different disc wind scenarios (see Figs. 3 and 4). The combination of a strong magnetic disc wind in the early phase and the onset of internal photoevaporation at a later stage, when the column density of the magnetic wind drops, can lead to rapid depletion of the inner disc region, causing the stellar accretion rate to decrease early on during the evolution. This behaviour is also reflected by the larger fraction of discs opening up a cavity in the strong disc wind scenarios (~57–63% and ~30–34% for a weak and a strong FUV field, respectively; see Tables 2 and 3).

In Tables 2 and 3, we show the mean fraction of initial mass lost through different processes. We see that external photoevaporation is by far the most important process in dispersing the disc, both for weak (~50%) and strong (80% to 90%) FUV field environments. Values are larger in the cases with Σ-dependent torques. This is the effect of a lesser amount of radial mass flow in this scenario (Figs. 3 and 4), leaving more gas in the outer disc, which is then more susceptible to external photoevaporation. The mass fraction removed by magnetic winds is dependent on the disc wind scenario. Strong disc winds lead to larger losses through this process, which can reach up to ~42% in weak and ~13% in strong FUV fields. Constant torque scenarios also lead to greater losses through disc winds, although the difference is lower than that among strong and weak disc winds. Internal photoevaporation is less effective except for the case of weak disc winds and Σ-dependent torques. The combination of low disc winds (from the weak disc winds scenario) and lower initial radial material redistribution (from Σ-dependent torques) lead to a lower column density of the magnetic wind so that internal photoevaporation is active earlier on.

Stellar accretion plays only a subordinate role for disc dispersal in our simulations. With strong disc winds, the losses are so large before the gas can reach the inner edge that the amount of the disc being accreted onto the star does not exceed on average 1 %. Weak disc winds allow for more mass to be accreted by the star (up to 17%).

Observational data from star forming regions Lupus (age ~2 Myr, Testi et al. 2022), Chamaeleon I (age ~3 Myr, Testi et al. 2022) and σ Orionis (age ~3–5 Myr, Ansdell et al. 2017) show a wide spread in accretion timescales Ṁ_acc/M_disc between approximately 0.1 and 10 Myr. Lupus and Chamaeleon I are known to have weak FUV field strengths⁵ and we thus compare them to our weak FUV field populations. We note that σ Orionis does contain an OB system (σ Ori) and its FUV fluxes are expected to be high (~8000 G₀ within 1 pc of σ Ori, according to Ansdell et al. 2017). Therefore, the strong FUV field populations were compared to data from σ Orionis. We note that the constituents of the star forming regions show a large spread in ages (according to Table 1 of Testi et al. 2022).

A good agreement with the observational data was obtained only for a weak magnetic disc wind. However, accretion rates beyond 10⁻⁸ M_⊙ yr⁻¹, as observed in the young star-forming regions Lupus and Chamaeleon I, were not reproduced by our model.

Fig. 6 Stellar accretion rate Ṁacc vs. gas disc mass Mdisc, analogous to Fig. 5, but for populations with a strong FUV field distribution (see Table 1). A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

3.3 Disc lifetime

As we have seen, the combination of stellar accretion, photoevaporation, and magnetically driven disc winds can lead to a rapid dispersal of the disc. We consider a disc to be dispersed when it becomes unobservable in the near-infrared (NIR). Following the dispersal condition of Kimura et al. (2016) the disc is considered unobservable in NIR when the optical depth τ drops below unity in the regions where the midplane temperature is 7_mid > 300K.

Figure 8 shows the time evolution of the fraction of stars possessing a disc. An exponential decay of the following form is typically used, (19)

to describe the evolution, with f₀ being the initial fraction of stars possessing a disc and τ_disc being the characteristic time scale⁶. We show fits from Mamajek et al. (2009); f₀ = 1 and τ_disc = 2.5 Myr, Richert et al. (2018); f₀ = 1 and τ_disc = 5 Myr and Michel et al. (2021); f₀ = 0.8 and τ_disc = 6.4 Myr. Furthermore, Richert et al. (2018) made use of magnetic pre-main sequence (PMS) models from Feiden (2016) for determining the cluster ages and we thus refer to it as Feiden16M. Both Mamajek et al. (2009) and Richert et al. (2018) have assumed that initially all stars possess discs, whereas Michel et al. (2021) assumed an initial fraction of f₀ = 0.8 to account for close binary systems (Kraus et al. 2012). The differences in characteristic time scales, τ_disc, emerge partially from different modelling approaches for the ages of the stellar clusters and partially from the choice of clusters considered. While the sample of Mamajek et al. (2009) also includes clusters containing O/B type stars (e.g. σOri, λOri, 25Ori), Michel et al. (2021) intentionally excluded them, since the presence of such stars implies a strong FUV field, which can have significant impact on the disc lifetime (see below and Fig. 7).

Our populations show very rapid dispersal for strong disc winds, regardless of the ambient FUV field distribution. Here, the lifetimes are below the values from observations. For weak disc winds, disc lifetimes are longer and show a dependence on the ambient FUV field distribution. For a strong FUV field, the fraction of stars with a protoplanetary discs roughly follow the data from Mamajek et al. (2009), whereas for a weak FUV field, the fraction of stars declines much more slowly, namely, in a way that is comparable to the data from Richert et al. (2018) and Michel et al. (2021). The dependency on the strength of the magnetic disc wind can be explained by the dispersal condition we used, where the evolution of the inner disc is the most important factor. When assuming strong disc winds, the decrease in surface density is much faster. When the wind gets weaker, internal photoevaporation sets in, preventing the inner disc from being replenished with gas, leading to rapid dispersal (Figs. 3 and 4). Again, we stress that in our model, the resulting disc lifetimes are a direct prediction of the model. This is in contrast to past approaches (e.g. Mordasini et al. 2009; Emsenhuber et al. 2021b), where the external photoevaporation rates were adjusted such that disc lifetimes were in agreement with observations.

As mentioned previously, the presence of nearby massive stars can have a large impact on the lifetime of circumstellar discs (e.g. Adams et al. 2004, 2006; Clarke 2007; Fang et al. 2012). Figure 7 shows the dependence of the lifetime on the ambient FUV field strength and stellar mass for two example disc wind scenarios of populations with strong FUV field strength distributions: strong disc wind with constant torque and weak disc wind with constant torque. Populations with weak FUV field distributions can be regarded as a subset (1 G₀ < ℱ_FUV < 10² G₀) of these simulations.

There is a strong dependence of the disc lifetime on the ambient FUV field strength for weak disc winds. Also the mass of the host star shows impact on the disc lifetime, with lifetime tending to increase with higher stellar mass.⁷ This is in contrast with the results of observations, which show that disc fractions decrease with increasing stellar mass (e.g. Bayo et al. 2012). However, Komaki et al. (2021) and Picogna et al. (2021) both suggested that this is a result of strong internal photoevaporation. For strong disc winds, the lifetime is only correlated for very high FUV field strengths (> 10³ G₀). This can explain why the disc fractions for strong disc winds in Fig. 8 depend only weakly on the FUV field distribution.

Fig. 7 Disc lifetime with regard to the ambient FUV field strength ℱFUV· Point sizes indicate the mass of the host star. We note that our grid only extends to 104 G0 and values beyond are treated as the respective boundary value, which is also seen in the results, given that beyond 104 G0, the lifetimes do not further decrease.

Fig. 8 Fraction of stars possessing a circumstellar disc as a function of time for both weak and strong FUV field distributions. Evolution tracks of the simulations are shown in comparison with exponential decay fits to observational data from Mamajek et al. (2009); Richert et al. (2018); Michel et al. (2021). Note: the fit from Richert et al. (2018) is denoted Feidenl6M, since it relies on magnetic PMS models from Feiden (2016) to determine the cluster ages.

3.4 Additional cases: Variable

Our simulations still fall short of the high accretors, especially in the case of a weak FUV field (Fig. 5). In order to investigate the dependency on the initial torque strength, we ran additional populations with distributed uniformly in log between 10^−5.5 and 10^−2.5 for two fixed FUV field strengths, 10 G₀ and 10³ G₀. The results are shown in Appendix A.

The tracks in Fig. A.2 show that high accretors can be reached by high initial torques. However, when comparing the simulations after a few million years with the observations, regions where Ṁ_acc > 10⁻⁸ M_⊙yr⁻¹ and M_disc/Ṁ_acc ≲ 0.1 Myr are almost never achieved. This has to do with the fact that discs beyond M_disc/Ṁ_acc ~ 0.1 Myr are dispersed very quickly. The simulation tracks show a strong correlation with the initial torque strength and the large spread in initial torque strength results in a large spread in accretion rates, similar to observations. The majority of simulations with a Σ-dependent torque stay at almost constant accretion rates and move horizontal in the log (M_disc)-log(Ṁ_acc) plane.

Regarding the disc lifetimes, the choice of the torque (i.e. constant or Σ-dependent) has less impact than the choice of the disc wind scenario (strong or weak). Discs exposed to a strong disc wind are short lived, which is in line with our previous findings.

3.5 Additional cases: Varying internal photoevaporation

Our model incorporates both MHD winds and photoevaporation. The interplay between these two types of outflows is currently still under debate (see Sect. 4.2). In order to take this into account, we ran the two limiting cases with ‘No EUV radiation shielding’ and ‘No internal photoevaporation’. We show the main results from Sect. 3.2 and Sect. 3.3 in Appendices B and C.

We find that internal photoevaporation only plays an important role at a late stage of the disc evolution by reducing the accretion rate and opening the inner cavity. In the absence of internal photoevaporation (Appendix C), discs do not open an inner cavity. Thus, the star remains accreting gas and, as a result, the NIR lifetimes are much longer. This can be seen also in the evolutionary paths in the M_disc – Ṁ_acc plane. Instead of having a gap forming, simulations follow lines of constant accretion or are dispersed outside-in. Also note that in the case of a constant magnetic field (Σ-dependent torque), simulations tend to have constant accretion towards the end, as a result of the increasing torque. However, it does not change the location where we find the simulations in the M_disc-Ṁ_acc plane and, in addition, the best correspondence is still found for weak disc winds with a strong torque. Thus, the presence of internal photoevaporation has mainly influence on the NIR lifetimes. This necessity of the interplay between MHD wind and internal photoevaporation to explain the inner disc lifetime has already been reported by Kunitomo et al. (2020).

4 Discussion

Our results give some first insights on the interplay of magnetically driven disc winds, magnetic braking, and internal and external photoevaporation in MRI inactive discs, as well as their influences on and their agreement with observables. However, there are several simplifications and caveats that have to be discussed.

4.1 Stellar accretion

Our model considers magnetospheric accretion driven by turbulence (MRI) and magnetic braking (Eq. 18). We find that stellar accretion makes up only a minor contribution to disc dispersal, which is in contrast to the common point of view (e.g. Ercolano & Pascucci 2017). We obtained upper limits on the stellar accretion rates of 10⁻¹⁰ M_⊙ yr⁻¹ to 10⁻⁸ M_⊙ yr⁻¹, dependent on the disc wind scenario, which is on the lower side of what we expect from observations. Especially the young star forming regions Lupus and Chamaeleon I show stellar accretion rates beyond 10⁻⁸ M_⊙ yr⁻¹. Higher magnetic torques could explain these high accretors, but it is still difficult to find such systems after a few million years. More massive and compact discs would also lead to higher accretion rates.

There may be other processes at work that could contribute to the accretion rates observed. Takasao et al. (2018) conducted 3D MHD simulations of magnetised accretion discs and found that a failed magnetic disc wind can drive fast accretion onto high latitudes, namely, so-called funnel-wall accretion. This new accretion process can coexist with the usual accretion processes, but is expected to drive much lower accretion rates. Thus, the discrepancy between our simulations and observed accretion rates above ≳10⁻⁸ M_⊙ yr⁻¹ persists.

We also note that some observations of transitional discs show relatively high accretion rates (~ 10⁻⁸ M_⊙ yr⁻¹, Owen 2016), despite the presence of an inner cavity. Wang & Goodman (2017) showed that such high accretion rates could, in principle, be obtained by a strong magnetised wind for moderately good coupling of the magnetic field. In our model, such high accretion rates can no longer be supported after opening the cavity.

4.2 Magnetic versus photoevaporative winds

In our model, we assume the presence of both magnetic and photoevaporative winds. By calculating the column density, we take into account the shielding of the disc from EUV radiation by the magnetic wind. As soon as the column density falls below a critical value, internal photoevaporation is added to the outflows. In reality, the interplay between MHD-driven and photoevaporative winds is expected to be much more complicated. Bai et al. (2016) found a magneto-thermal wind (also known as magneto-photoevaporation) which is a thermally assisted magnetic outflow. Wang et al. (2019) found that for magnetic winds, the mass loss rates are relatively robust with regard to high-energy photon luminosities and EUV photons can even reduce the mass-loss rates. At this point, the transition from MHD wind to photoevaporation driven outflow remains unsolved. The present opinion tend to give both processes important roles at different evolutionary stages (see e.g. recent reviews by Pascucci et al. 2022; Lesur et al. 2022).

Our simulations show that at a given time, the evolution of large parts of the disc are either dominated by magnetic winds or internal photoevaporation. This is shown by Fig. 9, which compares the mass-loss rates and strengths of different winds by location and time. As depicted in Fig. 9 in the early phase, magnetic winds dominate the mass-loss rates, whereas at a later stage, the internal photoevaporation gains importance. This aspect was already reported on and discussed by Kunitomo et al. (2020). Although external photoevaporation is strong initially, its mass-loss rates are very high and concentrated in the outer disc only, where they dominate magnetic winds by several orders of magnitude (as shown in Figs. 3 and 4). In Appendices B and C, we show the two extreme cases where we either add up MHD wind and internal photoevaporation or completely disable internal photoevaporation throughout the entire disc evolution. We see that the internal photoevaporation plays an important role in setting the disc life time, but it does not have big influence on other observables. We therefore do not expect our results to be very strongly affected by the interplay of magnetic and photoevaporative winds.

Fig. 9 Relative importance of different mass-loss processes for the exemplary case with weak DW, strong constant torque and ℱFUV = 10 G0 (bottom row of Fig. 3). Top: time evolution of mass-loss rates for different processes. Purple dashed line distinguishes MHD-wind dominated and internal PEW dominated period (analogous to Fig. 2c in Kunitomo et al. 2020). Bottom: ratios of Σ for different mass loss processes shown for radial extend and time. Blue areas are dominated by MHD-wind while yellow areas indicate PEW dominated. The hatched area highlight external PEW-dominated regions.

4.3 Dust evolution

Our model does not include dust evolution and we directly compare obtained gas disc masses with observed dust masses converted to gas mass, with f_D/G = 1/100. However, Sellek et al. (2020b) found that growth and radial drift of dust can lower the dust-to-gas ratio significantly over time. This implies that gas disc masses inferred from sub-mm observations stand as the lower limits, rather than actual disc masses. They further note that including dust evolution can also increase the scatter in accretion rates.

4.4 Disc dispersal

We essentially observed two different disc dispersal modes (i.e. outside-in and inside-out). Our simulations suggest that the dispersal mode is strongly correlated with the external FUV field strength. Furthermore, it depends on the disc wind scenario chosen, as a strong wind removes mass at a higher rate from the inner disc than a weak wind. Koepferl et al. (2013) classified over 1500 sources in nearby star-forming regions and found the inside-out clearing of discs to be the preferred way of disc dispersal. In our simulations, systems in a weak ℱ_FUV environment survive long enough, so that an inner cavity will eventually open up. As can be seen in Appendix C, the opening of an inner cavity is strongly dependent in the internal photoevaporation.

We would also like to point out that outside-in dispersal via external photoevaporation happens to be a very efficient process and disc lifetimes in strong FUV environments are often below 1 Myr, which would suggest that only very few disc are caught by observations during outside-in dispersal.

We also neglected X-ray as source of internal photoevaporation. While the EUV prescription we use has a integrated mass loss rate of about 10⁻¹⁰ M_⊙ yr⁻¹ (Fig. 9), X-ray internal photoevaporation can produce mass loss rates of the order of 10⁻⁸ M_⊙ yr⁻¹ (Owen et al. 2010, 2011, 2012). Including X-ray photoevaporation in our model would considerably increase mass loss rate by internal photoevaporation, which would cause a reduction of loss by MHD winds and stellar accretion. This would, in turn, make it even more difficult to match the observed stellar accretion rate and would lead to further reduced disc lifetimes.

4.5 Comparison with previous works

There are several models describing MHD wind-driven disc evolution in 1D (e.g. Armitage et al. 2013; Bai 2016; Suzuki et al. 2016; Tabone et al. 2022a). Our model closely follows Kunitomo et al. (2020), who combined MHD winds based on Suzuki et al. (2016) with internal EUV and X-ray photoevaporation, but with internal EUV and external FUV instead of internal X-ray photoevaporation, along with a consistent treatment of EUV shielding. We also find that disc dispersal is a result of both magnetic winds and photoevaporation, since without internal photoevaporation and considering weak external FUV fields, disc lifetimes are much longer than inferred from observations (Appendix C). This is in contrast to Armitage et al. (2013), Bai (2016), and Tabone et al. (2022b) who suggested that disc dispersal could be driven solely by magnetic disc winds.

The approach of running disc populations syntheses has gained popularity in recent years (e.g. Lodato et al. 2017; Mulders et al. 2017; Tabone et al. 2022b; Somigliana et al. 2022). However, models and approaches differ. Lodato et al. (2017) investigated the self-similar solution of the classical α-disc model (Lynden-Bell & Pringle 1974) for varying initial disc masses and viscous times. They were able to reproduce the correlation between accretion rate and disc mass for the Lupus star-forming region under the assumption that the efficiency of angular momentum transport is an increasing function of radius. Tabone et al. (2022b) managed to reproduce the Ṁ_acc – M_disc correlation and the spread around the mean trend in Lupus, as a result of MHD-wind driven disc evolution. They used a distribution of accretion timescales that was obtained by fitting observed disc fractions. This allowed them to obtain disc lifetimes that comply with observations without considering photoevaporation. We note that variations in accretion timescales is equivalent to variations in their α_DW (similar to our ) and naturally lead to a larger spread in accretion rates.

In our fiducial simulations, the results (i.e. spread in disc lifetimes and Ṁ_acc-M_disc evolution) are a direct consequence of the variations in initial conditions (i.e. stellar mass, disc properties and ambient FUV-field strength). In contrast to the models mentioned above, we do not assume variations in accretion or viscous timescales (i.e. variations in α). In Appendix A, we show additional cases, where we varied the initial torque strength , which introduced a large spread in observed accretion rates (see also Sect. 3.4 for a discussion).

5 Summary and conclusions

We constructed a global model for protoplanetary disc evolution including internal redistribution of angular momentum via a turbulent viscosity (α-disc) and angular momentum removal through magnetic braking, magnetically driven disc winds, and internal and external photoevaporation. We further developed a new, physically consistent (albeit simple) model for the shielding of EUV-driven photoevaporation by MHD winds. Taking into account the recent discussion of non-ideal MHD effects in the literature and suppressing MRI in large parts of the disc, we studied the impact of four different magnetic disc wind scenarios on MRI inactive discs for a wide range of initial conditions, resembling conditions found in young star-forming regions. We then compared observables such as stellar accretion rates, disc mass, disc dispersal mode, and disc lifetime with observational data.

We present the following conclusions based on our model:

1. We find that discs are primarily dispersed by the combination of outflows (magnetically driven disc winds, as well as internal and external photoevaporation), while only a minor fraction of the mass is accreted onto the star;

2. Weak wind-driven mass-loss rates (weak disc winds) are favoured over high mass-loss rates (strong winds) to reproduce the results of observations. Strong magnetic disc winds combined with internal photoevaporation lead to fast depletion of the inner disc region, the opening of an inner cavity and a swift decrease in stellar accretion rate;

3. A strong torque can support stellar accretion in MRI inactive discs up to Ṁ_acc ~ 10⁻⁸ M_⊙ yr⁻¹ or higher, as observed. Our comparison with observations therefore indicates that we need weak MHD-driven mass loss, but strong torques for agreement. This allows us to identify the subtype of MHD-wind models (weak disc wind + strong torque) that best reproduces the observations;

4. The inclusion of EUV internal photoevaporation only affects the late stage of the disc evolution, where it can lead to rapid disc dispersal. This has an important influence on disc lifetimes in our simulations;

5. We found that for weak disc winds, the disc lifetimes are strongly influenced by the ambient FUV field strength, driving external photoevaporation. Initial disc mass and stellar mass seem to play merely secondary roles. Discs under the influence of weak disc winds have lifetimes which are in line with observations for both weak- and strong ambient FUV field strengths. For strong disc winds, disc lifetimes are too short when compared with observations of NIR excess.

These findings support the recent shift towards magnetically driven accretion. However, there is further research underway towards improving our understanding the dynamics of the magnetic field evolution and more detailed modelling of the strength and evolution of the disc wind torque needed. Further investigation of the interplay between magnetic and photoevaporative winds is highly encouraged as both processes are expected to play important roles in disc evolution.

Future works will address planet formation in these wind-driven accretion discs and possible imprints on planet populations. Regarding external photoevaporation, we now have a model that is dependent on the FUV field strength, which is considered one of the most important quantities that affects disc evolution in star-forming clusters. This will enable us to also investigate the influence of the cluster environment on planetary formation, which is a hot topic in recent discussions (e.g. Winter et al. 2020; Adibekyan et al. 2021).

Acknowledgements

We want to thank Oliver Schib, Remo Burn, Thomas Haworth, and Ilaria Pascucci for constructive discussions. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF) under grants 51NF40_182901 and 51NF40_205606. The authors acknowledge the financial support of the SNSF. J.W. and C.M. acknowledge the support from the SNSF under grant 200021_204847 “PlanetsInTime”. We thank the anonymous referee for a thorough read and a constructive report that was very helpful for improving the manuscript.

Appendix A Variable

Here, we show additional cases with varying initial torque strength, , to assess its influence on disc evolution on a population level. It is varied as an initial condition, with values distributed uniformly in log between 10^−5.5 and 10^−2.5. We further evolved populations in a weak and a strong FUV field environment with values of 10G₀ and 10³G₀, respectively. We chose not to vary the FUV field alongside with in order to get a better understanding on how the torque strength influences the evolution. We show disc lifetimes in Table A.1 and evolutions for both FUV field strengths in Figs. A.2 and A.3. The characteristic numbers of the populations are shown in Tables A.1 and A.2.

The evolution tracks of the simulations show a clear separation by initial torque strength. Here, even high accretors are reached by simulation tracks with high initial torques. However, we still did not recover simulations at these high accretion rates after a few million years (i.e. see snapshots), especially beyond the 0.1 Myr line of constant accretion. Discs beyond this point are rapidly dispersed. The simulation tracks show a much larger spread in accretion rates, more similar to the observations. This comes to no surprise, as the has a direct influence on the accretion rate. Looking at the characteristic numbers of the 10G₀ populations in Table A.1, we also observe fewer discs opening an inner cavity for Σ-dependent torque cases (23% to 35%), compared to the constant torque cases (44% to 57%). The same can be seen in the 10³G₀ populations (Table A.2). The torques start initially with the same strength, namely, . As the surface density declines, the torques get stronger in the Σ-dependent cases, leading to higher accretion rates. This makes it less likely for internal photoevaporation opening up a gap.

The evolutions of disc fractions show similar features as our nominal cases in Fig. 8. Strong disc winds have short lifetimes, regardless of the ambient FUV field strength. Discs exposed to weak winds have slightly longer lifetimes and show a dependency on the ambient FUV field. Lifetimes of weak disc wind scenarios for the 10G₀ cases are in line with observations. We also note that for strong and weak disc winds, the disc fractions evolve in pairs, regardless of the torque scenario.

Fig. A.1 Fraction of stars possessing a circumstellar disc as a function of time, analogous to Fig. 8, but for populations with varying . Top panel shows the disc fraction evolution for a weak 10G0 field, while the bottom panel shows the same evolution for a 103G0 FUV field.

Fig. A.2 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 5, but for an external FUV field of 10G0 and varying . Evolution tracks of individual systems are coloured by . Snapshots at 2 Myr (green circles) and 3 Myr (orange circles) are displayed. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant are shown for 0.1, 1 and 10 Myr (thin dotted lines).

Fig. A.3 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 6, but for an external FUV field of 103G0 and varying . Evolution tracks of individual systems are coloured by . A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

Appendix B EUV photoevaporation always active (i.e. no shielding)

To investigate the effect of our prescription for the shielding of internal photoevaporation by MHD disc winds, we compare the nominal results with where internal photoevaporation is active throughout the whole evolution of the disc. We provide the disc lifetimes in Fig. B.1 and the mass and accretion rate evolution for the case with weak FUV field in Fig. B.2, with the quantitative results in Table B.1, while the same for the strong FUV environment are provided in Fig. B.3 and Table B.2.

The results are almost identical to those presented in the main text. The only noticeable differences are: slightly shorter disc lifetimes, a lower number of discs with non-zero stellar accretion rate, and a higher percentage of mass removal by internal photo-evaporation. Nevertheless, the main results are unaffected by the inclusion of internal photoevaporation through the disc lifetimes. We still obtain the result that weak disc winds in combination with a strong constant torque are necessary to reproduce stellar accretion rates comparable with observations and observed disc lifetime distributions. This confirms that our prescription for internal photoevaporation is sufficiently weak to avoid affecting disc evolution before the dispersal stage.

Fig. B.1 Fraction of stars possessing a circumstellar disc as a function of time, analogous to Fig. 8, but for populations with EUV photoevapo-ration always active (i.e. no shielding).

Fig. B.2 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 5, but without EUV radiation shielding. Exemplary cases are the same as in Figs. 3 and 4, but without EUV radiation shielding. Snapshots at 2 Myr (green circles) and 3 Myr (orange circles) are displayed. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant are shown for 0.1, 1 and 10 Myr (thin dotted lines).

Fig. B.3 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 6, but without EUV radiation shielding. Exemplary cases are the same as in Figs. 3 and 4 but without EUV radiation shielding. A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

Appendix C No internal photoevaporation

To asses how internal photoevaporation affect the dispersal of protoplanetary discs, we computed the other extreme case where internal photoevaporation has been removed altogether. We proceeded as detailed in the main text and Appendix B, namely, by analysing disc lifetimes and the evolution of disc masses and stellar accretion rates.

In Fig. C.1, we show the time evolution of the fraction of stars showing NIR signatures of a disc. When comparing to the simulations including internal photoevaporation it is evident that NIR lifetimes are generally longer without internal photoevaporation. This comes as no surprise as internal photoevaporation is removing mass from the disc which is then lacking to resupply the inner disc. In our simulations, discs are primarily dispersed by external photoevaporation, where as for the weak fUv field, MHD winds remove substantial fractions of gas. This combined with the longevity of the discs, we come to the conclusion that disc lifetimes cannot be reproduced by MHD winds only for the chosen initial conditions and disc wind scenarios. However, things may look different for stronger disc wind torques or radially dependent .

Analogously to Sect. 3.2, we show evolution tracks for disc populations in the plane for different disc wind scenarios and ambient FUV field strength distributions in Figs. C.2 and C.3. The corresponding characteristic numbers are listed in Tables C.1 and C.2.

The most striking difference is that the stellar accretion rates stay high during the evolution of the disc and do not drop while the disc is still massive. This is because the internal photoe-vaporation is largely responsible for the opening of the inner cavity – and it is the opening of the inner cavity that leads the stellar accretion rates to drop sharply. We note that inner cavities are only opened in the case of a strong wind paired with a Σ-dependent torque. The main difference between these populations and the ones shown in the main text is the sudden drop in accretion rate when it reaches ~ 10⁻¹¹ M_⊙/yr. In the absence of internal photoevaporation, no such drop is observed. However, the location of the synthetic discs in the plane in the snapshots does not change significantly in the absence of internal photoevaporation. As discussed above, disc lifetimes are strongly affected by the internal photoevaporation. As a consequence, some of the discs in the populations without internal photoevaporation exhibit a more extended disc dispersal stage, with relatively low masses and stellar accretion rates.

We conclude that in our model, internal photoevaporation (i.e. EUV photoevaporation only) only influences the late stage of protoplanetary disc evolution. It therefore has a big influence on the presence of inner cavities and disc lifetimes, but it does not have strong influence on the location where we find snapshots of our populations in the plane. Thus, our findings are not affected by the inclusion of internal photoevaporation.

Fig. C.1 Fraction of stars possessing a circumstellar disc as a function of time, analogous to Fig. 8, but for populations without internal photo-evaporation.

Fig. C.2 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 5, but without internal photoevaporation. Exemplary cases are the same as in Figs. 3 and 4 but without internal photoevaporation. Snapshots at 2 Myr (green circles) and 3 Myr (orange circles) are displayed. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant are shown for 0.1, 1 and 10 Myr (thin dotted lines).

Fig. C.3 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 6, but without internal photoevaporation. Exemplary cases are the same as in Figs. 3 and 4 but without internal photoevaporation. A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

References

All Tables

All Figures

	Fig. 1 Adopted distributions of the Monte Carlo variables for our disc populations: dust-to-gas ratio (panel a), initial dust masses (panel b), resulting gas mass (panel c), inner disc radius (panel d), FUV field strength for weak-FUV and strong-FUV environments (panel e), and stellar mass (panel f).
In the text

	Fig. 2 Initial relation of disc dust mass and characteristic radius for a synthetic population of 1000 discs, following Fig. 12 of Tobin et al. (2020).
In the text

Fig. 3 Exemplary disc evolution shown for initial conditions: Mdisc = 0.1 M⋆ with M⋆ = 0.25 M⊙, Rin = 0.04 AU and ℱFUV = 10 G0. We show temporal evolution of the surface density, radial mass flow rate and rate of change in surface density due to outflows (i.e. magnetic disc winds, internal- and external photoevaporation) for all disc wind scenarios. The individual contributions can be distinguished roughly when considering the outflows . External photoevaporation corresponds to the outermost peak while internal photoevaporation acts at ~ 0.3 AU to a few AU. MHD winds remove mass from the whole disc and increase in strength towards the inner disc. Dotted lines show the evolution at 104 yr, dashed lines at 5 × 104 yr. Coloured lines are spaced by 105 yr.

In the text

	Fig. 4 Exemplary disc evolution, analogous to Fig. 3 but for a strong ambient FUV field strength of 5000 G0.
In the text

Fig. 5 Stellar accretion rate Ṁacc vs. gas disc mass Mdisc shown for all disc wind scenarios (see Table 1) with a weak FUV field strength distribution. Evolution tracks of individual systems, coloured by ℱFUV, and snapshots at 2 Myr (green circles and 3 Myr (orange circles are displayed. Tracks of exemplary cases are highlighted with black dashed and dotted lines. We show the detailed evolution of those exemplary discs in Figs. 3 and 4. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant Mdisc/Ṁacc are shown for 0.1, 1 and 10 Myr (thin dotted lines).

In the text

Fig. 6 Stellar accretion rate Ṁacc vs. gas disc mass Mdisc, analogous to Fig. 5, but for populations with a strong FUV field distribution (see Table 1). A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

In the text

	Fig. 7 Disc lifetime with regard to the ambient FUV field strength ℱFUV· Point sizes indicate the mass of the host star. We note that our grid only extends to 104 G0 and values beyond are treated as the respective boundary value, which is also seen in the results, given that beyond 104 G0, the lifetimes do not further decrease.
In the text

Fig. 8 Fraction of stars possessing a circumstellar disc as a function of time for both weak and strong FUV field distributions. Evolution tracks of the simulations are shown in comparison with exponential decay fits to observational data from Mamajek et al. (2009); Richert et al. (2018); Michel et al. (2021). Note: the fit from Richert et al. (2018) is denoted Feidenl6M, since it relies on magnetic PMS models from Feiden (2016) to determine the cluster ages.

In the text

Fig. 9 Relative importance of different mass-loss processes for the exemplary case with weak DW, strong constant torque and ℱFUV = 10 G0 (bottom row of Fig. 3). Top: time evolution of mass-loss rates for different processes. Purple dashed line distinguishes MHD-wind dominated and internal PEW dominated period (analogous to Fig. 2c in Kunitomo et al. 2020). Bottom: ratios of Σ for different mass loss processes shown for radial extend and time. Blue areas are dominated by MHD-wind while yellow areas indicate PEW dominated. The hatched area highlight external PEW-dominated regions.

In the text

	Fig. A.1 Fraction of stars possessing a circumstellar disc as a function of time, analogous to Fig. 8, but for populations with varying . Top panel shows the disc fraction evolution for a weak 10G0 field, while the bottom panel shows the same evolution for a 103G0 FUV field.
In the text

Fig. A.2 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 5, but for an external FUV field of 10G0 and varying . Evolution tracks of individual systems are coloured by . Snapshots at 2 Myr (green circles) and 3 Myr (orange circles) are displayed. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant are shown for 0.1, 1 and 10 Myr (thin dotted lines).

In the text

Fig. A.3 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 6, but for an external FUV field of 103G0 and varying . Evolution tracks of individual systems are coloured by . A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

In the text

	Fig. B.1 Fraction of stars possessing a circumstellar disc as a function of time, analogous to Fig. 8, but for populations with EUV photoevapo-ration always active (i.e. no shielding).
In the text

Fig. B.2 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 5, but without EUV radiation shielding. Exemplary cases are the same as in Figs. 3 and 4, but without EUV radiation shielding. Snapshots at 2 Myr (green circles) and 3 Myr (orange circles) are displayed. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant are shown for 0.1, 1 and 10 Myr (thin dotted lines).

In the text

Fig. B.3 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 6, but without EUV radiation shielding. Exemplary cases are the same as in Figs. 3 and 4 but without EUV radiation shielding. A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

In the text

	Fig. C.1 Fraction of stars possessing a circumstellar disc as a function of time, analogous to Fig. 8, but for populations without internal photo-evaporation.
In the text

Fig. C.2 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 5, but without internal photoevaporation. Exemplary cases are the same as in Figs. 3 and 4 but without internal photoevaporation. Snapshots at 2 Myr (green circles) and 3 Myr (orange circles) are displayed. We compare our simulations with observed populations in Lupus and Chamaeleon I. Observational dust disc masses and stellar accretion rates are taken from Manara et al. (2019) and dust masses are converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Triangles denote upper limits on disc mass. Lines of constant are shown for 0.1, 1 and 10 Myr (thin dotted lines).

In the text

Fig. C.3 Stellar accretion rate vs. gas disc mass Mdisc, analogous to Fig. 6, but without internal photoevaporation. Exemplary cases are the same as in Figs. 3 and 4 but without internal photoevaporation. A snapshot for the individual systems is shown at 4 Myr (orange circles) in comparison with observational data from σ Orionis. Observational dust disc masses are taken from Ansdell et al. (2017) and converted to gas masses by assuming the standard dust-to-gas ratio of 0.01. Stellar accretion rates are taken from Rigliaco et al. (2011).

In the text