© ESO 2020

1 Introduction

Observations of star-forming regions have now provided a sufficiently large sample of protoplanetary discs to allow statistical analyses of their properties (Ansdell et al. 2016; Manara et al. 2017; Alcalá et al. 2017; Andrews et al. 2018). Disc properties such as gas and dust disc masses and their respective outer radii are important because they directly influence the formation of planets by setting the available mass reservoir and the time-span in which they can form. However, it is not trivial to interpret the distribution of disc properties in star-forming regions. For instance, the mass reservoir in observed protoplanetary discs does not appear to be large enough to be able to form the commonly observed exoplanetary systems (Najita & Kenyon 2014; Manara et al. 2018).

Star formation takes place in giant molecular clouds. These clouds fragment into dense cores that in turn collapse gravitationally to form stars (McKee & Ostriker 2007). Conservation of angular momentum causes some of the gas and dust in the cloud core to land in a disc surrounding the protostar. The young protostar grows in mass by accreting gas from the surrounding disc. This accretion is driven by a redistribution of angular momentum in the disc, although the physical process behind the angular momentum redistribution is uncertain (Turner et al. 2014).

In order to model the evolution of the dust component in protoplanetary discs, detailed models have been developed that include Takeuchi & Lin (2005); Takeuchi et al. (2005), Brauer et al. (2008), and Birnstiel et al. (2010, 2012). These works show that compact millimetre(mm)-sized dust particles are drained in less than 1 Myr, but that fluffy grains are able to remain for longer. Takeuchi et al. (2005) also found that in the presence of strong photoevaporation, the inner gas disc can be cleared out before the dust drains. The dust remaining in the outer disc then piles up near the inner edge of the gas disc. Takeuchi & Lin (2005) also showed that the rapid radial drift of dust weakens the emission from the disc at mm wavelengths, and that the dust needs to be in the size range between 1 mm and 1 cm for the emission at mm wavelengths to remain present for more than 1 Myr.

Surveys of protoplanetary discs provide insight into scaling relations between disc properties and the stellar mass. A super-linear scaling relation between the stellar mass and the dust disc mass was observed by Pascucci et al. (2016). In the gas accretion-rate-versus-stellar-mass scaling relation, indications of a broken power have been found (Manara et al. 2017; Alcalá et al. 2017). Manara et al. (2017) found $\dot{M}\propto M_{\star }^{4.31}$ for M_⋆ ≤ 0.3 M_⊙ and $\dot{M}\propto M_{\star }^{0.56}$ for M_⋆ ≥ 0.3 M_⊙, and Alcalá et al. (2017) found $\dot{M}\propto M_{\star }^{4.58}$ for M_⋆ ≤ 0.2 M_⊙ and $\dot{M}\propto M_{\star }^{1.37}$ for M_⋆ ≥ 0.2 M_⊙. Manara et al. (2017) also made a single power-law fit and found $\dot{M}\propto M_{\star }^{2.3}$.

Some progress has been made in trying to understand the scaling relations mentioned above. Dullemond et al. (2006) used a disc formation and evolution model to investigate the scaling between the gas accretion rate and the stellar mass. These latter authors found $\dot{M}\propto M_{\star }^{1.8}$, if the ratio between the rotational and gravitational energy of the cloud core and the viscous α parameter are kept the same for cloud cores of different masses. This agrees well with observations measuring the gas accretion rate and stellar mass scaling with a single power law (Muzerolle et al. 2003; Natta et al. 2004; Manara et al. 2016a; Alcalá et al. 2017), but is in contrast with recent indications of a broken power law scaling (Manara et al. 2017; Alcalá et al. 2017).

Using Monte Carlo methods, Alexander & Armitage (2009) modelled viscous disc evolution around 1000 protoplanetary discs including photoevaporation and gap opening due to the presence of a giant planet. With a maximum initial gas disc mass of ≈ 0.12 M_⊙ and a mean of 0.0045 M_⊙, the authors find median gas accretion rates between ≃0.5 × 10⁻⁸ M_⊙ yr⁻¹ at t ≃ 0.1 Myr and ≃10⁻¹⁰ M_⊙ yr⁻¹ at t ≃ 4 Myr. Disc lifetimes range between 2.3 and 10.7 Myr.

In contrast to these studies we aim to model the formation and evolution of a representative sample of protoplanetary discs in a star-forming region by evolving both the gas and dust component. We are especially interested in the evolution of the dust mass, as it sets the initial reservoir of mass that is available for planet formation and because it is an important observable in disc surveys.

This paper explores the formation and evolution of protoplanetary gas and dust discs. In Sect. 2 we describe the numerical model that we use for the formation and evolution of the protoplanetary disc. We do not use a singular istothermal sphere (Shu 1977; Hueso & Guillot 2005) to calculate the infall of mass onto the disc from the collapsing core, but instead a Bonnor-Ebert sphere (Bonnor 1956; Ebert 1957) to model the cloud core, as in Takahashi et al. (2013). We combine the infall with viscous evolution and disc gravitational instability in a one-dimensional numerical model to simulate disc evolution. At first, we look at the formation and evolution of a single disc and then we explore the effects of the angular momentum budget of the cloud core on disc evolution. We then use a Monte Carlo method to study a cluster of discs in a star-forming region that evolves over several million years. In Sect. 3, we present our main results. We find that during evolution of the disc the dust component goes through three main phases: the formation phase, a phase of slow radial dust drift due to high gas surface densities, and one of rapid radial dust drift once the gas surface density drops sufficiently. Our most important conclusion is that the low observed dust disc masses in star-forming regions (Ansdell et al. 2016 found ${\overline{M}}_{\text{d}}\approx 15M_{\oplus }$ for 1–3 Myr-old star-forming regions) is consistent with radial drift of dust pebbles. In Sect. 4 we focus on the observationally relevant scaling relationships. Of particular interest is the relation between the main observables: the stellar gas accretion rate and the dust disc mass. In Sect. 5 we discuss the results and compare them to observations and examine their implications on planet formation. Finally, in Sect. 6 we summarise our conclusions.

2 Model

2.1 Disc formation

We use a model for the collapse of the cloud core and the formation of the circumstellar disc similar to the one described in Takahashi et al. (2013). The molecular cloud core is modelled as a Bonnor-Ebert sphere (Bonnor 1956; Ebert 1957) with a constant temperature T_cc = 10 K throughout the entire cloud and a solid-body rotation rate Ω₀. In order to trigger the gravitational collapse of the cloud core, we follow Takahashi et al. (2013) and increase its density by a factor f = 1.4. For a discussion on the effect of using a different value of f we refer to Takahashi et al. (2013).

The cloud core collapses inside out, with the radius of the collapse front, r_cf, expanding outwards in time. Takahashi et al. (2013) gives the time it takes for the collapse to reach a radius r_cf in the cloud core as $$t_{\text{collapse}}=\sqrt{\frac{r_{\text{cf}}^3}{2GM(r_{\text{cf}})}}{\displaystyle {\int }_0^1\frac{\text{d}R}{\sqrt{\frac{1}{f}\ln R+\frac{1}{R}-1}}},$$(1)

where G is the gravitational constant, M(r_cf) is the mass inside the collapse front radius, and R = r_cf∕r_cc, where r_cc is the radius of the cloud core.

The rate at which gas will fall onto the disc and star, Ṁ_infall, is given by $${\dot{M}}_{\text{infall}}(t)=4\pi r_{\text{cf}}{(t)}^2\rho (r_{\text{cf}}(t))\frac{\text{d}{r}_{\text{cf}}}{\text{d}t},$$(2)

where ρ(r_cf) is the density in the cloud core at r_cf and dr_cf∕dt is the velocity at which the collapse front expands outwards. A derivation for an expression for d r_cf ∕dt can be found in Appendix A.

The gravitational collapse of the cloud conserves angular momentum. This causes the infalling material to land on an orbit in the protoplanetary disc where its Keplerian angular momentum is the same as the angular momentum it had in the cloud core. The outermost radius at which material lands at a given time is called the centrifugal radius, a_c, and is given by $$a_{\text{c}}(t)=\frac{{\Omega }_0^2{r}_{\text{cf}}{(t)}^4}{GM(r_{\text{cf}}(t))}.$$(3)

The infalling material is spread out across the disc within the centrifugal radius such that its angular momentum is conserved. The change in the surface density of the disc is given by the following equation from Hueso & Guillot (2005): $$S_{\text{g}}(r,t)=\frac{{\dot{M}}_{\text{infall}}}{8\pi a_{\text{c}}^2}{\left(\frac{r}{a_{\text{c}}}\right)}^{-3/2}{\left[1-{\left(\frac{r}{a_{\text{c}}}\right)}^{1/2}\right]}^{-1/2},$$(4)

where r is the orbital radius in the disc. For the dust component of the infall we assume the source term to be $$S_{\text{d}}(r,t)=0.01S_{\text{g}}(r,t),$$(5)

where we assumed the dust-to-gas ratio to be 1/100 to match interstellar medium values and to maintain a fixed initial metallicity in the disc.

2.2 Disc evolution

2.2.1 Gas disc

We evolve the protoplanetary gas disc by viscous evolution, where the disc expands outwards as matter is accreted onto the central star. The viscous evolution of Pringle (1981) is combined with the source term from Eq. (4), which gives a total gas surface density evolution of $$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}\left(\nu {\Sigma }_{\text{g}}{r}^{1/2}\right)\right]+S_{\text{g}}(r,t),$$(6)

where Σ_g is the gas surface density and ν is the turbulent viscosity. To model the turbulent viscosity we use the standard α parametrisation, $$\nu =\alpha c_{\text{s}}{H}_{\text{g}};$$(7)

(Shakura & Sunyaev 1973). Here, c_s is the sound speed in the disc and H_g = c_s∕Ω_k is the scale height of the gas disc. We assume that the total viscosity parameter is given by $$\alpha ={\alpha }_{\text{min}}+{\alpha }_{\text{Q}},$$(8)

where α_min = 10⁻² is the minimum value throughout the disc. This model parameter was chosen pragmatically to evolve the α-disc within the observed disc lifetimes of a few million years (Haisch et al. 2001; Fedele et al. 2010). However, recent numerical efforts have shown that disc midplanes are almost laminar (Bai & Stone 2013; Gressel et al. 2015) and point towards very low midplane alpha values (on the order of α_mid ≈ 10⁻⁴), in apparent agreement with protoplanetary disc scale heights (Pinte et al. 2016).

The effects of using a lower α_min are discussed in Appendix B. The parameter α_Q models the effect of gravitational instabilities. When the disc is young it can become very massive relative to the central star and become gravitationally unstable. Gravitational instability gives rise to gravitational torques in the disc that act as an increased disc viscosity. We model the gravitational torque following Zhu et al. (2010) and Takahashi et al. (2013) and use $${\alpha }_{\text{Q}}=A\exp \left(-B{Q}^4\right),$$(9)

with A = 1 and B = 1, and where Q is the Toomre parameter given by $$Q=\frac{c_{\text{s}}{\Omega }_{\text{k}}}{\pi G{\Sigma }_{\text{g}}}.$$(10)

These choices for A and B cause α_Q to become important relative to α_min when Q ≲ 1.4. Three-dimensional global simulations have shown this to be an appropriate limit for gravitational torques to become important (e.g. Boley et al. 2006). The precise choice of A and B has been demonstrated to have no significant effect on the long-term gas disc evolution (Takahashi et al. 2013).

2.2.2 Dust disc

We use a simplified model for the dust component of the disc. We assume that the dust is of constant size (r_p = 100 μm) throughout the disc. There are a number of motivations behind this choice. First, observations of the spectral index in envelopes around forming stars show that dust particles already as large as 100 μm to 1 mm in size could be present (Miotello et al. 2014; Agurto-Gangas et al. 2019). Likewise, observations of very young stars indicate the presence of mm-sized dust already in the formation phase of the disc (Harsono et al. 2018). Observations of polarised emission from protoplanetary discs also indicate that the maximum particle size is ≈ 100 μm if compact and spherical dust particles are assumed (Kataoka et al. 2016; Ohashi et al. 2018).

From laboratory experiments, it also appears that particle growth by coagulation in protoplanetary discs is inefficient past roughly mm sizes (Güttler et al. 2010; Zsom et al. 2010). This holds for both silicates and ice particles, because at temperatures below 180 K the surface energy of ice grains is similar to that of silicates (Musiolik & Wurm 2019).

Hence, our assumption of dust growth to 100 μm should accurately capture the drift of dust particles in the outer regions of the protoplanetary disc. Another simplification is that we do not include feedback of dust on gas, which is appropriate for low dust-to-gas ratios.

The aerodynamic coupling of the dust to the gas can be measured with the dimensionless Stokes number. In the Epstein regime, when the dust particles are smaller than the mean free path of the gas, and assuming vertical hydrostatic equilibrium in the disc, theStokes number is given by $${\tau }_{\text{s}}=\frac{\sqrt{2\pi }{\rho }_{\text{p}}{r}_{\text{p}}}{{\Sigma }_{\text{g}}},$$(11)

where ρ_p is the material density of the dust particles. We set ρ_p = 1.6 g cm⁻³ following Birnstiel et al. (2010). We also include the transition to the Stokes drag regime when the dust particles are larger than the mean free path of the disc, which can occur in the inner disc, as described in Youdin (2010).

The dust surface density is evolved by advecting it with the gas and by letting the dust drift due to gas drag. We ignore turbulent diffusion of dust as the midplane is assumed to be almost laminar. In the early stages of disc evolution the advection term dominates because the high gas surface densities result in low Stokes numbers. Eventually, the gas surface density drops and Stokes numbers increase. The dust begins to decouple from the gas and radial drift of dust becomes important. We model this as an additional advection term on the dust. To calculate the radial dust velocity we use the equation given by Nakagawa et al. (1986), $$v_{\text{r},\text{d}}=-\frac{2{\tau }_{\text{s}}}{1+{\tau }_{\text{s}}^2}\eta v_{\text{k}},$$(12)

where η is a measure of the radial gas pressure support and v_k is the Keplerian velocity. Here, η is given by $$\eta =-\frac{1}{2}{\left(\frac{H_{\text{g}}}{r}\right)}^2\frac{\partial \ln P}{\partial \ln r},$$(13)

where the last factor is the radial pressure gradient.

2.3 Disc temperature

For simplicity, we assume that the disc is locally isothermal, with the temperature given as in the Hayashi MMSN model (Hayashi 1981), but with the minimum value set by the cloud core temperature of T_cc = 10 K, $$T(r)=\{\begin{array}{ll}280{\left(\frac{L_{\star }}{L_{\odot }}\right)}^{1/4}{\left(\frac{r}{\text{au}}\right)}^{-1/2}\text{K},\hfill & \text{if}T>T_{\text{cc}}\hfill \\ T_{\text{cc}},\hfill & \text{else},\hfill \end{array}$$(14)

where L_⋆ is the luminosity of the central star. For low-mass stars that are a few million years old or less, the stellar luminosity scales with the stellar mass roughly as $$L_{\star }\propto M_{\star }^{1-2};$$(15)

(Baraffe et al. 2002). The stellar luminosity is therefore expected to increase as the star grows in mass over time. However, our model does not take into account the heating of the disc by the infalling material dissipating its energy in the disc. Therefore, in our model, evolving the stellar luminosity with the stellar mass results in an exceedingly cold disc during the initial stages of a discs life. Instead we assume the luminosity to be given by a star of the same mass as the cloud core after a few years of evolution, according to Eq. (15). For this study, we chose to use $$\frac{L_{\star }}{L_{\odot }}\propto {\left(\frac{M_{\star }}{M_{\odot }}\right)}^2.$$(16)

In our model we chose to ignore the effects of viscous heating, which is thought to contribute to the temperature structure of the disc at short orbital radii. However, Mori et al. (2019) used simulations of protoplanetary discs under non-ideal MHD conditions to show that in the inner disc, where accretion might be driven by magnetised disc winds and heating is caused by Joule dissipation rather than viscous dissipation, heating of the disc midplane by accretion is weak and stellar irradiation dominates the disc heating. Only at radii of ≲0.5 au can viscous heating become important. Outside of the viscously dominated region, stellar irradiation sets the disc temperature. Although, during the early phases of the disc when accretion rates are higher, and especially during accretion bursts, the region dominated by viscous heating can extend out to a few tens of astronomical units (Cieza et al. 2016). As we see below, in our model the radial drift of the dust occurs from the outside in. Therefore, the focus of our interest is the outer disc at radii of more than at least a few tens of astronomical units. Viscous heating can be neglected at these radii.

2.4 Simulations

In this section we describe the simulations that we run and explain the reasoning behind certain modelling choices. A few parameters are the same for all simulations. These are the temperature of the cloud core, which we assume to be T_cc = 10 K, the value of the minimum viscous α parameter α_min = 10⁻², and the size of the dust, which we assume to be r_p = 100 μm throughout the disc. We refer to Appendix B for the effect of using a lower α_min value and Appendix C for the effect of using a larger dust size.

In short, the simulations that we perform are: (a) a reference case for understanding the evolution of a single disc, (b) a set of runs where we look at the effect of the centrifugal radius and cloud core mass on disc evolution and (c) a set of Monte Carlo runs to model the evolution of a star-forming cluster.

2.4.1 Reference case

The first simulation we run is the reference case, which we use to understand how a single protoplanetary disc forms and evolves. This case consists of a cloud core mass of M_cc = 1 M_⊙. We set the angular momentum of the cloud core such that the final centrifugal radius is a_c = 50 au.

2.4.2 Changing the centrifugal radius and cloud core mass

After we have analysed the reference case we examine the effect of changing the final centrifugal radius of the disc. For this we use the centrifugal radii listed for the 1 M_⊙ case in Table 1. We chose an upper limit of 128 au for two reasons. One is that observations show that disc radii can be as high as 100–200 au (Ansdell et al. 2018). Observations also show that the centrifugal radius can be as high as 100–200 au even aroundyoung forming stars (Sakai et al. 2014, 2016). Additionally, we find that if the centrifugal radius becomes very large (≫ 100 au), the protoplanetary disc lifetime becomes much longer than the expected lifetimes of a few up to 10 Myr (Haisch et al. 2001; Fedele et al. 2010), assuming α_min = 10⁻².

After we have looked at the effect of the centrifugal radius on the disc evolution we also introduce a change in the cloud core mass and therefore also the resulting stellar mass. The other mass cases listed are in Table 1. We use the 1 M_⊙ case as the reference and then scale the centrifugal radius range for the other cases linearly with the cloud core mass. This is based on the fact that the ratio between the rotational and gravitational energy of molecular cloud cores is independent of the cloud core size (e.g. Goodman et al. 1993; Chen & Ostriker 2018). By setting this ratio to be the same for all cloud cores, one finds that the centrifugal radius of the protoplanetary disc scales linearly with the mass of a Bonnor-Ebert sphere. We note that from Eq. (3) one might expect the centrifugal radius to scale as the inverse of the cloud core mass. However, the radius of a Bonnor-Ebert sphere scales linearly with its mass, giving a linear scaling between the centrifugal radius and cloud core mass (see Appendix D for a more detailed description).

2.4.3 Monte Carlo cluster

Once we have explored the effects of both the centrifugal radius and the mass of the cloud core we move on to using a Monte Carlo method for selecting the final centrifugal radius and the cloud core mass. With this we aim to mimic the evolution of a star-formingcluster for several million years. A total of 100 young stars and their protoplanetary discs are simulated. We do not model any interactions between the different star–disc systems, for example through gravity or gas dispersion due to radiation from nearbystars. That is, each star–disc system is modelled completely independently from the others, but we analyse them as a group.

We pick cloud core masses from a stellar IMF under the assumption that the entire mass of the cloud core will eventually end up in the star. For this we use the IMF of Kroupa (2001) in the mass range 0.1–2.0 M_⊙. The final centrifugal radius is picked from a normal distribution with a mean of 40 and a standard deviation of 35. This value is then scaled with the mass of the cloud core. We refer to Appendix D for the motivation behind this scaling. We further put a lower and upper limit on the centrifugal radius of 2 au and 120 au. Due to the scaling of centrifugal radius with cloud core mass this lower limit results in the low-mass cases being biased towards discs that are large and massive relative to the star. The lower limit is constrained by computational costs and the upper limit is set because discs with overly large centrifugal radii live much longer than the expected protoplanetary disc lifetimes as discussed in Sect. 3.2.

Stellar clusters form from clumps inside giant molecular clouds (McKee & Ostriker 2007). Simulations suggest that the lifetimes of these clumps appear to range from 10⁵ to 10⁷ yr (Falceta-Gonçalves & Lazarian 2011). To mimic the stars forming at different times we give a time delay of between 0 and 1 Myr for the collapse of each cloud core from when the cluster begins to form.

Other distributions of the centrifugal radius within the same lower and upper limit were tested, including a normal distribution with different mean and standard deviation, and a uniform distribution, but no significant differences were found on the cluster evolution. We do find that a distribution that is weighted towards higher centrifugal radius cases gives a larger average dust disc mass. However, the age spread between the discs also significantly affects the average dust mass, and this is a parameter that can have a large spread.

Fig. 1 Snapshots of gas (solid lines) and dust (dashed lines) surface densities of the reference case. In the evolution phase, the inner disc is drained of gas as it is accreted by the star and the outer disc expands as gas is transported outwards to conserve angular momentum. In the first ≈1.0 Myr the dust surface density mimics that of the gas in most of the disc out to radii of about 100 au. Beyond this radius, radial drift of dust is already rapid, draining the dust disc. When the disc is between 2 and 2.5 Myr old, the dust disc is completely drained.

2.5 Method

We evolve the disc using a one-dimensional numerical model similar to that of Bath & Pringle (1981) and Booth et al. (2017). The computational domain is divided into 100 logarithmically spaced cells. The outer edge is 10 000 au and the inner edge is set such that it is 1/100 of the final centrifugal radius. Gas and dust that lands inside the inner edge of the domain is treated as being directly accreted onto the star.

3 Results

3.1 Reference case

First we investigate the reference case with parameters as explained in Sect. 2.4. We numerically evolve the gas and dust surface density. Therefore, we begin by examining the evolution of these two quantities. Figure 1 shows snapshots intime of the gas (solid lines) and dust (dashed lines).

We define the formation phase of the disc as the time period during which there is mass infalling from the collapsing cloud core onto the disc, that is, the time it takes for the cloud core to collapse. The collapse of a 1 M_⊙ cloud core at a temperature of 10 K lastsabout 1.8 × 10⁵ yr. This depends on the temperature of the cloud core, as it affects the size of the Bonnor-Ebert sphere.

The gas surface density in the disc during this formation phase changes both by the infalling material and due to viscous evolution, but the former dominates the mass evolution for most of this phase. Initially, this leads to high surface densities in the inner parts of the disc, almost reaching 10⁵ g cm⁻² at the inner edge of our simulation domain (0.5 au). Towards the end of the formation phase the disc becomes gravitationally unstable at radii from r ≈ 2 au to r ≈ 200 au. In these gravitationally unstable regions the viscosity increases because of gravitational torques according to the prescription in Eqs. (8) and (9). The increased viscosity leads to more efficient transport of gas compared to the stable regions.

Once the formation phase is over, the gas disc only evolves viscously. With constant particle size and density, the Stokes number is a function of the gas surface density only (assuming the Epstein drag regime, see Eq. (11)). Therefore, as long the gas surface density is sufficiently high, the Stokes number is kept low and dust drift becomes negligible. In this reference case it takes about 0.75 Myr for the gas surface density to decrease enough in the outer disc for the Stokes number to increase to the point where radial drift becomes significant. From this point, the dust disc is drained in about 1.75 Myr. This can be seen in the 2.5 Myr snapshot of Fig. 1.

We now explore how the stellar gas accretion rate evolves with the mass in the disc. This is shown in Fig. 2 for the reference case. There are three clear regimes for the dust disc and two for the gas disc, starting with the disc formation phase. During the first 0.1 Myr of this regime, both the stellar gas accretion rate and the dust-disc mass are increasing. In the remainder of the formation phase the gas accretion does not change significantly. It is during this later part of the formation phase that gravitational instability is active. The gas accretion rate peaks at Ṁ_g ≈ 5 × 10⁻⁶ M_⊙ yr⁻¹. The dust mass reaches a maximum mass of M_d ≈ 1.4 × 10³ M_⊕ and the gas mass reaches M_g ≈ 0.4 M_⊙. This far in the disc evolution, the dust-to-gas ratio remains close to Σ_d∕Σ_g ≈ 0.01 throughout the disc.

The second regime begins once the formation phase is over. During this regime the gas surface density is still high and as a consequence the Stokes numbers are low. Dust drift is therefore slow, leading to a very slow depletion of the dust disc. This regime last until about 1.3 Myr.

Subsequently, the dust clearing regime begins, which can be identified by the dust drift turn-over in Fig. 2. Now, the lower gas surface densities result in increased Stokes numbers of the particles and radial drift of dust can rapidly reduce the mass of the dust disc. At an orbital radius of 1 au, the Stokes number increases to τ_s ≈ 10⁻⁴ at t = 2 Myr. After 2.25 Myr the dust disc has been drained to less than 0.1 M_⊕.

In summary, we find that a single disc goes through three regimes in gas accretion rate versus dust disc mass. The first is the formation phase where both the accretion rate and the dust disc mass grow. After the formation phase is over, a phase of slow radial dust drift is entered. Eventually, the gas surface density drops and the Stokes number becomes large. Radial dust drift now becomes rapid and the dust disc is drained.

Fig. 2 Stellar gas accretion rate as a function of disc mass for the reference case. The black curve shows the dust disc mass and the red curve shows the gas disc mass. The black dots mark the age in years as labelled. There are three clear regimes in this relation. During the formation phase of the disc, the gas accretion rate and the disc dust mass are both increasing. Once the formation phase has ended, dust drift is slow and the gas accretion rate drops faster than the disc dust mass. Rapid draining of the dust disc mass begins around t ≈ 1.3 × 106 and lasts about 1 Myr.

3.2 Dependency on the centrifugal radius

In this section, we explore the effect of the centrifugal radius on the evolution of the protoplanetary disc. For a given mass of the cloud core, a change in the rotation frequency of the core, and thus the centrifugal radius, affect the evolutionary track of the disc in the phasespace of stellar gas accretion rate versus dust disc mass. This is shown in Fig. 3. The red and yellow lines show the evolutionary track of the highest and lowest centrifugal radius case, respectively. The different symbols represent different times as shown in the legend. For each symbol, the five points show the five final centrifugal radii cases of the 1 M_⊙ set in Table 1.

During the formation phase, the case with the smallest centrifugal radius (yellow line) has the highest gas accretion rate and the lowest dust mass of all cases. This trend holds until the end of the formation phase (t ≈1.8 × 10⁵ yr). As the discs evolve past the formation phase, a disc with a lower final centrifugal radius evolves faster. Therefore, it will reach the dust drift turn-over earlier and will drain its dust faster.

3.3 Dependency on cloud mass

We now investigate the effect of changing the mass of the cloud core, and therefore also the final stellar mass. We use the cases shown in Table 1. For a given cloud core temperature, the timescale for the formation phase of the disc increases nearly linearly with the mass of the cloud core.

Again, we look at the stellar gas accretion rate as a function of dust disc mass. This is shown in Fig. 4. Panels a–d show the discs at 0.25, 0.75, 1.5, and 3.0 Myr. Unsurprisingly, the more massive the cloud is, the more massive the disc becomes. However, this is not a strict relation, given our spread in centrifugal radii. Given a sufficiently high centrifugal radius from a lower-mass case, and a sufficiently low centrifugal radius from a higher-mass case, it is possible to find cases where the higher-mass core ends up with a lower-mass disc than the low-mass core. The reason for this is that with a high centrifugal radius, more mass is deposited onto the disc rather than directly onto the star. This can be seen in panels a and b of Fig. 4, where the 0.5 M_⊙ and 1.0 M_⊙ cases overlap in disc dust mass.

When the system is 0.75 Myr old (panel b of Fig. 4), we can see that the dust discs from the 0.1 M_⊙ core have been depleted of dust and that all the remaining discs line up, showing a trend of higher cloud core masses giving higher gas accretion rates and higher disc masses. As a general trend, the lower the beginning cloud core mass, the earlier the dust drift turn-over is reached. Because the discs from low-mass cloud cores reach the point of rapid radial drift earlier, theyare drained of dust faster. Therefore, at the 1.5 Myr mark (panel c of Fig. 4), the discs from the lower-mass cloud cores have much less massive dust discs. At 3 Myr (panel d) only the discs from the 2 M_⊙ still have a dust component remaining. The most massive of these two discs had a total mass of M_d ≈ 3 × 10³ M_⊕ at the end of the formation phase. This disc is drained to less than 10⁻¹ M_⊕ in 3.5 Myr.

To briefly summarise the present section and the previous one, we find that both gas and dust discs evolve faster for lower mass cloud cores (therefore those with lower resulting stellar mass) with smaller final centrifugal radii. As a result, the dust disc is drained faster around lower mass stars and in discs with lower final centrifugal radii.

Fig. 3 Stellar gas accretion rate and dust disc mass as a function of time, for discs with final centrifugal radii of 8, 16, 32, 64, and 128 au forming from 1 M⊙ cloud cores. The red and yellow lines show the tracks of the two cases with the highest and lowest final centrifugal radius, respectively. The formation phase ends at t ≈ 1.8 × 105 yr. During the formation phase a trend can be seen where systems with smaller centrifugal radii have lower dust disc masses and slightly higher gas accretion rates. Between 1.8 × 105 yr and 5.0 × 105 yr the gas accretion rates drop significantly, while dust disc masses decrease only slightly. Discs with a smaller final centrifugal radius experience a larger drop in accretion rate. A lower final centrifugal radius leads to faster evolution of the protoplanetary disc. At 2.5 × 106 yr, only the case with the highest final centrifugal radius remains, but the dust disc mass is only ≈ 10−2 M⊕.

3.4 Monte Carlo cluster

In this section, we continue to explore the relationship between the stellar gas accretion rate to disc mass, but now we use a Monte Carlo method to select the cloud core mass and the final centrifugal radius, with the aim of mimicking the evolution of a star-forming region for several million years. The parameter space we study is specified in Sect. 2.4.

Figure 5 shows the stellar gas accretion rate versus the disc gas and dust mass for the Monte Carlo cluster. The dust is shown in black circles and the gas in red squares. Panels a–d show the cluster at 1, 1.5, 2 and 3 Myr respectively. At 1 Myr, the majority of the discs have left the formation phase and are in the phase of slow dust drift. A few remain in the formation phase and are characterised by low disc masses and/or high accretion rates, and some of the lowest mass discs have passed the dust drift turn-over. At this point, the dust traces the gas very closely, at one-hundredth of the gas mass, since it is still tightly coupled to the gas in all but the least massive discs. At 1.5 Myr, all discs have left the formation phase (for the range of parameters used here this is inevitable: the most massive cores allowed take ≈ 0.4 Myr to collapse and the age spread between cores is only 1 Myr). Most discs are now near or at the dust drift turn-over. The lowest mass discs with the lowest centrifugal radii have passed the dust drift turn-over and the low-mass end of the dust track is now beginning to flatten out. When the cluster is 2 Myr old, nearly all discs have passed the dust drift turn-over and are in the regime of rapid radial dust drift, and many of the low-mass discs are devoid of dust. At 3 Myr we see that the majority of dust discs are drained, and of those that remain, only 4 out of 100 have a total mass of more than 1 M_⊕.

Once all the discs have left the formation phase, the gas component settles into a fairly narrow band in the accretion-rate–disc-mass space. This can be seen in panels b–d of Fig. 5. As the discs evolve forward in time, both the accretion rate and the disc masses decrease as the gas is accreted onto the star. The spread in the relation also decreases with time. That gas discs end up in such a tight space in the gas-accretion-rate-versus-disc-mass relation was also found theoretically by Lodato et al. (2017).

Fig. 4 Gas accretion rate as a function of dust disc mass for protoplanetary discs formed from cloud cores with masses and final centrifugal radii shown in Table 1. The yellow and blue lines trace the time evolution of the cases with the lowest and highest mass, and lowest and highest centrifugal radius, respectively. As a general trend, discs from low-mass cloud cores evolve faster than those from higher-mass cloud cores. However, the speed with which the disc evolves also depends on the final centrifugal radius. At 3 Myr (panel d), all dust discs except for the ones coming from 2 M⊙ cloud cores have been completely drained out.

4 Scaling relations

In this section we explore various scaling relations between the gas accretion rates, the disc, and the stellar masses for the reference case and the cluster of discs. We then compare the scaling relations in the cluster to observations. For the single disc and thecluster, the scaling relations are summarised in Tables 2 and 3, respectively.

4.1 The scaling of the mass accretion rate with the disc mass

Looking at the single star in the reference case, we find that during the formation phase the gas accretion rate scales with both gas and dust disc mass as $\dot{M}\propto M_{\text{d},\text{g}}^{0.5}$. After the formation phase but before the dust drift turn-over we find $\dot{M}\propto M_{\text{d}}^{4.0}$. After the turn-over we find $\dot{M}\propto M_{\text{d}}^{0.04}$. For the gas disc we find $\dot{M}\propto M_{\text{g}}^{3.0}$ in the post turn-over phase.

Now we examine the cluster of discs. In the gas accretion rate versus dust disc mass we do not find that a power law describes the data well at 2 Myr. At 3 Myr most of the discs are drained of dust, so we chose not to make a fit to the few remaining discs. A broken power law would likely provide a better fit than a single power law; one would represent discs that have not yet reached the dust drift turn-over, and the other would represent those that have passed it.

Observations report a near-linear scaling between the stellar accretion rate and the dust disc mass, $\dot{M}\propto M_{\text{d}}^{0.7}$ (Manara et al. 2016b), most similar to early times in our Monte Carlo run (≲ 1.5 Myr). These observations nevertheless also show a spread of several orders of magnitude in the gas accretion rate for a given dust mass. For an evolved cluster of discs, we find a much tighter relation with a spread that is only a factor of a few at most. The larger spread in observed discs could be due to a combination of different disc ages, intrinsic viscosity variations, and an intrinsic spread in dust sizes. We explore the effects of the viscosity and the dust size in Appendices B and C. Another difference is that we do not reach accretion rates as low as < 10⁻⁹ M_⊙ yr⁻¹; such low accretion rates are common in the sample of Manara et al. (2016b). However, our model does not include gas mass-loss mechanisms such as disc winds, disc photoevaporation, or accretion of gas onto planets. These could have the effect of decreasing the gas accretion rate for a given dust disc mass.

In the gas-accretion-rate-versus-gas-disc-mass relation in the cluster we find a near-linear scaling with $\dot{M}\propto M_{\text{g}}^{0.9}$ at 2 Myr and 3 Myr. This can be seen in the black dotted lines in panels c and b of Fig. 5. A common assumption made in observational studies is that the gas mass is 100 times the dust mass. However, at 2 and 3 Myr most or all of the discs in our model are in the phase of rapid radial dust drift where this assumption does not hold. Hence, a scaling relation measured for the dust mass does not necessarily translate directly onto the gas mass. Directly comparing the gas-accretion-rate-to-gas-disc-mass relation to observations is complicated because of the uncertainty in estimating the mass locked in H and He gas based on CO spectral line emission (Manara et al. 2016a).

Fig. 5 Gas accretion rate as a function of dust (black circles) and gas (red squares) disc mass for protoplanetary discs of the stellar cluster. The dashed black lines show power-law fits of the accretion rate as a function of the gas disc mass ($\dot{M}\propto M_{\text{g}}^{0.9}$). The discs are given a random individual age of between 0.0 and 1.0 Myr to mimic a spread in formation times. Panels from a to d: protoplanetary discs at the cluster ages of 1.0, 1.5, 2.0, and 3.0 Myr. For the first 1 Myr the dust is coupled well to the gas inmost discs, and dust disc masses evolve in a very similar way to the gas disc mass. Some of the lowest mass discs have reached the dust drift turn-over at this point. Beyond this time more and more of the dust discs reach the dust turn-over and begin to drain rapidly. At 1.5 Myr the lower mass discs have reached the turn-over point, as can be seen in the flattening of the low-mass end of the dust discs in panel b. At 2 Myr only the most massive dust discs are yet to reach the turn-over and at 3 Myr nearly all discs are drained of dust. Only 4 out of 100 discs still have a dust mass greater than 1 M⊕, although two of these have a mass of more than 100 Earth masses.

Fig. 6 Gas accretion rate as a function of stellar mass of the cluster of discs at 2 Myr (black) and 3 Myr (red). The blue and grey dashed lines are fits to the data at 2 and 3 Myr respectively. Both scaling relations have a slope of ≈ 1.4.

4.2 The scaling of the gas mass and stellar accretion rate with stellar mass.

The gas disc mass scales slightly steeper than linear with the stellar mass in the cluster. We find $M_{\text{g}}\propto M_{\star }^{1.4}$ at 2 Myr and 3 Myr. The gas accretion rate scales nearly linearly with stellar mass. At 2 and 3 Myr we find $\dot{M}\propto M_{\star }^{1.4}$. This is shown in Fig. 6. With a single power-law fit, a slope of 2.2 has been found. However, this relation appears to show a broken power-law scaling with a much shallower slope of 0.56 for M_⋆ > 0.3 M_⊙ (and a much steeper slope of 4.3 below 0.3 M_⊙) (Manara et al. 2017; Alcalá et al. 2017). A direct comparison is therefore difficult at this point. Nonetheless, Fig. 6 shows no broken scaling. However, at the high-mass end, the slope that we find (1.4) sits between thesingle power-law slope (2.2) and the high-mass end of the broken power-law slope (0.56) of Manara et al. (2017).

4.3 The scaling of the dust mass with the stellar mass

In terms of how the dust mass scales with the stellar mass, the situation is slightly more complicated. The dust disc mass as a function of stellar mass is shown in Fig. 7. At 2 Myr, we find $M_{\text{d}}\propto M_{\star }^{4.1}$ if all discs with a dust mass of less than 10⁻² M_⊕ are excluded from the fit. This is the criterion we have previously used to define when a disc is drained. If instead we exclude discs less massive than 10⁻¹ M_⊕, we find $M_{\text{d}}\propto M_{\star }^{3.8}$. Hence, the slope does not appear sensitive to the cut off we use.

If alternatively we look at when the cluster is 1.5 Myr old and apply the same criteria, we find $M_{\text{d}}\propto M_{\star }^{3.1}$ when excluding discs below 10⁻² M_⊕. When the cluster is 1 Myr old we find $M_{\text{d}}\propto M_{\star }^{1.4}$.

Comparing this to observations, the slope we find is steeper than expected. Pascucci et al. (2016) found $M_{\text{d}}\propto M_{\star }^{1.1-1.6}$ in the 1–2 Myr-old Taurus star-formingregion, $M_{\text{d}}\propto M_{\star }^{1.1-1.8}$ in the 1–3 Myr-old Lupus region, $M_{\text{d}}\propto M_{\star }^{1.3-1.9}$ in the 2–3 Myr-old Chamaeleon I, and finally $M_{\text{d}}\propto M_{\star }^{1.9-2.7}$ in the 10 Myr-old Upper Sco region. Similarly, Barenfeld et al. (2016) found a scaling in Upper Sco of $M_{\text{d}}\propto M_{\star }^{1.67\pm 0.37}$. Pascucci et al. (2016) suggest that the difference in slope between the first three regions is due to the lower sensitivity of the survey of Taurus, and that the scaling is similar for star-forming regions between 1 and 3 Myr. Only when our cluster is 1 Myr old does thescaling we find agree with observations.

Our findings and those of Pascucci et al. (2016) both suggest that the scaling between stellar mass and dust disc mass becomes steeper as the cluster becomes older. Although, in our model the steepening appears to occur significantly faster. Pascucciet al. (2016) suggested that the steepening of the relation can be explained when dust growth is fragmentation limited, as this leads to more rapid dust clearing around low-mass stars. Similarly, in our model we find that discs from lower-mass cloud cores, and consequently around lower-mass stars, are more quickly drained of dust. This leads to a steepening of the scaling relation. Our findings of lower-mass stars having dust discs that live shorter do however seem to contradict observations that report a larger fraction of discs around lower mass stars than around higher mass stars in Upper Sco and the 5–12 Myr old Collinder 69 cluster (Carpenter et al. 2006; Bayo et al. 2012).

Fig. 7 Dust disc mass as a function of the stellar mass at 1 Myr (black) and 2 Myr (red). The dashed blue and grey lines show fits to the data. The scaling relation steepens with time from a slope of 1.4 at 1 Myr and 4.1 at 2 Myr.

5 Discussion

5.1 Observed disc masses

In the Monte Carlo parameter exploration we find dust disc masses up to ≈ 2.3 × 10³ M_⊕ and gas disc masses up to ≈0.7 M_⊙. The most massive disc comes from the most massive cloud core with the highest angular momentum budget at the end of its formation phase, when the disc is 0.34 Myr old. At a cluster age of 1 Myr we find dust masses ranging from zero, i.e. a completely drained disc, to ≈1.3 × 10³ M_⊕, and gas masses between ≈3 × 10⁻³ M_⊙ and ≈ 0.4 M_⊙. Once the cluster is 2 Myr old, many dust discs have started to drain and several of the lowest mass discs are completely drained. At this age we find dust masses up to ≈1 × 10³ M_⊕ and gas masses are in the range ≈4 × 10⁻³ M_⊙ to ≈ 0.3 M_⊙.

Observational measurements from CO and HD put gas disc masses in the range of ≈ 10⁻⁵ M_⊙−10⁻¹ M_⊙ (Bergin et al. 2013; Ansdell et al. 2016; McClure et al. 2016; Miotello et al. 2017). However, it should be noted that gas masses measured from CO are very uncertain; see e.g. Fig. 3 in Manara et al. (2016b). This puts the lower limit of our gas masses a few orders of magnitude above the lowest observed values and our highest value a factor of a few higher than the observed sample.

We can now identify at what ages the mean dust mass in our cluster of discs matches the mean dust mass in these observed star-formingregions. In both the 1–3 Myr-old Lupus and the 1–2 Myr-old Taurus star-forming regions, dust disc masses range between ≈0.1 M_⊕ and ≈ 100 M_⊕, and the average mass is ≈15 M_⊕ (Ansdell et al. 2016). A smaller average of 5 M_⊕ was found in the 10 Myr-old Upper Sco star-forming region. The average dust disc mass as a function of time in our population synthesis model is shown in Fig. 8. When including all dust discs regardless of their mass we find that our cluster has an average dust mass of 15 M_⊕ at an age of 2.5 Myr, which is the average in Lupus and Taurus. This agrees well with the estimated ages of 1–3 Myr for Lupus and 1–2 Myr for Taurus. We find an average dust mass of 5 M_⊕ at 2.9 Myr. If we exclude discs with masses <10⁻² M_⊕ we find an average dust disc mass of 15 M_⊕ at 3.3 Myr, which again agrees well with observations. Figure 8 also shows that the mean dust mass declines from ≈ 200 M_⊕ at 1 Myr to ≈3 M_⊕ at ≈ 3 Myr. This relatively rapid decline does not agree well with the retention of dust in the 10 Myr-old Upper Sco region, which has a mean mass of 5 M_⊕. However, as with the scaling relations, this result is sensitive to the age spread given to the disc. An age spread of 2 Myr instead of 1 Myr gives a mean dust mass of15 M_⊕ at 3.5 Myr when all discs are included, 1 Myr later than when the age spread is 1 Myr.

There are some indications that protoplanetary dust disc masses may be underestimated. For example, Galván-Madrid et al. (2018) investigated the effect of self-obscuration in protoplanetary discs due to them being partially optically thick. These latter authors find that protoplanetary disc masses may be underestimated by a factor of ten or more. Recent work by Ballering & Eisner (2019) also showed that protoplanetary disc masses may be underestimated. These authors used radiative transfer models to fit the SED of Taurus discs and found that masses may be one to five times larger than previously thought. Another recent study has shown that dust scattering is capable of reducing the emission from an optically thick disc, making it appear optically thin, and therefore leading to an underestimation of its mass (Zhu et al. 2019). Furthermore, the previously mentioned discrepancy between the masses of protoplanetary discs and the masses of exoplanetary systems also hints at a possible underestimation of protoplanetary disc masses (Najita & Kenyon 2014; Manara et al. 2018).

Manara et al. (2018) proposed two solutions to the discrepancy between observed disc masses and exoplanet system masses. The first solution was that protoplanetary discs are continuously replenished with dust, meaning that the total mass present in the disc over time is larger than the observed masses. However, it is difficult to imagine dust replenishment without replenishing gas at the same time. Since removing gas is less efficient than removing dust, additional replenishment of gas would likely lead to exceedingly long gas disc lifetimes. The second suggestion by Manara et al. (2018) is that planets form early, when discs are more massive. This solution seems more plausible given the large dust masses we find in young discs.

In Lupus and Taurus, the dust disc masses range between 0.1 and 100 M_⊕. At 2 Myr, the majority of the dust discs that we find are within or close to this range, but the most massive dust disc we find is ≈ 1000 M_⊕. This is significantly higher than in Lupus and Taurus, implying that our model overestimates the angular momentum budget of some of the cloud cores because high-angular momentum budgets lead to larger and more massive discs. It has been suggested that magnetised cloud cores regulate the angular momentum budget such that even with quite different initial conditions, the protoplanetary discs form with similar physical conditions, such as their mass (Hennebelle et al. 2020). If this is the case, cloud cores with high angular momenta do not necessarily produce discs with equivalently large angular momenta. This could possibly prevent the formation of these very high-mass discs.

An alternative explanation for why we find such very high-mass discs could be that a yet-to-be found correlation exists between disc mass and formation time such that stars with large discs typically form before those with smaller discs. Indeed, observations of very young protoplanetary discs report the presence of dust discs as massive as 1000 M_⊕ (Tychoniec et al. 2018), suggesting that our discs that form with 1000 M_⊕ are not unreasonable. However, it is possible that such discs only form in the early stages of the cluster lifetime and are then significantly depleted of dust by the time the cluster becomes comparable in age to Lupus and Taurus. In contrast, in our model the formation time was randomly drawn between 0 and 1 Myr, allowing large discs to from relatively late.

One could also speculate that particle growth could be more efficient in these more massive and larger discs, which is possible because the growth of particles in the outer dilute regions of these discs is not well understood. More efficient growth creates larger dust particles that drift faster, resulting in a more rapid depletion of the dust disc. However, such a process would have to apply only to the largest discs.

Fig. 8 Mean dust disc mass as a function of cluster age. The black curve shows the mean mass when including all dust discs and the red curve shows the mean mass when dust discs with masses < 10−2 M⊕ are excluded.

5.2 Observed disc radii

In the literature, the outer radius of the dust disc is often calculated as the radius containing some percentage of the total observed flux density, for example 68 or 95% (see e.g. Huang et al. 2018; Facchini et al. 2019). Under the assumption of an optically thin dust disc, the flux density is directly proportional to dust mass as $$F_{\nu }=\frac{B\left({\overline{T}}_{\text{d}}\right){\kappa }_{\nu }}{D^2}{M}_{\text{d}},$$(17)

where $B\left(T_{\text{d}}\right)$ is the Planck function at the average dust temperature T_d, κ_ν is the frequency dependent opacity, and D is the distance to the source. The disc temperature is often taken as T_d = 20 (see e.g. Ansdell et al. 2016), although models where T_d depends on the stellar luminosity have also been used (Pascucci et al. 2016). In this work we simply define the disc radius as the distance from the central star that contains 95% of the total mass.

We first look at the reference case (Sect. 3.1). Figure 9 shows the stellar gas accretion rate as a function of the dust disc radius. This relation shows a similar behaviour to the stellar gas accretion rate versus dust disc mass shown in Fig. 2, but with a few notable differences. During the formation phase there are two regimes. One during the first half of the formation phase where both the accretion rate and the dust disc radius increase as the disc builds up, and one in the second half where the radius increases but the accretion rate remains almost constant. After the formation phase is over, the accretion rate begins to decrease, while the dust disc radius still increases because it is strongly coupled to the viscously evolving gas. At t ≈ 7.5 × 10⁵ yr, the dust in the outer parts of the disc begins to drift rapidly and the disc radius starts to decrease. Comparing this to Fig. 2, we can see that the dust drift turn-over takes place earlier when considering the dust disc radius instead of the dust disc mass. There is a lag between when dust starts drifting rapidly in the outer portions of the disc and when the mass of the disc starts to decrease.

We can also look at the stellar gas accretion rate versus dust disc radius for the cluster of stars. This is shown in Fig. 10. The black circles show the cluster at 1 Myr and the red circles show the cluster at 2 Myr. In general, a larger disc radius corresponds to a higher gas accretion rate. This is because larger discs typically come from more massive cloud cores. Thesecores form more massive discs with higher gas surface densities, which leads to higher gas accretion rates. Discs with very high accretion rates are still in their formation phase.

This range of values for dust disc radius line up well with observational measurements of the Lupus and Taurus clusters, that place dust disc radii between a few to a few hundred astronomical units (Ansdell et al. 2018; Long et al. 2019).

For a given disc, observations show that the gas disc has a radius that is typically 1.5–8 times larger than the dust disc radius (Najita & Bergin 2018). Figure 11 shows the gas disc radius as a function of the dust disc radius for the modelled cluster at 1 and 2 Myr in black and red circles respectively, and the data from Najita & Bergin (2018) are shown as purple crosses. At a cluster age of 1 Myr our results appear to agree well with these observations. These observed discs, drawn from different star-forming regions, have ages that have been estimated to be between 0.6 and 20.9 Myr. We show in Fig. 11 our results at 1 and 2 Myr. At a cluster age of 1 Myr our results appear to agree well with these observations. However, at 2 Myr the gas disc radii in our cluster are all larger than the observations by a factor of a few. One explanation for why we find much larger discs is that we do not include external photoevaporation in our model, which can play a significant role in truncating or clearing gas discs entirely (e.g. Winter et al. 2018; Nicholson et al. 2019; Concha-Ramírez et al. 2019). The very large disc radii we find can also indicate, as the very high disc masses we find also do, that the angular momentum budget of the cloud cores in our model are overestimated, since this leads to a larger centrifugal radius. During the collapse of magnetised cloud cores it is possible that the angular momentum is regulated such that discs form with similar initial radii, which could prevent the formation of discs with very large centrifugal radii (Hennebelle et al. 2020). An additional possibility is that the viscous expansion of our discs is too strong.

Fig. 9 Gas accretion rate as a function of the dust disc radius for the reference case. During the formation phase there is a first regime where the disc radius and accretion rate both increase. In the latter half of the formation phase the disc radius continues to increase, while the gas accretion rate remains almost constant. Once the formation phase is over the accretion rate drops while the dust disc radius still continues to increase. At t ≈ 7.5 × 105 yr the disc radius starts to decrease as the outer parts of the disc begin to drift rapidly. After t ≈ 2.5 × 106 yr the disc is drained of dust.

Fig. 10 Gas accretion rate as a function of the dust disc radius of the cluster of discs at 1 Myr (black) and 2 Myr (red). Larger disc radii generally correspond to higher accretion rates. As the discs evolve, both disc radii and accretion rate decrease.

Fig. 11 Gas disc radius as a function of dust disc radius for the cluster of discs at 1 Myr (black circles) and 2 Myr (red circles). The purple crosses are observational measurements taken from Najita & Bergin (2018). The blue line shows Rg = Rd. At 1 Myr our results agree well with observations, but at 2 Myr the gas disc radii that we find are a factor of two to three higherthan observations.

5.3 Stellar gas accretion rates

In this study we find that it is difficult to reach low accretion rates (≲ 10⁻⁹ M_⊙ yr⁻¹) within a fewmillion years for a large range of disc masses and centrifugal radii; such low accretion rates are nevertheless often observed (see e.g. Manara et al. 2016b). One way to attempt to remedy this is to decrease the value of the viscosity parameter α_min, because thislatter directly sets the gas accretion rate for a given surface density of gas. However, decreasing α_min also decreases the speed of the viscous spreading, which in turn keeps the gas surface density higher for longer, increasing the gas accretion rate. Decreasing α_min by an order of magnitude to α_min = 10⁻³ in the reference case, we find that the gas accretion rate at the end of the 5 Myr simulation is nearly the same as when α_min = 10⁻², but that the gas disc is roughly twice as massive. This is discussed in more detail in Appendix B.

If some of the accreting gas near the star is removed from the inner disc through some other physical mechanism, such as for example photoevaporation of the disc or disc winds, the measured stellar gas accretion rate drops. Indeed, work by Alexander & Armitage (2009) and Kimura et al. (2016) demonstrates that photoevaporation of the disc can already be important after the first million years of disc evolution when stellar gas accretion rates are in the range 10⁻⁹−10⁻⁸ M_⊙ yr⁻¹. The mechanism removing gas would have to be efficient enough to remove enough mass such that the accretion rate drops by an order of magnitude or more, down to less than 10⁻¹⁰ M_⊙ yr⁻¹, as seen for example in Manara et al. (2016b).

Considering the stellar gas accretion rates and gas masses found in the evolution of the stellar cluster, we generally find gas disc lifetimes that are very long. Crudely estimating the remaining gas accretion timescale as the ratio between the gas accretion rate and the gas disc mass at 3 Myr yields remaining lifetimes of between 3 and 6 Myr in the cluster of discs. However, since the gas accretion rate continues to decrease, the remaining lifetime of the disc is underestimated. Our model also lacks gas dissipation channels other than accretion onto the star. Such methods of dissipating gas from the disc include photoevaporation (Owen et al. 2010, 2011, 2012) and accretion onto planets. If photoevaporation were to occur at a rate of around 10⁻⁹ M_⊙ yr⁻¹ then this process would begin to be important in the disc with the lowest gas masses that we find at 3 Myr, reducing the remaining lifetime. In a model similar to ours, but including photoevaporation, typical gas disc lifetimes between 2 and 4 Myr have been found (Kimura et al. 2016).

Fig. 12 Contour plot of the dust-to-gas ratio over time and disc radius for the reference case. For the first ≈ 0.8 Myr the dust-to-gas ratio remains at the initial ISM value of 1/100. After this the dust-to-gas ratio begins to increase in the inner 100 au of the disc due to the pile up of dust caused by rapid radial dust drift. The increase is strongest in the inner ≈ 5 au. The pile up occurs because dust drift is faster in the outer parts of the disc compared to the inner parts.

5.4 Dust pile up

During the rapid dust drift phase we find that dust piles up, increasing the local dust-to-gas ratio defined as $$\varepsilon ={\Sigma }_{\text{d}}/{\Sigma }_{\text{g}},$$(18)

above the initial value of ɛ_in = 0.01 at orbital radii out to 100 au. The evolution of the dust-to-gas ratio is shown in Fig. 12, which is a contour plot where the colour indicates the dust-to-gas ratio as a function of time and orbital radius of the reference case. We can see that in the range ≈ 1−2 Myr, the dust-to-gas ratio increases by a factor of several in the inner 10 au of the disc, and by a factor of two out to almost 100 au. In the inner few astronomical units we reach dust-to-gas ratios as high as 0.045.

This pile up of dust occurs because drift rates for fixed-radius particles are such that rapid drift begins furthest out in the disc. This leads to an increase in the dust surface density within this inwards-expanding front of rapid dust drift. A similar pile-up for constant-sized dust was also reported by Youdin & Shu (2002). The pile up is driven by the same physical mechanism in both models. However, our models differ, because we use an evolving gas disc, which calculates the outer disc radius self-consistently and allows for viscous spreading of dust.

The increase in dust-to-gas ratio is favourable for the formation of planetesimals. The streaming instability path for forming planetesimals can become active if the dust-to-gas ratio is greater than some critical value (Carrera et al. 2015; Yang et al. 2017). This critical dust-to-gas ratio ranges from ɛ_c ≈ 0.035 at τ_s = 10⁻³ to a lower value of ɛ_c ≈ 0.014 at τ_s = 10⁻¹, after which ɛ_c increases again; see Eqs. (8) and (9) of Yang et al. (2017). Applying these criteria to the dust-to-gas ratio from the reference case seen in Fig. 12, we find the results shown in Fig. 13. Favourable conditions for streaming instability are present starting at orbital radii of ≈100 au at ≈0.8 Myr. This streaming instability-favourable region then expands inwards to ≈ 8 au at ≈1.5 Myr. After the region reaches the innermost orbital radius it quickly disappears as the dust is drained entirely from the disc. Therefore, after 1 Myr of evolution the assumption of a canonical dust-to-gas ratio of 0.01 no longer holds.

We can also analyse whether we have conditions for the continued growth of planetary embryos into full-sized planets. To investigate this we carry out a simple test where we place planetary embryos of mass M = 0.1 M_⊕ at orbital radii and times where we have streaming instability conditions. We then let these cores grow via pebble accretion in the 2D Hill accretion regime (Lambrechts & Johansen 2014). Planet migration is not included in this test. One core is placed at an orbital radius of r = 30 au at t = 1.2 Myr and the other at an orbital radius of r = 10 au at t = 1.5 Myr. We find that the core at r = 30 au grows to 10 M_⊕ at t ≈ 1.6 Myr. At t ≈ 2 Myr it has grown to 17 M_⊕. After this point the dust outside of the planets orbital radius is nearly completely drained and growth is halted. The core at r = 10 au is able to reach a mass of 3.2 M_⊕ before the dust disc is drained. From the time the planet started growing, it takes ≈ 0.5 Myr for the dust to drain at this orbital radius.

Tanaka & Tsukamoto (2019) explored the possibility of planets growing by pebble accretion during the early Class 0 and Class 1I YSO phase, which lasts about 0.5 Myr. These authors found that a planet embryo starting at an orbital radius of 10 au is able to grow to a mass of 10 M_⊕ during this early phase. While we do not find conditions for planetesimal formation this early in the reference case, which we take as a prerequisite for forming the planetary embryos, we find in Appendix C that using dust with a size of 1 mm gives conditions that are favourable of streaming instability already at 0.35 Myr. These findings, combined with those of Tanaka & Tsukamoto (2019) that planet formation could occur quickly during this phase and with the observed difference in dust disc masses and exoplanetary system masses of Najita & Kenyon (2014) and Manara et al. (2018), suggest that planet formation is rapid and possibly already takes place during the early stages of disc lifetimes, when dust discs are massive. Signs of grain growth to mm sizes already in discs of ≈ 10⁵ yr old have also been found, suggesting that planet formation takes place early in the lives of protoplanetary discs (Harsono et al. 2018).

Fig. 13 Contour plot of the ratio between the dust-to-gas ratio and the critical dust-to-gas ratio needed for planetesimal formation via the streaming instability in the reference case. The dust-to-gas ratio starts to exceed the critical value at ≈ 100 au at 0.8 Myr. Over the next ≈0.7 Myr this super-critical region expands inwards to an inner edge of ≈8 au. The region disappears as the dust is drained from the disc.

5.5 Metallicity

In our model, we assume that the initial dust-to-gas ratio (metallicity) of the cloud core is ɛ = 0.01 in order to match ISM values. This value is also quite close to the solar mass fraction of metals Z_⊙ = 0.0134 (Asplund et al. 2009).

In Appendix E we look at the effects of varying the initial dust-to-gas ratio in our model, using the metallicity in the solar neighbourhood and in nearby star-forming regions to set the distribution of the initial dust-to-gas ratio. Variations in this initial dust-to-gas ratio will have some effects, such as changing the initial dust disc mass and the evolved dust-to-gas ratio. This will in turn affect when streaming instability conditions are reached, with higher initial dust-to-gas ratios resulting in these conditions being met more easily. We find that introducing a scatter in the initial dust-to-gas ratio that is consistent with the scatter of metallicities in nearby star-forming regions introduces a small scatter of the dust disc masses, but the effect is very small.

6 Conclusions

Here, we investigate the formation and evolution of protoplanetary gas and dust discs using a 1D numerical model. We performed a parameter study to explore how the centrifugal radius of the disc and the mass of the cloud core affect the disc evolution. We then used population synthesis to model the evolution of a population of discs in a forming stellar cluster for several million years. We looked into how the dust disc mass evolves over time. A study of the scaling relation between disc and stellar properties was also carried out.

We find that the evolution of the dust component in a protoplanetary disc can be separated into three distinct phases, assuming dust has an almost constant size. The first is the formation phase during which both the stellar gas accretion rate and the disc mass increase. This phase lasts on the order of ≈ 10⁴− to × 10⁵ yr, depending on the mass of the cloud core. The second regime is the slow-dust-drift phase. In this phase the gas surface density is high, which keeps the Stokes number of the dust low and therefore radial dust drift is slow. This results in a very steep dependence of the gas accretion rate on the dust disc mass. The third regime begins once the gas surface density drops sufficiently that dust drift becomes rapid; at this point the dust disc mass decreases on a timescale that is much shorter than the timescale over which gas accretion decreases. This results in an almost flat scaling between gas accretion rate and dust disc mass. For dust particles of 100 μm in size, completely clearing the protoplanetary disc of dust takes between 2 and 4 Myr, depending on the mass of the cloud core from which the disc formed and on the angular momentum budget of the cloud core.

We find that when using viscous evolution alone, it is difficult to reach gas accretion rates that match observed values, namely of Ṁ < 10⁻⁹, for the range of cloud core angular momentum and mass we consider. This is probably due to our neglect of disc winds, photoevaporation, and gas-accreting planets. We also find that the most massive discs in our population synthesis model are significantly more massive and larger than the most massive and largest observed protoplanetary discs in star-forming regions of similar ages. We suspect this could be due to our model overestimating the occurrence rate of protoplanetary discs with large angular momentum budgets (Hennebelle et al. 2020). However, observations of very young discs show that such very massive discs are present (Tychoniec et al. 2018), indicating that the very massive discs we find are younger than the youngest disc in the star-forming regions of similar cluster ages.

When evolving a whole cluster of discs we find that the average dust disc masses agree well with observations, and that the low observed average dust disc mass is consistent with radial drift of pebbles. However, we also note that studies indicate that dust disc masses may be higher than previously thought (Galván-Madrid et al. 2018; Ballering & Eisner 2019; Manara et al. 2018).

For the cluster of discs we find a near-linear relation between the stellar gas accretion rate and the gas disc mass, with $\dot{M}\propto M_{\text{g}}^{0.9}$ between 2 and 3 Myr. We find that a simple power law does not accurately describe the dependence of the gas accretion rate on the dust disc mass. The age spread of discs in the cluster will affect how many discs are in the regime of either slow or fast radial dust drift and in the first of the two regimes the scaling is very steep while in the second it is very flat. Since both regimes are somewhat populated during the first few million years of evolution, a single power law does not fit well.

We also found scaling relations of the disc mass as a function of the stellar mass. Over time this scaling becomes steeper and takes the value $M_{\text{d}}\propto M_{\star }^{1.4-4.1}$ between 1 and 2 Myr. This steepening with time qualitatively agrees with observations (Pascucci et al. 2016), but the steepening in our model is faster than observations suggest.

We find in all our simulations that the drifting dust piles up in the inner regions of the disc. The pile up occurs when the disc is 1–2 Myr old. During the pile up, the dust-to-gas ratio increases in the inner 100 au of the disc, reaching values four times higher than the initial value (1/100). The dust-to-gas ratio becomes high enough that, according to the criteria from Yang et al. (2017), the streaming instability can become active at orbital radii between 10 au and 100 au, between 0.8 and 1.5 Myr.

In the present study, we use modelling of the formation and evolution of protoplanetary discs to show that for dust particles of constant size, the protoplanetary dust discs go through three distinct phases. The build up of the disc, the slow radial drift regime, and the fast radial drift regime. Using a population synthesis model we are able to show that the low observed dust disc masses in star-forming regions are consistent with the rapid radial drift of pebbles. Dust discs can be very massive when they are young and still agree with the low observed average masses at their later evolutionary stages. This implies that planets form rapidly during the early, massive phase of the dust disc evolution.

Acknowledgements

We want to thank the anonymous referee for helpful comments that helped to improve the original manuscript. We also want to thank Carlo Manara, Giovanni Rosotti and Per Bjerkeli for their helpful discussions, and Paul McMillan for his help. J.A. and M.L. thank the Wallenberg Academy Fellow (grant 2017.0287). A.J. thanks the Swedish Research Council (grant 2018-04867), the Knut and Alice Wallenberg Foundation (grant 2017.0287) and the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) for research support.

Appendix A Cloud core collapse

The velocity at which the collapse front expands outwards can be found from the equation of the collapse time, which Takahashi et al. (2013) found to be $$t(r_{\text{cf}})=\sqrt{\frac{r_{\text{cf}}^3}{2GM\left(r_{\text{cf}}\right)}}{\displaystyle {\int }_0^1\frac{\text{d}R}{\sqrt{\frac{1}{f}\ln R+\frac{1}{R}-1}}},$$(A.1)

where r_cf is the radius of the collapse front, G is the gravitational constant, M(r_cf) is all the mass contained within the collapse front radius, and R = r_cf∕r_cc where r_cc is the radius of the cloud core. Combining all the constants into one as $$\gamma =\sqrt{\frac{1}{2G}}{\displaystyle {\int }_0^1\frac{\text{d}R}{\sqrt{\frac{1}{f}\ln R+\frac{1}{R}-1}}},$$(A.2)

we can rewrite Eq. (A.1) as $${\displaystyle {\int }_0^{r_{\text{cf}}}4}\pi r^2\rho (r)\text{d}r=\frac{r_{\text{cf}}^3}{t^2}{\gamma }^2.$$(A.3)

Differentiating both sides of this equation with respect to r_cf gives $$4\pi r_{\text{cf}}^2\rho (r_{\text{cf}})=\frac{3r_{\text{cf}}^2{\gamma }^2}{t^2}-\frac{2r_{\text{cf}}^3{\gamma }^2}{t^3}\frac{\text{d}t}{\text{d}{r}_{\text{cf}}}.$$(A.4)

Solving this for dt∕dr_cf and inverting it gives $$\frac{\text{d}{r}_{\text{cf}}}{\text{d}t}=\frac{2r_{\text{cf}}{\gamma }^2}{3t{\gamma }^2-4\pi t^3\rho (r_{\text{cf}})}.$$(A.5)

Appendix B The effect of changing α_min

Increasing or decreasing the minimum value, α_min, that the viscosity parameter α can take affects the rate of viscous evolution. To find out the effect of changing α_min we run a case that is identical to the reference case, except that α_min = 0.001. Because the transport of gas is less efficient for lower viscosities, a lower value for α_min will lead to a disc with a smaller radius, and hence a higher gas surface density, if the disc mass is the same. The lower viscosity will decrease the gas accretion rate compared to the reference case, but the increased gas surface density will act to increase the gas accretion rate. Therefore, the combined effect of decreasing α_min on the gas accretion rate is not immediately clear. The matter becomes further complicated when gravitational instability is taken into account. The more compact disc from a low α_min case can become more gravitationally unstable due to the same mass being confined to a smaller region.

The result of using α_min = 0.001 is shown in Fig. B.1, which shows the gas accretion rate as a function of disc mass. During the formation phase, this appears to be very different from the reference case in Fig. 2. In the very early stages, the gas accretion rate and disc mass increase in a similar way as the reference case, but at a lower accretion rate. However, once the disc becomes gravitationally unstable the accretion rate increases very quickly. This dramatic change occurs because when α_min = 0.001 the disc does not viscously expand significantly. Therefore this disc becomes more gravitationally unstable than the reference case, and the effective α becomes much higher. The highest α value in the disc during the phase of gravitational instability is nearly two orders of magnitude higher when α_min = 0.001 compared to when α_min = 0.01.

Fig. B.1 Stellar gas accretion rate as a function of disc mass using αmin = 10−3. The greylines show the reference case dust mass (solid line) and gas (dashed lines). Decreasing αmin results in the disc becoming more gravitationally unstable. This causes the very large increase in the accretion rate in the second half of the formation phase. The accretion rate then drops rapidly until about 0.5 Myr. After this it decreases more slowly. The dust drift turn-over occurs later with lower αmin and the dust disc drains slower. The final gas accretion rate is similar for αmin = 10−3 and αmin = 10−2, but the disc masses of the αmin = 10−3 case are higher.

Once the formation phase is over, the gravitational instability phase also ends quickly and the gas accretion rate now drops very quickly at first, dropping four orders of magnitude in the first 3 × 10⁵ yr after the end of the formation phase. After this, the evolution of the accretion rate slows down. The dust drift turn-over is reached at 3.5 Myr, instead of at 1.3 Myr as for the reference case. There is also several tens of Earth masses left in the disc at 5 Myr compared to the reference case that was completely drained in 2.5 Myr. The gas mass is about twice as high at the end of the simulation in the α_min = 0.001 compared to the reference case.

The final gas accretion rate at 5 Myr is not significantly different when using α_min = 0.001 and α_min = 0.01. Both have accretion rates of Ṁ ≈ 10⁻⁸ M_⊙ yr⁻¹ at this point.

Appendix C The effect of changing the dust particle size

In this section we analyse the effect of changing the size of the dust particles. For this we evolve a system that is identical to the reference case, but using a dust size of r_p = 1 mm, compared to r_p = 100 μm in the reference case. Figure C.1 shows the stellar gas accretion rate as a function of disc mass. Increasing the dust size does not change the evolutionary path of the dust disc, but it speeds it up. Instead of draining in ≈ 2.5 Myr as for the reference case, it drains in ≲0.8 Myr. The onset of the rapid dust drift regime is earlier for larger dust because the Stokes numbers are higher. Since the rapid drift begins earlier, the gas accretion rate during the dust draining phase is higher.

Fig. C.1 Stellar gas accretion rate as a function of disc mass using a dust size of rp = 0.1 cm. The red curve shows the dust and the black curve shows the gas. The black dots mark points in time as labelled. The grey line shows the dust disc mass of the reference case. The rapid draining of the dust mass starts at about 0.5 Myr. Within 1 Myr the dust disc has been drained to the point where only 0.017 × 10−2 M⊕ of dust remains.

Fig. C.2 Contour plot of the ratio between the critical dust-to-gas ratio and the critical dust-to-gas ratio for the streaming instability. Using dust with a size of 1 mm the critical dust-to-gas ratio passed after ≈ 3 × 105 yr. The super-critical region remains until t ≈ 7 × 105 yr.

The dust accretion rate through the disc is also higher for 1 mm sized particles. This in turn leads to a larger increase in the dust-to-gas ratio during the phase of rapid dust drift than for the reference case. This increase also occurs earlier. The larger and earlier pile-up of dust makes this dust reach the critical dust-to-gas ratio for streaming instability earlier. This can be seen in Fig. C.2. A super-critical region exists from t ≈ 0.35 Myr to t ≈ 0.7 Myr at radii between r ≈ 2 au and r ≈ 70 au. For the reference case the super-critical region exists rom t ≈ 0.8 Myr to t ≈ 1.5 Myr at radii between r ≈ 8 au and r ≈ 100 au.

Appendix D a_c − M_cc relation

The final centrifugal radius is given by $$a_{\text{c}}=\frac{{\Omega }_0^2{r}_{\text{cc}}^4}{G{M}_{\text{cc}}},$$(D.1)

where Ω₀ is the rotation rate of the cloud core, r_cc the cloud core radius, and M_cc the cloud core mass. The ratio of rotational and gravitational energy in the cloud core is given by $$\beta =\frac{E_{\text{rot}}}{E_{\text{grav}}}=\frac{\frac{1}{2}I{\Omega }_0^2}{q\frac{G{M}_{\text{cc}}^2}{r_{\text{cc}}}},$$(D.2)

where I is the moment of inertia of the cloud core and q is a parameter related to the density profile of the cloud core, e.q. q = 3∕5 for a constant density sphere. β is thought to be independent of the size of a cloud core (Goodman et al. 1993; Chen & Ostriker 2018). We therefore assume β to be the same for all cloud cores. From Eq. (D.2) we find ${\Omega }_0^2$ as $${\Omega }_0^2\propto \frac{\beta \frac{M_{\text{cc}}^2}{r_{\text{cc}}}}{I}.$$(D.3)

The moment of inertia, I, of a sphere with varying density can be calculated as $$I=\frac{2}{3}M\frac{{\displaystyle {\int }_0^R{r}^4}\rho (r)\text{d}r}{{\displaystyle {\int }_0^R{r}^2}\rho (r)\text{d}r},$$(D.4)

where M is the mass of the sphere, ρ(r) its density, and R its radius. For a radial density profile that follows ρ(r) ∝ r^a where a≠ − 3 the above relation gives $$I\propto M{r}^2.$$(D.5)

Inserting this into Eq. (D.3) gives $${\Omega }_0^2\propto \frac{\beta \frac{M_{\text{cc}}^2}{r_{\text{cc}}}}{M_{\text{cc}}{r}_{\text{cc}}^2}=\frac{\beta M_{\text{cc}}}{r_{\text{cc}}^3}.$$(D.6)

If we insert this into Eq. (D.1), using that β is not proportional to the cloud core radius or mass we get $$a_{\text{c}}\propto r_{\text{cc}}.$$(D.7)

For a Bonnor-Ebert sphere of a given temperature the radius of the sphere is linearly proportional to the mass, giving $$a_{\text{c}}\propto M_{\text{cc}}.$$(D.8)

Appendix E The effect of adding a metallicity spread

Here we look at how increasing or decreasing the initial dust-to-gas ratio (i.e. that of the cloud core) affects the disc evolution. We begin by looking at two cases, one where the initial dust-to-gas ratio has been halved and another where it has been doubled, in order to understand how this affects the dust disc. We then repeat the population synthesis of a star-forming region as in the main study, but introduce a scatter in the initial dust-to-gas ratio. The amount of scatter we introduce is based on the scatter in metallicity of nearby star-forming regions. This is quite low, and therefore including variation in the initial dust-to-gas ratio only introduces a small scatter in dust disc masses into our results.

Figure E.1 shows the gas accretion rate as a function of disc mass for the two cases with halved and doubled initial dust accretion rate in blue and green respectively. The reference case is plotted in grey. From this we can see that halving or doubling the initial dust-to-gas ratio halves or doubles the dust disc mass, respectively. Since the gas disc is unaffected by the dust in out model, the dust drift turn-over is unchanged.

Fig. E.1 Gas accretion rate as a function of dust disc mass for three different initial dust-to-gas ratio cases. The grey curve shows the reference case, the blue curve shows the case with the initial dust-to-gas ratio half of the reference case, and the green curve shows the case with double the reference case’s initial dust-to-gas ratio. From this we can see that in this relation, increasing or decreasing the initial dust-to-gas ratio by some factor results in an equally large change in the dust disc mass. The dust drift turn-over point and the gas accretion rate are unaffected in our model.

An effect of increasing the initial dust-to-gas ratio is that the backreaction of dust on gas would become stronger. This effect can lead to self-induced dust traps which further increase the concentration of dust (Gonzalez et al. 2017). However, we do not model this backreaction so we would not see such an effect in our model. The dust-growth timescale in coagulation models is also shortened for larger dust-to-gas ratios (Birnstiel et al. 2012), but since we do not include dust growth in our model, it does not affect our findings.

To a good approximation, the metallicity in the solar neighbourhood has a mean that is slightly subsolar with a small scatter [Fe/H] = −0.1±0.1 dex (Schönrich & Binney 2009; Nieva & Przybilla 2008; Bragança et al. 2019). This value is similar to that of stars in nearby star-forming regions where the metallicity is solar or slightly subsolar. Reported values are: [Fe/H] = −0.01±0.04 dex in Orion (D’Orazi et al. 2009), [Fe/H] = −0.01±0.05 dex, in Taurus (D’Orazi et al. 2011), and [Fe/H] = −0.08±0.04 dex in Chameleon (Spina et al. 2014).

Fig. E.2 Gas accretion rate as a function of disc mass for a stellar cluster where the scatter in initial dust-to-gas ratio matches that of the scatter in metallicites in the solar neighbourhood and nearby star-forming regions. Comparing this to when the initial dust to gas ratio is constant, which is shown in grey, (see also Fig. 5) we can see that the two appear very similar.

With this in mind, we rerun the population synthesis model but now we pick the initial dust-to-gas ratio from a normal distribution with a mean of 0.01 and a standard deviation of 0.0025 (corresponding to 0.1 dex). The result of this is shown in Fig. E.2. Comparing this to the population synthesis without a scatter in the initial dust-to-gas ratio shown in grey symbols, the results are very similar. Because the scatter of metallicities we introduce is small, most cloud cores will have an initial dust-to-gas ratio within a factor of two from the mean of 0.01. In Sects. 3.2 and 3.3 we showed that variations in the centrifugal radius have a similar or larger effect on the dust disc mass for a given cloud core mass. Hence, the effect of changing the initial dust-to-gas ratio is indistinguishable from the effect of variations in the centrifugal radius when one only looks at the total mass of the dust disc.

References

All Tables

All Figures

Fig. 1 Snapshots of gas (solid lines) and dust (dashed lines) surface densities of the reference case. In the evolution phase, the inner disc is drained of gas as it is accreted by the star and the outer disc expands as gas is transported outwards to conserve angular momentum. In the first ≈1.0 Myr the dust surface density mimics that of the gas in most of the disc out to radii of about 100 au. Beyond this radius, radial drift of dust is already rapid, draining the dust disc. When the disc is between 2 and 2.5 Myr old, the dust disc is completely drained.

In the text

Fig. 2 Stellar gas accretion rate as a function of disc mass for the reference case. The black curve shows the dust disc mass and the red curve shows the gas disc mass. The black dots mark the age in years as labelled. There are three clear regimes in this relation. During the formation phase of the disc, the gas accretion rate and the disc dust mass are both increasing. Once the formation phase has ended, dust drift is slow and the gas accretion rate drops faster than the disc dust mass. Rapid draining of the dust disc mass begins around t ≈ 1.3 × 106 and lasts about 1 Myr.

In the text

Fig. 3 Stellar gas accretion rate and dust disc mass as a function of time, for discs with final centrifugal radii of 8, 16, 32, 64, and 128 au forming from 1 M⊙ cloud cores. The red and yellow lines show the tracks of the two cases with the highest and lowest final centrifugal radius, respectively. The formation phase ends at t ≈ 1.8 × 105 yr. During the formation phase a trend can be seen where systems with smaller centrifugal radii have lower dust disc masses and slightly higher gas accretion rates. Between 1.8 × 105 yr and 5.0 × 105 yr the gas accretion rates drop significantly, while dust disc masses decrease only slightly. Discs with a smaller final centrifugal radius experience a larger drop in accretion rate. A lower final centrifugal radius leads to faster evolution of the protoplanetary disc. At 2.5 × 106 yr, only the case with the highest final centrifugal radius remains, but the dust disc mass is only ≈ 10−2 M⊕.

In the text

Fig. 4 Gas accretion rate as a function of dust disc mass for protoplanetary discs formed from cloud cores with masses and final centrifugal radii shown in Table 1. The yellow and blue lines trace the time evolution of the cases with the lowest and highest mass, and lowest and highest centrifugal radius, respectively. As a general trend, discs from low-mass cloud cores evolve faster than those from higher-mass cloud cores. However, the speed with which the disc evolves also depends on the final centrifugal radius. At 3 Myr (panel d), all dust discs except for the ones coming from 2 M⊙ cloud cores have been completely drained out.

In the text

Fig. 5 Gas accretion rate as a function of dust (black circles) and gas (red squares) disc mass for protoplanetary discs of the stellar cluster. The dashed black lines show power-law fits of the accretion rate as a function of the gas disc mass ($\dot{M}\propto M_{\text{g}}^{0.9}$). The discs are given a random individual age of between 0.0 and 1.0 Myr to mimic a spread in formation times. Panels from a to d: protoplanetary discs at the cluster ages of 1.0, 1.5, 2.0, and 3.0 Myr. For the first 1 Myr the dust is coupled well to the gas inmost discs, and dust disc masses evolve in a very similar way to the gas disc mass. Some of the lowest mass discs have reached the dust drift turn-over at this point. Beyond this time more and more of the dust discs reach the dust turn-over and begin to drain rapidly. At 1.5 Myr the lower mass discs have reached the turn-over point, as can be seen in the flattening of the low-mass end of the dust discs in panel b. At 2 Myr only the most massive dust discs are yet to reach the turn-over and at 3 Myr nearly all discs are drained of dust. Only 4 out of 100 discs still have a dust mass greater than 1 M⊕, although two of these have a mass of more than 100 Earth masses.

In the text

	Fig. 6 Gas accretion rate as a function of stellar mass of the cluster of discs at 2 Myr (black) and 3 Myr (red). The blue and grey dashed lines are fits to the data at 2 and 3 Myr respectively. Both scaling relations have a slope of ≈ 1.4.
In the text

	Fig. 7 Dust disc mass as a function of the stellar mass at 1 Myr (black) and 2 Myr (red). The dashed blue and grey lines show fits to the data. The scaling relation steepens with time from a slope of 1.4 at 1 Myr and 4.1 at 2 Myr.
In the text

	Fig. 8 Mean dust disc mass as a function of cluster age. The black curve shows the mean mass when including all dust discs and the red curve shows the mean mass when dust discs with masses < 10−2 M⊕ are excluded.
In the text

Fig. 9 Gas accretion rate as a function of the dust disc radius for the reference case. During the formation phase there is a first regime where the disc radius and accretion rate both increase. In the latter half of the formation phase the disc radius continues to increase, while the gas accretion rate remains almost constant. Once the formation phase is over the accretion rate drops while the dust disc radius still continues to increase. At t ≈ 7.5 × 105 yr the disc radius starts to decrease as the outer parts of the disc begin to drift rapidly. After t ≈ 2.5 × 106 yr the disc is drained of dust.

In the text

	Fig. 10 Gas accretion rate as a function of the dust disc radius of the cluster of discs at 1 Myr (black) and 2 Myr (red). Larger disc radii generally correspond to higher accretion rates. As the discs evolve, both disc radii and accretion rate decrease.
In the text

	Fig. 11 Gas disc radius as a function of dust disc radius for the cluster of discs at 1 Myr (black circles) and 2 Myr (red circles). The purple crosses are observational measurements taken from Najita & Bergin (2018). The blue line shows Rg = Rd. At 1 Myr our results agree well with observations, but at 2 Myr the gas disc radii that we find are a factor of two to three higherthan observations.
In the text

Fig. 12 Contour plot of the dust-to-gas ratio over time and disc radius for the reference case. For the first ≈ 0.8 Myr the dust-to-gas ratio remains at the initial ISM value of 1/100. After this the dust-to-gas ratio begins to increase in the inner 100 au of the disc due to the pile up of dust caused by rapid radial dust drift. The increase is strongest in the inner ≈ 5 au. The pile up occurs because dust drift is faster in the outer parts of the disc compared to the inner parts.

In the text

Fig. 13 Contour plot of the ratio between the dust-to-gas ratio and the critical dust-to-gas ratio needed for planetesimal formation via the streaming instability in the reference case. The dust-to-gas ratio starts to exceed the critical value at ≈ 100 au at 0.8 Myr. Over the next ≈0.7 Myr this super-critical region expands inwards to an inner edge of ≈8 au. The region disappears as the dust is drained from the disc.

In the text

Fig. B.1 Stellar gas accretion rate as a function of disc mass using αmin = 10−3. The greylines show the reference case dust mass (solid line) and gas (dashed lines). Decreasing αmin results in the disc becoming more gravitationally unstable. This causes the very large increase in the accretion rate in the second half of the formation phase. The accretion rate then drops rapidly until about 0.5 Myr. After this it decreases more slowly. The dust drift turn-over occurs later with lower αmin and the dust disc drains slower. The final gas accretion rate is similar for αmin = 10−3 and αmin = 10−2, but the disc masses of the αmin = 10−3 case are higher.

In the text

Fig. C.1 Stellar gas accretion rate as a function of disc mass using a dust size of rp = 0.1 cm. The red curve shows the dust and the black curve shows the gas. The black dots mark points in time as labelled. The grey line shows the dust disc mass of the reference case. The rapid draining of the dust mass starts at about 0.5 Myr. Within 1 Myr the dust disc has been drained to the point where only 0.017 × 10−2 M⊕ of dust remains.

In the text

	Fig. C.2 Contour plot of the ratio between the critical dust-to-gas ratio and the critical dust-to-gas ratio for the streaming instability. Using dust with a size of 1 mm the critical dust-to-gas ratio passed after ≈ 3 × 105 yr. The super-critical region remains until t ≈ 7 × 105 yr.
In the text

Fig. E.1 Gas accretion rate as a function of dust disc mass for three different initial dust-to-gas ratio cases. The grey curve shows the reference case, the blue curve shows the case with the initial dust-to-gas ratio half of the reference case, and the green curve shows the case with double the reference case’s initial dust-to-gas ratio. From this we can see that in this relation, increasing or decreasing the initial dust-to-gas ratio by some factor results in an equally large change in the dust disc mass. The dust drift turn-over point and the gas accretion rate are unaffected in our model.

In the text

	Fig. E.2 Gas accretion rate as a function of disc mass for a stellar cluster where the scatter in initial dust-to-gas ratio matches that of the scatter in metallicites in the solar neighbourhood and nearby star-forming regions. Comparing this to when the initial dust to gas ratio is constant, which is shown in grey, (see also Fig. 5) we can see that the two appear very similar.
In the text