© A. J. Cridland et al. 2022

1 Introduction

Some form of planet formation undoubtedly begins early in the life cycle of stars. This statement is largely based on the fact that protoplanetary disks – the source of gas for giant planets – tend to last no longer than about 10 Myr before dissipating (Haisch et al. 2001). Moreover, recent high-resolution surveys of the dust in protoplanetary disks (van der Marel et al. 2013; Zhang et al. 2016; Isella et al. 2016; Fedele et al. 2018; van Terwisga et al. 2018; Andrews et al. 2018; Andrews 2020; Sheehan et al. 2020) have shown a ubiquity of substructure in extended disks that is often attributed to the presence of massive planets. This includes substructure in the young class I objects GY 91 (Sheehan & Eisner 2018), IRS 63 (Segura-Cox et al. 2020) and the well-known HL Tau disk (ALMA Partnership et al. 2015), all of which are no older than 500 kyr. Furthermore, the few known planets embedded in disks that have been directly observed (i.e. Keppler et al. 2018; Haffert et al. 2019) are relatively massive planets and show signs of effectively being at the end stage of their formation. This is because their inferred accretion rate is very low and is unlikely to contribute much more to their total mass¹.

This observational evidence suggests that the bulk of the planet formation process does not occur during the protoplanetary disk phase (typically called class II disks), but in an earlier phase when the disk is still surrounded by its nascent envelope of gas and dust: its protostellar stage (Class O/I). In the protostellar stage, the mass flux in the system is split between gas and dust falling onto the embedded disk from the envelope, and mass accretion onto the star through the disk (Hueso & Guillot 2005). Classes 0 and I are either differentiated by their observed properties such as bolometric temperature and submillimeter excess (Andre et al. 1993; Chen et al. 1995) or by physical evolution. In class 0, most of the mass is contained in the envelope, while class I systems contain most of their mass in the star and in the embedded disk (van Kempen et al. 2009).

Recent surveys by Tobin et al. (2020) and Tychoniec et al. (2020) studied the distribution of dust masses in class 0, I, and II objects in Orion and Perseus with ALMA² and the VLA³ and compared them to inferred exoplanet metallicities from Thorngren & Fortney (2019). They find that the average class II disk does not contain enough total dust mass to reproduce the observed population of giant exoplanets. Class 0 and I objects, on the other hand, do contain enough dust mass in their embedded disks to build a giant planet, and they require a rather realistic formation efficiency ranging between ~ 10 and 30% to do so. Thorngren & Fortney demonstrated that at the very least, the physical processes in the class O/I phase (at evolution times <150 kyr; Dunham et al. 2014) of protostellar evolution lead to the disappearance or conversion of dust mass prior to the class II phase (ages >500 kyr). In this work we investigate whether planetesimal formation, which represents the first step in the bottom-up planet formation model known as core accretion, can contribute to the observed loss of dust mass between class 0, class I, and class II objects.

Core accretion broadly describes the growth of a giant planet by first accumulating up to a few tens of Earth masses (M_⊕) of solid material before then accreting a significant amount of gas until it reaches its final mass, which is about that of Jupiter (=317 M_⊕) or higher. This process differs from the other popular model of planet formation, gravitational instability (GI), which might also be active in the embedded disk during the collapse phase of star formation. The general assumption of GI is that the resulting giant planet reflects the metallicity of its natal disk and does not generate a significant solid core. This assumption means that GI is incompatible with observed planetary populations (as seen in Thorngren & Fortney 2019). However, several numerical studies have shown that through gravitational and/or hydrodynamic interactions between the solid component of the disk and the gravitationally unstable clumps, planets with higher metallicity can be generated, and might include a solid core (Helled et al. 2006; Boley & Durisen 2010; Boley et al. 2011). While we recognise that embedded giant planets could have been generated through GI, we therefore focus this work on the production of planetesimals because they might represent the seeds of giant planets and certainly represent the seeds of asteroids and comets.

The preferred scenario of building planetesimals from dust grains is known as the streaming instability. This instability arises in the regions of a disk in which the dust is sufficiently dense relative to the gas to have a dynamical effect on the gas flow. This effect locally increases the gas flow velocity and reduces the radial drift speed of the dust. The instability is then fed with more dust from regions in which the radial drift rate remains faster, which further increases the dust density in the unstable region. For strictly radial perturbations, however, the instability is suppressed by strong thermal pressure gradients. Perturbations in either the vertical or azimuthal direction are therefore also needed for the instability to grow. These perturbations can lead to related instabilities such as the settling instability described in Squire & Hopkins (2018, 2020). This build-up of dust density further feeds the instability, which causes a runaway process in which the dust density rapidly increases until it becomes gravitationally unstable and collapses into asteroid-sized bound objects.

The pioneering work on the streaming instability was presented by Youdin & Goodman (2005). The authors showed that the optimal environment for the streaming instability to occur was dust grains that are moderately coupled to the gas. These types of grains are characterised by a Stokes number (St) that is about 0.1. St is the ratio of the stopping time of a dust grain to the relevant dynamic timescale of the system. This means that a grain with a very small Stokes number effectively follows the flow of the surrounding gas, while a grain with a larger St is subject to the drag force from the gas, but mainly evolves following other external forces such as gravity. In addition, the streaming instability requires a relatively high ratio of the midplane dust to gas mass; this reportedly lies between 0.2 and 3 for dust grains with St ~ 0.1. This high dust-to-gas ratio is the main limitation of the streaming instability, considering that the standard interstellar medium (ISM) ratio is 0.01. We here largely follow the prescription of Youdin & Goodman (2005), with small changes discussed below.

Recent numerical work by Li & Youdin (2021) has updated the estimates for optimal dust-to-gas mass ratios for a wide range of St. It is our intention to compute the collapse of a dusty molecular cloud, which is initialised with the typical ISM dust-to-gas mass ratio, ϵ_ISM = 0.01, to show that the mass ratio in some regions in the young embedded disk exceeds the critical dust-to-gas ratio and can drive the instability. In regions of the disk in which the threshold dust-to-gas ratio can be reached, the streaming instability could turn the available dust mass into planetesimals.

An early attempt to model the dusty collapse of a protostellar cloud was reported by Morfill et al. (1978), who pointed out that computing the evolution of the cloud’s ‘solid phase’ is the only way to understand the initial conditions that set the formation of the planets and meteorites. This work can thus help us to understand planet formation as a whole, but can be further extended to understand the formation and chemical properties of known comets and meteorites as well. Morfill et al. (1978) reported that the dust does indeed fall rapidly, forming a dusty embedded disk with enough dust sedimentation that further instabilities may lead to rapidly planetary growth and the generation of planets even before the Sun is formed. While this final point may exaggerate how quickly planets can form in protostellar systems, the idea that planet formation should occur while the young star is still in its embedded phase has not received much attention since this pioneering work.

The recent study of planetesimal formation in a young stellar system by Drążkowska & Dullemond (2018) assumed that the dust in the envelope was strictly monomer-sized (about 0.1 µm) and that the growth and subsequent enhancement of dust density was driven exclusively by processes in the embedded disk. Their method of enhancing the dust-to-gas ratio relied on the well-known traffic-jam effect (Pinilla et al. 2016) at the water-ice line of the embedded disk. This effect drives high dust-to-gas ratios because the removal of the water-ice layer from the grains makes them more susceptible to fragmentation. The maximum grain size is thus reduced by an order of magnitude (as seen in Birnstiel et al. 2012), which reduces the radial drift rate by a corresponding amount. Hence the radial drift rate leaving the ice line is an order of magnitude lower than the approaching rate, which drives a density enhancement that is ideal for the streaming instability. Drążkowska & Dullemond (2018) reported that only a minimal number of their final distribution of planetesimals were formed during the class I phase of their model, and the majority were still being formed in the class II phase.

The assumption of strictly monomer-sized dust grains in the envelope is not supported by continuum observations of young stellar systems. Preliminary studies in the youngest dense cores already shows evidence of grain growth up to ~1 µm sized grains through infrared core-shine (Pagani et al. 2010). In this process, background light is strongly scattered by larger grains in the core. In the slightly more evolved prestellar core, L1544, Chacón-Tanarro et al. (2019) found changes in the millimetre opacity towards the inner ~2000 AU, which is indicative of grain sizes up to 3–4 µm. Surveys of class 0 envelopes have found low dust emissivity indices that are consistent with millimeter-sized dust grains (Jørgensen et al. 2009; Kwon et al. 2009; Galametz et al. 2019). Class I objects also show similarly small emissivity indices that indicate dust grain growth up to between 100 µm and millimeter-sized grains (Agurto-Gangas et al. 2019). Furthermore, modelling efforts at submillimeter wavelengths require maximum grain sizes up to ~1 mm to explain the observed continuum emission profiles (Miotello et al. 2014), and polarized dust emission studies have shown evidence for grain growth up to at least 10 µm during the early stages of the collapse (t ~ 10⁵ yr; Valdivia et al. 2019).

Theoretical models can account for some of this grain growth within the lifetime of the collapse. For example, Guillet et al. (2020) reported hydrodynamical turbulence-driven grain growth of up to 1-10 µm in one free-fall time of their collapse simulation. Ormel et al. (2009) argued that dust growth up to a size of 100 µm is only possible if additional support mechanisms such as ambipolar diffusion delay the collapse to longer than a free-fall time, or if a significant ice mantle is present on the grains. Otherwise, dust growth is inefficient during the collapse phase of their models. Because the typically cold temperatures of the outer envelopes of pre- and protostellar objects and the fact that water ice makes up a large fraction of the overall oxygen abundance (Boogert et al. 2015; van Dishoeck et al. 2021), the latter of the above two requirements is very likely. Even when this requirement is met, however, it is still an open question whether grain growth within the envelope can produce the observed approximately millimeter grain sizes.

More detailed hydrodynamical studies of a dusty collapse, including Bate & Lorén-Aguilar (2017) and Lebreuilly et al. (2020), have demonstrated that dust grains larger than about 10 µm sufficiently decouple from the gas to dramatically alter their collapse trajectory. In particular, Bate & Lorén-Aguilar (2017) reported a hierarchy of dust dynamics that is dependent on the grain size. Monomer-sized grains are fully coupled to the gas and maintain nearly the same dust-to-gas ratio throughout the collapse. Grains of about 100 µm are sufficiently decoupled that they do fall faster than the pressure-supported gas, and they gather in the midplane of the system to form the first embedded dusty disk. Their largest grains (1 mm) are nearly completely decoupled and rapidly fall onto the midplane. Near the midplane, these large grains rapidly oscillate around z = 0, causing maximum dust-to-gas ratios that exceed 10⁵ for a short time during the embedded phase. Lebreuilly et al. (2020) similarly showed that millimeter-sized grains produce the highest dust enhancement in their magnetised collapse simulations. They additionally reported that magnetised models increased the amount of decoupling between gas and dust relative to purely hydrodynamic collapse models because of the additional magnetic pressure that was supplied to the gas. This additional source of pressure, however, does not drastically alter the grain-size-dependent trajectories found in purely hydrodynamic simulations.

These simulations demonstrate the key feature that we wish to test with the semi-analytical model we present here: that differences in the gas and dust trajectories caused by the decoupling of large grains lead to a build-up of dust density that is sufficiently large to induce the streaming instability and might even lead to the generation of the first generation of planetesimals. The simplicity of our semi-analytic approach allows us to extend the above studies by directly connecting our results to the physics of planetesimal formation and the global mass evolution of embedded disks. Our simulations focus on a single grain size, that is, we ignore the impact of dust grain growth during the collapse. We test a range of dust grain sizes to probe the effect of different levels of hydrodynamic coupling on the clumping and future collapse due to the streaming instability. The paper is organised as follows: in Sect. 2 we introduce our computational method and outline the physics on which we focus. In Sect. 3 we present the results of our collapse calculations, including estimates of the evolution of the local dust-to-gas ratio. In Sect. 4 we discuss the relevance of our results and estimate the mass conversion of dust mass to planetesimals, while in Sect. 5 we present our conclusions.

2 Methods

The key physics that we wish to capture with our model is shown in Fig. 1. The gas and dust are subject to different forces during their infall phase, which changes the way in which they collapse towards the midplane. While the gas is largely held in place by a hydrostatic balance between gravity and gas pressure (until the arrival of the rarefaction wave), the dust is not subject to gas pressure like this. Dust grains begin their collapse immediately, while the gas must wait for the loss of gas pressure.

While the dust moves through the gas, however, a headwind drains angular momentum and alters its trajectory with respect to pure projectile motion. This effect is especially strong when a dust grain enters the embedded (rotationally supported) gas disk, where it quickly begins to radially drift through the disk. Dust mass is thus removed from the midplane of the embedded disk (the optimal location for the generation of planetesimals) and is replenished by incoming dust that accretes onto and then settles down from the upper surface of the embedded disk.

2.1 Gas model

The backbone of this work is the semi-analytic collapse model of Visser et al. (2009). We outline some of the primary features of the model that are important for our current work. The model computes the collapse of an initially isothermal, spherical, and axisymmetric cloud of gas following the theoretical work of Shu (1977) and Terebey et al. (1984). The slow rotation and angular momentum conservation forces the infalling gas towards the equatorial plane, forming an embedded disk. Our model follows the work of Cassen & Moosman (1981) in handling the flux of gas onto the central star and the disk. The mass flux onto the disk modifies the usual diffusion equation to (Lynden-Bell & Pringle 1974; Dullemond et al. 2006) (1)

where Σ is the disk gas surface density, and S = 2Nρ_su_z is the mass flux from the envelope through both the bottom and top surfaces of the disk, thus the factor of 2. ρ_s is the volume mass density of the gas at the disk-envelope interface, and N is a normalisation factor that ensures an equilibrium between the mass accretion rates onto the disk, through the disk, and onto the central star. In the Visser et al. (2009) model, the interface between the disk and envelope is defined at the level at which the gas pressures in the disk and envelope are equal. The radial speed of the gas while in the disk in the direction of the cylindrical radius R is (2)

where (3)

is the viscosity under the standard a-prescription (Shakura & Sunyaev 1973), with α = 0.01. Equations (1) and (2) are only valid if the specific angular momentum of the inflowing gas is the same as in the disk where it is delivered. This is not always true (Lin & Pringle 1990), and our model contains an additional term in Eq. (2) to account for the redistribution of angular momentum of incoming gas based on Hueso & Guillot (2005). Equation (2) thus becomes (4)

where η_r = 0.002 is a numerical factor that reproduces the results of Hueso & Guillot (2005) best.

The disk-gas scale height H and disk sound speed c_s,d are linked through the Kepler rotation rate, (5)

where M(R, t) is the total mass inward of R at time t, including the host star such that (6)

The cloud initially has a solid-body rotation rate of Ω₀ and an initial gas volume density that varies with the spherical radius r (Shu 1977), (7)

where G is the gravitational constant, and c_s is the effective gas sound speed of the cloud. We set a maximum envelope radius r_env, so that the initial cloud mass is (8)

The structure and evolution of the volume density of the gas depends on the spherical radius r, polar angle θ, and time t through the dimensionless functions , y, and w (Cassen & Moosman 1981; Terebey et al. 1984), (9) (10) (11)

and the dimensionless variables x = r / c_st and τ = Ω₀t. The rotational speed u_ϕ is kept constant in this solution and assumes solid-body rotation at a rate of Ω₀. The values of , y, and w were derived numerically in Terebey et al. (1984)⁴. The cylindrical radius R, height above the disk midplane z, spherical radius r, and the polar angle θ are related through

Broadly speaking, the cloud collapses inside out, with a rarefaction wave (or collapse front) moving outward at the sound speed (Shu 1977). The wave front has the dimensionless position of x = 1 over all time. At a given location with x’ > 1, the gas remains in hydrostatic balance between its gravity and pressure from all gas with x < x’. Once x’ = 1, the shell loses its supporting gas pressure and begins to collapse towards the centre of the cloud. Because of a Taylor expansion in their derivation, when x’ = τ², the Terebey et al. (1984) solution breaks down and the model switches to using the solution derived by Cassen & Moosman(1981).

The Cassen & Moosman (1981) solution follows the trajectory of individual gas streamlines from the envelope to the disk rather than building a self-similar solution, as was done in Shu (1977) and Terebey et al. (1984). The evolution of the gas density, radial, polar, and rotational velocities is (12) (13) (14) (15)

respectively, where is the centrifugal radius, is the total accretion rate from the envelope onto the host star and disk, m₀ = 0.975 is a numerical factor, P₂ is the second-order Legendre polynomial, and θ₀ is the polar angle with which a given streamline enters the Cassen & Moosman (1981) solution. The outer shell of the envelope reaches the star and disk at a time of t_acc = M₀/Ṁ (= 252.5 kyr in our model) and signifies the end of the collapse phase of the protostar. Visser et al. (2009) pointed out, however, that this end point does not exactly correspond to other mature class II systems because the disks produced here tend to be warmer at t = t_acc than the average protoplanetary disk observed in the literature.

The gas collapse was computed for two different envelope initial conditions: the so-called reference case with Ω₀ = Ω_ref ≡ 10⁻¹³ s⁻¹, and the standard case with Ω₀ = Ω_std ≡ 10⁻¹⁴ s⁻¹ (cases 7 and 3, respectively, from Visser et al. 2009). Both simulations were run for the same initial cloud mass M₀ = M_⊙ and the same initial envelope size r_env = 6686 AU, which corresponds to a gas sound speed of 0.26 kms⁻¹ and a gas temperature of approximately 18 K. The main difference between these simulations is that the reference case results in a larger embedded disk, because the majority of the gas falls at any given time near or inward of the centrifugal radius R_c and R_c ~ Ω₀.

Fig. 1 Schematic view of the embedded disk within a collapsing envelope. Unlike in traditional protoplanetary disks, i.e. Class II objects, embedded disks (Class O/I) are continually fed with new gas and dust from the envelope (infall). The gas in the embedded disk is subject to turbulent flow and evolves viscously. The dust, on the other hand, tends to sink to the midplane of the embedded disk and radially drifts towards the host star. Because the gas, but not the dust, is held aloft by gas pressure, the gas and dust arriving at the embedded disk at any given time originate from different initial locations in the envelope; dust arrives from larger initial radii than gas. The grey scale broadly shows the total volume mass density of the gas and dust in the protostellar system.

2.2 Tracing individual dust trajectories

The model includes a postprocessing scheme to trace the trajectories of Lagrangian particles throughout the collapse. This method has been used to trace the collapse of individual parcels of gas in order to study their chemical evolution (Visser et al. 2009, 2011; Harsono et al. 2013; Drozdovskaya et al. 2014). The method was designed to sample the local gas density and temperature, the local gas velocity, and the dust temperature to determine the properties relevant to chemistry as well as the direction in which the parcel will evolve next. We used the fact that the parcel algorithm looks up the local gas density and velocity, and modified it to compute the collapse trajectory of dust in the envelope, which we outline below.

The primary difference between the trajectory of the dust and gas is that dust with sufficiently high Stis not directly impacted by the gas pressure. Dust parcels that are only moderately coupled to the gas, with St ~ 0.1-1, therefore begin their collapse at t = 0 rather than at x = 1. Because it takes a finite time for the collapse front to reach the other regions of the envelope, there is an immediate drift between the collapse trajectories of the gas and dust there. As the dust accelerates towards the centre of the cloud, differences between the gas and dust velocities cause drag that drains angular momentum from the dust, accelerating their collapse towards the equatorial plane and producing a dusty disk with an outer radius that is smaller than that of the gas disk.

The drift of the dust is proportional to the total velocity difference between the gas and dust (Weidenschilling 1977). In the Epstein regime, when the mean free path of the gas is longer than the size of the dust grain (s), the hydrodynamical coupling of the dust grains, known as its stopping time, is defined as (Whipple 1973) (16)

where ρ_s is the density of an individual grain (assumed here to be 3.6 g cm⁻³), and ρ is the volume density of the surrounding gas. The grain St is the dimensionless quantity of the stopping time, defined as St = Ωt_s. In the envelope, we assume that Ω = Ω₀, the solid-body rotation rate of the cloud, while in the disk Ω = Ω_k, we assume Kepler rotation frequency. When the mean free path of the gas is nearly the size of the dust grain, the drag on the dust enters the Stokes regime, but this occurs nowhere in the collapse. In the Epstein regime, the drag force is (Weidenschilling 1977) (17)

where ∆v = u – v is the difference between the gas (u, computed in Sect. 2.1) and dust (v) velocities. The acceleration of the dust due to gas drag is thus (18)

When the dust falls faster than the gas (∆v< 0), a headwind slows its descent, while in the opposite case (∆v > 0), the dust is swept up by the falling gas.

At the beginning of the trajectory calculation, 10⁵ dust parcels were placed randomly throughout the cloud with spherical radii between 100 AU < r < r_env and polar angles between 0 < θ < π / 2. Initially, each dust parcel began with v_R = v_z = 0 and . These parcels began to fall immediately because they were not held up by hydrostatic balance. Dust parcels whose initial position can be described with x > 1 therefore began to collapse before the gas that surrounds them. This setup at t = 0 assumes that the dust and gas are both distributed homogeneously throughout the core prior to the beginning of the collapse. Because of the differential collapse that we expect to see, it is likely too simple for both species to start with the same spatial distribution. In reality, the spatial distribution of the gas and dust are set by turbulent motions in the larger molecular cloud, and they may largely begin with a more compact dust distribution than the gas. Such a distribution would thus lead to a more pronounced differential collapse between the gas and dust than we simulate here.

The equations of motion for the dust parcels are (19) (20) (21)

where r is the magnitude of the spherical radius. We note the units of each equation above to recognise that v_R and v_z are in linear units, while v_ϕ is the rotational speed in angular units. While the simulation itself is in two dimensions (R,z), gas drag depends on the difference between the gas and dust velocities in all three dimensions. Rather than change the underlying code of Visser et al. (2009) to include a third dimension, we assumed that the coordinate system rotates with the rotating dust parcel at v_ϕ;. Therefore, ϕ = 0 for all time in the rotating coordinate system, and we included the Coriolis force (third term in Eq. (20)). The mass contributing to the gravitational acceleration, M(r), includes the cloud and disk mass orbiting inward of the dust parcel, as well as the host star.

We evolved the velocity over every time step following the explicit expression of Rosotti et al. (2016). The velocity in the next timesteps follows (22)

Expressing the velocity evolution in this way ensures that if the time step exceeds the stopping time of the grain, the algorithm remains well behaved. In this limit (Δt >> t_s), the dust parcel ‘forgets’ its previous velocity and evolves with the gas, but with a slight drift proportional to the difference between the gas and dust accelerations and the stopping time. In the opposite limit (Δt >> t_s), the time is too short for the dust drag to alter the trajectory of the dust, and it evolves as if it were in free fall.

The gas evolution, including the gas velocities (u), were precomputed using the model of Visser et al. (2009), that is, we did not self-consistently evolve the gas using a_gas, nor did we include any back-reaction on the gas evolution caused by the clumping of dust. For the purposes of Eq. (22), however, we computed the forces on the gas at the instantaneous position of the dust parcel. Before the collapse wave reached an annulus of gas with x > 1, we assumed that a_gas = 0, u_R = u_z = 0, and u_ϕ = Ω₀. When the dust parcel reached an annulus of gas with x < 1, we computed the acceleration of the gas (23)

where and are unit vectors in the direction of the cylindrical radius and rotation, and is the unit vector in the spherical radius direction. The pressure gradient ∇P and local gas density ρ were taken from the precomputed gas properties of the Visser et al. (2009) collapse model.

Equations (19)-(23) were implemented into a routine that is called by the DVODE integrator (Brown et al. 1989) to integrate the trajectories of the collapsing dust particles. Before the dust trajectories were computed, we computed the evolution of the gas collapse using the methods described in Sect. 2.1. With these precomputed collapse models, we then traced the trajectories of each dust parcel using the method described in this subsection. The routine that computes the dust particle speed uses the current location of the particle to estimate the current gas speed and density from the precomputed collapse model. When the dust particle was in the disk, and when the evolution of the surrounding gas was described by the Terebey et al. (1984) solution, we assumed that the gas rotated around the R = 0 axis with sub-Keplerian and solid-body rotation, respectively. When it was in a region in which the gas evolution can be described by the Cassen & Moosman (1981) solution, we assumed that the gas rotated as shown in Eq. (15). The routine uses the gas speeds as inputs to compute Δv and the current evolution of the dust speeds. A benefit of using DVODE is its adaptive time stepping, which efficiently handles early times when the dust motion is slow, and later times when the dust parcel drifts radially through the embedded disk; this causes the problem to become stiffer.

2.3 Postprocessing the simulations

We computed the trajectory of 100 000 dust parcels, and with these trajectories, we traced the mass evolution of the dusty envelope by tagging each dust parcel with a particular amount of initial dust mass. We then investigated where that mass was delivered. This was done as a postprocessing of the simulation, and we outline it here. The general evolution of a dust parcel is to rapidly fall from the envelope into the disk, after which it begins to radially drift into the central star. The initial collapse rate and later radial drift rate are all dependent on the size of the dust grain. This is discussed in more detail in the following section. For some of the dust parcels, the above algorithm forces very small time steps, which ultimately makes their further evolution intractable. This generally occurs for particles entering the disk at later times when the gas densities and radial drift rates are highest. We extrapolate the evolution of these particles based on their last successfully integrated step and on the assumption that they drift radially in the disk.

2.3.1 Trajectory extrapolation

The extrapolation step takes any dust parcel for which further integration becomes intractable, and takes its position R_last, Z_last and time t_last to estimate and their final position at t_acc. Because all of the parcels that grind to a halt can be found in the disk, the extrapolation was done using two assumptions: the first assumption was that the velocity in the cylindrical radius direction is ‘forgotten’ after one stopping time has elapsed. The radial speed was then replaced by the typical radial drift speed, defined as (Birnstiel et al. 2012) (24)

where Y = |d log P/d log R| is the pressure gradient in the gas, and v_θ = ΩR is the rotational speed in the disk (where Ω is the Kepler frequency). This was done smoothly, so that the extrapolated radial speed is (25)

7

where v_R(t_last) is the last recorded radial speed that the simulation computed for a particular dust parcel, and Δt is the amount of time since t_last where v_R was computed. In this way, the extrapolated radial position is (26)

The second assumption refers to the z direction, where a typical particle trajectory in the disk is a rapid approach to the midplane, then proceeding in a damped oscillation around z = 0. In the simulation itself, particles that cross z = 0 are mirrored in the z > 0 plane, and thus these oscillations are not observed from the typical code output. To account for this oscillatory motion, we extrapolated the z-position of the dust parcel as follows: (27)

where ϕ_last ensures that when Δt = 0, z = z_last, and Ω has the same definition as in Eq. (24). In this way, the extrapolated particles both oscillate and are damped towards the midplane on a dynamical timescale. The damping rate of the dust must also depend on the stopping time of the grain because of the role that turbulent mixing has on slowing the dust settling. Hence the smaller grains may be damped more slowly than the larger grains. However, because the original calculation of the dust collapse ignores the impact of diffusion (discussed below), we also neglected its effect during the extrapolation.

2.3.2 Dust mass assignment

To compute the dust-to-gas mass ratio throughout the embedded disk, we must track the movement of dust mass throughout the collapse. Each dust parcel was assigned a certain dust mass at the beginning of the collapse. This initial assignment was made as follows: first, a logarithmically spaced grid from 0.1 AU out to 10⁴ AU in both the cylindrical radius and height was set at t = 0. The average gas density in each of these cells was computed from the initial gas density distribution (Eq. (7)), and the initial dust density was set to ρ_dust,0 = ϵ_ISMρ_gas,0. The total dust mass in a given grid cell was computed based on the volume⁵ of the cell and its average density. This dust mass was then equally distributed to all dust parcels that were present in the grid cell. We assumed that the dust mass in each parcel is conserved during the collapse, so that no diffusion occurs in or out of the parcels.

In all future time steps, we used the same grid as defined before and located all dust parcels that entered a given grid cell. We summed their dust mass and computed the average dust density in that cell using this mass and the physical volume of the cell. Finally, we used the average gas density in the grid cell from the Visser et al. (2009) model and the dust density to compute the local dust-to-gas ratio in the cell. This was then repeated for each grid cell. Because a dust parcel must be in a given cell at a given time step in order to contribute a dust mass, some regions of the disk will show a remarkably low dust-to-gas ratio in cells without dust parcels. We do not claim that regions like this will exist in practice because diffusion, which we ignored, will play a role in spreading the dust mass more evenly through the disk.

2.3.3 Defining the available mass for the streaming instability

Youdin & Goodman (2005) showed that the optimal environment for the streaming instability is dust grains that are moderately coupled to the gas (St < 0.1) and have relatively high dust-to-gas ratios. Their optimal dust-to-gas mass ratios were 0.2 and 3, with a slight reduction in growth speed at ratios of exactly unity, which was caused by the instability wave speed reaching a minimum. Recent numerical work by Li & Youdin (2021) have updated the estimates for optimal dust-to-gas mass ratios for a wide range of St. They reported an updated relation between the grain St and threshold mass density ratio for efficient particle clumping, leading to an streaming instability of (28)

where ϵ = ρ_d/ρ_g is the dust-to-gas ratio (Li & Youdin 2021). Two regimes show different Stokes-dependences on the mass ratio threshold. Thus their fitted constants have the values A = B = 0 and C = log₁₀(2.5) when St < 0.015, and A = 0.48, B = 0.87, and C = -0.11 in the opposite case. In each cell that contains dust particles, we compared the computed dust-to-gas ratio to the threshold in Eq. (28), and in any cells in which the threshold was reached, we assumed that their dust was available for the streaming instability.

2.4 Caveats

One simplification of our model is the lack of turbulent diffusion acting on our collapsing dust parcels. As a result, we find that dust settling is very efficient in our model of an embedded disk. For the sake of setting the accretion rate throughout the embedded disk, we assumed that the turbulence is α = 0.01, which would imply a high degree of turbulence. Strong turbulence in the vertical direction would tend to act against dust grain settling, in particular for the smallest grains, which would spread the dust mass off of the midplane more smoothly throughout the whole disk. In the standard case of isotropic turbulence, it is generally assumed that the viscosity driving mass accretion throughout the disk and the turbulence in the vertical direction are described with the same a parameter. If the turbulence is not isotropic, or if the angular momentum transport responsible for disk accretion is driven by another physical process such as disk winds, then the vertical a and the disk viscosity a could be different.

The planetesimal formation model of Drążkowska & Dullemond (2018) favoured lower turbulence (α ~ 10⁻⁴) in the vertical direction and intermediate turbulence (α = 10⁻²-10⁻³) in the radial direction (i.e. responsible for disk accretion). Observational evidence also tends to favour lower vertical turbulence (e.g. Pinte et al. 2016) in young protostellar systems, while observational constraints on the role of turbulence in global disk accretion have yet to be confidently found. Our method of assuming a high α for the global evolution of the embedded gas disk while ignoring the impact of vertical turbulence on the evolution of the dust is a simplification, but is consistent with a disk model that drives angular momentum transport through magnetohydrdynamic disk winds, but otherwise has low intrinsic turbulence. Such a model is consistent with observational and theoretical constrains, but generally evolves differently in size and surface density distribution from what we assumed in our model (see e.g. Tabone et al. 2021).

We assumed here that the dust can be described by a single grain size. Young clouds, however, are better described by a distribution of dust grain sizes with a maximum grain size that is similar to thoses that we have modelled in this work. Furthermore, the dust mass distribution can change through coagulation and fragmentation, which greatly complicates the nature of this problem. Depending on the efficiency of grain growth, the mass of the larger grains could be maintained against radial drift for an extended period of time, which changes the situation that we investigate here. The large grains have the highest St and thus require lower dust-to-gas mass ratio thresholds to drive the streaming instability than the smaller grains (we recall Eq. (28)). This would allow that more regions in the midplane of the disk contribute to the production of planetesimals, potentially increasing its efficiency. Including a simple formulation of dust grain growth would be an interesting next step that is currently beyond the scope of this work.

Finally, our simulation ignores the impact on magnetic fields on the collapse trajectories of the gas and the dust. As previously mentioned, Lebreuilly et al. (2020) found that models that include magnetic fields strengthen the drift between collapsing dust and gas, mainly driven by the fact that magnetic fields provide an additional source of pressure support for the gas but not for the (neutrally charged) dust. Dust grains can also carry a charge by absorbing free electrons that are liberated from elemental and molecular sources through thermal and radiative effects (explored e.g. in Akimkin et al. 2015, for H II regions). They would therefore be dynamicaly affected by the large-scale magnetic fields in the envelope during their collapse. Furthermore, the magnetic fields include torques that impact the distribution of the angular momentum of the gas as it collapses. Here, choosing different initial angular momentum models partially accounts for this physical effect because the standard model results in smaller disks, as predicted in magnetised collapses and noted in Visser et al. (2009). The additional gas pressure found by Lebreuilly et al. (2020) might enhance the build-up of dust-to-gas ratio at earlier times because the collapse of the gas would be further delayed.

3 Results: Dusty collapse

Below we present the results of four collapse simulations based on the prescriptions above. Three of the simulations use initial conditions consistent with the standard case from Visser et al. (2009): initial solid-body rotation rate Ω_std = 10⁻¹⁴ s⁻¹, cloud mass M₀ = 1 M_⊙, and average sound speed c_s,0 = 0.26 km s⁻¹. This sound speed is equivalent to an average gas temperature of 18.5 K. The parameter that varies between these three simulations is the assumed size of the dust grain. These sizes are 0.34 µm (consistent with monomer sizes, known as the Std-sml model), 34 µm (consistent with Agurto-Gangas et al. 2019, known as the Std-med model), and 3.40 mm (consistent with Galametz et al. 2019, known as the Std-lrg model). Their respective sizes are consistent with average St at the beginning of the simulation of 0.0012, 0.12, and 12, respectively. We define the initial St in the envelope as (29)

where ρ_g,0 is the average density for the spherical cloud with a radius of 6686 AU, ρ_s = 3.6 g cm⁻³ is the average density of the grains, and is the timescale related to the collapse. This timescale is most relevant to the hydrodynamic properties of the dust because early on in their collapse, the speeds in the direction of the spherical radius tend to be higher than their rotational speeds. As the cloud collapses and the dust grains fall towards the centre, their local St decreases because the gas density tends to increase and we do not invoke any grain growth. When they reach the disk, however, the local dynamical timescale (Ω) increases because orbital rates are faster than radial speeds in a rotationally supported disk, while the gas volume density increases as well. For most of the grain models, the St increases when the grain reaches the disk.

The fourth simulation used the 34 µm dust grains and the same initial cloud conditions, but a higher solid-body rotation rate of Ω_ref = 10⁻¹³ s⁻¹, and is known as the Ref-med model. The standard and the reference cases are known as the spread-and infall-dominated cases in Drozdovskaya et al. (2014) because they reflect the primary method for which the embedded disk grows given our (relatively high) choice of α. The lower initial rotation rate and subsequently lower gas angular momentum of the standard case causes the embedded disk to start small (~ 10 AU) and primarily grow through viscous spreading. In the reference (infall-dominated) case, however, the embedded disk starts larger (~ 100 AU) and further grows through a combination of viscous spreading and subsequent envelope accretion (see Table 1 for final sizes). Because we assume that the dust grains begin their evolution with the same rotation rate as the gas, dust falling in the standard case will tend to fall into the disk closer to the host star than grains falling in the reference case.

We summarise the important physical parameters in Table 1. Estimates of the St are included for each of the dust grains initially (St_env) as well as at a future time when the grain has entered the embedded disk (St_disk). When the grains enter the disk, they encounter a higher gas density, but also much faster rotation speeds, which drives a higher St of each of the models.

3.1 Single dust trajectory

As a first step in our analysis, we present a set of trajectories of the same parcel of dust in each simulation. In this way, we directly compare the effects of the initial angular momentum and grain size on the overall collapse history of the dust.

In Fig. 2 we show the collapse trajectories for a parcel of dust beginning at the same location in each model. In the standard model, nearly all of the different grain sizes follow a similar path, but the rate of the collapse changes as a function of grain size: the point marking 120 kyr (nearly halfway to t_acc) is farther along the collapse trajectory of the larger grains than the smaller ones. Because the larger grains are less strongly coupled to the gas, they collapse more freely than smaller grains, which are partly kept aloft by the gas. When the grains enter the disk, their local St increases because the gas rotates faster and the gas densities are higher. Radial drift becomes a dominant source of further evolution in the embedded disk, and the larger grains tend to drift to smaller radii than the smaller ones.

The higher rotation rate of the reference case keeps the collapsing dust grains farther out in the envelope as it collapses. When in the disk, the particles do not radially drift as fast as the dust in the standard model. This is partially due to the higher outward flux of the gas as the disk is being built, and partially due to the high rotational support that arises from the higher rotation rate. As a result, the general conclusion that envelopes with high rotation rates produce larger embedded disks continues to apply to the dust as much as it does to the gas, even when we consider the impact of radial drift.

In Fig. 3 we show the evolution of the total angular momentum of the dust parcels from Fig. 2. For the first few 10s of kyr, the angular momentum is strictly conserved as the dust begins its collapse because the relative collapse speeds of the gas and dust are low initially. After a while, the dust reaches a region of the envelope with higher gas densities and has accelerated such that its speed relative to the gas is sufficiently high for it to begin to lose angular momentum via gas drag more efficiently. When the dust parcels reach the disk, they begin to interact with a gas of higher density that is rotationally supported. Hence the relative speeds between the disk gas and the dust are high, and the dust is efficiently spun up in the first few time steps after it crosses the boundary of the disk. As a result, the angular momentum of the dust particle suddenly spikes when it enters the disk, but then flattens out when its rotation speed is once again similar to the local gas rotation speed. While the particle is spun up, it begins to radially drift through the disk, losing angular momentum more slowly than during its collapse through the envelope. The evolution during the radial drift stage of the dust evolution can best be seen in the reference model (red line), while the dust parcels in the other models require the aforementioned extrapolation to reach a time of t_acc ~ 252 kyr.

During the short time frame in which the dust regains angular momentum after passing into the disk, the gas should also lose a similar amount of specific angular momentum as the dust gains. This back-reaction may in principle impact the flow of the gas in the embedded disk. However, the total dust mass represented by a given parcel in our model is much lower than the total mass in an annulus of gas from which the dust extracts its angular momentum. We therefore do not expect a significant change in the rotation rate of the gas caused by this momentum exchange.

Fig. 2 Single dust parcel trajectory starting from the same location in the envelope for all four models. Broadly speaking, time flows from the top right of the figure to the bottom left. The coloured point denotes the location of the parcel in each model at a collapse time of 120 kyr. The size of the dust particles in each collapse model denotes how quickly the dust collapses, with larger grains being farther along in their collapse than smaller grains. The different rotation rates of the clouds produces a small change in the envelope collapse, and the embedded disk grows more quickly, which causes a large change in the trajectory when the dust particle encounters the disk.

Fig. 3 Evolution of the specific angular momentum (total angular momentum per mass) of the dust parcels for each of the models. The dashed line highlights the time of 120 kyr that was featured in Fig. 2. Angular momentum is strictly conserved for the first few 10 s of kyr until the dust has fallen into a higher-density region of the cloud, where interactions between the gas and dust begin to efficiently drain angular momentum. When the dust reaches the disk, it enters a region with much higher angular momentum and is quickly spun up by the gas.

3.2 Embedded gas disk formation

In the following section, we report the build up of the embedded dust disk. For reference, we first show in Fig. 4 trajectories for the collapse of the gas particles in the standard model. The method for computing these trajectories is exactly the same as was developed by Visser et al. (2009). The colour-coding of the points shows the location in the initial envelope in which the gas particles originated. Brown points shows initial radii smaller than 1000 AU, dark orange shows initial radii between 1000 and 2000 AU, light orange shows radii between 2000 and 3000 AU, and so on. The four different sets of three panels show four different time steps in the collapse at three different length scales. The left panels shows the envelope length scale, the right panels shows the length scale relevant to the disk, and the middle panel shows an intermediate length scale. The time steps shown here are at t = 0.05, 0.25, 0.5, and 0.75 t_acc.

As the gas collapses, it is either lost to the host star; enters the disk, where it begins to evolve viscously; or is ejected by the outflow. As previously mentioned, the embedded disk in the standard model primarily grows through viscous spreading rather than direct infall. The particles coming from smaller initial radii in the envelope that remain in the disk and slowly expand their cylindrical radius show this. Gas parcels beginning farther out in the envelope take a longer time to reach the disk, and enter the disk at slightly higher heights and at larger radii before further expanding viscously. In the middle of each set of three panels, we show the development of the growing outflow cavity. The size and evolution of this outflow are handled very simply by parametrisation (see Visser et al. 2009, for details). In the trajectory calculation, the gas particles that enter the outflow are immediately removed.

Fig. 4 Collapse and disk build-up of the gas parcels in the standard model. Each set of three panels shows the collapse at 0.05,0.25,0.5, and 0.75 tacc over three length scales. The left panel shows the largest scale, the envelope scale, and the right panel shows the length scale relevant to the disk. The colour of each point denotes the initial radius range from which the parcel originates. Dark brown denotes an initial radius smaller than 1000 AU, orange a radius between 1000 and 2000 AU, and so on. The embedded disk builds up from gas that successively falls from higher heights as well as from the viscous spreading of the disk itself. For reference, the black line shows an estimate of the location of the rarefaction wave defined by x = 1.

Fig. 5 Collapse of dust parcels in the Std-med model. Top left set of three panels: very early time step in the collapse. The colour of each point is the same as in Fig. 4. Each set of three panels shows the same time step at three different lengths scales. The right panel shows the length scales relevant for the embedded disk. The second, third, and forth set represent the simulation at 0.25, 0.5, and 0.75 tacc respectively. As the disk builds, some of the dust drifts radially inward, some is advected outward with the growing gas disk, and missing dust is replenished by newly infalling dust. For comparison, the black line shows an estimate of the location of the gas rarefaction wave defined by x = 1.

3.3 Dust disk formation

3.3.1 Standard model

In Fig. 5 we show the collapse of the dust in the standard model with a dust grain size of 34 (called the Std-med model in Table 1) at t = 0.05, 0.25, 0.5, and 0.75 t_acc. The colour of the points represents the same colours as in Fig. 4 to show the variation in the collaps rate between the gas and dust.

By comparing Figs. 5—4, we confirm that the dust falls faster than the gas because it lacks pressure support. This is shown in the second set of panels (t = 0.25t_acc) in Figs. 4 and 5, where the small gas disk consists of a gas parcel originating from between 1000 and 2000 AU (dark orange), while the dust component of the disk originates from between 2000 and 3000 AU (light orange). In addition, we note the relative location of the coloured parcels with respect to the gas rarefaction wave (solid black line) in the left panel of the second set.

At later times, in the third set of panels (t = 0.5t_acc), the dust disk is made primarily of dust parcels originating from between 3000 and 5000 AU (off-white and light blue points), while the gas disk still contains gas parcels originating from between 1000 and 2000 AU. This reflects the impact of radial drift, which drives the inward evolution of dust particles in the embedded disk. This ensures that the bulk of the embedded disk consists of fresh dust grains falling from the envelope because the colour of the embedded disk changes with every time step. At intermediate time-steps (0.25 and 0.5 t_acc; second and third sets of three panels, respectively), the dust disk follows a similar extent as the gas disk, but at later times (0.75 t_acc; bottom set of three panels), due to radial drift and settling, the dust parcels are not as extended in height as are the gas parcels. Unlike in the gas simulation, only a few dust parcels remain in the envelope by 0.75 t_acc, and at late times, the mass flux onto the embedded disk is therefore dramatically reduced.

Some similar features are seen in the collapse calculations for the standard model with 3.4 mm grains (Std-lrg) and with the 0.34 µm grains (Std-sml); see Figs. A.1 and A.2, respectively. The dust grains in the Std-lrg model evolve very quickly and only develop a very small, compact dust disk before this disk radially drifts into the host star. The dust in the Std-sml model, on the other hand, follows the gas more closely than either of the Std-med and Std-lrg grain models because its coupling is stronger. As a result, a significant amount of dust still remains in the envelope at late times, which will help to maintain the dust mass in the Std-sml model. Its distribution also follows the extent of the gas disk more closely. The Std-sml model shows some parcels being captured into the outflow cavity. In our analysis, we generally ignored dust parcels that entered the outflow, assuming they are carried outward and not inward.

Fig. 6 Collapse of dust parcels in the Ref-med model. The colour and time steps are the same as in Figs. 4 and 5, but the x-axis scale changes in the right panels. While the overall collapse is over a similar timescale as in the Std-med model, the disk build-up and subsequent evolution in the disk is different. The gas disk is larger in the reference model, and dust grains fall predominately at the outward edge of the disk on average. Viscous spreading is still present in the embedded disk, but most of the dust evolution is inward via radial drift. At late times, the dust ring is therefore extended, and most of the dust is inward of the ring in the midplane or at the upper surface of the disk. For comparison, the black line shows an estimate of the location of the gas rarefaction wave defined by x = 1.

3.3.2 Reference model

In Fig. 6 we show the reference model (Ref-med) over the same time steps as in Fig. 5. Because this model uses the same dust grain sizes as the Std-med model, the collapse occurs over a similar time scale. The disk build-up, however, is different mainly because the size of the gas disk in the reference model is larger than in the standard model. Because of the size of the disk, we find that settling and radial drift are both slightly less efficient for the 34 in the reference model, as can be seen by comparing the colour of the dust parcels contributing to the embedded disk at intermediate times. At both 0.25 and 0.5 t_acc, there are dust parcels that originated between 1000 and 2000 AU and 2000 and 3000 AU, respectively, in the Ref-med embedded disk, while dust parcels originating from these radii have completely disappeared in the Std-med model. At 0.75 t_acc, the envelope has been depleted of dust, but the embedded disk remains large, with many dust grains still in the process of settling towards the midplane.

Overall, there may be a sweet spot for planetesimal formation with particular dust grain sizes given their collapse trajectories. Large grains fall very quickly, then settle and drift radially into the host star so that they may not build up enough dust density to spur the streaming instability. Small grains that are highly coupled to the gas stay largely extended and are much more slowly accreted onto the embedded disk. They also do not settle as efficiently in the midplane. Medium-sized grains are between the two extremes: they collapse more quickly than the small grains because they are only lightly decoupled with the gas, and they settle efficiently enough to possibly grow a density enhancement along the midplane. Furthermore, the radial drift rate is slower than that of the large grains, so that the density enhancements last longer in the embedded disk.

3.4 Total mass evolution

One of the important observables in protostellar systems is the total dust mass contained within the embedded disk. In this section we report the evolution of the total embedded dust and gas mass in our models, and compare the embedded dust mass to the observational results of Tychoniec et al. (2020). The total embedded dust mass is computed by summing the assigned mass (as explained in Sect. 2.3.2) of all particle that exist within the embedded disk surface at a given time. The embedded gas mass is computed by integrating the analytic surface density profile that describes the disk out to its outer radius.

Figure 7 presents the evolution of the gas and dust mass within the embedded disk. The black curves denote the temporal evolution of the embedded gas disk mass for the standard (solid line) and reference (dashed line) models. These masses are in units of Earth masses and scaled by the ISM dust-to-gas ratio to better compare the curves to the dust masses (coloured lines). Because the embedded disk in the reference model is physically larger (recall Table 1), its mass also grows higher than the embedded disk in the standard model.

We defined the dust mass by keeping track of all dust parcels that are found inward of the disk surface at a given time, and summing all of their dust masses together. Because the inner disk radius does not extend all the way to R = 0, dust that radially drifts past the inner boundary of the disk is considered lost to the host star. We find that the embedded dusty disk mass slowly increases over nearly the whole lifetime for models with grain sizes ≤34 µm. The largest grains grow to a maximum mass of nearly 100 M_⊕ at t ~ 0.5 t_acc before they rapidly drain due to radial drift. The smaller grains only begin to have a negative flux of dust mass at late times when their envelope has largely dissipated. The left vertical dashed line shows the location in which the embedded disk first appears in the simulation. We denote this time the emergence of our class 0 object.

The class I phase of our simulated systems is defined as the point at which the gas mass in the embedded disk is equal to the remaining gas mass in the envelope. This occurs just over 0.5 t_acc and is denoted by the right vertical dashed line. During this phase, the dust mass flux from the envelope in the Std-med and Std-sml models is sustained because grains of this size are more strongly coupled to the gas than in the Std-lrg model, which delays their collapse with respect to the larger grain sizes. Because the 34 µm grains are less strongly coupled than the 0.34 µm grains, the Std-med model halts its growth earlier than the Std-sml model. The grains in the reference model show a faster build-up than embedded disks in the standard models. Primarily because the gas disk grows larger and faster in the reference model, more dust grains are caught early in the collapse of the Ref-med model, leading to the higher-mass embedded disk. Furthermore, because the settling and radial drift are less efficient in the Ref-med model than in the standard models, the disk holds longer onto its dust, allowing for a larger build-up of mass.

Comparing the evolution of the dust mass to the gas mass in each of the models, we can already see that the reference model may struggle to create planetesimals because its global dust mass is always lower than the mass of the embedded disk gas, scaled by the ISM dust-to-gas ratio. The standard models, on the other hand, appear to favour the generation of planetesimals because their dust masses in the Std-med and Std-sml models are always higher than the scaled embedded disk gas mass. We investigate this further in the following sections.

We include in Fig. 7 the median dust mass (horizontal dashed grey line) and its confidence interval (hashes) inferred by Tychoniec et al. (2020) for a population of protostellar systems studied with ALMA and the VLA. The difference between their inferred median mass and the confidence intervals comes from the different dust opacities used in converting luminosity to mass. We discuss this in more detail in Sect. 4.2.1. We used the python package lifelines⁶ to compute the median dust masses and their confidence interval. This population includes both class 0 and class I, which we separated in order to compute the median dust mass for each protostellar type. The observational data highlight a particular feature that our simulations do not tend to reproduce: class 0 objects tend to have higher dust masses than the older class I objects. We discuss our interpretations of this in Sect. 4.

Fig. 7 Time evolution of the total dust mass inside of the embedded disk in each model over one tacc (~252 kyr). The grey bands show the median disk masses for ALMA-detected (right hatch) and VLA-detected (left hatch) protostars reported by Tychoniec et al. (2020) along with the confidence interval of the median. The ALMA-inferred masses are likely lower limits because the dust is optically thick, while the VLA-inferred masses are more likely to be optically thin and use realistic opacities assuming the existence of large grains. The vertical dashed lines show the time at which the embedded disk first appears (classifying the system as class 0) and when the disk gas mass exceeds the remaining gas mass in the envelope (classifying the system as class I). Class 0 and class I objects were both included by Tychoniec et al. (2020), and we have separated them to compute individual medians and confidence intervals. For comparison, we include the mass evolution of the disk gas in the standard (solid black) and reference (dashed black) models, scaled by the ISM dust-to-gas ratio.

3.5 Dust-to-gas ratio evolution

In Fig. 8 we show the evolution of the dust-to-gas ratio during the collapse of the Std-med model at the same time steps as in Fig. 5. The gas density was computed using the smooth solution of the Visser et al. (2009) model (rather than using the trajectories seen in Fig. 4), and the dust density was computed using the collapse trajectories shown in Fig. 5 that we discussed in Sect. 2.3.2. The colours in the figure denote the local dust-to-gas ratio in the simulation. More orange points show progressively more dust-rich regions of the system, and more blue points show more gas-rich regions.

The set of four panels shows the dust-to-gas ratio in the system at t = 0.05, 0.25, 0.5, and 0.75 t_acc in the top left, top right, bottom left, and bottom right positions of the figure, respectively. The dashed grey line shows the edge of the gas outflow, and the solid grey line shows the surface of the embedded disk, which is defined as the location in which the gas pressure in the disk model exceeds the gas pressure in the cloud. When the gas disk forms, it immediately acts as a trap for the falling dust parcels, and most of the dust is trapped near the outer and upper edge of the disk surface.

The dust in the apparent jammed region at the upper edge of the embedded disk is only temporarily in this region. The speeds of the dust parcels in this region are subsonic, with typical speeds up to a few 100s of cm s⁻¹ (compared to the typical sounds speed there of ~10⁴ cm s⁻¹). This means that while the inward velocity of the collapsing dust parcels is greatly slowed as they transition from the envelope to the disk, they return to accelerating towards the midplane after a brief stay at the top of the disk atmosphere. This brief stay could provide enough time for grain growth to occur at the top of the disk atmosphere. We investigate this in the discussion section below.

Because we ignored the impact of turbulent diffusion when we computed the dust density in the embedded disk, no dust mass is found in some regions. The three main dust components in our models are therefore the aforementioned dust trap at the disk surface, the dust in the midplane, and the still-falling dust in the envelope. In reality, we would expect that turbulent diffusion would mix the dust as it passes from the disk surface to the midplane, and would tend to suppress the settling of the smaller grains. We further discuss the role of turbulence in our model below.

Fig. 8 Dust-to-gas mass ratio in the Std-med collapse model. The colour scale shows thes log of the local dust-to-gas ratio, hence blue denote dust-to-gas ratios below the ISM, and orange denotes dust-to-gas ratios above the ISM. The dashed lined shows the extent of the outflow in the gas model, and the solid line shows the surface of the embedded disk. We assumed that planetesimal formation occurs largely along the midplane, which is highlighted by the grey arrow. We show the midplane dust-to-gas ratio below.

3.5.1 Midplane dust-to-gas ratio

The relevant portion of the model in which we assume that planetesimal formation occurs is near the midplane of the embedded disk. This region is small in Fig. 8 (shown by an arrow in the forth panel), and hence we highlight this region in Fig. 9. The colour of the solid and dashed lines denotes the particular time step in the simulation. The solid lines show the dust-to-gas ratio in the midplane of the embedded disk, and the dashed lines show the threshold dust-to-gas ratio needed for the streaming instability as computed in Eq. (28). At 0.05 t_acc, there is insufficient dust along the midplane to be shown in the figure. The general structure of the dust-to-gas ratio in the midplane of the disk includes a strong peak in dust density near the outer edge of the disk, driven by the recent accretion of dust, a relatively high density of dust between 1 and 10 AU driven by the flux of dust from the outer disk through radial drift, and settling from the upper edge. Finally, inward of 1 AU, the dust density slowly builds with time as it is fed by inward-drifting dust grains. The gas density in the innermost region in the midplane is slowly vacated by the growing inner cavity, and thus radial drift becomes inefficient there.

By 0.25 t_acc, there is already enough dust in the midplane of the disk to raise the dust-to-gas ratio above the ISM value by an order of magnitude outside of 1 AU. By 0.5 t_acc, the outer 10 AU of the embedded disk (20–30 AU) reaches a dust-to-gas ratio of unity, but it does not exceed the threshold needed for the streaming instability at that time. By 0.75 t_acc, the dust-to-gas ratio exceeds the threshold needed for the streaming instability in a region between ~0 and 40 AU.

We reiterate that generally speaking, the dust in a given time step at the top and outer edges of the disk has only recently been accreted. This is shown in Fig. 5 because there is little mixing between dust grains originating from different radii. Dust parcels originating from smaller radii have radially drifted inward before parcels from larger radii arrive. This ensures that the high dust-to-gas ratio regions at intermediate to late times are caused by the arrival of new dust parcels, rather than by the outward diffusion of dust parcels that entered the disk previously. This distinction is important because when dust parcels are in a region with a high dust-to-gas ratio, they would presumably begin to produce planetesimals, which would make them unavailable for further evolution. We do not currently take this sink into account because our interest here is to show that high dust-to-gas regions in the disk are possible given drift-driven deviations in collapse trajectories between gas and dust.

Fig. 9 Midplane dust-to-gas ratios for the Std-med collapse model at the same time steps as shown in Fig. 8. At the earliest times (solid blue curve), the dust disk has yet to form and does not appear in this plot. The black line shows the ISM dust-to-gas ratio, which represents the initial conditions in the envelope, and the coloured dashed lines show the streaming instability threshold at each time step. The hashed region is an estimate of the streaming instability threshold based on Youdin & Goodman (2005), and the dashed lines show the threshold based on Li & Youdin (2021). The threshold is met along the midplane of the embedded disk only at the last time step.

3.5.2 Grain size dependence

The dependence on the size of the dust grain is as expected, and the results of the Std-lrg and Std-sml models are presented in Figs. B.1 and B.3, respectively. The large grains evolve quickly and build up a dust-to-gas ratio no higher than a few factors of ten times the size of the ISM dust-to-gas ratio. Furthermore, Fig. B.2 shows that while the dust disk initially grows from t = 0.25 t_acc to t = 0.50 t_acc, by 0.75 t_acc,the dust disk has drifted radially inward to about 0.8 AU, effectively disappearing.

The very small grains, on the other hand, are more strongly coupled to the gas and hence show a much slower evolution than either the Std-lrg or Std-med models. The midplane dust-to-gas ratio therefore takes longer to develop, as shown in Fig. B.4, the radial drift has not reached an equilibrium with the accretion onto the disk, which causes the dust density between 1 and 10 AU to constantly grow over the collapse of the system. The inner disk (< 1 AU) also takes longer to collect dust. The slower radial drift speed causes the outer disk to become very enhanced in the dust, resulting in a dust-to-gas ratio of almost 10. This high mass ratio likely breaks our assumption that the dust does not have any back-reaction on the evolution of the gas, but for our simple experiment, we ignored the possibility of these effects.

3.5.3 Reference model

In Fig. 10 we show the evolution of the dust-to-gas ratio for the Ref-med model. As already pointed out, this model has a less efficient dust settling, which results in lower dust-to-gas ratios in the midplane. However, because the disk in the Ref-med model is large, more dust is captured by the growing embedded disk than in the Std-med model. The jamming near the disk edge, however, occurs at a lower density in the Ref-med model, and thus less dust is captured there. The dust-to-gas ratio in the midplane in Fig. 11 shows that the radial evolution of the dust is clearly different in this model than in the standard models. Unlike in the standard models, the dust density in the disk midplane grows nearly uniformly over all radii. Furthermore, while the gas disk in the reference model is much larger than in the standard model, most of the dust mass in the midplane extends to roughly the same radii (30–40 AU), as it does in the standard model. It therefore appears that while the amount of angular momentum initially in the cloud has a strong impact on the final size of the gas disk, the final size of the dust disk seems to be less sensitive to Ω₀.

Fig. 10 Same as Fig. 8, but for the reference model Ref-med. Settling is less efficient than in the standard model, and most of the dust grains remain near the upper surface of the disk. Regions just outward of the disk surface along the midplane show sufficiently high dust-to-gas ratios. These regions, however, do not survive to the last snapshot shown here. Planetesimal formation predominantly occurs along the midplane, highlighted by the grey arrow in the forth panel, and shown below.

Fig. 11 Midplane dust-to-gas ratios for the Ref-med collapse model at the same time steps as shown in Fig. 10. The lines and hashed regions have the same definition as in Fig. 9. The dust-to-gas ratio never exceeds the threshold for the onset of the streaming instability at any of the shown time steps.

4 Discussion

4.1 Available mass for planetesimal formation

Our main focus is to analyse whether the differential collapse between the dust and gas can lead to dust enhancements in the embedded disk that can generate planetesimals. To determine the quantity of dust mass that is available for the production of planetesimals, we computed the threshold dust-to-gas ratio for each of the grid cells that are near the midplane of our embedded disks (within z < 0.2 AU). We computed the St in the gas model for grains of each size and the corresponding threshold dust-to-gas ratio in every grid cell near the midplane in each model. Next, we determined every grid cell in which the dust-to-gas ratio exceeded the threshold ratio in Eq. (28) for z < 0.2 AU, and selected all dust parcels in the corresponding grid cells. Finally, we computed the total mass of all of the dust parcels contributing to the grid cell dust density to estimate the amount of mass that might be transformed into planetesimals through the streaming instability.

Figure 12 presents the total amount of dust mass in regions of the disk in which the clumping threshold ratio is reached. The Ref-med model shows the least efficient clumping, with the lowest amount of dust mass available for the streaming instability over nearly the whole simulation. While the gas disk is largest in the reference model, it shows inefficient dust settling compared to the standard models; most of the dust mass is kept off of the midplane. As a result, the dust-to-gas ratio does not rise sufficiently high throughout the disk in order drive the streaming instability.

The standard models show different available dust masses, dependent on their dust size. During the class 0 phase, only the Std-lrg model produces regions in its embedded disk in which planetesimals can be formed. The figure clearly shows that this reservoir would either be quickly used up to form planetesimals, or lost to radial drift, which is what occurs in our simulation. The large grains can produce about 7 M_⊕ of planetesimals in the class 0 phase at most before the dust mass is exhausted. The Std-sml and Std-med models, however, can both only start to generate planetesimals during the class I phase because of the amount of time needed for dust to build up in the embedded disk. During the class I phase, the Std-med model can produce 22 M_⊕ at most before it too loses its dust reservoir to radial drift. Only during the last ~50 kyr does the dust in the Std-sml model begin to produce planetesimals, and could proceed to do so as the system transitions from class I to II (at t > t_acc). This model has 35 M_⊕ at most available for the generation of planetesimals near the end of the class I phase.

The amount of planetesimals produced through the streaming instability, between 7 and 35 M_⊕, is confined to a relatively small ribbon in the embedded disk (recall Fig. 9). It is therefore possible that this population of planetesimals could lead to the first collapse of a large protoplanetary embryo. This mass range should be viewed as the lowest possible planetesimal mass in embedded disks. At the high-mass end of the above range, the protoplanetary embryo is on the order of the pebble isolation mass and would thus produce a pressure bump in the embedded disk that could contribute to dust trapping and further planetesimal formation (as explored e.g. in Lenz et al. 2019; Eriksson et al. 2020). At the low-mass end, the embryo would be sufficiently massive to efficiently accrete pebbles via pebble accretion (e.g. as discussed in Bitsch et al. 2015), which would increase its mass. This growing planet would eventually reach the pebble isolation mass and similarly cause a pressure bump and dust trapping in the embedded disk. While this initial generation of planetesimal mass is small with respect to the total mass in the embedded disk (recall Fig. 7), it could therefore lead to further efficient planetesimal formation and to the production of a system of planets.

Fig. 12 Temporal evolution of the mass available for the streaming instability. The streaming instability is most efficient for the 34 µm grains as the Std-med model transitions into the class I phase of its evolution.

Fig. 13 Evolution of the dust mass, while removing the mass available for the streaming instability. In doing so, we can more closely match the class I disk masses found in the observations (grey regions, same as in Fig. 7). The Ref-med model finds almost no generation of planetesimals because dust mass in the disk midplane is very low (see Fig. 10).

4.2 Interpreting the dust mass evolution of young systems

Figure 13 reproduces Fig. 7 (with solid coloured curves) and removes the mass that was available to the streaming instability (shown in Fig. 12), which results in the dashed coloured curves and replicates the production of planetesimals from the available dust mass. We assumed that the streaming instability is efficient enough to use all of the available mass within a few orbits, so that we can remove the available streaming instability mass from the total dust mass at any given time. The assumption that dust from previous time steps is no longer available in the embedded disk is partly supported by our previous observation that the dust parcels in the disk at a given time are largely accreted recently from the envelope. We therefoer assume that dust in the embedded disk is removed either by radial drift or by planetesimal formation. Figure 13 clearly shows that the production of planetesimals has only a small impact on the total dust in the embedded disk in our model.

The streaming instability-subtracted final (at t = t_acc) masses range between ~300 and 400 M_⊕ for the Std-med, Std-sml, and Ref-med models. This mass range represents the top ~20% of the observed class I disks in Tychoniec et al. (2020), which could imply that ourmodel represents systems that are atthe higherend of mass in the observed population. However, our initial envelope mass of 1 M_⊙ (with an additional 1% of dust mass) is near to the peak in the typical clump mass function (Motte et al. 1998; Testi & Sargent 1998; Könyves et al. 2010), which likely represents the distribution of initial masses of star-forming clumps. Hence it is more likely that our initial setup represents a typical cloud.

Our final stellar mass is roughly 0.7 M_⊙, representing a star formation efficiency of about 70%. This is much higher than the typical 30% that is seen when the clump mass function is compared to the to the initial mass function. This suggests that we are missing some inefficiency in the collapse phase of star formation, which is likely related to our simple treatment of the outflow. Whether dust is entrained in the outflow remains an open question, but given the missing inefficiencies in the evolution of the gas, it is likely that our final embedded dust disk masses are also on the high-mass side.

4.2.1 Are class 0 or class I embedded disks that are more massive?

One of the primary unknowns in interpreting continuum emission data as they pertain to the amount of dust mass in young stellar systems is the opacity law, which governs radiative transfer in the cloud. Tychoniec et al. (2020) used the opacity for their ALMA emission (at 1.3 mm) of Ansdell et al. (2016) in order to directly compare their results to available class II dust mass data. Their value of 2.3 cM² g⁻¹ is on the higher end of the possible dust opacities computed by Draine (2006). Furthermore, Ribas et al. (2020) showed through spectral energy distribution modelling that the dust masses observed by ALMA could be as high as three times the masses inferred by Tychoniec et al. (2020). In addition, Zhu et al. (2019) have argued that dust scattering, as an additional source of opacity, will generally cause underestimates of the continuum emission dust masses seen by ALMA. The computed masses of Tychoniec et al. (2020) from the ALMA continuum should thus be considered the minimum dust masses of young systems in their survey, which could improve the comparison between our models and observational data.

The dust masses inferred with VLA data use a different dust opacity model and observe at longer wavelengths. Their light is thus more likely to be optically thin and may more accurately represent the real spread in dust masses in class O/I systems. The observed range of class 0 dust masses span a similar mass range as we found for our modelled embedded disks when in their class 0 phase. During the class I phase, the modelled disks overpredict the dust mass compared to observations. The data inferred by ALMA and the VLA both report the same property of class O/I disks: class 0 disks tend to be more dust rich than class I disks.

Our model tends to predict that the embedded disk is more dust rich in its class I phase than in its class 0 phase because it takes time for the dust from the envelope to fall into the disk. Furthermore, dust that reaches the disk early can already radially drift into the host star. This appears to be a main component of dust evolution in embedded disks because successive time steps in Figs. 5 and 6 typically show that the disk is mainly populated by dust that has only recently arrived from farther out in the envelope; dust that had arrived from earlier time steps is no longer found in the disk. This appears to be independent of the initial angular momentum in the cloud because both the reference and standard models share this property. Because our class I disks have a higher mass than class 0, they tend to overestimate the dust mass in observed class I systems, except for the Std-lrg model, where the large grains quickly drift radially in the embedded disk. Unless dust growth is very fast in the first few 10s of kyr into the collapse, it seems unlikely that the observed class O/I dusty disk masses are dominated by millimeter-dust grains because of their weak retention in the embedded disk.

4.2.2 Impact of grain growth

As argued in the introduction, there is observational evidence to support a wide range of possible grain sizes in young stellar systems. Our results suggest that these populations of large grains are unlikely to be part of the dust grain distribution in the initial cloud because the bulk of the dust mass larger grains is rapidly lost to radial drift during the class I phase. On the other hand, the smaller grains take between 0.3 and 0.45 t_acc (~75–113 kyr) to build up the observed median dust masses in the standard model. To account for the observational evidence of larger grains in protostellar systems, some amount of grain growth would need to occur during the collapse.

As previously mentioned, the dust at the upper edge of the disk is jammed, which could lead to some grain growth that would change the overall mass evolution of the disk. The typical growth timescale through coagulation is (Birnstiel et al. 2012) (30)

To estimate the possible grain growth within the jamming, we considered the 0.5 t_acc snapshot of the Std-med model (bottom left panel of Fig. 8). Near 10 AU, the thickness of the jammed region is approximately 1 AU, and the dust-to-gas ratio is within a factor of a few higher than the typical ISM ratio. The amount of time that a dust parcel moving at 100 cm s⁻¹ would take to cross the jammed region at 10 AU is almost 5 kyr. If the dust-to-gas ratio at 100 AU near the disk edge is equal to the typical ISM ratio, then the growth timescale of the dust grains is 0.5 kyr at most. Because the disk is marginally more dust rich than in the ISM (slightly more orange in Fig. 8), the growth timescale is even shorter. Because the growth timescale is shorter than the crossing timescale, it is possible that some grain growth could occur in the jammed region.

These estimates are sensitive to the local dust-to-gas ratio and to the kinematics of the grain parcels moving through the upper regions of the disk. This simple experiment shows, however, that grain growth near the upper edge of the embedded disk could have important implications for the overall evolution of the dust. In particular, grain growth would tend to produce larger grains, which would contribute to the dust grain population simulated in the Std-lrg model. These larger grains can settle more quickly in the midplane, even in the presence of vertical turbulence and diffusion, and would provide a higher and more prolonged mass flux to the midplane than was found in the Std-lrg model. Because the streaming instability is more efficient for larger grains, we would then expect that more planetesimal mass is produced. We leave the study of this possible effect to future work. A simulation that self-consistency computes the dust grain size evolution in particular must be included in the collapse.

Even with this possible region of the efficient grain growth, the large dust grains in the envelope are not fully explained. To do this, grain growth in the embedded disk would need a second mechanism to cycle the larger grains outward into the envelope where they could be observed. Such a mechanism would also help to sustain the dust mass in protostellar systems, and is discussed below.

4.2.3 Sustaining the dust mass in class O/I disks

Models have suggested that in the context of our own Solar System that dust grains are cycled from the inner disk outward to larger radii. Such a mechanism would counterbalance radial drift and could help to accelerate the dust mass built up during the class 0 phase. These models centre on calcium-aluminium inclusion (CAI) formation and broadly argue that CAIs are formed in the high-temperature region of the young solar nebula before they are mixed outward by an X-wind (Shu et al. 1996, 1997, 2001), winds driven by magnetohydrodynamics and turbulence (Cuzzi & Hogan 2003; Cuzzi et al. 2003), or hydrodynamic turbulence (Desch et al. 2010). While these models generally focus on an older protostellar system than the ~50 kyr class 0 disk we discuss here, and often include a young Jupiter to set the distribution of CAIs (McKeegan et al. 1998; Desch et al. 2018; van Kooten et al. 2021), the physical mechanisms responsible for mixing CAIs throughout the young Solar System could sufficiently work against radial drift to enhance the dust mass in class 0 disks. The role of winds and the diffusive effect of turbulence on the radial evolution of the dust is not included in our model, but offers an interesting line of further study.

In addition to helping us understand the distribution of CAIs in the protosolar nebula, the mixing of dust grains from the inner part of the protostellar envelope could also help to explain the existence of larger grains that are observed at large scales of young systems. While larger grains have been inferred observationally (Galametz et al. 2019), the physical mechanism to explain their existence in the outer envelope remains unknown. If the cloud and its dust simply began collapsing from the surrounding ISM, the dust would be expected to be mainly monomer size. These large grains might be the result of early accretion of small grains onto a class 0 embedded disk or of limited grain growth within the disk, followed by outward recycling back into the envelope. Solving this part of the problem may require large-scale hydrodynamic simulations to demonstrate whether this recycling is possible. This is beyond the scope of this work.

If an additional mechanism could be included that slowed the global impact of radial drift, then our overestimation of the class I disk mass would be exaggerated. However, preserving dust mass in the disk, coupled with a certain amount of grain growth, will lead to a higher available dust mass for the streaming instability. This means that while the preserved mass in class I disks may be higher, the production rate of planetesimals would also be higher unless the dust mass is kept away from the midplane. This seems unlikely because as grains grow, they will also tend to sink to the midplane, ready to be converted into planetesimals. If the physical mechanisms that are left out of our simple experiment tend to preserve dust mass, then it is possible that the generation of the planetesimals can be efficient enough to reverse the tendency we see here: that class I embedded disks contain a higher dust mass than class 0 disks.

5 Conclusion

We modelled the collapse of a dusty protostellar envelope to study the accumulation of the dusty embedded disk. The variations in the collapse trajectory of dust grains, caused by hydrodynamic drag and the immediate collapse of dust, can generate regions in the embedded disk in which the local dust-to-gas ratio exceeds the typical 1% of the ISM. These regions of enhanced dust mass become susceptible to the streaming instability, which would lead to the generation of planetesimals very early in protostellar evolution. This process and the subsequent beginning of planet formation is most efficient in the class I phase for grain sizes ≤34 µm, but even earlier in the class 0 phase for larger grains. This generation of planetesimals can help to understand observations made in surveys of class O/I systems that show that class I disks tend to be less massive than the younger class 0 objects.

We have developed a simple model to account for the gas drag on a spherical cloud of dust embedded in the gas of a collapsing core. Our model uses the semi-analytic gas model of Visser et al. (2009) and computes the trajectory of 100000 dust parcels as they collapse. This model includes the effect of gas drag on the dust particles as they fall towards the host star, which changes their collapse trajectory with respect to the gas particles. Overall, the lack of gas pressure support in the dust grains causes them to fall much faster than the gas, forming dust enhancements in the midplane of the embedded disk. These enhancements survive for an extended period of time against radial drift in the midplane of the embedded disk because of ongoing mass flux from the collapsing envelope. This balance, however, tends to lean towards radial drift for the larger grains, which are eventually lost to the host star on a timescale related to the grain size. Through our simple experiment, we have found the results we list below.

1. All grain sizes form regions of their embedded disk in which planetesimals can be formed. The larger grains do this earlier.

2. Our models have between 7 and 35 M_⊕ available for the generation of planetesimals from the streaming instability during the class O/I phases.

3. The submicron-sized grains reach a total available mass of about 35 M_⊕ near the transition between class I and II.

4. The range of formed planetesimal masses is on the order of the shift in median masses between class 0 and class I objects as inferred by continuum emission surveys.

5. Our total dust masses, with the mass available for the streaming instability removed, are higher than the median dust masses in observed class I disks. They represent the top ~20% of observed disk masses.

The differences between our model and observations is likely linked to missing physics such as the outward radial mixing of dust, which is required to explain the distribution of CAIs in the Solar System. Furthermore, we have neglected turbulent diffusion, which can act to suppress mass clumping in the disk midplane (through the suppression of settling) and thus the streaming instability, but can also further counteract radial drift to slow the loss of dust mass. It is clear that while this semianalytic approach to the collapse of a dusty pre-stellar cloud is illustrative, more rigorous numerical methods are needed to fully explore the mass evolution of young disks.

The quantity of planetesimal mass that is generated in our models could represent the first step in the planet formation process. Specifically, because the region of the embedded disk in which planetesimals can be produced is restricted to roughly 10 AU, these planetesimal could hypothetically coagulate into the first planetary core early on in the class O/I phase. This young planetary core (with a core mass of ~20 M_⊕) would then produce the first pressure bump in the embedded disk gas, causing dust trapping, and the earliest substructures in embedded disks. These dust traps would also promote further planetesimal growth as they are a prime location for high dust-to-gas ratios, which can further generate new planet cores. Thus, our model provides a minimum mass in planetesimals in this phase of star formation.

Under this scenario, the formation of our own Solar System would proceed with the formation of the Jupiter core first, producing a pressure maximum where the planetesimals that form the Saturn core are collected. The Saturn core would then grow and produce its own pressure maxima, in which the next generation of planetesimals would form. This next generation of planetesimals would proceed to form the ice giants. The delayed formation of Saturn with respect to Jupiter would be in line with the assumptions of the Grand Tack model (Walsh et al. 2011), while the further delay in the formation of Uranus and Neptune would explain their lower masses in the pebble accretion paradigm (see e.g. Helled et al. 2020; Brouwers et al. 2021).

It should be clear from our work that the generation of planetesimals through a process such as the streaming instability can potentially occur during the collapse phase of protostellar formation. This therefore pushes the start time of the planet formation process back, possibly as far as into the class 0 phase. Such an early start suggests that the first planets to form around a young protostar may emerge nearly fully formed by the time their host star transitions from the class I to the class II phase, which is colloquially called the protoplanetary disk phase.

Acknowledgements

Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA). BT acknowledges support from the research programme Dutch Astrochemistry Network II with project number 614.001.751, which is financed by the Dutch Research Council (NWO). GR acknowledges support from the NWO (program number 016.Veni.192.233) and from an STFC Ernest Rutherford Fellowship (grant number ST/T003855/1).

Appendix A Collapse of the Std-lrg and Std-sml models

For brevity in the main text, we plot the collapse trajectories of the Std-lrg and Std-sml here. The colours of the points have the same definition as in figures 4-6. The grains in the Std-lrg model (figure A.1) evolve much faster than the small grains in the Std-sml model (figure A.2), as can be seen by comparing the colour of the dust parcels at any given time.

Fig. A.1 Same as in Figure 5, but for the Std-lrg model. The larger grains are less strongly coupled to the gas and collapse much more quickly. Once in the disk, their radial drift is also faster than in the Std-med model, and they rapidly deplete into the host star. Unlike in the Std-med model, grains originating from different initial radii ranges therefore do not coexist in the embedded disk with dust parcels of different colours.

Fig. A.2 Same as in Figure 5, but forthe Std-smlmodel. These smaller grains are more dynamically coupled to the gas and follow similar trajectories. The dust disk that they form is also more similar to the gas as radial drift and settling are both less efficient for this population of dust grains.

Appendix B Gas-to-dust ratio in Std-lrg and Std-sml models

Again we plot the results of the Std-lrg and Std-sml here for brevity. As discussed, it is clear from Figures B.2 and B.4 that the dust in the Std-lrg model evolves rapidly, but still produces regions in which dust exceeds a dust-to-gas ratio of ~0.5. This is enough to form planetesimals. The dust in the Std-sml dust model takes longer to build up in the disk midplane, but due to its slow radial drift in the disk, the dust-to-gas ratio can rapidly build up at later times, producing an exceedingly high dust-to-gas ratio near the end of the simulation.

Fig. B.1 Same as Figure 8, but for the Std-lrg model. These larger grains evolve faster than the smaller grains. In particular, in the high density of the embedded disk, they radially drift very efficiently into the host star.

Fig. B.2 Midplane dust-to-gas ratios for the Std-lrg collapse model at the same time steps as shown in Figure B.1. The lines and hashed regions have the same definition as in Fig. 9.

Fig. B.3 Same as Figure 8, but for the Std-sml model. Because of the small dust grains used in this model, the collapse and evolution in the disk are more strongly coupled to the gas. They therefore tend to evolve in a similar was as the gas, creating mainly ISM dust-to-gas ratios in the midplane of the disk. In the thin regions around the surface of the disk, the sudden pressure change captures the dust, similarly to typical dust traps in class II disks.

Fig. B.4 Midplane dust-to-gas ratios for the Std-lrg collapse model at the same time steps as shown in Figure B.3. The lines and hashed regions have the same definition as in Fig. 9.

References

All Tables

All Figures

Fig. 1 Schematic view of the embedded disk within a collapsing envelope. Unlike in traditional protoplanetary disks, i.e. Class II objects, embedded disks (Class O/I) are continually fed with new gas and dust from the envelope (infall). The gas in the embedded disk is subject to turbulent flow and evolves viscously. The dust, on the other hand, tends to sink to the midplane of the embedded disk and radially drifts towards the host star. Because the gas, but not the dust, is held aloft by gas pressure, the gas and dust arriving at the embedded disk at any given time originate from different initial locations in the envelope; dust arrives from larger initial radii than gas. The grey scale broadly shows the total volume mass density of the gas and dust in the protostellar system.

In the text

Fig. 2 Single dust parcel trajectory starting from the same location in the envelope for all four models. Broadly speaking, time flows from the top right of the figure to the bottom left. The coloured point denotes the location of the parcel in each model at a collapse time of 120 kyr. The size of the dust particles in each collapse model denotes how quickly the dust collapses, with larger grains being farther along in their collapse than smaller grains. The different rotation rates of the clouds produces a small change in the envelope collapse, and the embedded disk grows more quickly, which causes a large change in the trajectory when the dust particle encounters the disk.

In the text

Fig. 3 Evolution of the specific angular momentum (total angular momentum per mass) of the dust parcels for each of the models. The dashed line highlights the time of 120 kyr that was featured in Fig. 2. Angular momentum is strictly conserved for the first few 10 s of kyr until the dust has fallen into a higher-density region of the cloud, where interactions between the gas and dust begin to efficiently drain angular momentum. When the dust reaches the disk, it enters a region with much higher angular momentum and is quickly spun up by the gas.

In the text

Fig. 4 Collapse and disk build-up of the gas parcels in the standard model. Each set of three panels shows the collapse at 0.05,0.25,0.5, and 0.75 tacc over three length scales. The left panel shows the largest scale, the envelope scale, and the right panel shows the length scale relevant to the disk. The colour of each point denotes the initial radius range from which the parcel originates. Dark brown denotes an initial radius smaller than 1000 AU, orange a radius between 1000 and 2000 AU, and so on. The embedded disk builds up from gas that successively falls from higher heights as well as from the viscous spreading of the disk itself. For reference, the black line shows an estimate of the location of the rarefaction wave defined by x = 1.

In the text

Fig. 5 Collapse of dust parcels in the Std-med model. Top left set of three panels: very early time step in the collapse. The colour of each point is the same as in Fig. 4. Each set of three panels shows the same time step at three different lengths scales. The right panel shows the length scales relevant for the embedded disk. The second, third, and forth set represent the simulation at 0.25, 0.5, and 0.75 tacc respectively. As the disk builds, some of the dust drifts radially inward, some is advected outward with the growing gas disk, and missing dust is replenished by newly infalling dust. For comparison, the black line shows an estimate of the location of the gas rarefaction wave defined by x = 1.

In the text

Fig. 6 Collapse of dust parcels in the Ref-med model. The colour and time steps are the same as in Figs. 4 and 5, but the x-axis scale changes in the right panels. While the overall collapse is over a similar timescale as in the Std-med model, the disk build-up and subsequent evolution in the disk is different. The gas disk is larger in the reference model, and dust grains fall predominately at the outward edge of the disk on average. Viscous spreading is still present in the embedded disk, but most of the dust evolution is inward via radial drift. At late times, the dust ring is therefore extended, and most of the dust is inward of the ring in the midplane or at the upper surface of the disk. For comparison, the black line shows an estimate of the location of the gas rarefaction wave defined by x = 1.

In the text

Fig. 7 Time evolution of the total dust mass inside of the embedded disk in each model over one tacc (~252 kyr). The grey bands show the median disk masses for ALMA-detected (right hatch) and VLA-detected (left hatch) protostars reported by Tychoniec et al. (2020) along with the confidence interval of the median. The ALMA-inferred masses are likely lower limits because the dust is optically thick, while the VLA-inferred masses are more likely to be optically thin and use realistic opacities assuming the existence of large grains. The vertical dashed lines show the time at which the embedded disk first appears (classifying the system as class 0) and when the disk gas mass exceeds the remaining gas mass in the envelope (classifying the system as class I). Class 0 and class I objects were both included by Tychoniec et al. (2020), and we have separated them to compute individual medians and confidence intervals. For comparison, we include the mass evolution of the disk gas in the standard (solid black) and reference (dashed black) models, scaled by the ISM dust-to-gas ratio.

In the text

Fig. 8 Dust-to-gas mass ratio in the Std-med collapse model. The colour scale shows thes log of the local dust-to-gas ratio, hence blue denote dust-to-gas ratios below the ISM, and orange denotes dust-to-gas ratios above the ISM. The dashed lined shows the extent of the outflow in the gas model, and the solid line shows the surface of the embedded disk. We assumed that planetesimal formation occurs largely along the midplane, which is highlighted by the grey arrow. We show the midplane dust-to-gas ratio below.

In the text

Fig. 9 Midplane dust-to-gas ratios for the Std-med collapse model at the same time steps as shown in Fig. 8. At the earliest times (solid blue curve), the dust disk has yet to form and does not appear in this plot. The black line shows the ISM dust-to-gas ratio, which represents the initial conditions in the envelope, and the coloured dashed lines show the streaming instability threshold at each time step. The hashed region is an estimate of the streaming instability threshold based on Youdin & Goodman (2005), and the dashed lines show the threshold based on Li & Youdin (2021). The threshold is met along the midplane of the embedded disk only at the last time step.

In the text

Fig. 10 Same as Fig. 8, but for the reference model Ref-med. Settling is less efficient than in the standard model, and most of the dust grains remain near the upper surface of the disk. Regions just outward of the disk surface along the midplane show sufficiently high dust-to-gas ratios. These regions, however, do not survive to the last snapshot shown here. Planetesimal formation predominantly occurs along the midplane, highlighted by the grey arrow in the forth panel, and shown below.

In the text

	Fig. 11 Midplane dust-to-gas ratios for the Ref-med collapse model at the same time steps as shown in Fig. 10. The lines and hashed regions have the same definition as in Fig. 9. The dust-to-gas ratio never exceeds the threshold for the onset of the streaming instability at any of the shown time steps.
In the text

	Fig. 12 Temporal evolution of the mass available for the streaming instability. The streaming instability is most efficient for the 34 µm grains as the Std-med model transitions into the class I phase of its evolution.
In the text

	Fig. 13 Evolution of the dust mass, while removing the mass available for the streaming instability. In doing so, we can more closely match the class I disk masses found in the observations (grey regions, same as in Fig. 7). The Ref-med model finds almost no generation of planetesimals because dust mass in the disk midplane is very low (see Fig. 10).
In the text

Fig. A.1 Same as in Figure 5, but for the Std-lrg model. The larger grains are less strongly coupled to the gas and collapse much more quickly. Once in the disk, their radial drift is also faster than in the Std-med model, and they rapidly deplete into the host star. Unlike in the Std-med model, grains originating from different initial radii ranges therefore do not coexist in the embedded disk with dust parcels of different colours.

In the text

	Fig. A.2 Same as in Figure 5, but forthe Std-smlmodel. These smaller grains are more dynamically coupled to the gas and follow similar trajectories. The dust disk that they form is also more similar to the gas as radial drift and settling are both less efficient for this population of dust grains.
In the text

	Fig. B.1 Same as Figure 8, but for the Std-lrg model. These larger grains evolve faster than the smaller grains. In particular, in the high density of the embedded disk, they radially drift very efficiently into the host star.
In the text

	Fig. B.2 Midplane dust-to-gas ratios for the Std-lrg collapse model at the same time steps as shown in Figure B.1. The lines and hashed regions have the same definition as in Fig. 9.
In the text

Fig. B.3 Same as Figure 8, but for the Std-sml model. Because of the small dust grains used in this model, the collapse and evolution in the disk are more strongly coupled to the gas. They therefore tend to evolve in a similar was as the gas, creating mainly ISM dust-to-gas ratios in the midplane of the disk. In the thin regions around the surface of the disk, the sudden pressure change captures the dust, similarly to typical dust traps in class II disks.

In the text

	Fig. B.4 Midplane dust-to-gas ratios for the Std-lrg collapse model at the same time steps as shown in Figure B.3. The lines and hashed regions have the same definition as in Fig. 9.
In the text