Disc population synthesis: Decrease in the solid mass reservoir through pebble drift

Surveys of star-forming regions reveal that the dust mass of protoplanetary discs decreases by several orders of magnitude on a timescale of a few million years. This decrease in the mass budget of solids is likely due to the gas-drag-induced radial drift of mm-sized solids, called pebbles. However, quantifying the evolution of this dust component in young stellar clusters is difficult due to the inherent large spread in stellar masses and formation times. Therefore, we aim to model the collective evolution of a cluster to investigate the effectiveness of radial drift in clearing the discs of mm-sized particles. We use a protoplanetary disc model that numerically solves for disc formation, and the viscous evolution and photoevaporative clearing of the gas component, while also including the drift of particles limited in size by fragmentation. We find that discs are born with dust masses between 50 Earth masses and 1000 Earth masses, for stars with, respectively, masses between 0.1 solar masses and 1 solar masses. The majority of this initial dust reservoir is typically lost through drift before photoevaporation opens a gap in the gas disc for models both with and without strong X-ray-driven mass loss rates. We conclude that the decrease in time of the mass locked in fragmentation-limited pebbles is consistent with the evolution of dust masses and ages inferred from nearby star-forming regions when assuming viscous evolution rates corresponding to mean gas disc lifetimes between 3 Myr and 8 Myr.


Introduction
Surveys of young star-forming regions, such as the 1 Myr-old Perseus cluster, reveal typical dust disc masses that exceed 100 M C (Tychoniec et al. 2020;Tobin et al. 2020). Older starforming regions, on the other hand, often show vastly reduced dust masses by up to two orders of magnitude, which are albeit characterised by a large spread (Andrews et al. 2013;Ansdell et al. 2016;Pascucci et al. 2016;Eisner et al. 2018;Ruíz-Rodríguez et al. 2018;Ansdell et al. 2018;Cieza et al. 2019;Cazzoletti et al. 2019;Tobin et al. 2020;Grant et al. 2021). However, care should be taken in the interpretation of individual dust disc mass measurements, as these depend on model choices such as the chosen opacity prescription, particle size distribution, and assumptions of optically thin continuum emission (Liu 2019;Zhu et al. 2019).
The dominant process responsible for removing dust is not yet completely clear. Two likely candidates are the radial drift of dust pebbles and the formation of planetesimals and planets. Radial drift occurs due to gas drag from gas on sub-Keplerian orbits removing angular momentum from the pebbles, which causes them to drift inwards in the disc (Weidenschilling 1977). If dust particles are sufficiently large, this occurs in most regions of protoplanetary discs but can be halted at locations with a pressure bump where gas drag is not present (Brauer et al. 2008;Pascucci et al. 2016). In Appelgren et al. (2020), we studied the effect of radial drift in population synthesis models and found radial drift to be compatible with typical measurements of disc dust masses. However, this study only traced a single dust size and did not include gas disc clearing by photoevaporation. The effect of limiting dust growth by fragmentation and the impact of photoevaporation is the focus of the present population synthesis study.
Estimates of dust grain sizes in protoplanetary discs have found maximum sizes on the order of mm-cm, based on measured opacity slopes (Pérez et al. 2012(Pérez et al. , 2015Tazzari et al. 2016;Carrasco-González et al. 2019). However, if maximum dust sizes are estimated from the degree of polarisation in the emission from protoplanetary discs, typically lower maximum sizes of 0.10-0.15 mm are found (Kataoka et al. 2015(Kataoka et al. , 2016, although 1 mm is still possible for very low turbulence values (Ueda et al. 2021). Dust coagulation is a complex process whereby the collision of two particles can result in very different outcomes, including sticking, bouncing, fragmentation, mass transfer, cratering, and erosion, depending on the sizes and velocities of the two particles (Blum & Wurm 2008;Zsom et al. 2010;Okuzumi et al. 2011a,b;Blum 2018). The outcome of these coagulation processes sets the size, and therefore Stokes number, of the largest particles that characteristically hold a very large fraction of the total disc dust mass (Birnstiel et al. 2010). Importantly, the dominant particle size affects not only the speed of the radial drift of these particles but also their ability to form planetesimals via the streaming instability (Youdin & Goodman 2005;Johansen et al. 2007;Carrera et al. 2015;Yang et al. 2017;) and Article number, page 1 of 16 arXiv:2303.07297v3 [astro-ph.EP] 1 Aug 2023 A&A proofs: manuscript no. 45252corr the efficiency of pebble accretion (Ormel & Klahr 2010;Lambrechts & Johansen 2012. Protoplanetary discs are formed together with their host star when dense cores of giant molecular clouds collapse under their own gravity. The formation of the disc is thought to last a few hundred thousand years up to " 1 Myr (Hueso & Guillot 2005;McKee & Ostriker 2007). As the disc continues to evolve, gas accretion removes gas from the disc over a timescale of a few to several million years (Hartmann et al. 1998;Manara et al. 2022). The accretion of gas is driven by a redistribution of angular momentum in the discs. The mechanism behind the redistribution is not yet fully understood, but it is believed to be disc turbulence driven by the magnetorotational instability (MRI) (Balbus & Hawley 1998) or disc winds (Bai & Stone 2013), or a combination of the two. For a recent review on magneto-hydrodynamic disc processes, we refer to Lesur et al. (2022). Due to the poor understanding of what drives the turbulence, this is commonly parameterised by the well-known α-disc model (Shakura & Sunyaev 1973). During this time, gas accretion rates drop from as high as 10´6 M d /yr, for the youngest discs, to À 10´9 M d /yr for discs as old as a few million years (Manara et al. 2016). Simultaneously, the gas disc mass is reduced from " 0.1 M d to " 0.001 M d (Miotello et al. 2022). However, estimating gas disc masses from, for example, CO gas emission, is highly uncertain. Therefore, dust disc masses have been used to infer gas disc masses, under the crude assumption of a global dust-to-gas ratio of 0.01 (e.g. Manara et al. 2022). This value is typical of the interstellar medium, but no longer valid when considering the effects of pebble drift (Appelgren et al. 2020). Furthermore, dust disc masses themselves may be underestimated if the millimeter emission is optically thick.
It is important to include the correct range of stellar masses when trying to model the evolution of a cluster of discs. Lower mass stars are more commonly formed than solar-mass stars according to the stellar initial mass function (Salpeter 1955;Kroupa 2001;Chabrier 2003). The dominating presence of subsolar stars is of importance because the mass of dust discs has been found to scale with the stellar mass in nearby star-forming regions as M dust 9M 1.3´2.4 ‹ , although with a large scatter (Pascucci et al. 2016;Ansdell et al. 2017). Similarly, a trend has been found where higher stellar gas accretion rates relate to more massive dust discs, but again with a significant scatter (Manara et al. 2016).
The final dispersal of the gas component of protoplanetary discs is suggested to be driven by photoevaporative winds (e.g. Hollenbach et al. 1994;Shu et al. 1994;Clarke et al. 2001). A combination of extreme ultraviolet (EUV), X-ray, and far ultraviolet (FUV) radiation heats the gas at the disc surface to the point where the gas becomes gravitationally unbound and is lost from the disc. When the gas removal rate via photoevaporation and the gas accretion rate through the disc become comparable, the inner disc is no longer efficiently replenished by gas. This results in the efficient dispersal of the inner disc, creating an inner disc cavity that expands outwards with time.
Several prescriptions of photoevaporation aimed for use in population synthesis models have been developed. Owen et al. (2012) have provided a prescription for X-ray-driven photoevaporation, which was extended by an updated X-ray photoevaporation prescription (Picogna et al. 2019;Ercolano et al. 2021;Picogna et al. 2021). Komaki et al. (2021) have provided an alternative prescription for photoevaporation driven by EUV, FUV, and X-rays, using a method similar to that of Nakatani et al. (2018b). The importance of each of these drivers has been the subject of a number of studies. In particular, EUV and X-rays are more straightforward to model, whereas FUV is a bigger challenge because it is affected by dust evolution and chemistry (Gorti et al. 2015). Owen et al. (2012) concluded that X-rays are the dominant driver of photoevaporation. Picogna et al. (2019) similarly argued that X-rays dominate over EUV. In contrast, Nakatani et al. (2018b) found that FUV-driven photoevaporation is dominant. We explore the impact of both the Picogna et al. (2021) and the Komaki et al. (2021) models in the following.
In this paper, we first elaborate on the method used for the population synthesis model (Sect. 2). We then present the results of the study. We focus on the evolution of the disc dust mass and how its evolution is affected by the strength of the assumed disc viscosity and the two different photoevaporation prescriptions available in the literature (Sect. 3). This is followed by a discussion of our results, their stellar mass dependency, and the implications on the gas and dust radius evolution with time (Sect. 4). Finally, we present our conclusions (Sect. 5).

Model
In this paper, we present an update to the disc model used in Appelgren et al. (2020). As in the previous paper, we calculate the evolution of discs around stars in a newly formed stellar cluster using viscous gas evolution and radial drift of solids. We briefly review the used methodology here. In addition, here we evolve the size of the dust particles, while including the clearing of the gas disc by photoevaporation and using a more complete temperature treatment.

Disc formation
The formation of the protoplanetary disc is modelled from the collapse of an over-dense Bonnor-Ebert sphere, similar to the work in Hueso & Guillot (2005); Takahashi et al. (2013). The collapse occurs inside-out with the collapse front expanding outwards over time. The infall rate of gas mass onto the disc is given by: where r cf is the radius of the collapse front, ρpr cf q is the density of the cloud core at this radius and dr cf {dt is the velocity at which it expands outwards. A more detailed description of how these quantities are calculated can be found in Appendix A of Appelgren et al. (2020). The gas is distributed across the disc according to the following equation: The infalling material lands on the disc within the centrifugal radius R c , given by: where Ω 0 is the rotation rate of the cloud core and Mpr c q is the mass interior to the collapse front radius. We assume here that the angular momentum of the infalling material is conserved. This is a simplification that ignores the role of magnetic braking (e.g. Galli et al. 2006;Wurster et al. 2019) and possible anistropic accretion (e.g. Smith et al. 2011;Kuffmeier et al. 2017). However, a detailed modelling of the angular momentum of molecular cloud cores and their evolution during the collapse phase is beyond the scope of this paper.

Gas disc evolution
We evolved the disc viscously using the well-established alphadisc prescription of Shakura & Sunyaev (1973), where the viscosity of the disc is described by ν " αc 2 s {Ω. Here, c s is the sound speed and Ω is the rotation rate in the disc, while α sets the strength of the viscosity. We treat α as the value which evolves the disc on a timescale consistent with measured disc lifetimes. We do not treat it as a measurement of the true level of midplane turbulence stirring the dust midplane of the disc, which has been observationally inferred to be much lower in the disc (Pinte et al. 2016;Villenave et al. 2020Villenave et al. , 2022. Instead, we parameterise the latter as α t , which we also use to calculate the size of fragmentation-limited dust (see Sect. 2.3.2).
Using accretion rates of pre-main sequence stars, Hartmann et al. (1998) found α " 10´2. However, recent estimates of linewidth broadening by turbulent motions in outer disc regions have found values of α À 10´3 Flaherty et al. 2018Flaherty et al. , 2020. Moreover, measurements of the fraction of stars with discs in star-forming regions suggest that the discs dissipate with a characteristic timescale of 2.5-3 Myr (Haisch et al. 2001;Mamajek 2009). Such disc lifetimes of " 3 Myr would suggest high α-values to drive disc evolution. However, more recent measurements of the disc fraction finds that, in those starforming regions that are not subject to external photoevaporation, the characteristic timescale for disc dissipation is " 8 Myr (Michel et al. 2021;Pfalzner et al. 2022). The slower disc evolution (indicated by the longer disc dissipation timescale) would suggest lower values of α. We used a base value, which we labelled as α ν , of either 10´2 or 10´3 to explore both scenarios of faster and slower disc evolution. In practise, α is calculated as: which increases the disc viscosity when the disc becomes gravitationally unstable following Zhu et al. (2010). Here, Q is the Toomre parameter expressing the gravitational stability of the disc. The viscous evolution of the disc is calculated using the equation from Pringle (1981). The total surface density evolution of the disc is then given by: where 9 Σ g,PE is the mass loss due to internal photoevaporation. The way we calculate this is explained in Sect. 2.5 and Appendix A.

Dust evolution
The infall rate of dust onto the disc is the infall rate of gas multiplied by the assumed dust-to-gas ratio of the cloud core. We assume this to be uniformly Z " 0.01.

Dust drift
The evolution of the dust in the disc is governed by the radial gas flow and the radial drift of the dust particles. Radial drift occurs because the dust feels a slight headwind from the gas that orbits at a sub-Keplerian velocity (Weidenschilling 1977). This headwind removes angular momentum from the dust causing it to drift inwards in the disc. The speed at which the dust drifts is given by: Here, v k is the Keplerian orbital velocity and η is a measure of the radial pressure support in the disc, given by: The Stokes number, St, of the dust particles is a measure of how aerodynamically coupled particles are to the gas. The ratio H{r expresses the gas disc aspect ratio. To calculate the Stokes number we take into account the Epstein and the Stokes regimes, given by: Here, ρ ‚ is the material density of the dust particle which we set to 1.6 g/cm 3 (Birnstiel et al. 2012). The size of the dust grains, a p , is discussed further in Sect. 2.3.2. Because the dust is coupled to the gas, small grains simply move along with the gas in the disc. Taking this into account, the total radial velocity of the dust, v d , is given by: where v g is the radial velocity of the gas. The surface density of dust is then evolved with the following continuity equation:

Dust sizes
We followed the evolution of the largest grains that carry most of the mass in typical dust size distributions (Birnstiel et al. 2010). The dust grows from an initial monomer size of 1 µm to the fragmentation limit on a timescale given by: For grains that reach their maximum size limit by fragmenting, the maximum grain size is set via where f f " 0.37 is a fudge factor to better match detailed coagulation simulations (Birnstiel et al. 2012). We set the level of the midplane turbulence α t to 10´4, in line with the lower α values inferred from observations Pinte et al. (2016); Villenave et al. (2022). We note that we treat it as a separate parameter from α ν (Eq. 4). The fragmentation threshold velocity, u f , depends on the material properties of the grains. Laboratory experiments of silicate dust grain collisions suggest limits of about 1 m/s (Blum & Wurm 2008). Water ice was thought to be stickier than silicates and to have a fragmentation speed of about 10 m/s (Gundlach & Blum 2015). However, more recent studies have concluded that the fragmentation speed of water ice is likely to be lower. The stickiness of CO 2 -ice is also thought to be quite low, with fragmentation speeds of about 1 m/s (Musiolik et al. 2016a). Similar results are found for mixtures of water and CO 2 ices (Musiolik et al. 2016b). Therefore, the choice of 1 m/s as a global fragmentation velocity should be suitable. Particles that are not limited by fragmentation in their growth naturally approach their drift-limited sizes in our model. Their maximum size is then bound by: where f d " 0.55 is a fudge factor with similar purpose as f f and γ " | d ln P d ln r | (Birnstiel et al. 2012).

Disc temperature
The temperature of the disc was calculated similarly to the work in Kimura et al. (2016). We included irradiation and viscous heating. The temperature due to irradiation is given by (e.g. Menou & Goodman 2004): where L ‹ is the stellar luminosity. We assume that the luminosity of the star scales linearly as L ‹ " L d pM ‹ {M d q ( Baraffe et al. 2002;Liu et al. 2019). The heating rate due to irradiation is given by (e.g. Menou & Goodman 2004): where T core " 10 K is the temperature of the cloud core, τ " κΣ{2 is the optical depth, and κ is the opacity. The heating rate due to viscous heating is given by (e.g. Kimura et al. 2016): Heat is lost from the disc by radiative cooling (e.g. Menou & Goodman 2004): where T mid is the midplane temperature of the disc. The midplane temperature due to irradiation and viscous heating is given by the balance of the heating and cooling rates. We used an opacity prescription similar to that in Kimura et al. (2016), where the opacity is that of solid grains and given by κ " 0.052 cm 2 {g pT {Kq 0.738 (Zhu et al. 2009). This opacity scaling assumed a solar composition of dust grains. However, above a temperature of " 1500 K, most dust will have sublimated. This will cause the opacity to drop and the disc to cool until the dust solidifies again. Unless the much lower gas opacity is enough to drive the disc to higher temperatures, this results in a temperature plateau around 1500 K. This effect has been found in more complex disc temperature models (see e.g. Birnstiel et al. 2010;Schobert et al. 2019;. We therefore imposed a temperature limit of 1500 K in our model. Photoevaporation removes gas from the surface of the protoplanetary disc by either EUV, FUV, or X-rays heating the gas until it can escape the disc (Gorti et al. 2009(Gorti et al. , 2015Owen et al. 2010Owen et al. , 2011Owen et al. , 2012Wang & Goodman 2017;Nakatani et al. 2018a,b). In this work, we look at two different photoevapora- The two prescriptions differ in the physics that is included, and also in the methods used to derive the photoevaporation rates. To model the photoevaporation, Pic21 included X-rays and EUV radiation, although the X-ray component is the dominant contributor to the heating and mass loss rate (Picogna et al. 2019). The model from Kom21 includes X-rays, EUV, and FUV radiation, with the FUV radiation being the dominant source of heating. In this model, the X-ray contribution is negligible (as illustrated in Fig. 1), where mass loss rates are driven by the FUV radiation, even if the FUV fraction would be reduced by a factor of 100. We therefore chose to present both of these models to reflect the uncertainty in modeling photoevaporative mass loss (for a further discussion, see Sect. 4.4).

Photoevaporation
The equations and the parameters used to calculate the photoevaporation rates are given in Appendix A. To find the parameters for stellar masses that fall in between the ones given in either model, we linearly interpolated over the grid of masses. For stellar masses that fall outside the range of masses modelled in the prescriptions, we used the parameters for the lowest and highest mass given to avoid extrapolating. We did this because, in general, the parameters of either prescription are not monotonically increasing or decreasing with stellar mass. For the Pic21 model, this upper mass is 1 M d and stars above this mass will prove very uncommon in our sample; therefore, the choice of using the 1 M d cut off will be of little consequence. With the Kom21 model the lower mass is 0.3 M d and thus more stars will be affected since our stellar masses extend to À 0.1 M d , possibly leading to moderately overestimating the photoevaporation rate for very low-mass stars. , gas surface density (middle), and dust surface density (right) evolution for a disc forming from a 1 M d cloud core with a centrifugal radius of 10 au and photoevaporation according to Pic21. The dust is allowed to grow to the fragmentation limit. In the innermost part of the disc the temperature can reach 1500 K due to viscous heating and in the outer disc the temperature is set purely by irradiation. Photoevaporation begins to open up a gap in the gas disc around 3.75 Myr. By 4 Myr, the inner gas disc is completely dispersed, leaving behind an outer depleted disc. We find that 0.6 M Earth of the dust is retained outside this gap and is pushed outwards with time.

Population synthesis
The main cloud core parameters that are varied in the population synthesis are the mass and angular momentum. The masses of cores in molecular clouds are described by the so-called coremass-function (CMF). This is similar in form to the initial mass function of stars (IMF) (McKee & Ostriker 2007). Therefore, we sampled masses of the cloud cores from the Kroupa IMF (Kroupa 2001). We limited our sampling to masses between 0.1-2 M d .
To set the angular momentum of the cloud core, we set the maximum centrifugal radius, which is connected to the angular momentum of the cloud core through Eq. 3. This is the largest radius at which any material will land on the disc during its formation. We select centrifugal radii from a normal distribution with a standard deviation of 30 au and a mean scaled linearly with the cloud core mass as: as in Appelgren et al. (2020, specifically, their Appendix D). In order to limit the number of parameters varied in the population synthesis, we chose to use a fixed value for the temperature of the cloud core of T core " 10 K across all simulations. The stars in a young star-forming cluster do not all form at the same time. For instance, star formation in the star-forming region of Taurus is thought to have been ongoing for 1-2 Myr (Kenyon & Hartmann 1995). Across Orion, for discs in isolation, the spread in disc mass is likely related to an age spread (van Terwisga et al. 2022). Therefore, for the internal age spread of the discs, we used a value of 1 Myr.
When a protoplanetary disc forms, it will inherit the dust-togas ratio of its parent molecular cloud core. The metallicity of a clump in a giant molecular cloud which forms s stellar cluster is thought to be very homogeneous, due to a very short mixing timescale (Feng & Krumholz 2014). Hence, we did not vary the initial dust-to-gas ratio of the cloud cores and kept the nominal value of Z " 0.01 for all cloud cores in the population synthesis.
Finally, we ran two population synthesis copies with the same cloud core masses and centrifugal radii, but one using the Pic21 prescription and the other using the Kom21 prescription. The discs are evolved for 10 Myr using α ν " 10´2 (Sect. 3.2) and α ν " 10´3 (Sect. 3.3) for disc lifetimes of 2.5 Myr and 8 Myr (Mamajek 2009;Michel et al. 2021), respectively.

Evolution of a single disc
Before presenting the population synthesis results, we show the evolution of a single protoplanetary disc to illustrate the key physical processes at play. For this purpose, we initiated a model with a molecular cloud core with a mass of 1 M d and a maximum centrifugal radius of 10 au. We used an α-viscosity of α ν " 10´2, together with the Pic21-photoevaporation prescription. Other model parameters are the same as described in Sect. 2.6.
Viscous evolution dominates the gas disc for the first " 3.7 Myr, after which photoevaporation begins to quickly carve out a gap in the gas. A gap in the disc forms at an orbital distance of about 3-4 au. This is seen in Fig. 2, which shows the evolution of the temperature (left panel), gas surface density (middle panel), and dust surface density (right panel). Once the gap opens up, dust drift over this gap is prevented because the decreased pressure gradient at the gap edge (Fig. 2, panel b). For this nominal case, we find that 0.6 M Earth of dust was trapped outside the gap, which gradually is pushed outwards as the inner cavity expands.
The radially inwards flux of pebbles is highest in the inner disc, early in its evolution (Fig. 3, left panel). However, during the early times (À 0.35 Myr) dust also has an outwards motion in the outer disc region between 10-400 au, because the gas in this region is expanding viscously outwards and the dust is wellcoupled to it. The dust flux drops over time until a gap is opened up by photoevaporation. The dust remaining within the gap is quickly pushed inwards, whilst the dust retained outside the gap is slowly pushed outwards. This latter effect is represented by the thin sliver of dust shown in green in the top right of the left panel in Fig. 3.
The dust-to-gas ratio exceeds the initial value of 0.01 at two regions in time and space (Fig. 3, middle panel). It occurs when the disc is 0.2-0.4 Myr old at radii between 0.1-5 au, and at a later time, between 0.6-0.9 Myr, at radii between 0.1-1.5 au. The high dust-to-gas ratio makes these potential sites for planetesimal formation (Yang et al. 2017;Appelgren et al. 2020;Li & Youdin 2021) (see also Sect. 4.1.3).
Particles do not grow beyond cm in size, except at early times, within À 1 Myr (Fig. 3, right panel). Around 1 au there is a local minimum in the particle size distribution. This region is located where the disc transitions from the maximum tempera-Article number, page 5 of 16  Fig. 3: Dust flux (left, absolute values), dust-to-gas ratio (middle), and dust grain size (right) evolution of a disc forming from a 1 M d cloud core with a centrifugal radius of 10 au and photoevaporation according to Pic21. When the disc is young and the gas expands outwards by viscous evolution, the dust in the outer disc also moves outwards as it is strongly coupled to the gas. The red line marks the radius at which the dust transitions from moving inwards to moving outwards in the disc. Once photoevaporation opens a gap in the gas, the dust interior to this line is quickly drained. The dust remaining outside the gap is slowly pushed outwards as the gas at the outside edge of the gap photoevaporates. Two islands of a dust-to-gas ratio higher than the initial value of 0.01 appear. One at 0.2-0.4 Myr at radii between 0.1-5 au, and one at 0.6-0.9 Myr at radii from 0.1-1.5 au. The particle size generally is larger in the inner regions of the disc and at earlier times.
ture of 1500 K to lower temperatures at wider radii, where stellar irradiation is dominant in setting the disc temperature rather than viscous heating.

Population synthesis with α ν " 10´2
In this section, we present the results from the population synthesis model using α ν " 10´2. The value of α ν which best matches observations of protoplanetary discs lifetimes is still a matter of debate, as discussed in Sect. 2.2. We show the results for α ν " 10´3 in Sect 3.3.

Gas disc evolution
The evolution of the gas accretion rate over time for all discs in the population synthesis is shown in Fig. 4. Discs evolving under Pic21 photoevaporation (red curves) are dominated by viscous evolution for about 1-3 Myr before photoevaporation takes over. In the Kom21 prescription (black curves) discs undergo viscous evolution for longer times. Photoevaporation clears the inner parts of discs at the earliest after about 5-6 Myr and for some discs, photoevaporation never plays an important role during the 10 Myr that we simulated in this work. Newly formed discs have gas masses ranging from 0.03 M d to 0.5 M d . These gas discs are then depleted by accretion and photoevaporation (Fig. 5). For both the Kom21 and Pic21 photoevaporation prescriptions, the mass reservoir diminishes until approximately 1 % of the initial mass remains. This remainder of gas is found outside the photoevaporative gap at large orbital radii. Therefore, it takes excessively long to deplete by photoevaporation.

Dust disc evolution
The radial drift of pebbles reduces the dust mass of protoplanetary discs on a timescale of a few Myr. This is illustrated in the middle panel of Fig. 5, which shows the evolution of the dust mass as a function of disc age. Discs are formed from cloud cores with gas masses between " 0.1 M d and " 1 M d , and are born with dust disc masses between of 50 M C and 1000 M C respectively. The least massive and shortest lived discs are depleted of dust within less than 1 Myr, whilst the most massive and longest  Fig. 4: Gas accretion rate as a function of disc age for all disc in both photoevaporation prescriptions. Once photoevaporation opens up a gap in the gas disc, the gas within disc gap is quickly cleared by photoevaporation and viscous evolution. This results in a rapid reduction in the gas accretion rate. In this figure, we show the results without including the spread in the time of formation of the stars in a cluster (see Sect. 2.6).
lived discs last up to 7 Myr. This process is also seen in the evolution of the cumulative distribution of dust disc masses at cluster ages between 0.5 Myr and 5 Myr (as shown in Fig. 6). In the 0.5 Myr-old cluster, 73 % of discs have a disc dust mass higher than 10 M C . At 1 Myr, this has decreased to only 45 % of discs having over 10 M C of dust. From here on, we see a continued decrease in the available disc dust mass as more and more dust drifts onto the star. We do not find that photoevaporation is able to retain significant amounts of dust by gap opening. Using the Kom21 prescription, all discs retain ă 0.1 M C of dust. With the ) model, we find that 20 % of all discs are able to retain at least 0.1 M C of dust. However, among these discs, the median mass retained is only 0.3 M C . The median mass of these discs right at the end of their formation, which is when they are the most massive, is 367 M C . Hence, these discs retain less than 0.01 % of their maximum dust mass. For comparison, the Fig. 5: Evolution of the gas disc masses (left), dust disc masses (middle), and dust-to-gas ratios (right) as a function of disc age using the α ν " 10´2. Gas discs are formed with masses between 0.03 M d to 0.5 M d . After 10 Myr of evolution about 1 % of this initial mass remains. Dust disc masses are initially between 10 M C to 1000 M C . The shortest lived dust discs are depleted in " 1 Myr, whilst the longest lived last up to 3-8 Myr, depending on if either the Pic21 or Kom21 photoevaporation prescription is used. Discs are formed from material with a dust-to-gas ratio of Z " 0.01, but, due to rapid loss of dust in the earliest phases of disc formation, the global dust-to-gas ratio of the discs is always below this value.
observed median dust masses in protoplanetary disc surveys at " 2 Myr are 3.0 M C (Lupus), 3.2 M C (Cham I), 6.2 M C (Taurus), and 11.5 M C (Cham II) (Andrews et al. 2013;Ansdell et al. 2016;Pinte et al. 2016;Villenave et al. 2021). The retention of dust by gap opening can also be seen in Fig. 7, where the gas accretion rate as a function of dust disc mass is shown for each disc using both photoevaporation models. Discs start out with high gas accretion rates (" 10´6 M d {yr) and low disc dust masses. For the Pic21 model (red curves) some of these discs reach a phase where the gas accretion turns over at rates between 10´8 and 10´9 M d {yr, due to photoevaporation clearing the inner disc of gas, whilst retaining some dust. These are the discs where photoevaporation has trapped dust outside the gas gap. None of the discs using the Kom21 model (black lines) go through this phase, showing that all these discs lost all of their dust before photoevaporation opened up a gap.
The efficient depletion of dust by radial drift agrees with our previous results from Appelgren et al. (2020) where radial drift is the dominant process for removing dust from protoplanetary discs. Other mechanisms of dust removal, such as planetesimal formation and subsequent planet formation, would only further deplete the dust reservoir (further discussed in Sect. 4.1.3).
If radial drift would not be effective, the dust would evolve by advecting with the gas and therefore be lost at the same rate as the gas. The cumulative distribution of disc dust masses in this case is shown with the thin coloured lines in the left plot of Fig.  6. For this experiment, we calculated these dust masses by multiplying the gas mass of each disc with the initial dust-to-gas ratio and kept it fixed at a constant value of Z " 0.01 in our simulations (Sect. 2.6). In this case, the dust is lost very slowly and the profile of the cumulative distribution remains similar throughout the evolution of the cluster.
For comparison, we show a sample of observed dust disc masses in different star-forming regions in the right plot of Fig. 6.  Manara et al. (2022). To get the occurrence fraction with respect to the full population of stars in the cluster, we scaled the cumulative distributions of this sample according to Eq. B.1 (for more, see Appendix B), assuming a gas disc lifetime of 2.5 Myr (Mamajek 2009). We note that the model disc masses presented in this paper are the masses as given directly by the model, that is the integrated surface densities. We do not take into account optical depth effects of the dust to model the dust disc mass that would effectively be observable. The ob-servable dust disc mass would be expected to be somewhat lower than the model masses.
At a cluster age of 0.5 Myr and 1 Myr, the median dust disc mass of the model is " 30 M C (0.5 Myr) and " 5.5 M C (1 Myr). This can be compared to the Perseus sample, which has significantly higher median masses for class 0 objects (158 M C ) and class I objects (52 M C ) (Tychoniec et al. 2020). In these early stages of disc evolution, radial dust drift could be less efficient than we find here. However, determining the masses of young embedded discs remains challenging (Tobin et al. 2018).
For the more evolved discs we find that the distribution of observed discs aged between 1-3 Myr all lie in between the 1-3 Myr old cluster in our population synthesis. The observed disc dust masses lie closer to the 2-3 Myr older model clusters than the 1 Myr model clusters. For the oldest observed cluster, Upper Scorpius, only the oldest model cluster (5 Myr) with the Kom21 photoevaporation model lies somewhat close. Using the Pic21 photoevaporation prescription, too much dust is retained in too many discs. However, the fate of solids concentrated at the cavity edge is uncertain (see Sects. 4.1.4 and 4.4).
If we compare the observed sample to the expected dust mass for models without radial drift (thin coloured lines in Fig. 6) it is clear that the match is very poor, with the exception of the Perseus discs. This indicates that dust in discs that are a few Myr old has been lost at a higher rate than the gas. Therefore, the sometimes-used assumption that the gas mass in protoplanetary discs is 100 times their dust mass does not hold. We see this in the right panel of Fig. 5, which shows the evolution of the global dust-to-gas ratio of all discs as a function of their age. We note that, even though the cloud cores from which the discs form initially have dust-to-gas ratios of 0.01, efficient radial drift of dust very early in the disc evolution decreases the global dust-togas ratio already very early on in the disc's evolution (À 10 5 yr).
Although radial drift is able to qualitatively match the cumulative disc dust mass distributions well, there is one component of the observational sample which is not well reproduced: the observed clusters show that there is a fraction (" 10%) of discs that retain large amounts (ě 10M C ) of dust at late times (van der Marel & Mulders 2021). In Appendix C we show an experiment of halting dust drift in a subset of the most massive discs. We find that stopping dust drift in a fraction as low as 5 % of the discs can extend the high-mass tail of the CDF, bringing the model more in line with the observations. This supports the claim that " 5%´10% of discs have a reduced efficiency, or complete halting, of radial dust drift. This group of discs can represent a  Fig. 7: Gas accretion rate as a function of dust disc mass for the models using the Pic21 (black) and the Kom21 (red) photoevaporation prescriptions . The grey dots show the observational sample from Manara et al. (2022). Each line represents one disc's evolutionary track. Discs start out at low dust masses and high gas accretion rates (" 10´6 M d {yr). As they leave the formation phase, they pass through a phase where the gas accretion rate drops rapidly and the dust mass changes slowly. After this, the disc dust mass begins to quickly drain. If photoevaporation opens up a gap in the gas and clears the inner gas disc during this phase, the gas accretion rate drops rapidly and some dust is trapped. Otherwise, the gas accretion rate remains nearly unchanged while the dust drains.
population of discs in which radial drift is less efficient. A likely mechanism for preventing radial drift of pebbles is the presence of pressure bumps in the protoplanetary disc, possibly triggered by wide-orbit planets. We discuss this further in Sect. 4.1.2.
An important measure of disc evolution is the relation between the gas accretion rate onto the star and the dust disc mass. This is shown in Fig. 7, where the lines shown the model and the grey dots show observational measurements from (Manara et al. 2022). Here, we see that a large region of parameter space with relatively low accretion rates of less than about 10´8 M d /yr and dust disc masses of over 1 M C is populated by the observed sample. Our synthesis model with the Pic21 photoevaporation does not produce such dust-rich discs with low gas accretion rates. This region contains roughly half of the of observed sample. Using the Kom21 photoevaporation prescription, an even smaller region of the observed sample is reached. Only those discs with accretion rates of about 10´8 M d /yr or more are reproduced. The disagreement between model and observation in the gas accretion rate and dust disc relation could originate from an overestimation in the model of the gas accretion rate onto the host star. Lower gas accretion rates onto the star are possible if disc winds remove some gas from the inner disc before it is accreted onto the host star (Manara et al. 2022). Alternatively, dust could have been consistently retained in heavily gas-depleted discs, in apparent conflict with the evolution of the cumulative mass distribution (Fig. 6). Finally, our models could be consistent with the observations if the estimated mass flux onto the host star is somehow underestimated from the accretion observations. J. Appelgren , M. Lambrechts and N. van

Population synthesis with α ν " 10´3
Lowering the viscous α ν has a number of effects on the evolution of a protoplanetary disc. In general, the disc will evolve more slowly as angular momentum transport is less efficient. However, when the disc is very young, the increased viscosity due to gravitational instabilities can dominate over the base viscosity, making early evolution somewhat similar for different base α ν values (Eq. 4). The slower viscous evolution results in photoevaporation becoming dominant at later times and therefore longer gas disc lifetimes. Viscous heating is also less efficient, leading to a cooler disc, further slowing down the evolution.
The effect on the dust component is less straightforward. A slower gas evolution results in generally higher gas surface densities. This, in turn, results in lower Stokes numbers and slower dust drift (Eq. 6 and 8). However, the particle size in the fragmentation limit is linearly dependant on the gas surface density (Eq. 12). Therefore, Stokes numbers do not change significantly (for a fixed value of α t ). However, the dust is also advected with the gas. At α ν " 10´3 the radial gas velocity is slower than at α ν " 10´2, resulting in less dust being removed from the disc by pure advection with the gas. The dust disc is therefore effectively drained at a lower rate at α ν " 10´3 compared to α ν " 10´2.
We do not find that the retention of dust by photoevaporative gaps is significantly different at α ν " 10´3 compared to α ν " 10´2 (Fig. 8). The dust discs are longer lived by a few million years at α ν " 10´3 compared to α ν " 10´2, the gas discs also evolve slower and a photoevaporative gap opens up typically 12 Myr later. The fraction of discs which retain ą 0.1 M C of dust outside the photoevaporative gap is 20% at α ν " 10´2 and 15% at α ν " 10´3. Using α ν " 10´2 the median dust mass retained outside the photoevaporative gap is 0.3 M C , compared to 0.2 M C when using α ν " 10´3.
The observed cumulative distributions (right panel of Fig. 8) were (as before) scaled according to Eq. B.1, but with a disc lifetime of τ " 8 Myr, following (Michel et al. 2021). This lifetime is derived by excluding regions with high external photoevaporation. It should therefore be more representative of the observational sample we are comparing to, which are not strongly exposed to external photoevaporation (Michel et al. 2021).
For α ν " 10´3 and τ " 8 Myr, radial drift is overall consistent with the observed depletion of dust, as seen in the two panels of Fig. 8. Also, here we find that radial drift alone is not able to reproduce a small fraction of the population of discs in the observed sample, which retain large amounts of dust for more than several Myr -similar to what we found when assuming α ν " 10´2. This further supports the idea that there are some discs where radial drift may be less efficient.
More slowly evolving discs lead to lower accretion rates onto host stars, which are manifested in the 9 M g´Mdust relation. However, we find that the combined evolution of both quantities is scarcely different from the α ν " 10´2 case. There is a slightly larger scatter in gas accretion rates and a larger scatter in the disc dust mass at the end of the disc formation, as shown in Fig. 9.

Dependency on the stellar mass
The fraction of the dust disc mass removed in a given time varies with disc and stellar mass. The dependency of the solid mass budget on stellar mass, for different cluster ages, is explicitly shown in Fig. 10. We find that dust clearing by radial drift is more efficient around low-mass stars, as also found in previous works (Appelgren et al. 2020;Pinilla et al. 2020 Fig. 9: Gas accretion rate, assuming α ν " 10´3, as a function of dust disc mass for the models using the Pic21 photoevaporation prescription (black) and the Kom21 photoevaporation prescription (red).
vations point towards a linear or super-linear trend of dust mass with stellar mass, possibly steepening with time (Pascucci et al. 2016). We find steeper power law profiles, given by M dust 9M β ‹ with β " 2.4 at 2 Myr and β " 6.9 at 7.5 Myr, although in the latter case the data is poorly described by a single power law. When limiting, instead, the model sample to those discs with detectable masses above 0.1 M C , we can recover much more moderate power law dependencies, with β " 2.2 at 2 Myr and β " 1.6 at 7.5 Myr, in line with the observed exponents. These relations may further be influenced by the stellar-mass dependent occurrence of giant planets that can reduce pebble drift, as recently argued in Pinilla et al. (2020); van der Marel & Mulders (2021). This dependency on stellar mass is also reflected in the evolution over time of the profile of the cumulative disc dust mass distribution (Figs. 6 and 8). The cumulative distribution becomes less steep with time. If dust removal was more efficient around more massive stars with more massive discs, the profile would become steeper with time. Conversely, if the efficiency is the same across all disc and stellar masses, the profile would appear similar in shape to the gas tracer profiles in Figs. 6 and 8. In summary, radial drift not only removes dust on a timescale that agrees with observations, but it is also able to match the profile of the cumulative distribution well.

Role of pressure bumps
Disc surveys have revealed that many of the resolved protoplanetary discs show significant substructures, with very pronounced rings and gaps in the dust (e.g. Andrews et al. 2018). These structures are thought to originate from pressure bumps trapping dust (Dullemond et al. 2018;Pinilla et al. 2020Pinilla et al. , 2021. The high occurrence of substructures in observed discs suggests that they might be common in all discs. However, by design, such surveys often target the largest and brightest protoplanetary discs, which could be more prone to forming gap-opening giant planets (van der Marel & Mulders 2021).
The frequent occurrence of deep gaps and rings in the dust continuum of these large and long-lived discs may signal regions where pebble drift is delayed. The link between exoplanet statistics and disc substructures was investigated in detail by van der  Fig. 10: Dust mass as a function of stellar mass for the population synthesis model using α ν " 10´3 and the Kom21 photoevaporation prescription. In discs around lower mass stars, the dust component is removed faster than in discs around higher mass stars.
Marel & Mulders (2021). They found that the increased occurrence of extended, structured discs with stellar mass is similar to the increased occurrence of giant planets with stellar mass. This supports the idea that large scale gaps in protoplanetary discs originate from giant planets. Such pressure bumps can arise when growing planetary cores reach the pebble isolation mass and carve shallow gaps that are sufficient to halt pebble drift (Lambrechts & Johansen 2014;Rosotti et al. 2016;Bitsch et al. 2018). Other sources of pressure bumps that can slow down or halt pebble drift also exist, for example, changes in the level of turbulence at the transition from regions where magnetohydrodynamic instabilities are active (Lesur et al. 2022). It is thus plausible that the majority of discs are radial-drift dominated, but that in " 5%´15% of discs, radial drift is reduced by the presence of giant planets. This may be connected to the frequently inferred presence of giant planets interior to debris discs (Pearce et al. 2022), where the latter could have been the natural outcome of left-over pebbles clumping into planetesimals during disc clearing .

Role of planet formation
Radial drift is not the sole process by which pebbles are removed from protoplanetary discs, as we know that a fraction of them will form the planets and planetesimals that we observe in our own solar system, as well as in exoplanet systems and debris discs. In this work, we show that radial drift removes dust from protoplanetary discs on a timescale that agrees well with dust disc lifetimes measured from observations. This suggests that radial drift is the dominant dust-removal process. However, since we do not include planetesimal or planet formation in the model, we cannot dismiss their potential importance.
Planetesimal formation could remove a substantial fraction of dust particles present in protoplanetary discs, but its efficiency, as well as how that might vary across stellar and disc mass, is currently not well understood. Planetesimals are thought to form by the streaming instability (Youdin & Goodman 2005;Johansen et al. 2007), which is likely to happen in regions with a supersolar dust-to-gas ratio (Carrera et al. 2015;Yang et al. 2017;Li & Youdin 2021). Pressure bumps at the gap edges of giant planets are a promising location for this process, as dust piles up in these regions and large amounts of planetesimals could potentially form (Eriksson et al. 2020). Such planetesimals would then represent a second generation of late-formed planetesimals, as they are triggered by the presence of an already fully-formed giant planet.
Similarly, the formation of planets should also be responsible for the removal of some fraction of the dust in protoplanetary discs. To assess the importance of this process, a possible future extension of this work could include a exoplanet synthesis model as well (similarly to the study of e.g. Liu et al. 2019). In this context, it is interesting to note that the mass of solids present in typical observed exoplanet systems is comparable to that found in class II discs (Mulders et al. 2021). In the earlier disc stages, the mass budget in pebbles is larger by more than one order of magnitude, so up until the class II stage the decrease in the mass in pebbles is mainly driven by dust drift, rather than planet formation.

What happens to the dust retained outside the photoevaporative cavity
In our model, the dust retained outside the photoevaporative cavity eventually moves into the location of the pressure maxima where there is no radial drift. Since our model does not consider any additional treatment of this dust, this results in narrow longlived rings of dust that are slowly pushed outwards as the photoevaporative gap in the gas expands outwards. The high dust-togas ratio in these pressure bumps could make them favourable locations for the formation of planetesimals via the streaming instability (Carrera et al. 2021;Carrera & Simon 2022). However, the fate of the dust retained outside photoevaporative gaps is uncertain. For example, Owen & Kollmeier (2019) found that this dust is effectively lost in a process where the largest particles grow until they fragment by collisions, replenishing the reservoir of small dust particles that are subsequently removed by radiation pressure. If the dust is turned into planetesimals, this could lead to the formation of debris discs. However, the median mass of dust outside the photoevaporative gaps is only 0.3 M C for α ν " 10´2 and 0.2 M C for α ν " 10´3. While this dust mass is similar to that present in the known debris disc population, the planetesimal population that creates this dust in debris discs could be significantly higher (Holland et al. 2017;Krivov & Wyatt 2021).

Role of stellar binarity
Most solar mass stars reside in multiple systems, and the stellar multiplicity fraction increases with stellar mass (Offner et al. 2022). Observations have found that discs around binaries typically host dust masses below 50 M C , and that for binary separations below 100 au, discs have dust masses below 10 M C (Zurlo et al. 2020). Therefore, measuring the mean disc masses in young stellar clusters without taking into account the binarity of stars can lead to underestimations of the typical masses of Gas disc radius as a function of dust disc radius, using α ν " 10´3 and the Pic21 photoevaporation prescription for different cluster ages. Triangles show discs which have retained more than 0.1 M C of dust, but whose gas accretion rates are lower than 10´1 2 M d /yr. These are the discs, corresponding to 15% of the initial disc population, in which photoevaporation has opened up a gap, with a detectable amount of dust is retained outside the gap. Grey markers show observational measurements from Long et al. (2022). Grey lines show the R gas " R dust (solid) and R gas " 4R dust (dashed).
discs around single stars. A future extension of this work could include the effect of stellar binarity in driving faster disc dispersal (Kraus et al. 2012) through tidal effects and stellar encounters truncating disc radii in the cluster stage(Portegies Zwart 2016)

Transition discs
Transition discs with dust-depleted inner cavities and low accretion rates are suggested to to have been caused by the photoevaporation of the inner disc (Owen et al. 2011;van der Marel 2022). However, we find that it is difficult to form discs that both show signs of gas accretion onto the hosts star and massive extended dust disc at wider radii. Gas within the photoevaporation gap is lost on very short timescales (À 10 5 yr), as seen in Fig. 4, while only retaining small amounts of dust trapped at the outer gap edge.
As discussed earlier in ths paper, the fate of the retained dust is such that it might be lost on short timescales (Sect. 4.1.4). Therefore we do not find that photoevaporation, as formulated here, is a likely pathway for forming the majority of transition discs. This is in line with the low occurrence of discs with a cavity radius within 10 au (van der Marel et al. 2022) that would have been expected in classic photoevaporation models (Picogna et al. 2019).

Evolution of the outer gas and dust disc radius
Pebble drift could also, in principle, be identified through tracing the evolution of the outer dust disc radius with respect to the outer gas disc radius (Trapman et al. 2019;Toci et al. 2021). This is nevertheless challenging as both radii could be influenced by late accretion events of gas and dust onto the disc (Kuffmeier et al. 2017). Moreover, in our model, the outer gas edge continuously viscously expands outwards with time, while recent models have argued against this based on disc wind models (Manara et al. 2022). However, an outwardly expanding discs is seen in disc wind models that include the truncation of the outer disc (Yang & Bai 2021). Figure 11 illustrates the evolution of the gas radius as a function of the dust radius for the population synthesis model using α ν " 10´3 and the Pic21 photoevaporation prescription. Here, we calculate the disc radius as the radius within which either 90% of the disc gas or dust mass is present. If a disc has opened up a gap by photoevaporation and retained > 0.1 M C of dust outside the gap, the dust radius is taken as the outer radius where the dust is trapped and indicated with a triangle symbol.
For very young discs, the gas and dust radii are similar because the dust is tightly coupled to the gas. As the disc evolves, the dust decouples from the gas and the dust radius begins to decrease, while the gas radius still increases. In the approximation used here, which consists of only tracing the largest particle size, the radius of the dust disc shrinks rapidly once radial drift becomes efficient while gas radii continue to expand due to viscous evolution. Finally, the triangular symbols represent those discs where photoevaporation has opened up a gap in the gas, whilst retaining ą 0.1 M C of dust outside the gap. The dust disc radii also increase, as the dust is pushed outwards as the outer edge of photoevaporative gap expands.
A direct comparison of synthetic with observed dust and gas radii is not straightforward. This would require post-processing our results with radiative transfer in order to determine the observable outer dust edge and gas edge, the latter using a CO gas tracer (Toci et al. 2021). Nevertheless, for illustrative purposes, we show in Fig. 11 (in grey points) the dust and gas (CO) radii determined for the discs in the sample by Long et al. (2022). This uniformly analysed collection of discs, drawn mainly from Taurus, Lupus, Ophiuchus, and Upper Sco, spans stellar masses between 0.15-2 M d and ages between 0.5 to 20 Myr. Similarly to our results, all observed discs have dust radii located within their gas disc radii and only a few discs have severely drift-depleted discs with R gas ą 4R dust . Previous works have argued that radial drift would rapidly deplete pebble discs and result in discs with a gas-to-dust size ratio smaller than four (Toci et al. 2021). However, this is likely an effect of assuming larger fragmentation velocities (u f " 10 m/s), compared to the parameters used in the current study (u f " 1 m/s, see Sect. 2.3.2). This results in larger pebbles closer to the drift limit and an outer dust radius that moves inwards too rapidly.

Uncertainty in photoevaporation prescriptions
The disagreement in derived mass-loss rates in different photoevaporation models reflects the fact that the underlying physics driving photoevaporation is not entirely understood. Sellek et al. (2022) found that the chosen X-ray spectrum frequency matters significantly for the derived photoevaporation rates. This can explain why some studies (e.g. Wang & Goodman 2017) find very low mass loss rates for for the derived photoevaporation rates X-rays. Also, X-ray photoevaporation models, such as the one from Pic21 used here, may underestimate the cooling efficiency in the discs and therefore overestimate the photoevaporation rate Wang & Goodman (2017); Sellek et al. (2022). Since we cover both low and high photoevaporation efficiencies here, we believe our conclusion on their relative small effect on dust evolution are robust. Nevertheless, this remains a clear area for future work, especially as recent works point towards photoevaporation and magnetic disc wind losses to be tightly linked (Bai et al. 2016).
Since our model does not track the population of small (À 100 µm) dust grains and the one from Kom21 provides no scaling with the dust-to-gas ratio (a ratio of 0.01 is used), we did not model the effect that a decreasing or increasing amount of small dust present has on the FUV photoevaporation rate. A decrease in the available small dust grains reduces the opacity and thus increases the penetration depth of the photons, which can enhance the mass-loss rate. Simultaneously, fewer small dust grains also leads to less photoelectric heating, which reduces the mass-loss rate. Gorti et al. (2015) found that these two effects can to a large degree negate each other. However, at dust-to-gas ratios above " 0.03, Nakatani et al. (2018b) found that a decrease in the dust-to-gas ratio increases the FUV photoevaporation rate. At lower dust-to-gas ratios, the photoevaporation rate instead decreases, but the inclusion of hydrogen ionisation by X-rays lowers the reduction in the photoevaporation rate due to the large amount of free electrons allowing for rapid recombination onto the dust grains. Similarly, Nakatani et al. (2021) also concluded that at dust-to-gas ratios below 0.01 photoevaporation by FUV is reduced due to inefficient photoelectric heating. With these results in mind, coupled with the fact that the discs in our model rarely exceed the initial dust-to-gas ratio of 0.01 (see middle panel of Fig. 3), including the effect of small grains of the photoevaporation rate would likely result in a reduced FUV photoevaporation rate compared to our current model. Since we already find that FUV photoevaporation has a negligible effect on the dust disc evolution, our conclusions would not change.
We also did not take into account that small dust particles may also be entrained in photoevaporative winds (Owen et al. 2011;Hutchison et al. 2016b,a), but this only affects very small particles of dust, as the maximum particle size that can be effectively entrained in photoevaporative winds is " 1 µm (Hutchison & Clarke 2021;Booth & Clarke 2021). Therefore, this process does not impact significantly the mass budget of pebbles in the disc, but may be important for the determination of gas mass loss rates in the upper parts of protoplanetary discs Wang & Goodman (2017).
Finally, in this project, we did not include external photoevaporation. This choice is appropriate for the well studied nearby stellar clusters considered here. Therefore, we did not include it in our model. Nevertheless, protoplanetary discs in more massive clusters with O and B stars would be subject to external photoevaporation from the strong FUV background radiation of said stars. This process may significantly shorten the disc lifetimes (Winter & Haworth 2022).

Conclusions
1. We found that the depletion of protoplanetary dust by radial drift of pebbles is compatible with the depletion trend seen in observed protoplanetary discs. We explored dust depletion in synthetic clusters discs with disc depletion timescales of 2.5 Myr or 8 Myr, corresponding to, respectively, a high and low viscous alpha (α ν " 10´2 and α ν " 10´3). Both cases broadly reproduce the cumulative loss off pebbles from protoplanetary discs.
2. Our synthesis model is consistent with the evolution of the cumulative dust disc mass distribution with time. However, a fraction of about 5-10% of observed discs show signs of reduced radial drift compared to the model, possibly triggered by wide-orbit planet formation. Other relations that are reproduced well are the moderate decrease of dust radius with gas radius with time and the trend between disc mass and stellar mass. However, the synthesis model does not recover a population of massive dust discs with low stellar accretion rates. The latter may point to our model not including inner disc winds, which could reduce the mass accretion rate onto the host star, or the possibility that the observations underestimate stellar accretion rates. 3. Discs forming with a global dust-to-gas ratio of 0.01 do not retain this proportion of dust. Efficient radial drift, even during the build-up of the discs, lowers the disc's global dust-togas ratio below that of its birth environment. Hence, gas disc masses cannot be estimated by multiplying the dust mass by the ISM-ratio of 100, even without considering dust optical depth effects. 4. Photoevaporation plays a minor role on the evolution of the disc dust mass. Strong X-ray photoevaporation is able to retain some dust outside the photoevaporative gaps, but only a very small fraction (À0.01 %) of the initial dust mass of the disc (50´1000 M C ). Photoevaporation may aid in reaching the low gas accretion rates observed (À 10´9 M d /yr) in protoplanetary discs (Sellek et al. 2020). However, this depends on the assumed photoevaporation efficiency, and is only achieved for maximally efficient X-ray photoevaporation prescriptions. The photoevaporation prescription of Kom21 is based on models which includes photoevaporation due to EUV, FUV, and X-rays. The rate of photoevaporation is given by Eq. (A.3): log˜9 Σ g,PE prq 1 g s´1 cm´2¸" c 5 x 5`c 4 x 4`c 3 x 3`c 2 x 2`c 1 x`c 0 , where x is given by x " logˆr r g˙, (A.4) and r g is the the gravitational radius given by: The photoevaporation prescription of Kom21 is in the range 0.1r g ď r ď 20r g , as this is where photoevaporation is thought to occur (Liffman 2003). We therefore only include photoevaporation in the same range. The parameters c 5,...,0 in Eq. (A.3) for different stellar masses are given in Table A.3. X-ray and FUV luminosities for the same stellar masses are given in Table A.4. To approximate the evolution of the photoevaporation rate as the star grows more massive we linearly interpolate between these stellar masses and parameter values.

Appendix B: Disc fractions
The dust present in protoplanetary discs will result in an excess of infrared emission at mm wavelengths. As a disc in a cluster of young stars evolves, the dust component is depleted and the infrared excess is lost. To determine the fraction of stars that have discs with mm emission, often only stars with discs that show an infrared excess are used as targets. A number of attempts have been made at fitting a function describing how the infrared disc fraction evolves over time given a typical disc lifetime (e.g. Haisch et al. 2001;Mamajek 2009;Ribas et al. 2014;Michel et al. 2021). The fitting function for the disc infrared fraction often used is where τ is the typical disc lifetime, and t is the age of the cluster. Disc lifetimes were long thought to be a few Myr, for example, τ " 2.5 Myr Mamajek (2009). However, these lifetimes were derived using samples of young star clusters where external disc photoevaporation is both present and not. External photoevaporation can significantly shorten disc lifetimes. A more recent