Towards a population synthesis of discs and planets II. Confronting disc models and observations at the population level

,


Introduction
Protoplanetary discs are the birthplace of planets (the 'nebular hypothesis' of Kant and Laplace).Discs serve as a source of gas and solids from which the planets accrete.Planet-disc interactions lead to planetary migration.To model planetary formation, it is therefore essential to have disc characteristics that are as close as possible to those observed to provide the highest possible fidelity.
Disc observations are not an entirely new subject of research.Disc masses (e.g.Beckwith & Sargent 1996;Andrews et al. 2009) and lifetimes (e.g.Haisch et al. 2001;Mamajek 2009;Fedele et al. 2010;Kraus et al. 2012;Ribas et al. 2014) have been observed for over two decades.However, there have been many new results concerning protoplanetary discs in the last several years, including the mass and physical extent of early discs (Tychoniec et al. 2018(Tychoniec et al. , 2020;;Tobin et al. 2020) and at later times (Hendler et al. 2020).
Nevertheless, some aspects of disc evolution are not captured by observations, such as the process that leads to transport of material.These are usually taken to be turbulent viscosity Tables 3 and 4 are only available in electronic form at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/cgi-bin/qcat?J/A+A/ generated by the magnetorotational instability and magnetically driven disc winds (Suzuki & Inutsuka 2009;Suzuki et al. 2016).The strength of the turbulent viscosity has not yet been properly determined and is usually parametrised using a factor α (Shakura & Sunyaev 1973).
There are indirect methods to estimate the value of α.Ultraviolet(UV)-excess measurements of the accretion luminosity were used to derive the accretion rate onto the star for the Chamaeleon I (Manara et al. 2016a(Manara et al. , 2017) ) and Lupus (Alcalá et al. 2014(Alcalá et al. , 2017) ) star forming regions.These measurements coupled to the disc masses for the same regions of Cha I (Pascucci et al. 2016) and Lupus (Ansdell et al. 2016) provide a relationship between mass and accretion rate (Manara et al. 2016b).Together, these can be used to calibrate numerical models (Manara et al. 2019) and to provide an estimate of the mass flux onto the disc (Mulders et al. 2017;Sellek et al. 2020).A second method to estimate α is to compare dust and gas emission, either using spatially resolved observations of disc substructures from ALMA (Andrews et al. 2018a) such as pressure bumps (Dullemond et al. 2018) or from the overall disc sizes (Toci et al. 2021).
Mass loss does not occur only due to accretion onto the star.For instance, observations also point towards protoplanetary disc dispersal occurring from the inside out and on relatively short timescales (Ercolano et al. 2011(Ercolano et al. , 2015;;Koepferl et al. 2013).
This suggests there is an additional mechanism removing gas close to the star, with one possibility being internal photoevaporation.Coupled with the findings that young stars emit a larger fraction of their flux in UV (e.g.Gómez de Castro 2009) and X-rays (e.g.Preibisch et al. 1996Preibisch et al. , 2005;;Feigelson & Montmerle 1999;Favata & Micela 2003), it is proposed that extreme UV (EUV; Hollenbach et al. 1994;Clarke et al. 2001) and/or Xrays (e.g.Alexander et al. 2004;Ercolano et al. 2008Ercolano et al. , 2009) ) are responsible for this mass loss.Using hydrodynamical simulations, it is possible to predict the mass-loss rate as a function of disc properties and stellar luminosity (Owen et al. 2012;Picogna et al. 2019Picogna et al. , 2021;;Ercolano et al. 2021).
Nevertheless, irradiation from the host star is not the only mass-loss mechanism: most stars are born in clusters where many stars form concurrently.Consequently, protoplanetary discs are exposed to a larger ambient radiation field than mature stars.This leaves an additional mechanism for mass removal by external photoevaporation (e.g.Matsuyama et al. 2003;Winter & Haworth 2022).This is supported by observational findings that discs near massive stars have lower masses than others (Ansdell et al. 2017;van Terwisga et al. 2020) and that clusters with a low ambient radiation field have longer disc lifetimes (Michel et al. 2021).As for external photoevaporation, hydrodynamical simulations were performed to predict mass-loss rate (Haworth et al. 2018) as function of disc properties and ambient flux.Together with simulations of cluster evolution (e.g.Adams et al. 2006;Qiao et al. 2022), this enables us to determine the mass-loss rate over an entire disc population.
All these observations and theoretical predictions put a lot of constraints on protoplanetary disc evolution, as the number of free parameters is limited.Whether or not the combination of initial disc properties and predicted accretion and mass-loss rates can be used to reproduce the distribution of, for instance, disc lifetimes remains to be determined.Previous studies in this direction usually consider only one type of photoevaporation, either internal (Gorti et al. 2009;Owen et al. 2011;Kunitomo et al. 2021) or external (Kunitomo et al. 2020).Burn et al. (2022, hereafter Paper I) introduced a relatively simple 1D radial disc model that is capable of consistently evolving gas, dust, pebbles, and planetesimals.In addition, this model is capable of predicting how the modelled disc would be observed by current instrumentation, such as ALMA (see Birnstiel et al. 2018).Further, the light computational requirements of that model make it possible to perform many such evolutions in order to study the effects of initial disc properties.
Our goals are twofold: first, we aim to determine whether we can understand the general picture of protoplanetary discs set by the observations and predictions highlighted above.Second, we want to find the combinations of disc properties that best reproduce the various observations by performing disc populations synthesis, such as Rosotti et al. (2017) and Somigliana et al. (2022).This should then serve as initial conditions for future planetary population syntheses, such as in Mordasini et al. (2009) or Emsenhuber et al. (2021b).
To fit the best parameters, many calculations need to be performed, each involving the evolution of a population of protoplanetary discs.To alleviate the computational requirements of this procedure, we use machine learning to fit neural networks that can reproduce the result of the underlying model with limited resources (Cambioni et al. 2019a(Cambioni et al. , 2021;;Emsenhuber et al. 2020).This 'surrogate model' can then be used as the forward model in the fitting procedure (Cambioni et al. 2019b).
In this work, we aim to find initial conditions for the disc evolution calculations that best match observations.For this pur-pose, we first compute two series of calculations using the model presented in Paper I (Sect.4.1).These data are then used to fit several surrogate models that hold the necessary outcomes for comparison with observations (Sect.4.2).Using these surrogate models, we study the effect of the photoevaporation prescriptions (Sect.4.3) and find initial conditions that best match the observational constraints discussed in Sect. 3 as a whole (Sect.5).A study dedicated to this last aspect using a Bayesian approach instead will be presented in Burn et al. (in prep., hereafter Paper III).

Methods
The disc evolution model is based on the Bern global model of planetary formation and evolution (e.g.Alibert et al. 2004Alibert et al. , 2005;;Mordasini et al. 2009;Fortier et al. 2013;Voelkel et al. 2020;Emsenhuber et al. 2021a) where planet formation has been turned off to retain only the disc part.Paper I presented an updated version of the coupled gas and solids model that includes proper modelling of the disc dispersal stage.As the model was extensively described in Paper I, we only provide a brief overview of the physical processes included in the model.

Gas disc
The gas disc is modelled by an azimuthally averaged 1D radial structure.Its evolution is obtained by solving the advectiondiffusion equation (Lüst 1952;Lynden-Bell & Pringle 1974) where is the surface density and ν = αc s H the viscosity (parametrised using the α prescription of Shakura & Sunyaev 1973), with c s and H being the sound speed and scale height of the disc, and Σint and Σext the sink terms due to internal and external photoevaporation, respectively.To compute the vertical structure of the disc (and with this ρ G , H, and c s ), we proceed as in Paper I and use the vertically integrated approach of Nakamoto & Nakagawa (1994), including stellar irradiation (Hueso & Guillot 2005) from an evolving stellar luminosity computed from Baraffe et al. (2015).

Internal photoevaporation
Internal photoevaporation is modelled assuming X-ray-driven mass loss.This prescription requires one parameter that is not obtained from elsewhere in the model, the stellar X-ray luminosity L X .This luminosity is converted into a total mass-loss rate ṀX , and then into a profile Σint using fits to hydrodynamical simulations performed by Picogna et al. (2019), Ercolano et al. (2021), andPicogna et al. (2021), as described in Paper I.

External photoevaporation
The mass-loss rate due to external photoevaporation is obtained from the FRIED grid (Haworth et al. 2018).Interpolation in the grid requires the stellar mass M , the current disc mass M G , its outer radius r out , and the ambient far-UV (FUV) field F. All but the latter parameter can be computed consistently from the disc structure.The grid spans values of the ambient field F between 10 and 10 4 G 0 , where G 0 = 1.6 × 10 −3 erg s −1 cm −2 approximately represents the interstellar value (Habing 1968).The total mass-loss rate is converted into a profile Σext assuming mass is lost in the outermost 10 % where the gas disc is present at a given time (Paper I).
The FRIED grid, in its current state, presents two shortcomings that we adapt here: (1) the lack of data for ambient fluxes below 10 G 0 and (2) a floor evaporation rate of 10 −10 M yr −1 .Both items lead to a significant external photoevaporation rate under any circumstances, which makes it very difficult to disentangle the effects of external photoevaporation from the rest (including internal photoevaporation).To remedy these problems, we make one addition and one change to the FRIED prescription.The change is to take the lower boundary of the external photoevaporation rate down to 1 × 10 −15 M yr −1 , which represents a negligible mass-loss rate.For this, we remove the floor value of 10 −10 M yr −1 from the value returned from the interpolation in the grid and ensure that the resulting value is at least 1 × 10 −15 M yr −1 .The change we make is to extend the domain down to 1 G 0 to be able to study low-ambient-field cases.In the region below 10 G 0 , we perform a linear interpolation between the value returned from the grid at that boundary and a fixed value of 1 × 10 −15 M yr −1 at 1 G 0 .

Solids disc model
The solid component of the disc is modelled using the twopopulation model of Birnstiel et al. (2012).This approximates the full size distribution using only its two extremes: the smaller a 0 well coupled to the gas (which can be seen as dust) and the larger, rapidly drifting a 1 (which can be seen as pebbles).The smaller size a 0 = 0.1 µm is fixed while the larger size is constrained by various limits.The fragmentation limit is given by where f F = 0.37 is a factor fitted to hydrodynamical simulations of Birnstiel et al. (2010) for the typical size of pebbles, ρ s = 1.675 g cm −3 is the bulk density, and v frag the fragmentation velocity.In the drift limit, the large size is given by where v K is Keplerian velocity, Σ 0 the surface density of dust only, and ζ is an efficiency parameter of the drift (to account for the fact that drift is more limited in discs with features such gaps created by planets; e.g.Zormpas et al. 2022, while we only study smooth discs).
The surface density of solids Σ D is divided into the two components with The factor f m = 0.75 when growth is fragmentation-limited and f m = 0.97 when drift-limited.Gas drag onto both dust and pebbles is assumed to be in the Epstein regime.The radial velocity of solids is made of two components, coupling to the radial gas flow and headwind (Nakagawa et al. 1986; Paper I), where u G,red is the reduced radial gas velocity according to Gárate et al. (2020) and Paper I, St is the Stokes number, = ρ G /(ρ G + ρ D ), and ρ D is the midplane dust density.Here, we introduce a drift efficiency ζ to parametrise mechanisms that reduce the headwind-induced drift velocity of dust.
In particular, it is possible to use this approach to represent the effect that radial substructures have on the drift of solids without modelling them in full detail.The mass-averaged radial velocity is then given by ū = (1 − f m ) u 0 + f m u 1 .As in the gas disc, time evolution is provided by an advection-diffusion equation, where D G is the gas diffusion coefficient.
The terms Σphoto and Σrad are sink terms due to dust being entrained by photoevaporative winds (e.g.Facchini et al. 2016;Franz et al. 2020) and ejected due to radiation pressure, respectively; they are both described in Paper I. In contrast to Paper I, we allow planetesimals to form.This is parametrised using the term which follows the prescription of Lenz et al. (2019) as implemented and described in detail in Voelkel et al. (2020).Here ε is a parameter that specifies the conversion efficiency into planetesimals over a length scale of d = 5H (Dittrich et al. 2013), ūdr is the drift component of the mass-averaged radial velocity, and ṀD the relative mass flux of dust and pebbles through the gas.

Conversion into observed disc masses
For consistency with disc mass observations, millimetre (mm) emission from dust and pebbles is computed from the disc surface density and temperature profiles.The method is similar to that of Birnstiel et al. (2018) and will be discussed in more detail in Paper III.The calculation is performed for a wavelength of λ = 0.89 mm to reproduce ALMA observations.The flux is converted back into a mass using a simple prescription assuming T = 20 K and the corresponding opacity.
For comparison, we also provide the unbiased disc masses of gas and solids.To be presented alongside disc gas masses, solid masses are multiplied by a factor 100, which is typically used as a gas-to-dust ratio in this context.

Model parameters and initial conditions
The evolution model requires several initial conditions and parameters for evolution.These are: the mass of the central star M ; the mass of the gas disc M G ; the initial dust-to-gas ratio f D/G ; the power-law index for the initial profile β; the inner edge of the disc r in ; the characteristic radius of the disc r 1 ; the turbulent viscosity parameter α; the planetesimals formation efficiency ε; the fragmentation velocity v frag ; the efficiency of drift ζ; the stellar X-ray luminosity L X ; and the ambient UV field strength F.
The initial surface density profile of the gas disc is set as (Veras & Armitage 2004) where Σ ini is the surface density at r 0 = 5.2 au, the reference distance.The conversion between that and the total mass is obtained with The initial solid profile of the disc Σ D (t = 0) is set by multiplying the initial gas profile Σ G (t = 0) by the dust-to-gas ratio f D/G .
In the remainder of this work, we do not provide all parameters as such.For instance, the inner edge is parametrised by its period P in , which we convert into distance by means of Kepler's third law.Also, we generally set the initial disc mass by its solid content M D .The ratio between the initial solids and gas masses is readily given by the dust-to-gas ratio, such that f D/G = M D /M G .

Simulation list
To generate the list of the simulations to be performed, we selected the Latin hypercube sampling (LHS) method (e.g.McKay et al. 1979).By dividing each dimension into n intervals and then selecting one random sample from each interval, LHS ensures that the entire range of possible values for each parameter is sampled with a uniform probability.Additional criteria are required to avoid correlation between selected values of different parameters (to disentangle their effects) and to ensure that the entire space is well sampled (to avoid locations with no results).To build the grid, we use the pyDOE2 Python package with the minmax setting.Each generated grid contains values in the [0,1] range with uniform probability.
These values have to be mapped into the range to be studied.For our main grids, we outline these in Table 1.The selection was made to encompass the needs of this and future works, as well as the limitations of the model.For instance, the stellar mass M is taken in steps of 0.1 M to lie on the stellar evolution tracks of Baraffe et al. (2015), and the limits of the ambient UV field strength F match those of the FRIED grid (Haworth et al. 2018) with the extrapolation for low field values from Paper I. The gas mass of the disc is given in terms of the stellar mass to roughly follow the scaling of Burn et al. (2021).The dust-to-gas ratio was selected to span the possible stellar metallicities, with a reference stellar metallicity f D/G, = 0.0149 (Lodders 2003).The range power-law index was selected to cover the possible values of Andrews et al. (2009).The lower boundary of the period at the inner edge P in corresponds to 0.7 d, which is nearly the maximum value of the stellar radii in the models of Baraffe et al. (2015).The fragmentation velocity was chosen to encompass the previously assumed value of ∼10 m s −1 (e.g.Dr ążkowska & Alibert 2017, and references therein), and more current values of ∼1 m s −1 , as ice was not found to be more sticky than silicates in recent experiments (Gundlach et al. 2018;Musiolik & Wurm 2019;Steinpilz et al. 2019).The range of planetesimal-formation efficiencies was selected to be able to study low efficiencies where only a small fraction of the mass of solids is converted into planetesimals and to cover the case ε = 0.05, which forms a sufficient amount of planetesimals.

Machine learning
Surrogate models of disc evolution are obtained by means of a neural network.These neural networks are trained, validated, and tested using the scikit-learn Python package (Pedregosa et al. 2011).scikit-learn uses cross-validation to train and validate the neural network with five passes.This means that the combined training and validation set is divided into five equalsized batches, and five successive training and validation steps are performed, each using four of the five batches for training and one batch for validation.The neural networks are fitted using either the L-BFGS-B algorithm (Byrd et al. 1995), which is part of the SciPy package (Virtanen et al. 2020) or the ADAM method (Kingma & Ba 2014).

Observational constraints
To compute disc populations that are comparable to observations, we must first describe the constraints on their initial properties and their outcomes.These are then used to set the initial conditions and the comparison point for the outcomes.

Stellar mass
The stellar initial mass function (IMF) has been determined (e.g.Chabrier 2003), and so it could in principle be used to reproduce the stellar population.However, the stellar mass functions for different star-forming regions deviate from the IMF.In the case of Taurus, Luhman (2000) found a peak around 0.6 to 1 M , while for the Orion Nebula Cluster (ONC), Da Rio et al. (2012) found that the best log-normal fit has a mean at log 10 (M /M ) = −0.45(corresponding to M = 0.35 M ), using the stellar evolution model of Baraffe et al. (1998).Corroborating this, the sample of Flaischlen et al. (2021), which is based on that of Manara et al. (2012), has a stellar mass distribution peaking around 0.4 M : a simple log-normal fit to that data gives log 10 (µ/M ) = −0.481and a narrower standard deviation of log 10 (σ) = 0.2383.As the main aim of this work is to compare our model with disc lifetimes, their mass, and the stellar accretion rate of nearby star-forming regions, we chose to follow the stellar mass function of Da Rio et al. ( 2012), with a mean of log 10 (µ/M ) = −0.45 and standard deviation log 10 (σ) = 0.44.This should offer a distribution that is representative of both nearby clusters in general and of stars for which observations of disc masses and stellar accretion rates are available.
Our model uses the stellar evolution tracks of Baraffe et al. (2015) to obtain the luminosity for disc irradiation.These are only defined for mass increments of 0.1 M from 0.1 to 1.4 M .To properly track stellar luminosities, we restricted ourselves to stellar masses that match these values.

Initial dust mass
Initial dust masses can be obtained from works targeting the youngest stars known to date, such as Tychoniec et al. (2018), Williams et al. (2019), or Tobin et al. (2020).Emsenhuber et al. (2021b) fitted the masses of the Class 0 discs of Tychoniec et al. (2018), which gave log 10 (µ/M ) = −3.49and σ = 0.35 dex, taking out the conversion from dust mass to gas mass using the  Preibisch et al. (2005) 1.44 ± 0.10 30.37 ± 0.06 0.65 Güdel et al. (2007) 1.52 ± 0.12 30.31 ± 0.06 0.54 Parameters a and b are those of the fit log L X /erg s −1 = a × log (M /M ) + b.The 'Sca.' column provides the scatter of the residuals from the fit.
standard factor of 100 that was used there.These disc masses were used for a population of stars with masses of 1 M while the populations around lower-mass stars of Burn et al. (2021) scaled the disc masses proportionally to the stellar masses.However, a complication arises from the fact that the mass of the central body is not known for the objects observed by Tychoniec et al. (2018) and Tobin et al. (2020).To properly convert the absolute masses into disc-to-star mass ratios, as we do in this work, we must assume a reference stellar mass M ref .The method of Emsenhuber et al. (2021b) and Burn et al. (2021) was equivalent to setting M ref = 1 M .We use this as our default conversion factor, although consistency with the stellar mass distribution discussed in the previous section would call for a lower value of M ref .We explore different values of this factor later in this work.

Sizes
Protoplanetary disc sizes have been found to be correlated with their mass.Andrews et al. (2010) found that discs in the Ophiuchus star-forming region have M D ∝ r (1.6±0.3) 1 and β = 0.9 ± 0.2.More recent studies, such as those of Tripathi et al. (2017) plus a residual scatter of the order of 0.1 dex.

Dust-to-gas ratio
The ratio between the initial masses of the gas and dust discs is given by f D/G .We select this parameter as in Emsenhuber et al. (2021b), that is, we assume it is the same as the stellar metallicity (Gáspár et al. 2016).Thus, we can use the relation et al. 2001), where f D/G, = 0.0149 is the primordial solar value (Lodders 2003).The distribution of metallicity is chosen to be that of the CORALIE RV search sample (Santos et al. 2005), which was modelled as a Gaussian with a mean of −0.02 and a standard deviation of 0.22.
To avoid extreme values, we restrict the parameter to within −0.6 < [Fe/H] < 0.5.

Stellar X-ray luminosity
A couple of surveys have been performed to determine the Xray luminosities of young stars, the relevant results of which are provided in Table 2.One is the Chandra Orion Ultradeep Project (COUP; Getman et al. 2005;Preibisch et al. 2005), which covers stellar masses M between 0.5 and 0.9 M .The survey found a stellar-mass dependency of L X ∝ M (1.44±0.10) .Another survey, the XMM-Newton Extended Survey of Taurus (XEST; Güdel et al. 2007), found that L X varies with stellar mass as L X ∝ M (1.54±0.12), which we used to correct for the stellar-mass effect and recompute the inherent scatter.The two surveys have similar stellar mass dependence, meaning that using one or the other to set the stellar X-ray luminosities should not affect the outcomes in any significant manner.For this work, we compute L X using a log-normal distribution with the parameters selected following XEST (Güdel et al. 2007), as the stellar mass dependence is consistent with the prescription used to compute the X-ray photoevaporation profiles in Picogna et al. (2021).

Ambient FUV field strength
The external photoevaporation prescription of Haworth et al. (2018) requires the stellar mass, disc mass, outer radius, and ambient FUV field strength.The first three can be readily obtained consistently from the simulation, but the latter, F, needs to be specified.
Most stars are formed in stellar clusters (e.g.Lada & Lada 2003), which result in high stellar densities.To retrieve the ambient FUV relevant during the lifetime of protoplanetary discs, we use the simulation of Adams et al. (2006).The authors determined that F is well described by a log-normal distribution with a median close to 10 3.25 G 0 , where G 0 = 1.6 × 10 −3 erg s −1 cm −2 is nearly the interstellar FUV field (Habing 1968).

Inner edge of the gas disc
The location of the inner edge of the gas disc is most relevant for the location of the close-in planets (such as hot Jupiters).As we are mostly interested in warm giants further away than the inner edge, this parameter is of less importance in this work.We chose this parameter in the same way as in Emsenhuber et al. (2021b), that is, by assuming that the disc is truncated at the corotation radius of the star.For the distribution of stellar rotation periods, we follow the results of Venuti et al. (2017).This gives a log-normal distribution with a median period of log 10 (µ/d) = 0.67617866 and deviation σ = 0.305 673 3 dex.
For comparison, the distribution of initial rotation periods used by Johnstone et al. (2021) has a median of log 10 (µ/d) = 0.5181 and a standard deviation of σ = 0.3236 dex.The median rotation period here is smaller here (3.3 d) than the 4.7 d value of Venuti et al. (2017) but not by a large amount, while the deviations are similar.The exact choice should therefore not affect the results significantly.

Full model
To generate the training, validation, and testing data for the surrogate models, we generated two sets of simulations.The first set contains 100 000 models that are used for the combined training and validation steps, while the second set contains 20 000 models and is used for the testing step.The values of the first set are provided in Table 3 while the values of the second set are provided in Table 4.Both tables are available at the CDS and have the same format; they contain the following columns: columns 1 to 12 are the initial conditions in the same order and units that are given in Table 1.Column 13 gives the disc lifetime according to when the mass becomes lower than 10 −6 M or when the surface density is lower than 1 × 10 −3 g cm −2 inside 100 au (or 30 au for M = 0.1 M or 0.2 M ) and 1 × 10 −2 g cm −2 outside that (this second criterion on the surface density is to avoid excessively long-lived discs when photoevaporation rates, particularly external ones, are low).Column 14 gives the lifetime using the minimum value of the criterion of column 13 and the observability criterion in the near-infrared (NIR) from Kimura et al. (2016).Columns 15-19 give the following outcomes at 1 × 10 5 yr: stellar accretion rate log 10 Ṁ /M yr −1 , the true gas mass log 10 (M G /M ), the true solids mass log 10 (M D /M ), the observed mass (Sect.2.3) log 10 (M obs /M ) and the radius encompassing 68 % of the flux log 10 (r 68 /au).Columns 20-24 repeat the same information, but at 2 × 10 6 yr.
Two epochs (1 × 10 5 yr and 2 × 10 6 yr) were selected to be compatible with the observations we are comparing to.The first epoch is for comparison with early discs, such as their initial masses.Its selection is a trade-off between two items: on the one hand, we would like to have the data as early as possible, while on the other hand, we need to wait until the initial dust growth has taken place.From the analysis of individual discs, we found that 1 × 10 5 yr represents a good compromise in that sense.The second epoch is for comparison with the star-forming regions of Lupus and Cha I.As the stars in these regions are between 1 and 3 Myr old, we take the results at 2 Myr, as in Manara et al. (2019).
Our results indicate that the two criteria for disc dispersal produce nearly identical results.In only about 10 % of the cases, the NIR criterion predicts a lower disc lifetime than the criterion based on the mass, and the difference remains small when this occurs (we do not check for the reverse, as calculations stop when the mass criterion is reached).These results are consistent with the findings of Kunitomo et al. (2021).As a consequence, hereafter we only report the disc lifetimes based on the NIR criterion of Kimura et al. (2016).Also, we stop the calculation at 100 Myr in any case.This affects some long-lived discs with minimal photoevaporation and accretion.In such cases, the lifetime based on the mass criterion is not reported while that based on the NIR emission is.

Performance of the surrogate model
We asses the performance of the surrogate models in terms of the best regression (obtained using ordinary least squares), the Pearson correlation coefficient R 2 , and the RMS of the differences between the predicted and target lifetimes (the square root of the mean square error).These were computed on the testing set (Table 4) that the surrogate model has not seen before.The hyper-parameters and results for all surrogate models that are part of this work are presented in Table 5.For three of them, we also show correlation plots in Fig. 1.In all cases, the fitting procedure was performed on the logarithm (base 10) of the values, and so all the reported performances are given in these units.
Concerning the different surrogate models, the ones for the disc lifetimes and for the stellar accretion rates provide the best performance.The ones that are based on the dust disc, namely masses and radii, show a lower performance, especially at 2 Myr.We note that these values are for each single prediction; they represent the level of additional uncertainty for the parameter studies (Sect.4.3 and Appendix A) while for the population studies (Sect.5) these errors can average out and result in an even better global accuracy.
The neural networks predicting the disc masses, their radii, and stellar accretion rates were fitted only on the discs that had not vanished at the time.This means that they are supported by a lower number of points than the ones predicting the lifetimes.This also implies that these surrogate models are only constrained in the region of the parameter space where lifetimes are larger than the time of the analysis.Thus, in the remainder of this work we only provide disc masses and stellar accretion rates for discs that have not yet dispersed.

Effects of photoevaporation
We began our investigations using the surrogate model, performing a parameter study of the effects of the photoevaporation prescriptions on disc lifetimes.For this purpose, we generated two maps, one for internal photoevaporation and one for external photoevaporation, which vary the stellar mass and the controlling parameter of each photoevaporation prescription.In each case, the value of the parameter controlling the other photoevaporation prescription was set at the minimum of the studied range in order to avoid cross effects.We assumed typical values of the remaining parameters: disc-to-star gas mass ratio M G /M = 0.1, dust-to-gas ratio f D/G = 0.0149, power-law index β = 0.9, period at the inner edge P in = 10 d, characteristic radius r 1 computed according to Eq. ( 9), viscosity parameter α = 1 × 10 −3 , fragmentation velocity v frag = 2 m s −1 , planetesimal formation efficiency ε = 1 × 10 −3 , and drift efficiency ζ = 1.

Internal photoevaporation
The resulting map for internal photoevaporation is provided in Fig. 2.Here we observe that the surrogate model predicts several sharp transitions of disc lifetime with stellar mass.The most evident are those between 0.2 and 0.3 M and between 0.8 and 0.9 M where disc lifetime increases.There are other transitions between 0.4 and 0.5 M and between 0.6 and 0.7 M where disc lifetime decreases, but only for large X-ray luminosities (L X > 10 30 erg s −1 ).These transitions match the switch from one photoevaporative profile to another, which are marked by the dashed white lines.This indicates that the profile of surfacedensity loss has a strong effect on disc lifetime and not only the total mass-loss rate, which gradually changes between each stellar mass.Also, the further out the location of the peak of internal photoevaporation (which is for the profiles of 0.3 M and 1.0 M stars; see top panel of Figure 7 of Picogna et al. 2021), the longer the disc lifetimes are in general.We find that this effect is due to a larger inner region where material is not evaporated at all and can only be dispersed by viscous accretion.In this case, the observed disc lifetime is set by the dispersal timescale of the inner disc, which is given by the viscous timescale at the outer radius of the inner disc.
As discussed in Sect.3.5, the X-ray luminosity is correlated with stellar mass.To highlight this, we show in addition the median stellar X-ray luminosity as a function of stellar mass from Güdel et al. (2007) with the green dashed curve.To determine the expected relationship between disc lifetime and stellar mass, one needs to follow this curve rather than a horizontal line on the plot.We see that internal photoevaporation leads to a limited change in disc lifetime with stellar mass.This is because more massive stars lead to stronger mass-loss rates (owing to a corresponding increase in stellar X-ray luminosity), which compensates for the increase in disc mass (as we assume disc mass to be proportional to stellar mass).This is shown by the black line that traces disc lifetimes of 3 Myr (a typical value in obser-

External photoevaporation
The resulting map for external photoevaporation is shown in Fig. 3. Unlike internal photoevaporation, the prescription for external photoevaporation provides for gradual changes of lifetime with stellar mass.However, these changes lead to a larger dependency of disc lifetime on stellar mass than what is expected from internal photoevaporation.This is illustrated by the black line, which tracks a 3 Myr lifetime, as in the map for internal photoevaporation.Its position in terms of ambient UV field strength varies across the entire parameter range studied here, from less that 2 G 0 for M = 0.1 M to about 4 × 10 3 G 0 for M = 1 M ; it becomes independent of stellar mass for M > 1 M and fluxes above ∼10 3 G 0 .
While the trend of reduced disc lifetimes in regions with strong ambient UV fields (e.g.Michel et al. 2021) is reproduced, the general behaviour of correlated disc lifetimes with stellar mass for a given ambient UV field strength is problematic for several reasons.First, this general behaviour is inconsistent with observations that suggest disc lifetimes are independent of, or slightly decreasing with, increasing stellar mass (Carpenter et al. 2006;Kennedy & Kenyon 2009;Bayo et al. 2012;Ribas et al. 2015).There are several possibilities to remedy this, although they are unlikely.To obtain a behaviour similar to observations, the mass loss would need to be correlated with stellar mass (Komaki et al. 2021), which in turn would require that the ambient field be correlated with stellar mass.However, the ambient field is usually dominated by the few most massive stars (Adams et al. 2006), which means that it depends more on the cluster as a whole than on the mass of the star in question.Another avenue is that disc masses scale to a lesser extent with stellar mass than assumed here.However, this would not yield the expected behaviour of stellar accretion rates with stellar masses (e.g.Hartmann et al. 2016;Alcalá et al. 2017;Flaischlen et al. 2021).We find that the FRIED grid prescription that we use in this work produces incompatible results that show at most a dependence of the disc lifetime on stellar mass.The second concern is that further lifetime analyses will be strongly affected by the selection of the stellar masses, in contrast to internal photoevaporation where this dependency is weaker.

Parameter sensitivity
The sensitivity of disc lifetimes, disc masses, and stellar accretion rates at 2 Myr is studied in detail in Appendix A. These results can be summarised as follows: all the outcomes are insensitive to the power-law index β and the inner edge of the gas disc r in .The characteristic radius r 1 and viscosity α control the viscous timescale of the disc, and therefore the stellar accretion rate.Disc lifetimes and observed dust masses are more strongly affected by the viscosity α than by the characteristic radius r 1 .
The twopop model parameters only affect the observed dust masses.Observed dust masses are less affected by the dust-togas ratio f D/G than by the initial mass of the gas disc M G , except for discs close to dispersal.The fragmentation velocity v frag and the drift efficiency strongly affect the observed dust masses, but only for v frag 200 cm s −1 , while the planetesimal formation efficiency ε only has a limited effect for values close to the maximum we study, namely of 0.1.
These results narrow down the parameter space that we explore in the remainder of this work.First, we keep the values of 2. Map of disc lifetimes as a function of stellar mass and X-ray luminosity, which is the main driver of internal photoevaporation, according to the surrogate model described in Sect.4.2.External photoevaporation was set to its minimum value (F = 1 G 0 ).Other parameters were selected as M G /M = 0.1, f D/G = 0.0149, β = 0.9, P in = 10 d, r 1 according to Eq. ( 9), α = 1 × 10 −3 , v frag = 2 m s −1 , ε = 1 × 10 −3 , and ζ = 1.The green dashed line represents the dependency of L X on M from Güdel et al. (2007) and the solid black line shows the location of a 3 Myr lifetime (a typical value).The results are discussed in Sect.4.3.1.the power-law index β and the inner edge of the gas disc r in as described in Sect.3, because they are of negligible importance.We then only use the initial mass of the gas disc M G to control disc masses, not the dust-to-gas ratio f D/G ; as the latter is in most cases of lower importance and well constrained by observations.Also, we keep the planetesimal formation efficiency to the minimum value of ε = 1 × 10 −3 , the fragmentation velocity to v frag = 200 cm s −1 , and the drift efficiency to ζ = 1.

Disc populations
We now compare synthetic disc populations with observations.For this, we proceed as follows: we draw 10 000 random discs whose initial conditions follow given distributions.The outcomes of each disc are obtained by means of the different surrogate models.For the analysis, we first compare the cumulative distribution of disc lifetimes so that it can be compared to the fraction of stars that have a protoplanetary disc for stellar clusters with different ages, as in Haisch et al. (2001).The second analysis is to compare observed disc masses and stellar accretion rates with the data of Manara et al. (2019).Here, we use the data at 2 Myr.Further, we only use discs whose lifetime, as determined by the surrogate model from the previous analysis, is larger than the time of analysis in order to avoid being in the region where the surrogate model is not supported by any underlying data.3. Map of disc lifetime as functions of stellar mass and ambient UV field strength, this latter being the main driver of external photoevaporation according to the surrogate model described in Sect.4.2.Internal photoevaporation is set to its minimum value (L X = 1 × 10 28 erg s −1 ).Other parameters were selected as M G /M = 0.1, f D/G = 0.0149, β = 0.9, P in = 10 d, r 1 according to Eq. ( 9  M D / f D/G β 0.9 log 10 (P in /d) N(0.67617866, 0.3056733 2 ) r 1 /au 70(M D /100 M ⊕ ) 0.25 × 10 N(0,0.1 2 ) log 10 (α) U(−3.5, −3.0)

Canonical
To determine if all the processes that are predicted from theory are able to reproduce disc observations, we compute a population of discs whose properties are as close as possible to observations from early discs.The only parameter that has some freedom is α.Here, we selected to draw log 10 (α) with a uniform probability of between −3.5 and −3.This was decided as a compromise between disc lifetime and stellar accretion rate, as we discuss below.For the other parameters, their distributions were selected as described in the discussion of Sect.3; for convenience, these are summarised in Table 6.
The resulting distribution of disc lifetimes is shown with the blue curve in Fig. 4. It becomes immediately apparent that the synthetic lifetimes are too short overall in comparison with observed discs.The median lifetime of the synthetic discs is 0.42 Myr.Disc lifetime depends on the assumed distribution of α, which we chose such that it results in the largest stellar accretion rates at 2 Myr.A histogram of stellar accretion rate versus observed dust mass for the observed discs in the Lupus and Chamaeleon I star-forming regions is shown in the top-left panel of Fig. 5.Only the synthetic discs that live beyond 2 Myr contribute to this histogram.The few remaining discs have low masses, as the discs are close to being dissipated.Using larger values of α, for instance between roughly 10 −3 and 10 −2 as was proposed by Mulders et al. (2017), would have resulted in even shorter disc lifetimes.This means that there would have been no discs that would live long enough to produce stellar accretion at 2 Myr.Conversely, selecting a distribution of α with even lower values would allow disc lifetimes to be matched by observations, but this would result in even lower stellar accretion rates, which would be in tension with the results of Dullemond et al. (2018), who concluded that α ≥ 1 × 10 −4 from the sizes of disc substructures.
The discrepancy between our modelled lifetimes and observations arises from the strong mass-loss rates predicted for internal and external photoevaporation, as we discuss in Sect.6.As disc lifetime depends on stellar mass (especially for external photoevaporation; Sect.4.3), this analysis is affected by the assumed stellar mass distribution.Had we selected larger stel-lar masses, the lifetimes would better match observations.However, selecting stellar masses of around 1 M , which would lead to a fairly good match including both photoevaporation prescriptions, is not representative of the star-forming regions that we are comparing to (Sect.3.1).

Effect of photoevaporation on disc populations
The mismatch in disc lifetimes and disc masses that we obtained is in contrast with other studies, such as that of Mulders et al. (2017), who did not include photoevaporation, Kunitomo et al. (2020), who included only internal photoevaporation, and Weder et al. (in prep.), who used an EUV-only internal photoevaporation prescription with much lower mass-loss rates.Further, as already discussed, photoevaporation leads to considerable shortening of disc lifetimes (Sect.4.3).Therefore, here we investigate how the photoevaporation prescriptions affect the evolution of disc populations and determine whether they are responsible for the mismatch.To this end, we created two additional populations, each time neglecting one of the photoevaporation process by taken the corresponding controlling parameter to the minimum value.For the population with internal photoevaporation only, the ambient flux has been set to F = 1G 0 , which corresponds to Ṁ = 1 × 10 −15 M yr −1 (Sect.2.1.2).The results are shown as 'Int.only' in Figs. 4 and 5.For the population with external photoevaporation only, we set L X = 10 28 erg s −1 , the results of which are shown as 'Ext.only' in the same figures.
In each population, the disc lifetimes strongly increase compared to the canonical case.However, the distributions are different: internal photoevaporation leads to a relatively narrow distribution around 2 to 3 Myr, while with external photoevaporation disc lifetimes are more spread out, including very short-lived discs.This leads to more discs being shown in the accretion versus mass diagram.
The population with only external photoevaporation has low disc masses combined with a narrow range of stellar accretion between 10 −10 and 10 −9 M yr −1 .This is because the mass loss occurs in the outer disc, which reduces its size, limiting the area from which the dust emits (as dust is also lost where there is no longer gas present).At the same time, the mass-accretion rate is weakly affected by the mass loss in the outer disc; again because external photoevaporation affects only the outer disc.
Conversely, using internal photoevaporation leads to a larger spread in stellar accretion rate.As internal photoevaporation removes material relatively close to the star, it competes with stellar accretion to some extent (Somigliana et al. 2020).This, coupled with the spread of X-ray luminosities, leads to a spread in accretion rate (Owen et al. 2011).Thus, internal photoevaporation is needed to reproduce the observed spread in stellar accretion rate (Manara et al. 2016b;Testi et al. 2022).We further note that neither population is able to reproduce the discs with an accretion rate larger than 10 −8 M yr −1 .
In addition, internal photoevaporation leaves discs with an inner cavity that have low accretion rates but where the outer disc is out of reach of internal photoevaporation and therefore takes a long time to dissipate; these are what Owen et al. (2011) refer to as relic discs.It is possible to avoid this situation to a large extent by having more compact discs initially, which leave nearly all of their mass within reach of internal photoevaporation.
We conclude that previous studies managed to reproduce disc lifetimes because they used only one photoevaporation mechanism as the main loss mechanism.However, when both are accounted for, the combined mass-loss rate is so large that discs are very short lived.We discuss the implications of this and possible remedies in Sect.6.

Towards a best match
Despite all the differences between models and observations highlighted so far, we now try to find a set of initial conditions that is able to better match disc evolution characteristics.To this end, we investigate how the initial conditions can be modified from our canonical values provided in Table 6.
The initial disc-to-star mass ratios were taken from Tychoniec et al. (2018) and Tobin et al. (2020) assuming they were measured on stars of M ref = 1 M .However, this assumption is inconsistent with the stellar mass distribution we selected, which has a median value of M = 0.35 M (Sect.3.1).Were we to select a lower reference stellar mass, such as M ref = 0.33 M , this would increase the disc masses.At the same time, selecting a 6.Kernel density estimate for two populations from Figs. 4 and 5, both for their initial conditions (dashed lines) and retrieved disc masses at 100 kyr (solid lines).The initial mass distribution of the best-match population with M ref = 0.33 M (in red) is consistent with the fit by Emsenhuber et al. (2021b) to the data of Tychoniec et al. (2018), which is shown as the dashed black line, while the retrieved masses at 100 kyr are compatible with the non-multiple discs of Tobin et al. (2020).
reference stellar mass that is similar to the median value from our initial conditions results in an agreement between the disc masses in observations and our initial conditions, as shown with the dashed black and red lines in Fig. 6, respectively.In addition, to account for observational biases, we also want to check the observed dust masses after a short evolution time of 1 × 10 5 yr, which we provide with the solid lines in Fig. 6.The results show that the discs in the nominal populations have low masses compared to the non-multiple discs measured by Tobin et al. (2020), while a population with more massive initial discs has slightly larger masses.Thus, the larger initial disc masses lead to a reasonable match with observations as a whole.
We point out that the disc-to-star mass ratio in the population with M ref = 0.33 M is about twice (2.1 times) that of Emsenhuber et al. (2021b) and Burn et al. (2021), rather than a factor three as one might assume from the change of M ref from 1 M to 0.33 M .This is due to an inconsistency in the selection of the disc masses in the previous work: there, the gas masses were taken as a Monte Carlo variable that were converted from dust observations using a dust-to-gas ratio of 1 %.However, the initial mass of solids in the model was recomputed from the gas mass using the same distribution as in this work (Sect.3.4), which has a median value of 1.42 % ( f D/G, = 1.49% with [Fe/H] = −0.02).In contrast, we use the solid disc mass as a Monte Carlo variable here.
Another possibility is that early protoplanetary discs are not as extended as what is suggested by the findings of Tobin et al. (2020).More compact discs are less susceptible to external pho-toevaporation, as there is less surface exposed to ambient radiation and they are more tightly bound to the star.At the same time, a more compact disc means that more material is concentrated in the region where internal photoevaporation is most efficient, which allows the stellar accretion rate to remain larger.Magnetohydrodynamics models of protoplanetary disc formation by Hennebelle et al. (2016), Lebreuilly et al. (2021), and Lee et al. (2021) could favour such a possibility.
Finally, nearby star-forming regions have low masses, which results in a low ambient UV field strength F; for instance, the value in Lupus is F ≈ 4G 0 (Cleeves et al. 2016).Our nominal distribution of ambient UV fields overestimates mass losses due to external photoevaporation.Therefore, we set F = 1G 0 , which results in negligible mass losses due to external photoevaporation ( Ṁ = 1 × 10 −15 M yr −1 ) and only internal photoevaporation remains.Using internal photoevaporation confers the further advantage that a wider range of stellar accretion rates is obtained.
Below, we investigate whether more massive and compact discs are able to improve the match with observations.The effects of the parameters mentioned above are discussed in Appendix B. From these, we find that the following modifications to the initial conditions given in Table 6 are best able to reproduce disc lifetimes and the accretion rate-mass relationship: -A decrease in the reference stellar mass to M ref = 0.33 M , which corresponds to an increase in the disc mass by a factor of three compared to the nominal population; r 1 = 2/3 × 70 au (M D /100 M ⊕ ) 0.25 , a factor 2/3 compared to Eq. ( 9); and only using internal photoevaporation (F = 1G 0 ).The population using these distributions is shown as 'Best match' in Figs. 4 and 5.
While the overall stellar mass distribution we chose is representative of the observed stars, our visual optimisation approach to reproducing the observed disc masses and accretion rates is independent of the exact stellar mass dependency of the observables.We will improve on this with a Bayesian framework in Paper III to optimise the initial parameters when reproducing a set of observations in four-dimensional space made up of stellar mass, disc mass, disc radius, and accretion rate.
Compared to the population with only internal photoevaporation (the other population that is closest in terms of initial conditions), we can see several differences.First, the larger disc masses result in an increased median lifetime.About 2.2 % of the discs have a lifetime of greater than 10 Myr, representing a 100-fold increase.Then, the combination of larger initial disc mass and smaller extent results in a certain number of cases with a stellar accretion rate of higher than 10 −8 M yr −1 , which was not previously seen.The smaller extent of the disc also causes less discs to be relics (towards the bottom of Fig. 5).Here, the highest concentration of discs is found near the best-fit value of Manara et al. (2016b, the pink dashed line in Fig. 5) with a similar number on either side.We are still failing to reproduce the discs with the large accretion rates.However, part of these discs with large accretion rates could be due to binaries (Zagaria et al. 2022), which we do not model.To corroborate this, the largest stellar accretion rates are biased to larger-mass stars (the mean stellar mass for systems with a stellar accretion rate higher than 3 × 10 −9 M yr −1 is 0.74 M versus 0.43 M for the general population), which at the same time are more likely to be in binary systems (e.g.Duchêne & Kraus 2013).The mismatch should therefore not be the source of significant concern.We find that this combination of parameters is able to reproduce the disc mass-accretion rate relationship and provide a reasonable match to most observations.

Discussion
It is difficult to reconcile the results of our synthetic disc populations with observations.We find that it is particularly hard to obtain discs with characteristic lifetimes of 2 to 3 Myr according to Mamajek (2009) or Fedele et al. (2010), and even less lifetimes of 5 to 10 Myr following Pfalzner et al. (2022).Here, we assumed that the initial mass is constrained by the dust-mass measurements of Tychoniec et al. (2018, M D /M ∼ 10 −3 ) , and dust-to-gas ratios similar to the solar initial abundance, namely f D/G = 1.49% (Lodders 2003), combined with the predictions of internal and external photoevaporation models.To illustrate this issue, in Fig. 7 we provide the evolution of the gas mass still present in the disc and removed by the processes modelled here.Here we choose typical values for the initial conditions, with a star of mass M = 0.5 M , a gas disc with a mass of M G /M = 0.1, a power-law index of β = 0.9, an inner radius of r in = 0.1 au, and a characteristic radius of r 1 = 87.8au, which was computed from Eq. ( 9) assuming a dust-to-gas ratio of f D/G = 0.0149.The ambient UV field strength F = 10 G 0 was set at the lower boundary of the computed grid while the X-ray luminosity L X = 7.02 × 10 29 erg s −1 follows the best fit of the XEST survey (Güdel et al. 2007, Table 2) for the given stellar mass.
Figure 7 shows that photoevaporation (both internal and external) is responsible for the loss of nearly all the gas; only 6 % of the gas is accreted onto the star.The final disc lifetime of 1.38 Myr is short compared to the characteristic lifetimes from observations.We therefore have a problem with the mass bud-get, which could be resolved by (i) larger initial disc masses, (ii) mass replenishment after disc formation, or (iii) lower mass-loss rates by photoevaporation.
Stellar accretion already plays a minor role in the mass budget and cannot be reduced further, in order to remain consistent with observed stellar accretion rates.However, radial mass transport could be the result of magnetically driven disc winds (Suzuki & Inutsuka 2009;Suzuki et al. 2010) rather than pure viscous dissipation as we assumed here.This would add another mass-loss channel, which would further exacerbate the problem (though it could shield stellar radiation, as discussed below).
We experiment with larger initial disc masses in this work.However, even with gas masses on the order of 10 % of the stellar mass, disc lifetimes are not sufficiently long.Therefore, increasing the masses even more would be required, but this increase would bring another series of problems.For one, such large discs are likely gravitationally unstable and produce spiral density waves.Second, such large discs would lead to strong gas-driven migration, which hinders planet formation (Nayakshin et al. 2022).Also, discs around 10 Myr-old stars HD163296 and TW Hya have at least 10 % of the stellar mass (Powell et al. 2019) and it is unclear what their initial mass would have been for them to remain so large at their age.We therefore do not find that massive initial discs would be able to solve the conundrum.
Discs do not form instantly.Rather, they grow from gas falling from the envelope.This allows for disc masses to increase during the early stages of disc evolution, which is at odds with the models described here.Modelling the infall stage would allow us to have longer-lived discs without them being very massive early on.Accretion can persist for several million years (Throop & Bally 2008), providing replenishment even at late times.The complex morphology of the gas disc around RU Lup (a ∼0.5 Myr old star) could be an outcome of this process (Huang et al. 2020).The long-lived discs could also be secondgeneration discs, formed following accretion from the molecular cloud (Kuffmeier et al. 2020) or the disruption of a planet (Nayakshin et al. 2020).While none of these individual processes are sufficient to explain the presence of massive discs, they should still be explored if they can explain certain characteristics of the overall disc population.
Finally, the photoevaporation rates used here could be overestimated.Lower mass-loss rates would increase the disc lifetimes and masses at later times.An argument in favour of this hypothesis is that young discs can be shielded from the radiation of both their host and nearby stars.The launching of magnetically driven disc winds occurs inside the location where internal photoevaporation is effective and would thus prevent EUV and X-ray photoevaporation during the early stages of disc evolution (Pascucci et al. 2022).Similarly, early discs are likely embedded, preventing the radiation of nearby stars from reaching the disc.Both effects would reduce the photoevaporation rates during the early stage of disc evolution compared to what we model here.

Summary and conclusion
In this work, we investigate whether protoplanetary disc observations can be reconciled with theoretical predictions of processes such as viscous accretion and photoevaporation (both internal and external).We first compute two sets of simulations that we use to fit neural networks (Fig. 1).With these neural networks, we can perform parameter studies and compute the outcomes of synthetic disc populations with limited computational resources.
We first compare how internal and external photoevaporation affect disc lifetime as a function of stellar mass.We find that because of a direct link between stellar mass and X-ray luminosity (e.g.Preibisch et al. 2005;Güdel et al. 2007), which means mass-loss rate due to internal photoevaporation, discs around more massive stars are not significantly longer-lived than those around low-mass stars (Fig. 2).Conversely, external photoevaporation leads to a strong positive correlation between disc lifetime and stellar mass (Fig. 3), because gas is more bound for higher-mass stars.This positive correlation is at odds with observations that find that disc lifetimes are either independent of or anticorrelated with stellar mass (Carpenter et al. 2006;Kennedy & Kenyon 2009;Bayo et al. 2012;Ribas et al. 2015;Pfalzner et al. 2022) and should be investigated in the future.
Turning to protoplanetary disc populations, we find that accounting for both internal and external photoevaporation according to theoretical predictions leads to disc lifetimes (Fig. 4) that are much too short.Discs whose initial mass is 10 % of the stellar mass are dispersed in roughly 1 Myr under the combined effects of internal and external photoevaporation (Fig. 7).
Despite the dissimilarities, a reasonably good match to the disc properties of the Lupus and Chaemeleon I low-mass starforming regions is obtained starting with more massive discs of smaller sizes, and with only internal photoevaporation.This is valid for clusters with low ambient field strengths, such as Lupus (4 G 0 ; Cleeves et al. 2016), where losses due to external photoevaporation are low.The larger masses and smaller sizes are needed to improve the match in stellar accretion rates and observed masses.The corresponding initial conditions and model parameters can be used to study planetary formation in similar environments.A more robust comparison with observations is performed in Paper III.
However, initial disc masses cannot be arbitrarily increased or discs would become gravitationally unstable.Instead, we suggest that future studies should include the modelling of the initial stages of disc formation, including the presence of an envelope.This envelope would allow discs to be replenished after their initial formation and provide shielding from UV radiation from nearby stars.Also, magnetically driven disc winds would shield UV and X-ray radiation from the central star.This would provide a reduction of losses by both internal and external photoevaporation.Both effects allow for longer lifetimes and larger masses at later times without the need for extremely large masses at earlier times.
We decided here to use observed dust masses as the main comparison point.However, this is not the only possible avenue.For instance, disc radii could be less susceptible to the degeneracy caused by regions that are optically thick (Pascucci et al. 2022).One likely difficulty would be in accounting for the large disc radii and sustained stellar accretion rate.Already with the comparatively small discs that we find to best match the disc mass-stellar accretion rate relationship, we are not able to reproduce the largest observed accretion rates.Having larger discs would lead to a reduction of the stellar accretion rates for a given disc mass (Fig. B.1); which would lead to a mismatch with observations of stellar accretion rates.
In this work, we assume that the gas discs evolve viscously.However, simulations of disc evolution that account for nonideal magnetohydrodynamical (MHD) effects find that the magnetorotational instability (MRI), which is the likely mechanism generating the turbulence, is largely suppressed (e.g.Bai & Stone 2013;Lesur et al. 2014).Instead, it has also been proposed that the evolution is driven by magnetically driven disc winds (Suzuki & Inutsuka 2009).Magnetically driven disc wind prescriptions, such as those of Suzuki et al. (2016) or Tabone et al. (2022), include several model possibilities, which can be narrowed down by performing a similar comparison to that presented here (Weder et al. in prep.).Once such a model is properly coupled with internal photoevaporation to account for shielding, a similar study to that presented here is possible.

Fig. 1 .
Fig. 1.Performance of three surrogate models based on the comparison of the predicted and actual values of the testing set.The insert values show the best regression (ordinary least squares), the Pearson correlation coefficient R 2 , and the RMS of the differences between each predicted and actual value.

Fig. 4 .
Fig. 4. Cumulative distribution of disc lifetimes for a population with canonical parameter distribution (see text).Two exponential decays following Fedele et al. (2010) with a characteristic time of 2.3 Myr (accretion) and 3 Myr (infrared excess) and the results of Ribas et al. (2014) are shown as well.

Fig. 5 .
Fig. 5. Histogram for stellar accretion rate vs. disc mass at 2 Myr and the same synthetic disc populations shown in Fig. 4. The observed data from the Lupus (in red) and Chamaeleon I (orange) star forming regions from Manara et al. (2019) are shown for comparison.Two disc dispersal timescales τ = ṀG /M G of 1 Myr and 10 Myr (assuming that the gas disc mass M G is 100 times the observed dust mass) and the best fit to the data from Manara et al. (2016b) are shown as well.

Fig. 7 .
Fig. 7. Evolution of the relative gas mass: in the disc (blue), accreted onto the star (green), and lost by internal (orange) or external (red) photoevaporation until disc dispersal at 1.38 Myr.This represents a typical disc, with M = 0.5 M , M G /M = 0.1, β = 0.9, R in = 0.1 au, r 1 = 87.8au, α = 1 × 10 −3 , L X = 7.02 × 10 29 erg s −1 , and F = 10 G 0 .The parameters for the solid disc are irrelevant, as this shows only the gas component, except that we used f D/G = 0.0149 to compute the solid mass needed to set the characteristic radius r 1 .

Fig. A. 1 .
Fig. A.1.Gas disc lifetime (left), observed dust masses at 2 Myr (centre), and stellar accretion rates at 2 Myr (right) as functions of the power-law index of the initial profile β and the inner edge r in (top) or viscosity parameter α and characteristic radius r 1 (bottom) In the bottom row, the white region is where disc lifetimes are less than 2 Myr and the values cannot be constrained by the model.

Fig. A. 2 .
Fig. A.2. Observed mass at 2 Myr as function of the parameters of the twopop model: gas-disc mass and dust-to-gas ratio (left), planetesimal formation efficiency and fragmentation velocity (centre), and drift efficiency and fragmentation velocity (right).
Fig. B.1.2D histograms for stellar accretion rate versus disc mass at 2 Myr, showing the effects of the different parameters, as given above each panel.The top row shows the effect of increasing the initial disc mass M D .The middle row shows the effect of reducing the disc characteristic radius r 1 .The bottom row shows the effect of the disc viscosity parameter α.The observed data from the Lupus (in red) and Chamaeleon I (orange) star-forming regions from Manara et al. (2019) are shown for comparison.

Table 1 .
Parameter range for the main simulation grids.

Table 2 .
Best-fit parameters for stellar X-ray luminosity.

Table 5 .
Hyper parameters and performance of the surrogate models.R 2 Test.R 2 Val.MSE Test.MSE Lifetime L-BFGS logistic 25, 50, 45 2.592 × 10 −3 0.99465 0.99436 vations), which is consistently lower than the median L X by a factor of a few.

Table 6 .
Random distributions for the canonical population