© O. Schib et al. 2021

1 Introduction

Planet formation takes place in protoplanetary discs around young stars. Two prevalent formation mechanisms have been proposed. The first is core accretion (CA) (Safronov 1972 and Pollack et al. 1996), where planetesimals, formed from dust, coagulate to form rocky cores and continue to become terrestrial planets or, by accretion of gas from the disc, gas giant planets. The second is gravitational instability (GI) (Kuiper 1951; Cameron 1978; Boss 1997, 1998), where entire regions of protoplanetary discs collapse under their own gravity (fragmentation) to form bound clumps consisting mainly of gas (fragments). Most if not all protoplanetary discs are self-gravitating early in their lives. In this phase, angular momentum transport is dominated by global instabilities (Sect. 2.5.1) or spiral arms, and the discs may fragment (Harsono et al. 2011; Vorobyov & Basu 2005 and Nixon et al. 2018, especially the introduction and references therein). This has a profound effect on properties such as the temperature in the disc during this time – and possibly much later as well. In particular, if protoplanets form during this early phase, their properties will be strongly influenced by the conditions in the disc. In some cases, the disc may fragment and produce bound clumps that could survive to become gaseous planets. Manara et al. (2018) compare masses of observed exoplanets to measured disc masses and found that exoplanetary system masses are comparable to or higher than the most massive discs. They conclude that dust disc masses are underestimated, or that planets form very rapidly, or that discs are being continuously replenished from the environment. Including the formation phase of the discs, when discs are most massive, could explain this apparent lack of planet formation material.

The formation of protoplanetary discs is coupled to the formation of stars: the collapse of molecular cloud cores (MCCs hereafter) leads to the formation of a protostar at the centre, surrounded by a disc of gas (and dust), or possibly multiple protostars. Gravitational instability typically begins at early times in the disc’s evolution (few 10³ to 10⁴ yr), when the MCC collapse is still ongoing (Lodato 2008). This means that if a fragment forms, the material removed from the disc can be replenished by infalling matter from the cloud, which may lead to further fragmentation.

Many planet formation studies begin from a given disc profile and do not consider any infall of envelope material. Initial disc masses are often chosen around a value obtained from the minimum mass solar nebula (MMSN) (Weidenschilling 1977 and Hayashi 1981) of (~ 0.01 M_⊙) or a few times this value (Mordasini 2018; Mutter et al. 2017; Coleman et al. 2017; Hopkins 2016; Lines et al. 2015 and Simon et al. 2015). Some studies do consider infall but use very simple models such as the classic Shu collapse model (Shu 1977) or the Bonnor–Ebert sphere (Ebert 1955 and Bonnor 1956): Hueso & Guillot (2005), Jin & Li (2014), Xiao & Jin (2015) and Kimura et al. (2016). Nixon et al. (2018) make a strong case against using the MMSN prescription and instead introduce the “maximum mass solar nebula”. They argue that the most realistic time to start planet formation simulations is when the disc reaches maximum mass. This typically corresponds to the end of the infall phase. However, as we discuss below, many discs become self-gravitating long before this stage. This typically leads to the transfer of significant amounts of mass and angular momentum in the disc, by global instabilities (Sect. 2.5.1), spiral arms (Sect. 2.5.2), or gravitationally bound clumps, which significantly changes the conditions in the disc. Furthermore, the properties of the host star are also influenced by the accretion of disc material and/or clumps onto the star.

There are several other aspects in the formation process of protoplanetary discs that are often neglected in planet formation studies: first, star formation is a highly turbulent process and it is not self-evident that protoplanetary discs can beadequately described by simple power-law density profiles B18. Second, the infalling material reaches the disc at sub-Keplerian velocity and therefore changes the angular momentum budget of the disc (Hueso & Guillot 2005). Third, the infalling material collides with the disc at some point along its trajectory. Adding the disc material at a distance from the star that corresponds to the streamlines intersecting the disc’s midplane, as done in several studies, is incorrect (Visser & Dullemond 2010). At last, collapsing MCCs are characterised by strong outflows from the inner region of the forming disc (Hartmann & Kenyon 1996). This can change the infall streamlines and prevent matter from falling near or directly onto the star as predicted for example by the Shu collapse model.

In this work we investigated the influence of the infall prescription on the statistical properties of the protoplanetary disc at the end of the infall phase, on their lifetimes, and on the number and properties of fragments formed. We performed a population synthesis of discs. Our model includes the formation of the discs and their evolution until their dispersal. We chose the infall times (duration of the infall phase) such that the resulting population of stars agrees with the initial mass function (IMF) of Chabrier (2007). Our inferred distributions of disc properties can be compared to observations.

The paper is organised as follows: in Sect. 2 we describe the model used in this study. Section 3details the parameter space investigated and the initial conditions used. In Sect. 4 we describe the formation and evolution of a system as an example. In Sect. 5 we present our results, which we compare to other studies and observations in Sect. 6. Section 7 contains a discussion about the influence of our assumptions. Our conclusions are summarised in Sect. 8.

2 Model description

Here we present the numerical model used in this work. After a brief overview we give some details about our treatment of the disc’s evolution, temperature model, viscosity, inner edge, stellar evolution, disc dispersal, and fragmentation.

The core of our simulation consists of a rotationally symmetric disc of gas described by the vertically integrated surface density Σ. The evolution of the surface density is calculated by solving the diffusion equation with an effective viscosity (see Sect. 2.1). The disc is truncated at the inner edge. Matter evolving across this truncation radius is assumed to accrete onto the central star of mass M_*≡ M_*(t), with a fraction considered as outflow (see Sect. 2.6).

For simplicity, in the following we refer to the young stellar object at the centre of the disc as “star”, even though not all such objects qualify for this term in a strict sense. In the early stages of the simulations the central objects are not on the pre-main sequence, and some objects never reach hydrogen fusion.

Our simulations start at a time of 1 kyr after the formation of the central object with an initial disc profile obtained from the data in Bate (2018, hereafter B18), as detailed in Sect. 3. The disc evolves viscously and is fed by infalling material from the MCC as described in Sect. 2.10.

As the disc evolves, the conditions for fragmentation (Sect. 2.9) are continuously monitored. When they are fulfilled, we expect the formation of a bound clump of disc material by gravitational collapse. Therefore, the corresponding clump’s mass is removed from the disc. At this stage of our project, we did not model the clump as an evolving protoplanet. Instead, we followed a similar approach as in Kratter et al. (2008) and added the clump’s mass to the star. The modelling of the clump would entail a number of complex physical processes related to the clump’s evolution such as mass accretion, mass loss, migration and gap formation. These processes are of great importance to the question whether some of these clumps can survive to become giant planets and we hope to address this topic in future research. In this study we focused on the conditions in the disc during and after clump formation, and we therefore assumeed that the clumps migrate inwards on a short time scale and accrete on the central star. Based on the results of hydrodynamic simulations we expect this to be a likely outcome of fragmentation (Baruteau et al. 2011 and Malik et al. 2015). Therefore, the material removed from the disc by fragmentation is simply added to the star after a type-I migration time scale τ₀ (Baruteau et al. 2016): $$\begin{array}{ll}{\tau }_0\hfill & =1.3\text{Myr}\hfill \\ \hfill & \times {\left(\frac{h}{0.05}\right)}^2{\left(\frac{\Sigma }{200\text{g}{\text{cm}}^{-2}}\right)}^{-1}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{3/2}{\left(\frac{r_{\text{p}}}{5\text{au}}\right)}^{-1/2}{\left(\frac{M_{\text{p}}}{M_{\oplus }}\right)}^{-1},\hfill \end{array}$$(1)

with h the disc’s aspect ratio, r_p and M_p the clump’s semi-major axis and mass, respectively. We considered 10% of the clump’s mass as outflow (see also Sect. 2.6). The central star is evolved according to the evolution tracks of Yorke & Bodenheimer (2008) as described in Sect. 2.7. We calculateed the disc’s vertical temperature structure by taking into account a number of physical effects (see Sect. 2.3).

2.1 Evolution of the protoplanetary disc

The evolution of the disc’s surface density Σ is described by the diffusion equation (Lüst 1952 and Lynden-Bell & Pringle 1974): $$\frac{\partial \Sigma }{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}\left(\nu \Sigma r^{1/2}\right)\right]+S.$$(2)

Here, r denotes the distance from the central star. S ≡ S(r, t) is a source/sink term: $$S\equiv S(r,t)=S_{\text{inf}}(r,t)-S_{\text{int}}(r,t)-S_{\text{ext}}(r,t)-S_{\text{frag}}(r,t).$$(3)

S_inf describes the infall from the MCC (see Sect. 2.10), S_int and S_ext the rates of internal and external photoevaporation, respectively (Sect. 2.4) and S_frag the removal of mass due to fragmentation (Sect. 2.9). We solve Eq. (2) using the implicit donor cell advection-diffusion scheme from Birnstiel et al. (2010).

2.2 Auto-gravitation

Since at early stages the discs are massive, their self-gravity cannot be neglected. We followed Hueso & Guillot (2005) in assuming a vertical disc structure that is isothermal ($p=\rho c_{\text{s}}^2$) and in hydrostatic equilibrium, and considered the vertical component of the star’s gravity as well as the local gravitation of the disc: $$\frac{1}{\rho }\frac{\text{d}p}{\text{d}z}=-\frac{G{M}_{*}}{R^3}z-4\pi G\Sigma ,$$(4)

where ρ ≡ ρ(r, z) and p ≡ p(r, z) are the density and pressure in the disc, respectively, G is the gravitational constant. The angular frequency of the disc becomes: $$\Omega (r,t)={[\frac{G{M}_{*}}{r^3}+\frac{1}{r}\frac{\text{d}{V}_{\text{d}}}{\text{d}r}]}^{1/2}.$$(5)

The Keplerian angular frequency therefore receives a modification through the disc’s gravitational potential V_d (Hueso & Guillot 2005): $$V_{\text{d}}(r)={\displaystyle {\int }_{R_{*}}^{\infty }\text{}}-G\frac{4\text{K}[-\frac{4r/r_1}{{(r/r_1-1)}^2}]}{\left|(r/r_1-1)\right|}\Sigma (r_1)\text{d}{r}_1.$$(6)

M_* and R_* denote the stellar mass and radius, respectively. K is the elliptic integral of the first kind.

The solution to Eq. (4) is: $$\rho (r,z)={\rho }_0(r)\exp (-(\frac{\left|z\right|}{H_0}+{(\frac{z}{H_1})}^2)),$$(7)

where ρ₀(r) ≡ ρ(r, 0) is the density in the midplane, H₀ represents the influence of the auto-gravitation of the disc and H₁ the gravitation of the star¹: $$H_0=\frac{c_{\text{s}}^2}{4\pi G\Sigma }$$(8) $$H_1=\frac{\sqrt{2}{c}_{\text{s}}}{\sqrt{\frac{\text{G}{M}_{*}}{r^3}}}$$(9)

with c_s the isothermal sound speed: $$c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu \text{u}}},$$(10)

k_B is the Boltzmann constant, T the midplane temperature, μ the mean molecular weight which we set to 2.3 in this study, and u the atomic mass unit.

For the relationship between ρ(r, z) and Σ(r) the following needs to hold: $$\Sigma (r)={\displaystyle {\int }_{-\infty }^{\infty }\text{}}\rho (r,z)\text{d}z.$$(11)

Requiring $${\rho }_0(r)=\frac{\Sigma (r)}{\sqrt{2\pi }H},$$(12)

we get for the vertical scale height H: $$H=\frac{H_1}{\sqrt{2}}\exp (\frac{H_1^2}{4H_0^2})(1-\text{erf}(\frac{H_1}{2H_0})).$$(13)

Equation (2) was derived assuming a Keplerian orbit. However, with the definition of Ω in Eq. (5) we introduce an angular frequency that is not Keplerian. Therefore, our treatment of the disc’s auto-gravitation on its evolution is not completely self-consistent and should be regarded as an approximation (as in Hueso & Guillot 2005).

2.3 Temperature model

Based on the vertical structure of the disc (Eq. (7)) we considered the following physical processes when determining the disc’s midplane temperature: viscous heating, irradiation from the remaining envelope of the MCC, shock heating from gas infalling on the disc’s surface, irradiation due to the star’s intrinsic luminosity and irradiation due to shock heating from accretion of material onto the star.

Following Nakamoto & Nakagawa (1994) in considering an optically thick as well as an optically thin regime we have an energy balance at the disc’s surface: $$\sigma T_{\text{S}}^4=\frac{1}{2}\left(1+\frac{1}{2\Sigma {\kappa }_{\text{p}}}\right)({\dot{E}}_{\nu }+{\dot{E}}_{\text{s}})+\sigma T_{\text{env}}^4+\sigma T_{\text{irr}}^4,$$(14)

with T_S the surface temperature and we obtain the following expression for the midplane temperature T_mid: $$\begin{array}{ll}\sigma T_{\text{mid}}^4=\hfill & \frac{1}{2}\left[\left(\frac{3\Sigma {\kappa }_{\text{R}}}{8}+\frac{1}{2\Sigma {\kappa }_{\text{P}}}\right){\dot{E}}_{\nu }+\left(1+\frac{1}{2\Sigma {\kappa }_{\text{P}}}\right){\dot{E}}_{\text{s}}\right]\hfill \\ \hfill & +\sigma T_{\text{env}}^4\hfill \\ \hfill & +\sigma T_{\text{irr}}^4.\hfill \end{array}$$(15)

Here, κ_R = κ_R(ρ₀, T_mid) and κ_P = κ_P(ρ₀, T_mid) denote the Rosseland and Planck opacity evaluated at the midplane temperature, respectively, ${\dot{E}}_{\nu }(T_{\text{mid}})=\Sigma \nu (T_{\text{mid}})\left(r\frac{\text{d}\Omega }{\text{d}r}\right)$ the energy dissipation rate due to viscous transport (we do not assume Keplerian frequencies), ${\dot{E}}_{\text{s}}=S_{\text{inf}}{(r\Omega )}^2/2$ the shock heating due to infalling material (Kimura et al. 2016) (K16 from now on). The expression for the infall source term S_inf is given in Sect. 2.10. T_env (set to 10 K in this study) is the background temperature in the envelope, and T_irr the irradiation temperature from the central star. T_irr contains thecontribution of the star’s intrinsic luminosity as well as the luminosity from accreted material (see Sect. 2.7). We used the irradiation model from Hueso & Guillot (2005). Equation (15) is solved numerically using Brent’s method to obtain the midplane temperature T_mid. We used gas opacities from Malygin et al. (2014) and dust opacities from Semenov et al. (2003) to calculate κ_R and κ_P, as in Marleau et al. (2017, 2019)². For the dust opacities we assumed dust grains made of “normal silicates” (NRM model in Semenov et al. 2003) as homogeneous spheres.

2.4 Photoevaporation

The discs are subject to both external and internal photoevaporation. We used an adaptation of the model for external FUV photoevaporation by massive stars from Matsuyama et al. (2003). Our model for internal photoevaporation includes EUV irradiation by the host star. We closely followed Clarke et al. (2001). The sink terms S_ext (r, t) and S_int (r, t) are given in Appendix A.

2.5 Viscosity

We used the α prescription of viscosity (Shakura & Sunyaev 1973): $$\nu =\alpha \frac{c_{\text{s}}^2}{\Omega }.$$(16)

In the calculation of the α-parameter we considered gravitational torques (α_G) as well as torques due to other (“background”, α_bg) processes (see Sect. 2.5.3) and calculate: $$\alpha =\max ({\alpha }_{\text{G}},{\alpha }_{\text{bg}}),$$(17)

where α_G = α_d during infall phase and α_G = α_GI afterwards. The contributions are explained in the following.

2.5.1 Global transport of angular momentum

At the beginning of our simulations, the disc is fed by infalling material and its mass is comparable to that of the host star and therefore the centre of mass can be shifted away from the central star. In this regime the disc can become globally unstable (Harsono et al. 2011). The resulting indirect gravitational potential leads to very efficient transport of mass and angular momentum. Kratter et al. (2010) develop a parametrisation that can be used to simulate this regime. They perform a numerical parameter study of rapidly accreting, gravitationally unstable discs and propose a characterisation of such systems by two dimensionless parameters ξ and Γ. The thermal parameter ξ is defined as $$\xi =\frac{{\dot{M}}_{\text{in}}G}{c_{\text{s}}^3}$$(18)

and relates the infall rate of material from the MCC, Ṁ_in, to the disc’s sound speed c_s (evaluated at the infall location). The rotational parameter Γ, $$\Gamma =\frac{{\dot{M}}_{\text{in}}}{M_{*\text{d}}{\Omega }_{\text{in}}},$$(19)

where M_*d is the combined mass of disc and star, and Ω_in the angular frequency in the disc at the infall location, compares the system’s infall time scale to its orbital time scale. Kratter et al. (2010) use this parametrisation to derive an effective Shakura–Sunyaev α to describe the transport of angular momentum in the globally unstable phase: $${\alpha }_{\text{d}}=\frac{1}{18{(2-k_{\Sigma })}^2{(1+l_j)}^2}\frac{{\xi }^{2/3}}{{\Gamma }^{1/3}}.$$(20)

Here, k_Σ is the power-law index of the disc’s surface density, l_j a parameter related to the infalling material’s angular momentum. We followed Kratter et al. (2010) in setting k_Σ = 3∕2 and l_j = 1. We thus set α_G = α_d during the infall phase.

2.5.2 Transport by spiral arms

At the end of the infall phase, when the disc mass M_d and disc-to-star mass ratio q are at their maximum and then start to decrease, the above description is no longer applicable (ξ = 0). We therefore followed K16 in using the approach of Zhu et al. (2010) to describe the effect of spiral arms on the angular momentum transport and set $${\alpha }_{\text{GI}}=\exp \left(-Q_{\text{Toomre}}^4\right)$$(21)

globally. We useed the minimum value of Q_Toomre (defined in Sect. 2.9.1) across the entire disc. There are also local prescriptions of this transport mechanism, such as in Kratter et al. (2008), so we may be overestimating α_GI. However, using a local prescription has an negligible influence on the evolution of the disc mass, and we wanted to avoid an underestimate of the angular momentum transport given the large extent of spiral arms seen in hydrodynamic simulations. As a result, we set α_G = α_GI after the infall phase.

2.5.3 Background α

Like K16 we used a high value of α_bg = 10⁻² in most of the simulations (see Sect. 3). α_bg is used to describe the torque due to any mechanism other than the global instability of the disc or spiral arms described above. For an overview of possible mechanisms, see Turner et al. (2014), see also Deng et al. (2020). The effect of MHD wind-driven accretion (magnetohydrodynamic disc winds; Suzuki & Inutsuka 2009; Bai & Stone 2013 and Gressel et al. 2015) will be addressed infuture work.

2.6 Innerdisc edge

The disc is truncated at the inner edge. This inner truncation radius is fixed at 0.05 au. Mass flowing across this radius is accreted on the central star, with 10% considered lost as outflow (Vorobyov 2010). This is a crude treatment of the inner region of the disc. However, since we were mostly interested in the regions at of tens of au in the disc, a detailed description of the innermost region is not necessary. We note, though, that the behaviour of the disc at its inner edge may have an influence on the disc’s observed lifetime.

2.7 Stellar model

During the simulation the star accretes a considerable amount of mass and its properties change significantly. We computed the stellar radius and its photospheric temperature at each time step, given its age and mass, according to pre-calculated tables from Yorke & Bodenheimer (2008). Furthermore, we considered heating of the disc by accretion of disc material onto the star, where the total stellar luminosity is given by $$L_{\text{tot}}=L_{\text{int}}+\frac{G{M}_{*}\stackrel{.}{M_{*}}}{2R_{*}},$$(22)

where L_int is the intrinsic luminosity from the stellar photosphere, M_* the mass, R_* the radius, and $\stackrel{.}{M_{*}}$ the accretion rate onto the star (K16 and Vorobyov et al. 2018). This assumes that one half of the potential energy of the gas relative to the star is dissipated in the accretion disc and only the remaining half at the surface of the star.

2.8 Disc dispersal condition and disc lifetime

We followedK16 in defining the condition for the dispersal of the disc based on the optical depth in the near-infrared (NIR) emitting region. This region is taken to be the location where the mid-plane temperature is above 300 K. The time of dispersal t_NIR is then defined as the moment when the vertical optical depth drops below unity in the NIR region.

K16 then go on to define the disc lifetime as: $$t_{\text{life}}=t_{\text{NIR}}-t_{\text{pms}},$$(23)

where t_pms is the starting time of the pre-main sequence phase given by the condition that t_grow ≥ 3t_evap, where t_evap = M_disc∕Ṁ_X and $t_{\text{grow}}=M_{*}/\stackrel{.}{M_{*}}$, ${\dot{M}}_{\text{X}}=1.62\times {10}^{-8}{(M_{*}/M_{\odot })}^{-1.57}{M}_{\odot }{\text{yr}}^{-1}$ are the disc’s evaporation and growth time scale, respectively. Thus, their lifetimes are reduce with respect to t_NIR. They argue that the evolutionary tracks of pre-main sequence stars are used to determine the age of young clusters observationally, and the age is therefore defined as the time after the pre-main sequence phase starts, not as the time after the collapse of the MCC. This is indeed how disc lifetimes are typically determined (Haisch et al. 2001). However, their definition of t_pms delays the start of the pre-main sequence to very late stages of the disc’s evolution, > 5 Myr in some cases, which is too long. Therefore, we used t_NIR as the disc’s lifetime in Sect. 5. We discuss the influence of this reduction in Sect. 6.

2.9 Fragmentation

2.9.1 Two regimes of fragmentation

The main condition for fragmentation is given by the Toomre Q parameter (Toomre 1964): $$Q_{\text{Toomre}}=\frac{c_s\kappa }{\pi G\Sigma },$$(24)

where κ denotes the epicyclic frequency. Our discs are not Keplerian and therefore κ is not equal to Ω. An axisymmetric razor-thin disc has exponentially growing modes when Q_Toomre < 1.

The conditions for disc fragmentation have been studied extensively (see Kratter & Lodato 2016 for a review). Our discs evolve in time while being fed by infalling material continuously. If matter reaches the disc at a rate faster than it can be transported away viscously, the disc will fragment independently of cooling (Boley 2009). The ξ and Γ parameters introduced in Sect. 2.5.1 can be used to judge if this is the case. The condition is given in Kratter & Lodato (2016): $$\Gamma <\frac{{\xi }^{2.5}}{850}.$$(25)

This condition is typically satisfied at the beginning of the disc’s evolution, when infall rates are high. We refer to this phase as the “infall-dominated regime”. Since we assumeed a short infall phase at constant infall rate, the condition from Eq. (25) is satisfied during most of the infall phase in almost all systems in our study. We assumed that the disc fragments when Q_Toomre < 1 in the infall-dominated regime, (i.e. when Eq. (25) is satisfied). We note that the exact Q value for fragmentation is uncertain, but is often taken as unity.

If the condition from Eq. (25) is no longer satisfied, the system transitions into the “cooling-dominated regime”. There, the disc can only fragment if it cools efficiently. This is stated in the Gammie criterion (Gammie 2001): $$\beta \equiv t_{\text{cool}}\Omega \lesssim {\beta }_c.$$(26)

In the cooling-dominated regime we therefore assumed that the disc fragments when Q_Toomre < 1 and the condition from Eq. (26) is satisfied (Kimura & Tsuribe 2012; Armitage 2007 and Baehr et al. 2017). In Eq. (26) we have (Mordasini et al. 2012): $$t_{\text{cool}}=\frac{3\gamma \Sigma c_{\text{s}}^2}{32(\gamma -1)\sigma T^4}{\tau }_{\text{eff}}$$(27)

with an effective optical depth τ_eff = κ_RΣ∕2 + 2∕(κ_RΣ), and β_c ≈ 3 the critical cooling parameter (Deng et al. 2017 and Baehr et al. 2017). The adiabatic index γ was fixed at 1.45. In reality it depends on temperature and on the ortho-to-para ratio of molecular hydrogen. The latter can vary substantially (D’Angelo & Bodenheimer 2013) and its precise value in protoplanetary discs is currently unknown (Boley et al. 2007). This may have an important effect on fragmentation and planet migration and should be studied in more detail in the future.

2.9.2 Initial fragment mass

When the conditions for fragmentation are satisfied, a region of the disc can collapse under its own gravity and form one or several bound clumps with a given initial mass. The number of clumps formed corresponds to the dominant azimuthal wave number m. Dong et al. (2015) find m ~ M_*∕M_disc. We therefore removed m times the initial fragment mass on a free-fall time scale ${\tau }_{\text{ff}}={\sqrt{3\pi /(32G{\rho }_{\text{c}})}}^{-1}$, where we substituted ρ_c with the density at the midplane.

There is currently no agreement about the value of the initial fragment mass. A rough estimate is the Toomre mass (Nelson 2006) $$M_{\text{T}}=\frac{\pi c_{\text{s}}^4}{G^2\Sigma }.$$(28)

Forgan & Rice (2011) estimate the local Jeans mass inside the spiral structure of a self-gravitating disc (using our convention for H, H = c_s∕Ω in absence of auto-gravitation) by³ $$M_{\text{J},\text{FR}}=\frac{4}{3}\frac{2^{1/4}{\pi }^{11/4}{c}_s^3{H}^{1/2}}{G^{3/2}{\Sigma }^{1/2}\sqrt{1+{\beta }^{-1}}}.$$(29)

Boley (2009) calculate a different initial fragment mass of $$M_{\text{F}}=\frac{1.6c_s^3}{G\Omega }.$$(30)

In a Keplerian disc with Q_Toomre = 1 it holds that $$M_{\text{F}}\approx 0.16M_{\text{T}}\approx 0.04M_{\text{J},\text{FR}}.$$

The result from Boley (2009) is confirmed by their SPH simulations as well as in Tamburello et al. (2015) (albeit in a different context) for a large parameter space. Therefore, we used M_F as the default value, but also investigate the other extreme M_J,FR.

2.10 Infall from the molecular cloud core

The simulations of B18 infer disc-to-star mass ratios of order unity shortly after the disc formation. Since these conditions are favourable for fragmentation, our simulations start as early as possible.

2.10.1 Disc formation

The results from B18 show one feature in particular: the process of disc formation is chaotic. Accretion from the environment may be interrupted and restart later. This may lead to inner and outer discs that are misaligned. Interaction between stars may lead to multiple systems or let multiple systems become unbound. Discs may be truncated or stripped by such processes. Due to the nature of our model, which includes a rotationally symmetric disc around a single star, we cannot reproduce this diversity. However, while B18’s simulation is a population synthesis of “protostellar” discs, our focus lies on “protoplanetary” discs. After all, for planet formation to proceed, favourable conditions will need to be in place for some time. We therefore concentrate on discs that survive the chaotic early phase described in B18 and evolve into “well-behaved” discs that can be reasonably well described by rotationally symmetric structures (although they could develop spiral arms). Of course the distinction between protostellar and protoplanetary discs is by no means sharp and we cannot exclude the possibility of planet formation in discs that do not meet our criteria of “well-behaved”. With the choice of our sample and the 1D framework we made some strong but necessary assumptions to study planet formation statistically.

We initialised our simulations 1 kyr after the formation of the protostar. During the first ~10³ to ~10⁵ yr the system’smass increases, typically by more than an order of magnitude, through infall from the MCC. We assumed that discs are formed in a short period of high and constant infall rates. This is in good qualitative agreement with B18 (see their Fig. 11) in the sense that most systems are dominated by very high infall rates at early times, followed by a rather sudden drop. The procedure for setting up the initial conditions is explained in Sect. 3.

2.10.2 Infall location

An important, yet difficult question is the location where the infalling material is added to the disc. Simple analytic models for the collapse of MCCs (for example Shu 1977) calculate the trajectories for the infalling material and therefore also their intersection with the disc’s midplane. But this treatment is too simple. First, the material does not land in the midplane but interacts with the disc higher up, which requires a treatment of the disc’s boundary. Second, the infalling material carries specific angular momentum different from the disc’s. Thus, a description of how the two components mix is required.

A treatment of these processes (albeit based on an obsolescent infall model) can be found in Visser & Dullemond (2010). Since in our model angular momentum and trajectories of the infalling material are unknown, we resorted to a simpler approach. At early times in the disc’s evolution, the infall rates are very high. Therefore, the disc’s mean specific angular momentum is comparable to that of the material incoming from the MCC. Thus we simply added the infalling material near the location characterised by a specific angular momentum equal to the mean specific angular momentum of the entire disc. Assuming a Keplerian disc with a surface density profile of the form $\Sigma ={\Sigma }_0{\left(\frac{r}{r_0}\right)}^{-1}$ (consistent with B18), this radius R_i is well-defined and it turns out (see Appendix B for a derivation): $$R_{\text{i}}\approx 0.7 R_{63.2},$$(31)

where R_63.2 is the disc radius we take from B18 (defined as the radius containing 63.2% of the disc mass). In practice, the material from the MCC is added to the disc using an infall source term of Gaussian shape centred at R_i, with a standard deviation σ_i = R_i∕3: $$S_{\text{inf}}(r,t)=\frac{{\dot{M}}_{\text{in}}}{2\sqrt{2}{\pi }^{3/2}{R}_{\text{i}}{\sigma }_{\text{i}}}\exp \left[-\left(\frac{r-R_{\text{i}}(t)}{\sqrt{2}{\sigma }_{\text{i}}}\right)\right].$$(32)

This is a compromise between having a source term that is too concentrated in radius and one that is very wide. We note that, aside extremes, the choice of this width does not influence the disc’s fragmentation and evolution significantly. R_i is assumed to grow at a constant rate b_disc: R_i ≡ R_i(t) = R_i,k + (t − 1 kyr)b_disc, where R_i,k is the initial value of R_i at t = 1 kyr. R_i,k and b_disc are obtained from the results of B18. For systems where we tried to mimic the effect of magnetic fields on the infall, R_i is kept constant. How we model this MHD-collapse, and how the initial conditions are determined in detail, is explained in the following section.

3 Investigated parameter space and initial conditions

We performed five runs to investigate the influence of different physical parameters on the inferred discs. Each run consists of 10 000 systems. The runs are summarised in Table 1.

RUN-1 is the default run in this study. It is based on the initial conditions obtained from hydrodynamical simulations where we included heating from infalling material of the MCC on the disc as well as accretion heating from disc material that is accreted onto the star. These heating mechanisms are also active in RUN-2 and RUN-3, but not in RUN-4 and RUN-5. RUN-2 is a modification of the first run where we aimed to mimic the effect of magnetohydrodynamics (MHD). Hennebelle et al. (2016) calculate the early disc radius r_H16 resulting from a magnetised collapse: $$r_{\text{H16}}=18\text{au} {\left(\frac{A}{0.1\text{s}}\right)}^{2/9}{\left(\frac{B_z}{0.1\text{G}}\right)}^{-4/9}{\left(\frac{M_{*\text{d}}}{0.1M_{\odot }}\right)}^{1/3}.$$(33)

Here, A is a measure for the ambipolar diffusivity⁴ and B_z the magnetic field in the inner part of the core. We used this expression, setting A = 0.1 s and B_z = 0.1 G to determine R_i for RUN-2. RUN-2 therefore differs from RUN-1 only in the infall location, which we chose very close to the star (1.7–9 au) and constant in time throughout the simulation. The locations are simply chosen such that the disc radius agrees with the analytic result from Eq. (33) at the end of the infall phase. In practice, finding the “correct” radii is an iterative process, but it works well (see the middle left panel in Fig. 5 in the Results section). The precise infall locations we used are given in Appendix D. RUN-3 is equal to RUN-1 with the exception of a lower assumed α_bg. In RUN-4 we turned off both infall and accretion heating. RUN-5 is similar to RUN-4 except for the choice of the initial fragment mass, where we used the larger M_J,FR (see Sect. 2.9.2).

Our model can cover a large range in stellar mass. We used a range of final (at the end of the simulation) stellar masses from 0.05 M_⊙ to 5 M_⊙ and divided this interval into 100 logarithmically spaced bins. We computed the evolution of 100 systems in each of thesebins to give a total of 10 000 systems per run. The initial conditions are chosen such that the resulting population of stars has a mass distribution in agreement with the stellar initial mass function (see following sct.). The final stellar mass of a system with a given set of initial conditions is a priori unknown. We therefore make an estimate, as we discuss in the following section.

Each simulation is initiated at 1 kyr by setting up a disc with an initial mass M_disc,i and an initial radius R_disc,i along with a protostar of mass M_*,i. We chose a power-law surface density profile with exponent −1, as expected by B18. The precise shape of the initial profile has little importance, since it is changed by infall and viscous evolution immediately. In order to choose suitable initial values, we selected 35 systems from the single star sample in B18 and constructed probability distributions for M_disc,i, M_*,i and R_disc,i. The (constant) infall rate is set as a probability distribution around Ṁ_in = (M_*,10k + M_disc,10k − M_*,i − M_disc,i)∕(9 kyr), where M_*,10k and M_disc,10k are the stellar mass and disc mass after 10 kyr, respectively. The time of formation of a protostar therefore agrees with that in B18 by construction. A similar procedure was used to set up the initial disc radius. We chose a sub-sample of 20 among the 35 previously chosen discs (some do not have well-defined radii during the first 10 kyr) to generate a probability distribution for R_i,k. We also used these discs to gauge the rate b_disc at which the infall radius expands in time (see Sect. 2.10.2). We assumed that protostellar mass, disc mass and infall rate are correlated, and that disc radius and its expansion rated are correlated. Mass and radius are assumed to be uncorrelated. The specific discs we selected, along with histograms showing the distributions of the initial parameters, can be found in Appendix E.

At this stage, there remains one important unknown: the infall duration t_infall. We could not obtain this quantity from B18, since their simulations do not run long enough for all systems to reach the end of the main accretion phase. Therefore, we chose a different approach⁵. At the end of the simulations, when the discs have dispersed, we expect the resulting distribution of stellar masses to obey the IMF. We used the IMF from Chabrier (2007) for individual stars:

Fig. 1 Stellar masses at the end of the simulation for all runs (see Table 1), compared to the Chabrier (2007) IMF.

$$\frac{\text{d}n}{\text{d}\log m}=\{\begin{array}{ll}0.093\times \exp \left\{-\frac{1}{2}\frac{{(\log \tilde{m}-\log 0.2)}^2}{{0.55}^2}\right\},\hfill & m\le 1M_{\odot }\hfill \\ 0.041{\tilde{m}}^{-1.35},\hfill & m\ge 1M_{\odot },\hfill \end{array}$$(34)

where dn∕dlogm denotes the stellar number density in pc⁻³ per logarithmic interval of mass and $\tilde{m}\equiv m/(1 M_{\odot })$. We chose the distribution for t_infall in such a way that in each run the resulting distribution of final stellar masses in our simulation agrees with this observed distribution. Figure 1 shows this works quite well for all our five runs. The distributions for t_infall can be found in Appendix E.

4 Formation and evolution of an example system

In this section, we demonstrate how a star-disc system forms and evolves by describing the time-evolution of a specific system (System 6410 from RUN-1, see Table 1). This is a typical system forming an ~ 1 M_⊙ star. The system is initialised with a protostar of mass 0.09 M_⊙ and a disc ofmass 0.02 M_⊙ with a radius of 41 au. The top left panel of Fig. 2 shows the evolution of the stellar mass, disc mass, and cumulative mass removed from the disc by fragmentation. The middle left panel depicts as a function of time the infall rate from the MCC, the accretion of disc material on the star, as well as the fraction of this material considered to flow out of the system. The bottom left panel shows the stellar luminosity and its two contributors, the stars intrinsic luminosity and the luminosity due to accretion of material onto the star. The right top and right middle panels show the surface density and temperature, respectively, as a function of radius at selected times. The contour plot in the right bottom panel color-codes the radial distribution of the surface density as a function of time. The fragmentation criteria are also shown. Fragmentation is forbidden because of the Gammie cooling criterion (Eq. (26)) in regions enclosed by the dash-dotted lines (inside of about 20 au for t > 0.1 Myr). The Gammie criterion thus never limits fragmentation in this simulation. We note that the condition for the infall-dominated regime of fragmentation (Eq. (25)) is always satisfied during the infall phase and is therefore not shown in the figure.

The evolution of the system can be divided into four phases. The first is the infall phase. In this example system, it lasts ~ 30 kyr. In this phase, the system is dominated by high infall rates from the MCC and fast transport of angular momentum due to the global instability of the disc (Sect. 2.5.1). The ratio of disc mass to stellar mass, q, is high (q~ 1). The disc is globally unstable and fragments 16 times in the outer region. Both star and disc grow in mass by a factor of about ten. The second phase is short (≈2 kyr) and barely seen in Fig. 2 as a sudden drop in disc mass, accretion rate and luminosity. Temperatures in the disc decrease quickly due to the reduced heating after infall has ceased. Consequently, the disc fragments another 15 times. q drops by ~ 10 % due to fast transport of angular momentum. A total of 0.07 M_⊙ of matter is removed from the disc by fragmentation. During the infall phase, of the order of one fragment per kyr is formed, and afterwards it is roughly ten times more. Both values are in the range of the number of fragments per time expected from hydrodynamic simulations of sufficiently high resolution (see Boley 2009 and Szulágyi et al. 2017 for an example with and without infall, respectively). We discuss the comparison to other studies further in Sect. 6.4. Fragments are added to the star in one time step (~ 1 yr during this phase). This is much longer than an orbital period at the truncation radius. Since we do not model the accretion process of the clumps onto the star, the accretion rates and luminosities in Fig. 2 are smoothed between ~ 20 and 40 kyr. The accretion of a fragment leads to a strong increase in stellar luminosity for a short amount of time. This may be observed as episodic accretion (Audard et al. 2014) and may affect the thermodynamics of the disc for a short time (Cieza et al. 2016). Our model cannot predict the duration or luminosity of these outbursts. However, the interval between outbursts is typically a few 10² to 10³ yr. The third phase lasts approximately 100 kyr and is characterised by a fast redistribution of angular momentum by spiral arms. The disc-to-star mass ratio is reduced to ~ 0.6. In the last phase the viscous evolution of the disc progresses much more slowly in the absence of gravitational instabilities until the disc disperses after 14 Myr. Figures showing the time-evolution of two more example systems, a low-mass system from RUN-2 and a high-mass system from RUN-5, can be found in Appendix C. Animations of all three systems are available online.

Fig. 2 Time evolution of System 6410 from RUN-1. Top left panel: stellar mass, disc mass and, cumulative mass removed from the disc by fragmentation. Middle left panel: accretion and outflow rates. Bottom left panel: stellar luminosity. Top right panel: surface density at different times. Middle right panel: midplane temperature at different times. Bottom left panel: stellar luminosity (accretion and intrinsic). Bottom right panel: contour plot of the surface density with fragmentation criteria (see text). The black vertical lines denote, in order of increasing time: tinfall ≈ 30 kyr, tpms ≈3.7 Myr and tNIR ≈ 14 Myr (see Sect. 4).

5 Results

We now consider all runs set up in Table 1 and study the statistics of their outcomes. Table 2 gives a summary of the disc properties at the end of the infall phase and the lifetimes of all runs. Mean and standard deviation of the respective quantities are shown. We discuss these together with the figures that follow. The systems’ properties at the end of infall are of particular interest, since they represent the initial conditions for planet formation models, where the collapse of the MCC is not modelled.

The results are organised as follows: in Sect. 5.1, we analyse the disc properties at the end of the infall phase as well as the disc lifetimes for RUN-1 and RUN-2. The same quantities are discussed in Sect. 5.2 for RUN-3, RUN-4 and RUN-5. The following subsection (Sect. 5.3) presents the dependency of these quantities on the final stellar mass. We then discuss the potential fragmentation for all runs (Sect. 5.4).

Figures 3 and 4 present the end-of-infall distributions of disc mass, disc radius, and disc-to-star mass ratio q_infall = M_disc,infall∕M_*,infall, as well as the disc lifetimes for the entire runs. Thus, they are a combination of the simulations performed in 100 mass bins (see Sect. 3). The different relative contributions of the stellar mass bins are considered and the slight deviations from the IMF (Fig. 1) are corrected. The results should therefore be representative of (unbiased) observationsof actual star forming regions, to the degree that the assumed IMF is representative of the stellar population in the region.

5.1 RUN-1 and RUN-2: “Hydro versus MHD”

The disc properties at the end of the infall phase for the “hydro” and the “MHD” runs and their inferred lifetimes are shown in Fig. 3. The mean disc mass in RUN-1 is higher than that in RUN-2 by almost a factor of three (top left panel of Fig. 3). This is expected since the infalling mass is deposited very close to the star throughout the simulation in RUN-2. Thetop left panel of Fig. 3 also shows observed Class 0 disc masses from Tychoniec et al. (2018; see Sect. 6.1 for a discussion).

The large difference between RUN-1 and RUN-2 is also visible in the disc radii (top right panel of Fig. 3), where amuch wider distribution is seen and the discs are found to be larger by almost a factor of six in RUN-1 compared to RUN-2. We used the same definition for the disc’s radius as B18 (the radius containing 63.2 % of the disc’s mass).

The distributions of the ratios of disc masses to stellar masses at the end of the infall phase (bottom left panel) showthe most striking difference between the two runs. In RUN-1 the distribution is centred around q_infall ≈ 1.0 and very wide, while the mean value is roughly a factor of three lower in RUN-2, with a very narrow distribution. A consequence of the high values of q_infall is that the stellar mass at this time (M_*,infall) differs significantly from the final stellar mass M_*,final.

The shape of the distribution of disc lifetimes (bottom right panel of Fig. 3) is very similar in RUN-1 compared to RUN-2, but shifted to longer lifetimes by almost 50%: a large fraction of the mass is transported to the star during the infall phase in RUN-2 early on. The mean lifetimes of 7.3 Myr (RUN-1) and 4.5 Myr (RUN-2) are long, see Sect. 6. The quantity displayed, t_NIR, corresponds to the disc’s total lifetime as explained in Sect. 2.8 with no reduction due to the start of the pre-main sequence phase applied. Figures of t_life (with the reduction applied) can be found in Appendix F.

5.2 RUN-3, RUN-4, and RUN-5: accretion heating and fragmentation

Here we show the same figures as in Sect. 5.1, but this time for RUN-3, RUN-4, and RUN-5. RUN-1 is also shown for comparison.

There isalmost no difference in the distribution of disc masses between RUN-1 and RUN-3 (top left panel of Fig. 4). This is expected: the infall phase is dominated by the global instability of the discs, the lower background viscosity parameter α_bg is only important later in the disc’s evolution. The disc mass distributions of RUN-4 and RUN-5 are similar in shape to that of RUN-1, but shifted to successively lower masses. Turning off accretion heating due to infalling matter from the MCC (infall heating) as well as accretion of disc material onto the star (stellar accretion heating) leads to lower temperatures in the disc and promotes fragmentation. In our model, fragmentation leads to the removal of mass from the disc to the star. This effect becomes more pronounced in RUN-5, where much more mass is removed from the disc each time the conditions for fragmentation are satisfied (see also Sect. 5.4). We note that the effect of stellar accretion heating is much more important than that of infall heating.

A similar effect on the disc radii is seen in the top right panel of Fig. 4. It is less pronounced than the effect on the disc masses.

The most obvious differences between the runs that are based on hydrodynamic simulations (RUN-1, RUN-3, RUN-4 and RUN-5) are seen in the distribution of q_infall (bottom left panel of Fig. 4). When the disc fragments more often due to lower temperatures in the disc, the disc mass is reduced, while the stellar mass grows. So q_infall exhibits the “double effect”. The distribution of q_infall is shifted to lower values in RUN-4, and even more so in RUN-5, compared to RUN-1 and RUN-3.

The distributions of lifetimes differ negligibly between RUN-1, RUN-4 and RUN-5 (bottom right panel of Fig. 4). Fragmentation can remove substantial amounts of mass from the disc. However, it leaves the inner disc (important for the determination of the lifetimes) mostly unchanged over long time scales. Lifetimes are longer by a factor of ~ 5 in RUN-3. This is expected: during the vast majority of the disc’s life, its evolution is dominated by the choice of α_bg.

Fig. 3 Distributions of the disc properties at the end of infall and lifetimes of the discs for RUN-1 and RUN-2. Top left: disc masses, top right: disc radii including the observational result of Tychoniec et al. (2018) (Ty18), bottom left: disc-to-star mass ratio qinfall, bottom right: disc lifetimes. All stellar masses are included in this figure.

5.3 Dependency on the final stellar mass

Here we show the same quantities as in the preceding section, but this time as a function of the final stellar mass M_*,final because this quantity is more easily observed than the stellar mass at the end of the infall phase. Figure 5 shows the mean values of disc mass, disc radius and disc-to-star mass ratio for each of the 100 mass bins. The spread in these parameters is also shown: the shaded region corresponds to ± one standard deviation. RUN-1 and RUN-2 are depicted in the left, RUN-1, RUN-3, RUN-4 and RUN-5 in the right panels of the figure. Histograms for specific mass bins can be found in Appendix G.

The dependency of M_disc,infall on M_*,final is roughly linear and very similar in shape for all runs. The spread is highest at low stellar masses. The reason is that these systems are more strongly dominated by the initial conditions: a larger fraction of the total mass is already in the system at the beginning of the simulation.

The discs’ radii, R_disc,infall, scale linearly with M_*,final to good approximation. The differences between the runs are very small, with exception of RUN-2, as expected. The middle left panel of Fig. 5 also contains an analytic estimate for the early disc radius (see Eq. (33) and Sect. 6.2 for a discussion).

The ratio M_disc,infall∕M_*,infall varies strongly with M_*,final (third row in Fig. 5). For all runs but RUN-2, it rises from around 0.5 at the lower end of the mass range considered to some maximum value at around 0.2 M_⊙, then decreases again to ≈0.5. The maximum value reached depends on the cooling and fragmentation mechanisms that are at play. When fragmentation is inhibited by stellar accretion heating, the mean values of q_infall reach 1.1 (RUN-1 and 3) while they stay below 0.8 in most systems when heating was turned off (RUN-4). They are further reduced when a larger initial fragment mass is used (RUN-5). This shows that the amount of mass that is removed from the disc by fragmentation is limited by the initial fragment mass in runs 1-4. The influence of the heating mechanisms discussed above is, however, more important than that of the initial fragment mass.

The behaviour of q_infall is different in RUN-2 compared to the other runs. Here q_infall depends only weakly on the final stellar mass. It is typically around 0.4 at the lowest stellar masses and decreases to a nearly constant ~0.25 at stellar masses ≳0.2 M_⊙.

The discs’ lifetimes as a function of final stellar mass are depicted in the bottom row of Fig. 5. Lifetimes exhibit a weak positive correlation with final stellar mass. We noted earlier, that runs 1, 4 and 5 show very little difference in disc lifetime. This remains true when looking at the dependence on stellar mass. The lifetimes in RUN-2 are roughly 40% shorter than those in RUN-1, due to the lower post-infall masses. Using a lower α_bg produces lifetimes that are around five longer except at the lowest masses. We compare the disc lifetimes to observations in Sect. 6.3.

Fig. 4 Properties at the end of the infall phase and lifetimes of the discs for RUN-1, RUN-3, RUN-4 and RUN-5. Top left: masses. Top right: radii. Bottom left: disc-to-star mass ratio qinfall. Bottom right: lifetimes.

5.4 Fragmentation

Here we look at the fragmentation behaviour of the discs. An overview is given in Table 3, where the global fragmentation properties of all runs are shown. Some of the simulated discs fragment, producing anywhere between ~ 1 and ~ 10³ fragments, while others do not fragment at all. The formation of several hundred fragments in a single disc has not been reported in the literature to our knowledge. Such large numbers of fragments are likely an overestimate related to our assumptions. We discuss this in Sect. 6. The fraction of fragmenting discs is given in the first column of Table 3. The following columns show the mean number of fragments, the mean initial mass of the fragments and the mean fragmentation location, respectively. The values in these three columns only take into account the discs that do fragment. For instance, if in one run every other disc fragments, with five fragments each, the mean number of fragments is 5.

Figure 6 displays the fraction of discs that fragment as a function of final stellar mass. Of the 10 000 discs in RUN-2, not a single fragments. This demonstrates the enormous importance of magnetic fields and magnetised collapse for the question whether protostellar and protoplanetary discs can be gravitationally unstable, and thus whether giant planets may form by GI. Figure 6 shows the background viscosity has very little influence on the fraction of discs fragmenting. There is only a slight difference at low stellar masses where only few discs fragment in RUN-1 and RUN-3. Our statistics are therefore not reliable in this region of parameter space and we do not deem this difference significant. There is however a significant difference when stellar accretion heating is turned off. In this case, the fraction of discs that fragment is increased by 50% as seen in Table 3. Stellar accretion heating heats the discs keeping them gravitationally stable, thus inhibiting fragmentation. Figure 6 reveals this happens mainly in systems with final stellar masses < 1 M_⊙.

For the remainder of this section, we concentrate on the discs that do fragment. RUN-2 is not included in the following figures.

5.4.1 Global fragmentation properties

Here, we focus on the fragmentation properties of the different runs as a whole. We discuss how many fragments are formed, their initial mass and the location in the disc at which the fragments form. We will refer to this location as “fragmentation location”. The top left panel in Fig. 7, shows the distribution of the number of fragments in fragmenting discs. RUN-1 and RUN-3 show a very similar behaviour. This is due to the very similar progress of the disc evolution during the infall phase, where fragmentation predominantly happens. Turning off the heating mechanisms gives rise to a strong increase in the number of fragments (RUN-4). The number of fragments is more than quadrupled. Using the local jeans mass criterion for the initial fragment mass reduces the number of fragments by around a factor of ≈ 8 (RUN-5 compared to RUN-4).

The mass distribution of the fragments (middle left panel of Fig. 7) is very similar in RUN-1 and RUN-3, for the same reason as above. The fragments formed in RUN-5 are much more massive initially, as expected from the different criterion. The mean mass is, however, not a factor of 25 larger than that in RUN4, as one might naively expect from the ratio of M_J,FR to M_f (see Sect. 2.9.2). Instead, the difference is only around a factor of ten. This is because the initial fragment mass is calculated self-consistently from the local conditions in the disc each time the conditions for fragmentation are satisfied. To illustrate this, let us consider two systems with identical initial conditions. The initial mass for the first fragment is indeed a factor of ~ 25 apart when the two different criteria are used. However, the subsequent evolution of the system is different and the initial fragment mass becomes self-limited also with the Jeans mass criterion.

The locations in the disc, where fragmentation happens, is depicted in the bottom left panel of Fig. 7. The main difference here is between the runs with stellar accretion heating (RUN-1, RUN-3) and the runs without (RUN-4, RUN-5): the former fragment predominantly at radii < 100 au, the latter at radii ≳100 au. This shows again the importance of this heating mechanism. It also gives an explanation for the difference in mass discussed above: initial fragment masses are different (and typically lower) closer to the star.

Fig. 5 End-of-infall properties and lifetimes of the discs as a function of the final stellar mass M*,final. Left column: RUN-1 and RUN-2. Right column: RUN-1, RUN-3, RUN-4 and RUN-5. We show, from top to bottom, masses, radii, disc-to-star mass ratios qinfall at the end of the infall phase, and lifetimes.

Fig. 6 Fraction of discs that fragment versus final stellar mass for all runs. The “MHD” run does not produce any fragments.

5.4.2 Fragmentation properties as a function of final stellar mass

The fraction of fragmenting discs around stars ≲0.1M_⊙ is low, in RUN-1 and RUN-3 in particular. (see Fig. 6). Therefore, the statistics in the corresponding mass bins is not very robust. This should be kept in mind when looking at the following figures.

The right column of Fig. 7 shows the mean number of fragments per disc, the mean initial fragment mass and the mean fragmentation location as a function of M_*,final, respectively. Masses and locations were averaged twice. For example, first, the mean fragment mass in one system (producing multiple fragments) was computed. Then, the mean of these was calculated in each mass bin. As seen in the figure, all three quantities increase with increasing final stellar mass. The dependency of the number of fragments n_fragment is strongest: $n_{\text{fragment}}~M_{\text{*},\text{final}}^{5/4}$ (top panel). For the initial fragment mass, the dependency is weak, $M_{\text{fragment},\text{init}}~M_{\text{*},\text{final}}^{1/2})$ (middle panel) and for the fragmentation radius we have: R_{fragment,init} ~ M_*,final (bottom panel).

6 Comparison to previous studies and to observations

Clearly, it isinteresting to compare our results to other theoretical studies as well as to observations. However, as we are going to elaborate in this section, such a comparison is difficult. This is due to limitations in existing studies (including our own) as well as observational uncertainties and biases. The population of stars resulting from our simulations agrees with the IMF by construction, as discussed (see Fig. 1). In the following we concentrate on the masses, radii and lifetimes of the discs as well as on their fragmentation.

6.1 Disc masses

The most obvious study to compare our disc masses to is B18. However, the results cannot be directly compared since we only used a sample to obtain initial conditions for our simulations. Nevertheless, we compare the results to identify the differences.

The topleft panel of Fig. 8 shows the disc masses at the end of the infall phase from RUN-1, together with the max. disc mass in two sets of data from B18. The first set comprises all the 183 “discs” in the simulation, the second includes only the 35 discs we chose to construct our initial distributions from (see Sect. 3). The representation of RUN-1 in Fig. 8 contains exactly the same data as the top left panel of Fig. 3, but with a different binning more appropriate to the other data displayed.

In our simulations, the end of the infall phase is a well defined point in time: infall rates are constant at the beginning, and zero later. The disc mass is always maximum (or very close to) at the end of this phase. This simple picture does not hold for the simulations in B18. There, infall rates are not constant and systems can stop accreting and restart later. Furthermore, some systems are still accreting at the end of the simulation, so the corresponding disc masses are still increasing by an unknown amount.

The mean value of the maximum disc masses of the full set shown in Fig. 8 is 0.08 M_⊙, a little less than a third of that from RUN-1. In the other set, containing only the 35 selected discs, it is 0.23 M_⊙, or about 80% of that from RUN-1. Given the difficulties explained, and considering we are comparing results from a 3D SPH calculation to those from an axis-symmetric 1D viscous evolution model, the agreement seems very reasonable.

Tychoniec et al. (2018) perform an observational analysis of dust emission as part of the VLA Nascent Disk and Multiplicity (VANDAM) survey. They present estimates for the distribution of the protostellar disc masses in Perseus, both for Class 0 and Class I phases. We compare their Class 0 data to our RUN-1 and RUN-2 in the top left panel of Fig. 3. Tychoniec et al. (2018) give a mean disc mass in Class 0 of 0.075 M_⊙, about 70% of our RUN-2 or a quarter of our RUN-1.

The similarity of the observed disc masses with our RUN-2 should not be overstated. First, our data shows the maximum masses the discs reach during their life, while the VANDAM survey probes discs with some (unknown) distribution of Class-0 ages. Second, our sample of discs is one that covers a large range of stellar masses, without any observational bias applied. As discussed in Sect. 5.3, the distribution of disc masses strongly depends on that of the stellar masses. The masses of embedded protostars are very difficult to measure, and are not given in Tychoniec et al. (2018). If, for instance, systems with final stellar masses larger than ≈ 0.3 M_⊙ are favoured in the survey, this is indicative of a distribution of disc masses more massive than that of our RUN-2.

Fig. 7 Fragmentation properties of all runs. Left column, from top to bottom: distribution of the number of fragments, mean initial fragment mass, and mean fragmentation location. Right column: same properties as a function of final stellar mass.

6.2 Disc radii

In RUN-2 (“MHD”), we chose an infall location close to the star and constant in time. This is an attempt to study the influence of magnetic fields on the formation and evolution of protoplanetary discs. As a consistency check, we show our disc radii together with the analytic expression from Hennebelle et al. (2016) in the second panel in the left column of Fig. 5. Hennebelle et al. (2016) compare their analytic estimate for the early disc radius to a large number of 3D, non-ideal MHD collapse simulations. They find agreement within a factor of two. This range is also shown (black shaded region). The disc radii from RUN-2, including uncertainty, lie inside this region across almost the entire range of stellar masses. Indeed, we chose the infall location specifically to produce this agreement. Nevertheless, the outcome is not self-evident. During the infall phase, the viscosity in the disc is very high and the disc radii grow very quickly even when the infall location does not progress outwards.

Tobin et al. (2020) (hereafter To20) performed a multi-wavelength survey of hundreds of protostars. They use dust continuum emission to measure Class 0 dust disc radii. In the right panel of their Fig. 11, they compare the radii of discs around non-multiple protostars to non-binary systems from B18. We plot their data along with the end-of-infall radii from RUN-1 and RUN-2 in the top right panel of Fig. 8. We also show two sets of data extracted from the online material of B18. The dash-dotted line depicts the temporal mean of the full sample, while the dotted line depicts the maximum radius the discs from our sample (see Sect. 3) reach during the simulation.

The observed disc radii lie in the middle between RUN-1 and RUN-2. The two representations of radii from B18 span a huge range of disc radii and seem compatible both with the observed radii and those from our RUN-1, but not with RUN-2. However, we would like to stress again how difficult this comparison is. The data from our runs comes from a very specific point in time, while the observed data has an unknown distribution of underlying ages. Additionally, the observations are from dust continuum emission while our data (along with that from B18) results from the simulation of pure gas discs. The distributions of gas and dust may ormay not be similar at early times. Therefore, we do not conclude discs stemming from a magnetised collapse are too small (these, too, grow quickly in time), or those from hydrodynamic collapse too large. An educated conclusion on this matter would necessitate the knowledge of the relationship between the distribution of dust of different sizes and that of the gas. Furthermore the ages of the observations would have to be well constrained. The uncertainties in the observations are also large. Nevertheless, disc radii somewhere between those from our RUN-1 and RUN-2 seem the most plausible. A promising path to learn from the comparison of observational data with theoretical models is to include the evolution of heavy elements in models such as our own and use that to predict observables.

Fig. 8 Comparisons of our results to previous studies and observations. Comparison of, top left: maximum disc masses from RUN-1 with those from the simulations of B18. Top right: disc radii at the end of the infall phase with observed Tobin et al. (2020, To20) and simulated B18 Class 0 radii. Bottom: disc fractions as a function of time with fits from Richert et al. (2018, Ri18); left: tNIR, right: tlife (reduced, see Sect. 2.8). RUN-4 and RUN-5 are indistinguishable from RUN-1 in the bottom right panel and omitted.

6.3 Disc lifetimes

There is another theoretical study relating collapsing MCCs to disc lifetimes: Li & Xiao (2016) investigate the dependence of mass, angular velocity and temperature of MCCs in solid rotation on disc lifetimes across a large parameter space. They assume the mass distribution of the MCCs (not that of the resulting stars) is given by the IMF by Kroupa (2001) and find a characteristic time of 3.7 Myr for the exponential decay of the disc fraction. They use a lower value for α_bg (1.5 × 10⁻³), however their collapse model is based on the Shu collapse and produces much less massive discs than we see in our work.

Richert et al. (2018) performe a study of the longevity of inner dust discs using data from 69 young clusters. The work is based on data from two projects: “Massive Young Star-Forming Complex Study in Infrared and X-ray” (MYStIX, Feigelson et al. 2013) and “Star Formation in Nearby Clouds” (SFiNCs, Getman et al. 2017). They determine the ages of the clusters by combining the empirical Age_JX method (described in Getman et al. 2014) with different PMS evolutionary models. In particular, they use the stellar evolutionary models from Siess et al. (2000) and the “magnetic” model from Feiden (2016) which considers magnetic inhibition of convection. The corresponding ages are called Age_JX-Siess00 and Age-Feiden16M, respectively. We compare the disc fraction as a function of time from our work with their results. The bottom left panel of Fig. 8 compares the disc fractions (based on our criterion for the disc lifetime, see Sect. 5.3) to two of their fits. We note that these fits are based on data points that only go from ≈ 0.5 − 4Myr. In said range our disc fractions are much higher. Furthermore, our disc fractions cannot be reasonably fit by an exponential decay of the form exp (−t∕τ) for some characteristic time τ.

We also show the disc fractions as a function of time based on t_life (see Sect. 5.3) in the bottom right panel of Fig. 8. The results from RUN-1 are (likely coincidentally) almost identical to the Richert et al. (2018)-fit to their combined MYStIX+SFiNCs sample based on Age-Feiden16M. This neither indicates that the Age-Feiden16M disc fractions are more accurate than those based on, for example, Age_JX-Siess00, nor that the reduction proposed by K16 is correct, nor that our disc fractions are correct. However, it does stress one important point: The choice of time zero may be crucial when comparing the lifetimes of protoplanetary discs from simulations with those from observations. We conclude that the mean total lifetime of protoplanetary discs may well be a factor of two to three higher than the oft-cited 1–3 Myr. Values for the background viscosity coefficient much lower than α_bg = 10⁻², as in our RUN-3, are clearly not compatible with observations (though we note that some discs with lifetimes of several 10 Myr have been observed, for example Lee et al. 2020). There has been a debate in the literature about the dependency of disc lifetimes on the host stellar mass. Some authors find longer disc lifetimes around low-mass stars (for example Ribas et al. 2015), which would be in disagreement with our results. Richert et al. (2018) find no evidence for such a dependency. The choice of the photoevaporation model and/or the strength of the evaporation rates (we use a constant value) also has a strong influence on the disc lifetimes. An example of such a model is Picogna et al. (2019). We will study the influence of the photoevaporation model on our results in future work.

6.4 Fragmentation

The numbers of fragments we find in RUN-1 and RUN-3 are substantial: 24 when averaged across the whole population and up to several hundred around stars ≳ 1 M_⊙. The number even approaches ~1000 in some cases when stellar accretion heating is ignored (RUN-4). We expect this to be an overestimate. The reason is that we ignored further gas accretion of the clumps once they have formed. If gas accretion were allowed, the surface density in the discs would be reduced and further fragmentation inhibited. It is, however, not guaranteed that mass accretion reduces the number of fragments by orders of magnitude in all cases. First, accretion competes with migration and/or gap formation, which means that fragments may migrate away from the gravitationally unstable region of the disc more quickly than they accrete (Müller et al. 2018). Second, fragments may also be tidally disrupted (Boley et al. 2010 and Nayakshin 2010), returning some or all of their mass back to the disc. The number of fragments is lower by almost an order of magnitude in RUN-5 compared to RUN-4 because we used a much higher initial fragment mass (see Sect. 5.4.1). A possible interpretation of the difference between the two initial fragment masses is as follows: while M_F is a measure of the fragment mass precisely at the time when the conditions for fragmentation are first satisfied, M_J,FR can be thought of as the fragment mass after accreting disc mass for a few orbits. Realistically, no more than a few dozen fragments around stars with ≳1 M_⊙ are expected. Nevertheless, this estimate is still higher than what is found in some high-resolution hydrodynamic simulations that typically find ~1−10 fragments (for example Hall et al. 2017). Such simulations, however, are so computationally expensive that they can only be run for a few kyr (~ 1 orbit at 100 au). Therefore, this could be one “burst” of fragmentation. Our discs are evolved for much longer and are replenished in mass by infall for several 10 kyr in the case of more massive systems, so many fragmentation “bursts” are possible. Another difficulty with studying fragmentation in hydrodynamic simulations is that fragmentation is only seen at sufficient numeric resolution (André Oliva & Kuiper 2020). B18 find only 10 out of the studied 183 systems fragment. However, while the number of SPH particles in their simulation is very high (3.5 × 10⁷), the simulation covers an entire star cluster, and the number of SPH particles in the discs is clearly insufficient to accurately study fragmentation. André Oliva & Kuiper (2020) perform simulations of a forming massive star using a grid-based self-gravity-radiation-hydrodynamics simulations. Around 26 million grid cells are used to study the disc physics. This represents one of the highest-resolution simulations on disc fragmentation performed so far. In their higher-resolution runs, several dozen fragments are found around massive stars. At lower resolution, fragmentation is suppressed.

7 Discussion

In this study, we made a number of simplifying assumptions: Protoplanetary discs are axisymmetric, discs consist exclusively of gas (except for the opacity in the temperature calculation), there is one main infall phase, and finally when a disc fragments, the initial fragment mass is removed from the disc and added to the star. Furthermore, the initial conditions were generated based on only a fraction of the 183 systems from B18.

Assuming an axisymmetric disc means we cannot reproduce per se the spiral structure and global instability of the disc. While we could self-consistently determine the location where fragmentation happens, we could not do this for the number of fragments. Furthermore, we could not model the stochastic nature of fragmentation (Paardekooper 2012). Instead, we relied on parametrisations (as described in Sect. 2.5 and 2.9). As a result, individual simulations cannot reproduce the 3D behaviour of discs accurately at early times. However, using parametrisations of the 2D and/or 3D results, we performed a parameter study, and these assumptions still allow us to study the formation and evolution of protoplanetary discs statistically and systematically. At later times, the discs’ evolution is no longer dominated by gravitational instabilities and the discs are expected to be closer to rotational symmetry as seen in many observations. The influence of our assumption concerning the evolution of the fragments has already been discussed in Sect. 6.4. In choosing a sample of “well behaved” discs from B18, we may have missed some discs that, despite going through a very chaotic formation process, may still form planetary systems. This limitation should be kept in mind for future studies of planet formation. It can be overcome when more advanced hydrodynamic studies of disc formation become available. In summary, our assumptions do not prevent reasonable estimates of our main results: the distributions of disc masses and radii at the end of the infall phase, numbers and properties of fragments, and disc lifetimes.

8 Summary and conclusions

We studied the formation and evolution of protoplanetary discs from their emergence in a collapsing MCC to their dispersal after a few million years. We performed five runs using different assumptions, each consisting of 10 000 systems. The systems were initialised by means of six parameters: initial stellar mass, initial disc mass, (constant) infall rate, disc radius, and the rate at which the latter increases with time. Distributions for these parameters were obtained from a selection of systems from the hydrodynamic disc population synthesis by Bate (2018).

We included the influence of the disc’s auto-gravitation on the angular frequency and scale height as well as transport of angular momentum by global instability, self-gravity and turbulent viscosity. Mass loss by photoevaporation was also included. Furthermore, we performed one run in which we investigated the effect magnetic fields might have on the infall and disc formation (RUN-2 “MHD”).

When comparing our default run (RUN-1) with RUN-2, we find that this produces very different discs. RUN-2 produces systems that typically have disc radii a factor of six smaller, disc masses a factor of 2.6 lower, and disc-to-star mass ratios that are a factor of 2.8 lower, with much narrower distributions. Systems in RUN-2 also have disc lifetimes ≈ 40% lower. Furthermore, none of the systems in RUN-2 fragment.

Fragmentation is, however, common in all other runs. We find that approximately half of the systems fragment when shock heating from disc material accreting on the star (accretion heating) heats the disc. If this effect is neglected, even three quarters of the systems fragment. Clearly, based on our current results we cannot make any robust statement about the feasibility of planet formation via GI yet. If at least some discs form in a way similar to our runs based on hydrodynamical simulations, many of them must fragment and would provide numerous bound clumps. It would then seem conceivable that at least some of these could survive and grow to form giant planets similar, for example, to the ones observed in the HR8799 (Marois et al. 2008) or GJ 3512 (Morales et al. 2019) system. We plan to study in more detail what happens to fragments once they have formed in future work. In summary we find:

If protoplanetary discs are formed in agreement with hydrodynamic simulations B18 as discussed in our work,

• 1.

  they are massive early on: (0.29 ± 0.06) M_⊙ or 1 ± 0.3 times their host stellar mass;

• 2.

  in systems with a final stellar mass of ~1 M_⊙ they are even more massive: 0.7 M_⊙ at the end of the infall phase;

• 3.

  they remain massive (a few tenths of their host stellar mass) for a large fraction of their lives;

• 4.

  they are large: (200 ± 100) au at the end of the infall phase

• 5.

  their total lifetimes are long: (7.3 ± 0.5) Myr, despite choosing a high value of 10⁻² for the background viscosity α-parameter;

• 6.

  half of the systems fragment. Neglecting stellar accretion heating raises this fraction to three quarters;

• 7.

  fragmentation potentially removes a lot of mass from the disc, either through numerous fragments, through a high initial fragment mass, or (not modelled) through subsequent accretion; and

• 8.

  the final stellar mass has a strong influence on fragmentation, with more massive systems fragmenting more often and producing more fragments of higher mass.

If discs are instead formed in a way expected by a magnetised collapse (Hennebelle et al. 2016),

• 9.

  they are almost a factor of three less massive, though their masses are still substantial: 0.1 M_⊙ at the end of the infall phase;

• 10.

  in systems with a final stellar mass of ~1 M_⊙ they have even a mean mass of ~0.2 M_⊙;

• 11.

  they are a factor of six smaller: (36 ± 7) au at the end of the infall phase;

• 12.

  their lifetimes are ~40% shorter: (4.5 ± 0.1) Myr

• 13.

  fragmentation is suppressed completely.

The comparison with masses and radii from observed young discs favours systems that are somewhere between the two extreme cases “hydro” and “MHD” studied here. In all cases the discs are found to be massive for a relatively long time. This could explain the apparent lack of planet-forming material (e.g. Manara et al. 2018). The discrepancy with observed lifetimes may in part be explained by the choice of “time zero”.

Movies

Access here

Access here

Access here

Acknowledgements

We thank Lucio Mayer, Eduard Vorobyov, Ken Rice, Til Birnstiel and Matthew Bate for the insightful discussions. We also thank the anonymous referee for valuable comments. Calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern. This work has been carriedout within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. O.S., C.M. and G.-D.M. acknowledge the support from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. R.H. acknowledges support from SNSF grant 200020_188460. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1).

Appendix A Photoevaporation

Here we give the expressions for the sink terms used in our implementation of photoevaporation (see Sect. 2.4).

External photoevaporation: $$S_{\text{ext}}(r,t)=\{\begin{array}{ll}{S}_{\text{wind}}\left(1-\frac{1}{1+s{m}_{\text{ext}}^{20}}\right),\hfill & \text{if }r>0.1{\beta }_{\text{M}}{r}_{\text{mI}}\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$(A.1)

where $$s{m}_{\text{ext}}=\frac{r}{{\beta }_{\text{M}}{r}_{\text{g},\text{ext}}},$$

$r_{\text{g},\text{ext}}(t)=G{M}_{*}(t)/c_{\text{s},\text{ext}}^2$ the gravitational radius. We use a smoothing term for numerical reasons. We choose c_s,ext to be 2.5 km s⁻¹, β_M = 0.14 and use S_wind = 2.8 × 10⁻⁸ g cm⁻¹ yr⁻¹ for all systems. This value corresponds to an evaporation rate of 10⁻⁸ M_⊙ yr⁻¹ if the disc extends to 1000 au and is comparable to previous studies (e.g. Armitage et al. 2003 and Mordasini et al. 2009).

Internal photoevaporation: $$S_{\text{int}}(r,t)=2c_{\text{s},\text{int}}{d}_r\text{u}s{m}_{\text{int}},$$(A.2)

where c_s,int = 11.1 km s⁻¹, $$d_{\text{r}}=1.8\times {10}^4{(M_{*}(t)/M_{\odot })}^{-0.25}{\left(\frac{r_{\text{g},\text{int}}}{{10}^{14}}\right)}^{-1.5}{\left(\frac{r}{r_{\text{g},\text{int}}}\right)}^{-2.5}.$$

sm_int is a smoothing factor for which it holds: $$s{m}_{\text{int}}=\{\begin{array}{ll}1-{\left(1+{\left(\frac{r}{0.14r_{\text{g},\text{int}}}\right)}^{20}\right)}^{-1},\hfill & \text{if }r>r_{\text{int}}\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$(A.3)

where $r_{\text{int}}=0.07r_{\text{g},\text{int}},r_{\text{g},\text{int}}(t)=G{M}_{*}(t)/c_{\text{s},\text{int}}^2$.

Appendix B Infall location

Here we derive the relationship between R_i and R_62.3 used in Sect. 2.10.2.

Consider a Keplerian disc with a surface density profile of the form ${\Sigma }_{\text{d}}(r)={\Sigma }_0{\left(\frac{r}{r_0}\right)}^{-1}$, for some constants Σ₀ and r₀. The gas at radius r will have specific angular momentum j: $$j=r^2{\Omega }_{\text{K}},$$(B.1)

where ${\Omega }_{\text{K}}=\sqrt{\frac{G{M}_{*}}{r^3}}$ is the Keplerian angular frequency. We can therefore calculate the mean specific angular momentum ⟨j_d ⟩ of this disc: $$\begin{array}{ll}\langle j_{\text{d}}\rangle =\hfill & \frac{2\pi }{M_{\text{d},\text{tot}}}{\displaystyle {\int }_{R_{\text{in}}}^{R_{\text{out}}}\text{}}rj{\Sigma }_{\text{d}}(r)\text{d}r\hfill \\ \hfill & \cong \frac{4\pi {\Sigma }_0\sqrt{G{M}_{*}}}{3M_{\text{d},\text{tot}}}{R}_{\text{out}}^{3/2}\hfill \\ \hfill & \cong 2/3\sqrt{G{M}_{*}{R}_{\text{out}}},\hfill \end{array}$$(B.2)

where R_in and R_out are the disc’s inner and outer edge, respectively. We used R_in ≪ R_out and the disc’s total mass M_d,tot≅2πr₀Σ₀R_out above.

If all the disc material were concentrated at one radius, this would be the radius characterised by a specific angular momentum ⟨j_d ⟩. We choose this radius to be our infall radius R_i. Combining Eqs. (B.1) and (B.2) yields: $$R_{\text{i}}=\frac{4}{9}{R}_{\text{out}}.$$(B.3)

Now we need to link R_out to R_63.2. We note that for the mass as a function of the radius in the disc it holds: $$M_{\text{d}}(r)\cong 2\pi {\Sigma }_0{r}_{0}r$$(B.4)

for r ≫ R_in. For the mass contained inside of R_63.2 it must therefore hold: $$\text{nnbegin}aligned2\pi {\Sigma }_0{r}_0{R}_{63.2}\text{nn \&}\cong M(R_{63.2})\equiv 0.632 M_{\text{d},\text{tot}}\text{nn \&}\cong 0.632\times 2\pi {\Sigma }_0{r}_0{R}_{\text{out}},\text{nnend}aligned$$(B.5)

which leads to: $$R_{\text{out}}\cong R_{63.2}/0.632$$(B.6)

Plugging this result into Eq. (B.1) we finally find: $$R_{\text{i}}\cong 0.7 R_{63.2}.$$(B.7)

Appendix C More examples of system evolution

Here we show two additional examples for the temporal evolution of star-and-disc systems (compare to Sect. 4). Figure C.1 shows the evolution of system 0004 from RUN-2 (“MHD”), an example of a very low mass system. The simulation starts with a 3.8 × 10⁻² M_⊙ star and a 1.5 × 10⁻² M_⊙ disc (top left panel). The stellar mass reaches 0.05 M_⊙ at t_NIR (≈ 0.95 Myr): a brown dwarf at the lower-mass end considered in our simulations. Nevertheless, the surface density (hence also the temperature, right top and middle panels) is already high (≫ 1000 g cm⁻²) at early stages. This is a consequence of the disc being very compact (~10 au initially). The stellar luminosity is dominated by accretion during the first 100 kyr and decreases steadily after the end of the infall phase (as does the accretion rate, left bottom and middle panels). Q_Toomre never drops much below 2 (right bottom panel), as a result the disc does not fragment and self-gravity plays a minor role in the system’s evolution.

Figure C.2 displays the evolution of system 9992 from RUN-5 (MJ), an example of a massive system. Initially it consists of a ≈ 0.09 M_⊙ star and a ≈ 0.03 M_⊙ disc. At the end of the simulation (t_NIR ≈ 19 Myr) the stellar mass is 5 M_⊙, the upper limit in our study (top left panel). The MCC is feeding the disc with material for 200 kyr. The disc fragments 72 times during the infall phase, but not afterwards. This is because in RUN-5 no accretion heating is considered, which facilitates fragmentation during and prevents a drop in temperature after the infall phase. Furthermore, M_J,FR is used as the initial fragment mass in this run. Therefore, fragmentation is not limited by the initial fragment mass. Approximately 2.1 M_⊙ of mass is removed from the disc by fragmentation. This is substantial: almost one third of the total system mass. Like in the example from Sect. 4, the Gammie criterion never limits fragmentation (right bottom panel). The accretion rate of disc material onto the star (left middle panel) shows an interesting oscillating feature between 4 and 10 kyr. This is not caused by the accretion of clumps onto the star (not included in the figure) but by a decrease in accretion rate that starts when the conditions for fragmentation are satisfied somewhere in the disc. A substantial amount of mass (≈ 10 M_♃ or around 10% of M_disc at this time) is removed from the disc during a free-fall time (see Sect. 2.9.2). This reduces the accretion from disc to star. Surface densities and temperatures (right top and middle panels) are significantly higher in this system than in the first example (Sect. 4) due to the higher stellar mass and disc mass.

Fig. C.1 Time evolution of system 0004 from RUN-2 “MHD”. Top left panel: stellar mass, disc mass and, cumulative mass removed from the disc by fragmentation. Middle left panel: accretion and outflow rates. Bottom left panel: stellar luminosity. Top right panel: surface density at different times. Middle right panel: midplane temperature at different times. Bottom left panel: stellar luminosity (accretion and intrinsic). Bottom right panel: contour plot of the surface density with fragmentation criteria (see main text). The black vertical lines denote, in order of increasing time: tinfall ≈ 1.02 kyr, tpms ≈0.75 Myr and tNIR ≈ 0.95 Myr (see Sect. 4).

Fig. C.2 Time evolution of system 9992 from RUN-5. Top left panel: stellar mass, disc mass and, cumulative mass removed from the disc by fragmentation. Middle left panel: accretion and outflow rates. Bottom left panel: stellar luminosity. Top right panel: surface density at different times. Middle right panel: midplane temperature at different times. Bottom left panel: stellar luminosity (no contribution from accretion considered in this run). Bottom right panel: contour plot of the surface density with fragmentation criteria (see main text). The black vertical lines denote, in order of increasing time: tinfall ≈ 215 kyr, tpms ≈425 kyr and tNIR ≈ 19.1 Myr (see Sect. 4).

Appendix D Infall radii in RUN-2

Table D.1 gives the infall radii used in each stellar mass bin for RUN-2 (see Sect. 3).

Appendix E Initial conditions

The 35 systems from B18 we used to obtain the probability distributions for M_*,i, M_disc,i and Ṁ_in are given inTable E.2. The 20 systems used for the distributions of R_disc,i and the expansion rate of R_i are given in Table E.1. The numbers given in these tables correspond to the protostar indices used in B18.

Figure E.1 shows kernel density estimates for M_*,i, M_disc,i, Ṁ_in, R_disc,i, b_disc, and t_infall. In order to create initial conditions for our runs, we did not use these estimates directly. Instead we used two multi-variate distributions, one combining M_*,i, M_disc,i, and Ṁ_in, the other one combining R_disc,i and b_disc.

Fig. E.1 Kernel density estimate for the initial distributions (neglecting correlation). Top left: protostellar mass. Top right: disc mass. Middle left: infall rate. Middle right: disc radius. Bottom left: expansion rate of the infall radius. Bottom right: tinfall.

Appendix F Reduced disc lifetimes

Here we show our results for the disc lifetimes when the reduction for the start of the PMS phase was applied (see Sects. 2.8 and 5). Table F.1 gives the global mean lifetimes for all runs, Fig. F.1 shows the corresponding distributions. The effect of the reduction is to decrease the disc lifetimes by t_PMS, which depends on the accretion history of the system as explained in Sect. 2.8. This leads to a stronger reduction for the runs based on hydrodynamic initial conditions (the PMS phase is reached later due to the higher disc masses and, hence, higher accretion rates). It is visible in the table as well as in the left panels of the figure. RUN-1 and RUN-2 are closer together although lifetimes in RUN-1 are still higher. The effect on runs 1, 4 and 5 is very similar, so the lifetimes in these runs are still almost the same, as seen both in the table as well as in the right panels of the figure.

Fig. F.1 Disc lifetimes with a reduction for the start of the PMS phase. Top: distribution of disc lifetimes, reduced by tpms (see Sect. 2.8). RUN-1 and RUN-2 are shown in the left, RUN-1, RUN-3, RUN-4, and RUN-5 in the right panel. RUN-1, RUN-4, and RUN-5 are almost indistinguishable. Bottom: same runs, tlife as a function of final stellar mass.

Appendix G Results in specific mass bins

Here we show the results from Sects. 5.1–5.3, the disc masses and disc radii at the end of the infall phase, as well as the disc lifetimes, in three specific mass bins. These are: 0.48 −0.5 M_⊙, 1.00 −1.04 M_⊙ and 1.44 −1.51 M_⊙. We kept the same x-axis in all three bins for comparison. Furthermore, some of the distributions are very narrow (see for example t_NIR, compare to the bottom left panel of Fig. 5). Sometimes all values fall into the same bin, hence the “spiky” nature of the figures shown here.

Fig. G.1 Distributions of disc masses (left column) and disc radii (right column) at the end of infall for specific mass bins. From top to bottom: 0.5 M⊙, 1 M⊙, and 1.5 M⊙.

Fig. G.2 Distribution of tNIR (left column) and tlife (right column) for specific mass bins. From top to bottom: 0.5 M⊙, 1 M⊙, and 1.5 M⊙.

References

All Tables

All Figures

	Fig. 1 Stellar masses at the end of the simulation for all runs (see Table 1), compared to the Chabrier (2007) IMF.
In the text

Fig. 2 Time evolution of System 6410 from RUN-1. Top left panel: stellar mass, disc mass and, cumulative mass removed from the disc by fragmentation. Middle left panel: accretion and outflow rates. Bottom left panel: stellar luminosity. Top right panel: surface density at different times. Middle right panel: midplane temperature at different times. Bottom left panel: stellar luminosity (accretion and intrinsic). Bottom right panel: contour plot of the surface density with fragmentation criteria (see text). The black vertical lines denote, in order of increasing time: tinfall ≈ 30 kyr, tpms ≈3.7 Myr and tNIR ≈ 14 Myr (see Sect. 4).

In the text

	Fig. 3 Distributions of the disc properties at the end of infall and lifetimes of the discs for RUN-1 and RUN-2. Top left: disc masses, top right: disc radii including the observational result of Tychoniec et al. (2018) (Ty18), bottom left: disc-to-star mass ratio qinfall, bottom right: disc lifetimes. All stellar masses are included in this figure.
In the text

	Fig. 4 Properties at the end of the infall phase and lifetimes of the discs for RUN-1, RUN-3, RUN-4 and RUN-5. Top left: masses. Top right: radii. Bottom left: disc-to-star mass ratio qinfall. Bottom right: lifetimes.
In the text

	Fig. 5 End-of-infall properties and lifetimes of the discs as a function of the final stellar mass M*,final. Left column: RUN-1 and RUN-2. Right column: RUN-1, RUN-3, RUN-4 and RUN-5. We show, from top to bottom, masses, radii, disc-to-star mass ratios qinfall at the end of the infall phase, and lifetimes.
In the text

	Fig. 6 Fraction of discs that fragment versus final stellar mass for all runs. The “MHD” run does not produce any fragments.
In the text

	Fig. 7 Fragmentation properties of all runs. Left column, from top to bottom: distribution of the number of fragments, mean initial fragment mass, and mean fragmentation location. Right column: same properties as a function of final stellar mass.
In the text

Fig. 8 Comparisons of our results to previous studies and observations. Comparison of, top left: maximum disc masses from RUN-1 with those from the simulations of B18. Top right: disc radii at the end of the infall phase with observed Tobin et al. (2020, To20) and simulated B18 Class 0 radii. Bottom: disc fractions as a function of time with fits from Richert et al. (2018, Ri18); left: tNIR, right: tlife (reduced, see Sect. 2.8). RUN-4 and RUN-5 are indistinguishable from RUN-1 in the bottom right panel and omitted.

In the text

Fig. C.1 Time evolution of system 0004 from RUN-2 “MHD”. Top left panel: stellar mass, disc mass and, cumulative mass removed from the disc by fragmentation. Middle left panel: accretion and outflow rates. Bottom left panel: stellar luminosity. Top right panel: surface density at different times. Middle right panel: midplane temperature at different times. Bottom left panel: stellar luminosity (accretion and intrinsic). Bottom right panel: contour plot of the surface density with fragmentation criteria (see main text). The black vertical lines denote, in order of increasing time: tinfall ≈ 1.02 kyr, tpms ≈0.75 Myr and tNIR ≈ 0.95 Myr (see Sect. 4).

In the text

Fig. C.2 Time evolution of system 9992 from RUN-5. Top left panel: stellar mass, disc mass and, cumulative mass removed from the disc by fragmentation. Middle left panel: accretion and outflow rates. Bottom left panel: stellar luminosity. Top right panel: surface density at different times. Middle right panel: midplane temperature at different times. Bottom left panel: stellar luminosity (no contribution from accretion considered in this run). Bottom right panel: contour plot of the surface density with fragmentation criteria (see main text). The black vertical lines denote, in order of increasing time: tinfall ≈ 215 kyr, tpms ≈425 kyr and tNIR ≈ 19.1 Myr (see Sect. 4).

In the text

	Fig. E.1 Kernel density estimate for the initial distributions (neglecting correlation). Top left: protostellar mass. Top right: disc mass. Middle left: infall rate. Middle right: disc radius. Bottom left: expansion rate of the infall radius. Bottom right: tinfall.
In the text

	Fig. F.1 Disc lifetimes with a reduction for the start of the PMS phase. Top: distribution of disc lifetimes, reduced by tpms (see Sect. 2.8). RUN-1 and RUN-2 are shown in the left, RUN-1, RUN-3, RUN-4, and RUN-5 in the right panel. RUN-1, RUN-4, and RUN-5 are almost indistinguishable. Bottom: same runs, tlife as a function of final stellar mass.
In the text

	Fig. G.1 Distributions of disc masses (left column) and disc radii (right column) at the end of infall for specific mass bins. From top to bottom: 0.5 M⊙, 1 M⊙, and 1.5 M⊙.
In the text

	Fig. G.2 Distribution of tNIR (left column) and tlife (right column) for specific mass bins. From top to bottom: 0.5 M⊙, 1 M⊙, and 1.5 M⊙.
In the text