Open Access
Issue
A&A
Volume 673, May 2023
Article Number A139
Number of page(s) 15
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202245252
Published online 23 May 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Surveys of young star-forming regions, such as the 1 Myr-old Perseus cluster, reveal typical dust disc masses that exceed 100 M (Tychoniec et al. 2020; Tobin et al. 2020). Older star-forming 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, 2018; Pascucci et al. 2016; Eisner et al. 2018; Ruíz-Rodríguez 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 photoevapora-tion. 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, 2015; Tazzari 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, 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, cra-tering, 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 mass dust (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 planetesi-mals via the streaming instability (Youdin & Goodman 2005; Johansen et al. 2007; Carrera et al. 2015; Yang et al. 2017; Li et al. 2021) and the efficiency of pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014).

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 yr−1, for the youngest discs, to ≲ 10−9 M yr−1 for discs as old as a few million years (Manara et al. 2016). Simultaneously, the gas disc mass is reduced from ~0.1 M to ~0.001 M (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 , 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 photoe-vaporation, which was extended by an updated X-ray photo-evaporation prescription (Picogna et al. 2019, 2021; Ercolano 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).

2 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.

2.1 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) and 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: (1)

where rcf is the radius of the collapse front, ρ(rcf) is the density of the cloud core at this radius and drcf /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: (2)

The infalling material lands on the disc within the centrifugal radius Rc, given by: (3)

where Ω0 is the rotation rate of the cloud core and M(rc) 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.

2.2 Gas disc evolution

We evolved the disc viscously using the well-established α-disc prescription of Shakura & Sunyaev (1973), where the viscosity of the disc is described by . Here, cs 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 mid-plane 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. 2020, 2022). Instead, we parame-terise 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 line-width broadening by turbulent motions in outer disc regions have found values of α ≲ 10−3 (Teague et al. 2018; Flaherty et al. 2018, 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 recentmea-surements of the disc fraction finds that, in those star-forming 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 αv, of either 10−2 or 10−3 to explore both scenarios of faster and slower disc evolution. In practise, α is calculated as: (4)

which increases the disc viscosity when the disc becomes grav-itationally 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: (5)

where is the mass loss due to internal photoevapora-tion. The way we calculate this is explained in Sect. 2.5 and Appendix A.

2.3 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.

2.3.1 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: (6)

Here, υk is the Keplerian orbital velocity and η is a measure of the radial pressure support in the disc, given by: (7)

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: (8)

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, ap, 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, υd, is given by: (9)

where υg is the radial velocity of the gas. The surface density of dust is then evolved with the following continuity equation: (10)

2.3.2 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: (11)

For grains that reach their maximum size limit by fragmenting, the maximum grain size is set via (12)

where ff = 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) and Villenave et al. (2022). We note that we treat it as a separate parameter from αv (Eq. (4)). The fragmentation threshold velocity, uf, depends on the material properties of the grains. Laboratory experiments of silicate dust grain collisions suggest limits of about 1 m s−1 (Blum & Wurm 2008). Water ice was thought to be stickier than silicates and to have a fragmentation speed of about 10 m s−1 (Gundlach & Blum 2015). However, more recent studies have concluded that the fragmentation speed of water ice is likely to be lower. The stickiness of CO2-ice is also thought to be quite low, with fragmentation speeds of about 1ms−1 (Musiolik et al. 2016a). Similar results are found for mixtures of water and CO2 ices (Musiolik et al. 2016b). Therefore, the choice of 1 m s−1 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: (13)

where fd = 0.55 is a fudge factor with similar purpose as ff and (Birnstiel et al. 2012).

2.4 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): (14)

where L is the stellar luminosity. We assume that the luminosity of the star scales linearly as L = L(M/M, Baraffe et al. 2002; Liu et al. 2019).

The heating rate due to irradiation is given by (e.g. Menou & Goodman 2004): (15)

where Tcore = 10 K is the temperature of the cloud core, τ = κΣ/2 Is the optical depth, and k Is the opacity.

The heating rate due to viscous heating is given by (e.g. Kimura et al. 2016): (16)

Heat is lost from the disc by radiative cooling (e.g. Menou & Goodman 2004): (17)

where Tmid is the midplane temperature of the disc. The mid-plane 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 cm2 g−1 (T/K)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; Li et al. 2021). We therefore imposed a temperature limit of 1500 K in our model.

thumbnail Fig. 1

Photoevaporation rates for a 1 M-star using identical X-ray luminosity for both the Pic21 (blue) and Kom21 (yellow) prescriptions. The two dashed lines show the Kom21 models where the FUV luminosity has been reduced by a factor of 10 and of 100.

2.5 Photoevaporation

Photoevaporation removes gas from the surface of the protoplan-etary disc by either EUV, FUV, or X-rays heating the gas until it can escape the disc (Gorti et al. 2009, 2015; Owen et al. 2010, 2011, 2012; Wang & Goodman 2017; Nakatani et al. 2018a,b). In this work, we look at two different photoevaporation prescriptions: one by Picogna et al. (2021, hereafter Pic21) and one by Komaki et al. (2021, hereafter Kom21).

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).

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 monoton-ically increasing or decreasing with stellar mass. For the Pic21 model, this upper mass is 1M and stars above this mass will prove very uncommon in our sample; therefore, the choice of using the 1M cut off will be of little consequence. With the Kom21 model the lower mass is 0.3 M and thus more stars will be affected since our stellar masses extend to ≲0.1 M, possibly leading to moderately overestimating the photoevaporation rate for very low-mass stars.

2.6 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 core-mass-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 and 2 M.

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: (18)

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 Tcore = 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.5Myr and 8 Myr (Mamajek 2009; Michel et al. 2021), respectively.

thumbnail Fig. 2

Temperature (left), gas surface density (middle), and dust surface density (right) evolution for a disc forming from a 1 M 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 of the dust is retained outside this gap and is pushed outwards with time.

3 Results

3.1 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 1M 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.7Myr, 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–4au. 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.6MEarth 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.35Myr) dust also has an outwards motion in the outer disc region between 10 and 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 and 5 au, and at a later time, between 0.6 and 0.9 Myr, at radii between 0.1 and 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 temperature of 1500 K to lower temperatures at wider radii, where stellar irradiation is dominant in setting the disc temperature rather than viscous heating.

3.2 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 αv = 10−3 in Sect. 3.3.

thumbnail 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 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 and 5 au, and one at 0.6–0.9 Myr at radii from 0.1 to 1.5 au. The particle size generally is larger in the inner regions of the disc and at earlier times.

thumbnail Fig. 4

Gas accretion rate as a function of disc age for all discs in both photoevaporation prescriptions. Once photoevaporation opens up a gap in the gas disc, the gas within disc gap is quickly cleared by photoevapo-ration 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).

3.2.1 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 to 0.5 M. These gas discs are then depleted by accretion and photoevaporation (Fig. 5). For both the Kom21 and Pic21 pho-toevaporation 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.

3.2.2 Dust disc evolution

The radial drift of pebbles reduces the dust mass of protoplan-etary 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 and ~1 M, and are born with dust disc masses between of 50 M and 1000 M respectively. The least massive and shortest lived discs are depleted of dust within less than 1 Myr, whilst the most massive and longest 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.5Myr-old cluster, 73% of discs have a disc dust mass higher than 10 M. At 1 Myr, this has decreased to only 45% of discs having over 10 M 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 of dust. With the (Picogna et al. 2021) model, we find that 20% of all discs are able to retain at least 0.1 M of dust. However, among these discs, the median mass retained is only 0.3 M. The median mass of these discs right at the end of their formation, which is when they are the most massive, is 367 M. Hence, these discs retain less than 0.01% of their maximum dust mass. For comparison, the observed median dust masses in protoplanetary disc surveys at ~2Myr are 3.0 M (Lupus), 3.2 M (Cham I), 6.2 M (Taurus), and 11.5 M (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 yr−1) 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 yr−1, 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. The observed sample is taken from Ansdell et al. (2016); Barenfeld et al. (2016); Pascucci et al. (2016); Tychoniec et al. (2020); Villenave et al. (2021) and 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 observable 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 (0.5 Myr) and ~5.5 M (1 Myr). This can be compared to the Perseus sample, which has significantly higher median masses for class 0 objects (158 M) and class I objects (52 M, 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 and 3 Myr all lie in between the 1 and 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 Kom2l photoevaporation model lies somewhat close. Using the Pic2l 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 (≲105 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) 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 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 Yr−1 and dust disc masses of over 1 M Is populated by the observed sample. Our synthesis model with the Pic2l 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 Kom2l photoevaporation prescription, an even smaller region of the observed sample is reached. Only those discs with accretion rates of about 10−8 M Yr−1 or more are reproduced.

The disagreement between model and observation in the gas accretion rate and dust disc relation could originate from an over-estimation 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.

thumbnail 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 αv = 10−2. Gas discs are formed with masses between 0.03 M to 0.5 M. After 10 Myr of evolution about 1% of this initial mass remains. Dust disc masses are initially between 10 M to 1000 M. 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 Pic2l or Kom2l 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.

thumbnail Fig. 6

Cumulative disc dust mass distributions for the model (left plot) using using α = 10−2. Here we show the Pic2l (solid) and the Kom2l (dashed) photoevaporation prescriptions. A sample of observed disc masses are shown in the right plot. The two photoevaporation prescriptions are very similar at 0.5 and 1 Myr and begin to differ at the low-mass end at 2 Myr. At 4–5 Myr the Pic2l retains dust in some discs, whereas the Kom2l has very little left. We note that the 3, 4, and 5 Myr lines of the Pic2l are on top of each other. The thin coloured lines in the left plot show the gas mass multiplied by the initial dust-to-gas ratio. This shows how the disc dust mass would evolve without any radial drift. The observed sample is colour coded according to the cluster’s oldest inferred age.

thumbnail Fig. 7

Gas accretion rate as a function of dust disc mass for the models using the Pic2l (black) and the Kom2l (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 yr−1). 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.

3.3 Population synthesis with αv = 10−3

Lowering the viscous αv 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 αv 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 (Eqs. (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 αv = 10−3 the radial gas velocity is slower than at αv = 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 αv = 10−3 compared to αv = 10−2.

We do not find that the retention of dust by photoevapo-rative gaps is significantly different at αv = 10−3 compared to αv = 10−2 (Fig. 8). The dust discs are longer lived by a few million years at αv = 10−3 compared to αv = 10−2, the gas discs also evolve slower and a photoevaporative gap opens up typically 1–2 Myr later. The fraction of discs which retain >0.1 M 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, compared to 0.2 M⊕ when using αν = 10−3.

The observed cumulative distributions (right panel of Fig. 8) were (as before) scaled according to (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 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.

thumbnail Fig. 8

Cumulative distributions of disc dust masses from the population synthesis model using α = 10−3 (left) and observed discs (right). The observed sample is taken from Tychoniec et al. (2020, Perseus) and Manara et al. (2022). The distributions of the observed sample are scaled according to the disc fractions from Michel et al. (2021), taking age of each region as the mid-point of the given age range.

4 Discussion

4.1 Interpretation of disc dust masses

4.1.1 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). Observations 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 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, 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) and 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.

thumbnail 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).

thumbnail 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.

4.1.2 Role of pressure bumps

Disc surveys have revealed that many of the resolved protoplan-etary 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. 2020, 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 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 magneto-hydrodynamic 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 (Carrera et al. 2017).

4.1.3 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 super-solar 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 planetesi-mals 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.

4.1.4 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 long-lived rings of dust that are slowly pushed outwards as the photoe-vaporative gap in the gas expands outwards. The high dust-to-gas 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 for αy = 10−2 and 0.2 M for αy = 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).

4.1.5 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, and that for binary separations below 100 au, discs have dust masses below 10 M (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 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).

4.2 Transition discs

Transition discs with dust-depleted inner cavities and low accretion rates are suggested to to have been caused by the photoevap-oration of the inner disc (Owen et al. 2011; van der Marel 2023). 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 (≲105 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).

4.3 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 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 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 and ages between 0.5 and 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 Rgas > 4Rdust. 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 (uf = 10 m s−1), compared to the parameters used in the current study (uf = 1 m s−1, 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.

thumbnail Fig. 11

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 of dust, but whose gas accretion rates are lower than 10−12 M Yr−1. 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 Rgas = Rdust (solid) and Rgas = 4Rdust (dashed).

4.4 Uncertainty in photoevaporation prescriptions

The disagreement in derived mass-loss rates in different photo-evaporation 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 photoevapo-ration 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 photoevap-oration rate Wang & Goodman (2017) and 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 photo-evaporation 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-togas ratios above ~0.03, Nakatani et al. (2018b) found that a decrease in the dust-to-gas ratio increases the FUV photoevap-oration 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 photoevapo-ration 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 photoe-vaporation. 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 photo-evaporation from the strong FUV background radiation of said stars. This process may significantly shorten the disc lifetimes (Winter & Haworth 2022).

5 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-to-gas 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). Photoevaporation may aid in reaching the low gas accretion rates observed (≲10−9 M yr−1) in pro-toplanetary discs (Sellek et al. 2020). However, this depends on the assumed photoevaporation efficiency, and is only achieved for maximally efficient X-ray photoevaporation prescriptions.

Acknowledgements

We thank Anders Johansen for his comments on the manuscript. We also thank Andrew Sellek for helpful discussions and input on the photoevaporation models. We also thank the anonymous referee for their comments that helped improve this paper. J.A. acknowledges the Swedish Research Council grant (2018-04867, PI A. Johansen). M.L. acknowledges the Wallenberg Academy Fellow grant (2017.0287, PI A. Johansen) and ERC starting grant 101041466-EXODOSS.

Appendix A Photoevaporation prescriptions

Appendix A.1 Picogna et al. (2021) photoevaporation

Pic21 provides a photoevaporation prescription derived from models of X-ray photoevaporation for different stellar masses using the photoevaporation model of Picogna et al. (2019); Ercolano et al. (2021).

Table A.1

Parameters of Eq. (A.1) and Eq. (A.2) for different stellar masses.

The photoevaporation rate is given by (A.1)

where (r) is the mass loss rate due to photoevaporation at a given disc radius. It is given by (A.2)

where (LX) is the total mass loss rate due to X-ray photoevaporation. The parameters a,…, g for different stellar masses are given by Table A.1. The X-ray luminosities for each stellar mass is shown in Table A.2.

Table A.2

Stellar masses and corresponding X-ray luminosities used in the Pic21 photoevaporation model.

Appendix A.2 Komaki et al. (2021) photoevaporation

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): (A.3)

where x is given by (A.4)

and rg is the gravitational radius given by: (A.5)

The photoevaporation prescription of Kom21 is in the range 0.1rgr ≤ 20rg, as this is where photoevaporation is thought to occur (Liffman 2003). We therefore only include photoevapo-ration in the same range.

Table A.3

Parameters of Eq. (A.3) for different stellar masses.

The parameters c5,…,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.

Table A.4

Stellar masses and corresponding FUV and X-ray luminosities used in the Kom21 photoevaporation model.

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 (B.1)

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 photoevap-oration can significantly shorten disc lifetimes. A more recent study of clusters without external photoevaporation found that discs in such environments have typical lifetimes of about 8 Myr (Michel et al. 2021).

The ages of the star-forming regions in the observational sample we made our comparison to were taken from Manara et al. (2022). We assumed that the of a cluster age is the midpoint of the given age range. These ages are then used in Eq. B.1 to calculate the disc fraction of a star forming region. The resulting disc fraction is then used to scale the cumulative distribution functions of dust disc masses.

Appendix C Halting dust drift in high-mass discs

Here, we illustrate how the halting of dust transport inwards in the disc, in a fraction of the most massive discs, can affect the cumulative dust disc mass distributions. Stopping (or slowing down) the inwards transport of dust could be caused by pressure bumps in the gas disc, which are discussed in Sect. 4.1.2. In order to approximate dust trapping in pressure bumps we halted the pebble drift between the end of the disc formation stage and 2 Myr in 17% of the 30% most massive discs (corresponding to 5% of the total disc population). We show our results in Fig. C.1 for a cluster with αν = 10−2 and using the Kom21 photoevaporation prescription.

Halting the dust evolution in these 5% of discs naturally results in a cumulative distribution of the dust disc mass that extends to higher masses, as seen in Fig. C.1. The thick lines shows the simulation including 5% of discs with halted dust evolution, and the thin lines shows the standard model. This experiment illustrates how halting pebble drift in a fraction of the most massive discs can help explain the presence of very massive dust discs seen in evolved star forming regions, whilst the majority of discs agree well with being drift-dominated.

thumbnail Fig. C.1

Cumulative dust mass distributions of the simulations using αν = 10−2 and the Kom21 photoevaporation prescription with pebble drift (thin lines) and with pebble drift halted in 5% of discs (thick lines). Halting pebble drift extends the high-mass tail of the CDF, bringing it more in line with what is seen in observations.

References

  1. Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129 [Google Scholar]
  2. Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46 [Google Scholar]
  4. Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240 [Google Scholar]
  5. Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
  6. Appelgren, J., Lambrechts, M., & Johansen, A. 2020, A&A, 638, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76 [Google Scholar]
  8. Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [Google Scholar]
  9. Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
  10. Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563 [CrossRef] [EDP Sciences] [Google Scholar]
  11. Barenfeld, S. A., Carpenter, J. M., Ricci, L., & Isella, A. 2016, ApJ, 827, 142 [Google Scholar]
  12. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bitsch, B., Morbidelli, A., Johansen, A., et al. 2018, A&A, 612, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Blum, J. 2018, Space Sci. Rev., 214, 52 [NASA ADS] [CrossRef] [Google Scholar]
  16. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [Google Scholar]
  17. Booth, R. A., & Clarke, C. J. 2021, MNRAS, 502, 1569 [NASA ADS] [CrossRef] [Google Scholar]
  18. Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71 [Google Scholar]
  20. Carrera, D., & Simon, J. B. 2022, ApJ, 933, L10 [NASA ADS] [CrossRef] [Google Scholar]
  21. Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [Google Scholar]
  23. Carrera, D., Simon, J. B., Li, R., Kretke, K. A., & Klahr, H. 2021, AJ, 161, 96 [Google Scholar]
  24. Cazzoletti, P., Manara, C. F., Liu, H. B., et al. 2019, A&A, 626, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
  26. Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698 [Google Scholar]
  27. Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 [NASA ADS] [CrossRef] [Google Scholar]
  28. Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46 [NASA ADS] [CrossRef] [Google Scholar]
  29. Eisner, J. A., Arce, H. G., Ballering, N. P., et al. 2018, ApJ, 860, 77 [CrossRef] [Google Scholar]
  30. Ercolano, B., Picogna, G., Monsch, K., Drake, J. J., & Preibisch, T. 2021, MNRAS, 508, 1675 [NASA ADS] [CrossRef] [Google Scholar]
  31. Eriksson, L. E. J., Johansen, A., & Liu, B. 2020, A&A, 635, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Feng, Y., & Krumholz, M. R. 2014, Nature, 513, 523 [NASA ADS] [CrossRef] [Google Scholar]
  33. Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117 [CrossRef] [Google Scholar]
  34. Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109 [NASA ADS] [CrossRef] [Google Scholar]
  35. Galli, D., Lizano, S., Shu, F. H., & Allen, A. 2006, ApJ, 647, 374 [NASA ADS] [CrossRef] [Google Scholar]
  36. Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237 [Google Scholar]
  37. Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29 [NASA ADS] [CrossRef] [Google Scholar]
  38. Grant, S. L., Espaillat, C. C., Wendeborn, J., et al. 2021, ApJ, 913, 123 [NASA ADS] [CrossRef] [Google Scholar]
  39. Gundlach, B., & Blum, J. 2015, ApJ, 798, 34 [Google Scholar]
  40. Haisch, K. E., Jr., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
  41. Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [Google Scholar]
  42. Holland, W. S., Matthews, B. C., Kennedy, G. M., et al. 2017, MNRAS, 470, 3606 [NASA ADS] [CrossRef] [Google Scholar]
  43. Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Hutchison, M. A., & Clarke, C. J. 2021, MNRAS, 501, 1127 [Google Scholar]
  46. Hutchison, M. A., Laibe, G., & Maddison, S. T. 2016a, MNRAS, 463, 2725 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hutchison, M. A., Price, D. J., Laibe, G., & Maddison, S. T. 2016b, MNRAS, 461, 742 [NASA ADS] [CrossRef] [Google Scholar]
  48. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  49. Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78 [Google Scholar]
  50. Kataoka, A., Tsukagoshi, T., Momose, M., et al. 2016, ApJ, 831, L12 [Google Scholar]
  51. Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117 [Google Scholar]
  52. Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016, MNRAS, 461, 2257 [NASA ADS] [CrossRef] [Google Scholar]
  53. Komaki, A., Nakatani, R., & Yoshida, N. 2021, ApJ, 910, 51 [NASA ADS] [CrossRef] [Google Scholar]
  54. Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, ApJ, 745, 19 [Google Scholar]
  55. Krivov, A. V., & Wyatt, M. C. 2021, MNRAS, 500, 718 [Google Scholar]
  56. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7 [NASA ADS] [CrossRef] [Google Scholar]
  58. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Lesur, G., Ercolano, B., Flock, M., et al. 2022, Protostars and Planets VII, eds. S.-I. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [Google Scholar]
  61. Li, R., & Youdin, A. N. 2021, ApJ, 919, 107 [NASA ADS] [CrossRef] [Google Scholar]
  62. Li, M., Huang, S., Zhu, Z., Petaev, M. I., & Steffen, J. H. 2021, MNRAS, 503, 5254 [NASA ADS] [CrossRef] [Google Scholar]
  63. Liffman, K. 2003, PASA, 20, 337 [NASA ADS] [CrossRef] [Google Scholar]
  64. Liu, H. B. 2019, ApJ, 877, L22 [NASA ADS] [CrossRef] [Google Scholar]
  65. Liu, B., Lambrechts, M., Johansen, A., & Liu, F. 2019, A&A, 632, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Long, F., Andrews, S. M., Rosotti, G., et al. 2022, ApJ, 931, 6 [NASA ADS] [CrossRef] [Google Scholar]
  67. Mamajek, E. E. 2009, in American Institute of Physics Conf. Ser., 1158, 3 [NASA ADS] [Google Scholar]
  68. Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2022, Protostars and Planets VII, eds. S.-I. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [Google Scholar]
  70. McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 [Google Scholar]
  71. Menou, K., & Goodman, J. 2004, ApJ, 606, 520 [NASA ADS] [CrossRef] [Google Scholar]
  72. Michel, A., van der Marel, N., & Matthews, B. C. 2021, ApJ, 921, 72 [CrossRef] [Google Scholar]
  73. Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. I., & Kataoka, A. 2022, Protostars and Planets VII, eds. S.-I. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [Google Scholar]
  74. Mulders, G. D., Pascucci, I., Ciesla, F. J., & Fernandes, R. B. 2021, ApJ, 920, 66 [NASA ADS] [CrossRef] [Google Scholar]
  75. Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016a, ApJ, 818, 16 [Google Scholar]
  76. Musiolik, G., Teiser, J., Jankowski, T., & Wurm, G. 2016b, ApJ, 827, 63 [NASA ADS] [CrossRef] [Google Scholar]
  77. Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018a, ApJ, 857, 57 [NASA ADS] [CrossRef] [Google Scholar]
  78. Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018b, ApJ, 865, 75 [Google Scholar]
  79. Nakatani, R., Kobayashi, H., Kuiper, R., Nomura, H., & Aikawa, Y. 2021, ApJ, 915, 90 [NASA ADS] [CrossRef] [Google Scholar]
  80. Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2022, Protostars and Planets VII, eds. S.-I. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura [Google Scholar]
  81. Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M.-A. 2011a, ApJ, 731, 95 [NASA ADS] [CrossRef] [Google Scholar]
  82. Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M.-A. 2011b, ApJ, 731, 96 [NASA ADS] [CrossRef] [Google Scholar]
  83. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Owen, J. E., & Kollmeier, J. A. 2019, MNRAS, 487, 3702 [Google Scholar]
  85. Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 [Google Scholar]
  86. Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [Google Scholar]
  87. Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880 [Google Scholar]
  88. Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125 [Google Scholar]
  89. Pearce, T. D., Launhardt, R., Ostermann, R., et al. 2022, A&A, 659, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Pérez, L. M., Carpenter, J. M., Chandler, C. J., et al. 2012, ApJ, 760, L17 [CrossRef] [Google Scholar]
  91. Pérez, L. M., Chandler, C. J., Isella, A., et al. 2015, ApJ, 813, 41 [CrossRef] [Google Scholar]
  92. Pfalzner, S., Dehghani, S., & Michel, A. 2022, ApJ, 939, L10 [NASA ADS] [CrossRef] [Google Scholar]
  93. Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691 [Google Scholar]
  94. Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611 [NASA ADS] [CrossRef] [Google Scholar]
  95. Pinilla, P., Pascucci, I., & Marino, S. 2020, A&A, 635, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Pinilla, P., Lenz, C. T., & Stammler, S. M. 2021, A&A, 645, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [Google Scholar]
  98. Portegies Zwart, S. F. 2016, MNRAS, 457, 313 [CrossRef] [Google Scholar]
  99. Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
  100. Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790 [Google Scholar]
  102. Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2018, MNRAS, 478, 3674 [Google Scholar]
  103. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  104. Schobert, B. N., Peeters, A. G., & Rath, F. 2019, ApJ, 881, 56 [NASA ADS] [CrossRef] [Google Scholar]
  105. Sellek, A. D., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 498, 2845 [NASA ADS] [CrossRef] [Google Scholar]
  106. Sellek, A. D., Clarke, C. J., & Ercolano, B. 2022, MNRAS, 514, 535 [NASA ADS] [CrossRef] [Google Scholar]
  107. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  108. Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 [Google Scholar]
  109. Smith, R. J., Glover, S. C. O., Bonnell, I. A., Clark, P. C., & Klessen, R. S. 2011, MNRAS, 411, 1354 [NASA ADS] [CrossRef] [Google Scholar]
  110. Takahashi, S. Z., Inutsuka, S.-I., & Machida, M. N. 2013, ApJ, 770, 71 [NASA ADS] [CrossRef] [Google Scholar]
  111. Tazzari, M., Testi, L., Ercolano, B., et al. 2016, A&A, 588, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133 [NASA ADS] [CrossRef] [Google Scholar]
  113. Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43 [CrossRef] [Google Scholar]
  114. Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130 [CrossRef] [Google Scholar]
  115. Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818 [NASA ADS] [CrossRef] [Google Scholar]
  116. Trapman, L., Facchini, S., Hogerheijde, M. R., van Dishoeck, E. F., & Bruderer, S. 2019, A&A, 629, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  117. Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Ueda, T., Kataoka, A., Zhang, S., et al. 2021, ApJ, 913, 117 [Google Scholar]
  119. van der Marel, N. 2023, EPJ+, 138, 225 [Google Scholar]
  120. van der Marel, N., & Mulders, G. D. 2021, AJ, 162, 28 [NASA ADS] [CrossRef] [Google Scholar]
  121. van der Marel, N., Williams, J. P., Picogna, G., et al. 2022, A&A, submitted, [arXiv:2204.08225] [Google Scholar]
  122. van Terwisga, S. E., Hacar, A., van Dishoeck, E. F., Oonk, R., & Portegies Zwart, S. 2022, A&A, 661, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164 [EDP Sciences] [Google Scholar]
  124. Villenave, M., Ménard, F., Dent, W. R. F., et al. 2021, A&A, 653, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11 [NASA ADS] [CrossRef] [Google Scholar]
  126. Wang, L., & Goodman, J. 2017, ApJ, 847, 11 [Google Scholar]
  127. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
  128. Winter, A. J., & Haworth, T. J. 2022, Eur. Phys. J. Plus, 137, 1132 [NASA ADS] [CrossRef] [Google Scholar]
  129. Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719 [NASA ADS] [CrossRef] [Google Scholar]
  130. Yang, H., & Bai, X.-N. 2021, ApJ, 922, 201 [CrossRef] [Google Scholar]
  131. Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  132. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [Google Scholar]
  133. Zhu, Z., Hartmann, L., & Gammie, C. 2009, ApJ, 694, 1045 [Google Scholar]
  134. Zhu, Z., Hartmann, L., Gammie, C. F., et al. 2010, ApJ, 713, 1134 [Google Scholar]
  135. Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18 [NASA ADS] [CrossRef] [Google Scholar]
  136. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  137. Zurlo, A., Cieza, L. A., Pérez, S., et al. 2020, MNRAS, 496, 5089 [Google Scholar]

All Tables

Table A.1

Parameters of Eq. (A.1) and Eq. (A.2) for different stellar masses.

Table A.2

Stellar masses and corresponding X-ray luminosities used in the Pic21 photoevaporation model.

Table A.3

Parameters of Eq. (A.3) for different stellar masses.

Table A.4

Stellar masses and corresponding FUV and X-ray luminosities used in the Kom21 photoevaporation model.

All Figures

thumbnail Fig. 1

Photoevaporation rates for a 1 M-star using identical X-ray luminosity for both the Pic21 (blue) and Kom21 (yellow) prescriptions. The two dashed lines show the Kom21 models where the FUV luminosity has been reduced by a factor of 10 and of 100.

In the text
thumbnail Fig. 2

Temperature (left), gas surface density (middle), and dust surface density (right) evolution for a disc forming from a 1 M 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 of the dust is retained outside this gap and is pushed outwards with time.

In the text
thumbnail 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 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 and 5 au, and one at 0.6–0.9 Myr at radii from 0.1 to 1.5 au. The particle size generally is larger in the inner regions of the disc and at earlier times.

In the text
thumbnail Fig. 4

Gas accretion rate as a function of disc age for all discs in both photoevaporation prescriptions. Once photoevaporation opens up a gap in the gas disc, the gas within disc gap is quickly cleared by photoevapo-ration 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).

In the text
thumbnail 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 αv = 10−2. Gas discs are formed with masses between 0.03 M to 0.5 M. After 10 Myr of evolution about 1% of this initial mass remains. Dust disc masses are initially between 10 M to 1000 M. 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 Pic2l or Kom2l 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.

In the text
thumbnail Fig. 6

Cumulative disc dust mass distributions for the model (left plot) using using α = 10−2. Here we show the Pic2l (solid) and the Kom2l (dashed) photoevaporation prescriptions. A sample of observed disc masses are shown in the right plot. The two photoevaporation prescriptions are very similar at 0.5 and 1 Myr and begin to differ at the low-mass end at 2 Myr. At 4–5 Myr the Pic2l retains dust in some discs, whereas the Kom2l has very little left. We note that the 3, 4, and 5 Myr lines of the Pic2l are on top of each other. The thin coloured lines in the left plot show the gas mass multiplied by the initial dust-to-gas ratio. This shows how the disc dust mass would evolve without any radial drift. The observed sample is colour coded according to the cluster’s oldest inferred age.

In the text
thumbnail Fig. 7

Gas accretion rate as a function of dust disc mass for the models using the Pic2l (black) and the Kom2l (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 yr−1). 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.

In the text
thumbnail Fig. 8

Cumulative distributions of disc dust masses from the population synthesis model using α = 10−3 (left) and observed discs (right). The observed sample is taken from Tychoniec et al. (2020, Perseus) and Manara et al. (2022). The distributions of the observed sample are scaled according to the disc fractions from Michel et al. (2021), taking age of each region as the mid-point of the given age range.

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

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

In the text
thumbnail Fig. 11

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 of dust, but whose gas accretion rates are lower than 10−12 M Yr−1. 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 Rgas = Rdust (solid) and Rgas = 4Rdust (dashed).

In the text
thumbnail Fig. C.1

Cumulative dust mass distributions of the simulations using αν = 10−2 and the Kom21 photoevaporation prescription with pebble drift (thin lines) and with pebble drift halted in 5% of discs (thick lines). Halting pebble drift extends the high-mass tail of the CDF, bringing it more in line with what is seen in observations.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.