© The Authors 2023

1 Introduction

Understanding planetesimal formation within protoplanetary disks is an important first step towards understanding planet formation. The distribution of these first solid bodies drives the locations where planetary embryos, which eventually form fully-fledged planets, grow (e.g. Chambers 2001; Walsh et al. 2011).

Observations of protoplanetary dust disks reveal two distinct properties: they are large and long-lasting. Their sizes range from 10 to 500 au with typical sizes of ~30 au (Tripathi et al. 2017; Andrews et al. 2018; Hendler et al. 2020), and they have lifetimes of millions of years (e.g. Barenfeld et al. 2017; Ruíz-Rodríguez et al. 2018). Because the disk formation occurs on much shorter timescales (of the order of 100 thousand years), dust is not continuously supplied to the system. It, therefore, needs to be preserved at large heliocentric distances for millions of years after disk formation.

The Solar System provides a set of additional constraints on the properties and evolution of the protosolar disk. However, it is unknown a priori whether they are common to most protoplanetary disks or specific to our own.

The existence and the properties of comets suggest that the protosolar disk was typical in terms of radial extension and lifetime. In fact, comets are thought to have formed at distances between 20 and 40 au from the Sun (Nesvorný et al. 2017; Nesvorný 2018), and cold classical Kuiper belt objects are thought to have formed in situ up to a distance of 45 au (Nesvorný et al. 2022). Additionally, comets likely formed late (e.g. Nakashima et al. 2015; Nimmo et al. 2018; Neumann et al. 2018), several million years after the formation of the first solids, the so-called calcium–aluminium-rich inclusions (CAIs; Amelin et al. 2010; Connelly et al. 2012). A late formation is needed to avoid any significant radiogenic heating, which would result in the loss of highly volatile ices such as CO₂ and CO (e.g. Eberhardt et al. 1987; Morse et al. 2015; Gasc et al. 2017). The presence of these highly volatile species in very large comets (~100 km) such as Hale-Bopp or Bernardinelli–Bernstein (Capria et al. 2000; Kelley et al. 2022) confirms that comets remained cold not because of their small sizes but rather because they formed late, at a time when most short-lived radioactive elements (e.g.²⁶ Al) had already decayed. Also, radioactive heating would have increased the bulk density of large objects to a degree inconsistent with the low density of icy bodies such as Trojans and Kuiper-belt objects (between 300 and 1500 km m⁻³; Preusker et al. 2017; Groussin et al. 2019; Berthier et al. 2020; Spencer et al. 2020), further supporting late formation.

We have additional evidence for a long-lasting protosolar disk. The meteoritic record contains samples from both differentiated and un-differentiated parent bodies. The latter formed significantly later – up to 5 million yr after CAI formation (Nimmo et al. 2018). Therefore, ample evidence suggests that our Solar System formed from an extended and long-lived protoplanetary disk. Because we focus in this work on the first generation of planetesimals, and the problem of long-lasting disks is an issue in itself, our first requirement for a good model of the Solar System disk is a large size.

Focusing on the first generation of planetesimals, the differentiated parent bodies of iron meteorites, we find that they can be divided into two isotopically distinct groups akin to carbonaceous chondrites (CCs) and non-carbonaceous chondrites (NCs; Warren 2011; Kruijer et al. 2017). Thus, they are usually referred to as CC- and NC-iron meteorites, respectively. Both groups of iron meteorites formed essentially simultaneously in the disk (Spitzer et al. 2021). Because they formed simultaneously, they must have formed at distinctly different locations in the disk, which may have different disk compositions. Therefore, our second requirement for a good model of the Solar System disk is that it produces planetesimals at two distinct locations.

Finally, the oldest Solar System solids, CAIs, are thought to have formed as high-temperature condensates very close (a few tenths of an au) to the proto-Sun (Scott & Krot 2003). The age of CAIs sets what is usually considered time zero for Solar System formation (see the review by Chaussidon & Liu 2015). Their age is 4567.30 ± 0.16 million yr according to Pb-Pb dating (Jacobsen et al. 2008; Connelly et al. 2012; Bouvier & Wadhwa 2010). Recent work argues for a revised age for CAIs of 4568.7 Myr (Piralla et al. 2023; Desch et al. 2023). The duration of CAI formation appears to be very short, from ~100 kyr (Connelly et al. 2012) to as low as ~10 kyr (Jacobsen et al. 2008). Importantly, the abundance of CAIs is significantly higher in CCs than NCs (Scott & Krot 2003), the latter of which are thought to have formed closer to the Sun than the former (Warren 2011). Furthermore, CAIs have even been found in comets (Brownlee et al. 2006; Zolensky et al. 2006), which descend from planetesimals formed the farthest away from the Sun. Therefore, even though CAIs were formed close to the Sun, the planetesimals that formed the farthest away are more enriched with them. This implies that these high-temperature condensates were efficiently transported to the outer disk, which therefore became enriched with CAIs, while the inner disk remained depleted in CAIs. The fact that the isotopic compositions of differentiated (early) and undifferentiated (late) planetesimals overlap within the CC and NC reservoirs, respectively (Kruijer et al. 2017), indicates that this division of a CAI-rich outer and CAI-depleted inner disk was already present at the time when the parent bodies of the iron meteorites formed. It has been proposed that CAIs were transported ballistically to the outer disk via magnetised winds (Shu et al. 2001). But modern simulations reveal that only particles much smaller than observed CAIs can be efficiently transported this way (Rodenkirch & Dullemond 2022). Thus, the radial transport of CAIs during the outward spreading of the disk (Jacquet et al. 2011; Pignatale et al. 2018) remains the best option.

In summary, for our Solar System, a disk formation and evolution scenario must satisfy at least the following three criteria: (1) it must develop an extended disk of gas and dust (up to 45 au for the dust); (2) in at least two distinct locations in the disk, the dust/gas ratio must be able to increase sufficiently to produce planetesimals and explain the early formation of NC- and CC-iron meteorite parent bodies; and (3) particles that condensed at high temperatures (i.e. CAIs) must be able to reach large heliocentric distances, that is to say, be transported from the star’s proximity to large distances.

In this work we try to build such a scenario. In Sect. 2 we describe the key processes in the formation of the disk, the evolution of its gas and dust components, and planetesimal formation. Then we describe the disk model we use (Sect. 3) before discussing the model setup (Sect. 4). In particular, we describe how four assumptions affect the ability to match (or not) the Solar System constraints. These are (1) the centrifugal radius, R_C; (2) the initial viscosity of the disk, α₀; (3) the infall timescale of material onto the disk, T_infall; and (4) the effect of a temperature-dependent fragmentation threshold for icy particles. Our results are presented in Sect. 5.

We show that an initial rapid expansion - signifying an inflationary disk phase – can result in large dust disks, forming planetesimals at two locations in the disk and transporting CAIs to the outer disk. We also show that disks forming from clouds with large angular momenta, which readily solves the problem of dust-disk sizes by directly delivering material at large distances, are unable to form planetesimals at two distinct locations and do not allow the transport of CAIs into the outer disk.

2 Key processes in disk evolution and planetesimal formation

We start by discussing key processes in the formation and evolution of the disk and planetesimal accretion, focusing on the unknowns we parametrise and test in our models.

2.1 Accretion of material into a protoplanetary disk

Whether protoplanetary disks are ‘born’ big (i.e. form from the outside in) or ‘grow up’ to be big (i.e. grow from the inside out) depends on the angular momentum of the infalling material. Thus, the angular momentum of the pre-stellar cloud determines where material falls into the disk. The larger the angular momentum of the material, the larger the distance at which it falls into the disk. The radius in the disk where the angular momentum of the infalling material is equal to the angular momentum of the Keplerian disk is called the centrifugal radius, R_C. If, for example, the pre-stellar cloud has a constant angular speed throughout, then shells of material closer to the centre collapse first and, having a small specific angular momentum, will fall very close to the protostar. More distant shells fall into the disk later and, having larger specific angular momentums, fall farther away from the star. Therefore, R_C increases with time for a prestellar cloud with a constant angular frequency. Depending on the pre-stellar cloud, the centrifugal radius can be as large as 100 au (e.g. Shu 1977; Hueso & Guillot 2005; Pignatale et al. 2018).

However, it is also possible the material falls continuously close to the star because of magnetic braking, which removes a significant amount of the angular momentum of the infalling material (Lee et al. 2021, magnetically braked material flows along the disk surface towards the protostar, sketched in Fig. 19). The formation of such small disks is observed in some magneto-hydrodynamic (MHD) simulations of the gravitational collapse of pre-stellar clouds (e.g. Machida & Basu 2019; Vaytet et al. 2018; Machida & Matsumoto 2011). These disks can then spread radially due to viscous evolution.

Current cloud collapse simulations do not yet provide a firm prescription on how a disk forms and where it collects the material falling from the molecular cloud. Thus, in the following we test different idealised recipes to identify which best fits the constraints enumerated in the introduction.

Observations suggest that the timescale of accretion of material into the disk is of the order of 10⁵ yr, with large uncertainties (Larson 1969; Vaytet et al. 2018; Wurster et al. 2021, 2022), so the infall timescale can be considered a free parameter within an order of magnitude. Late accretion through streamers is sometimes observed (Tobin et al. 2010; Yen et al. 2019; Pineda et al. 2020) but, given the stochastic nature of this process, we did not include it in our investigations.

The viscosity plays a key role in the evolution of the disk and its spreading away from R_C. There is a big discussion in the literature on the actual viscosity of protoplanetary disks, but it concerns isolated accretion disks. As long as the disk is accreting material from the molecular cloud, it is expected to be subject to strong Raynold stresses, which act as an effective viscosity (Kuznetsova et al. 2022). Thus, it seems legitimate to assume that a disk that is still accreting mass has a viscosity proportional to the mass infall rate, but the proportionality factor is poorly constrained, and therefore we consider different values in our study.

2.2 Motion of dust particles within the disk

For disks forming with a small R_C where, for example, the material never falls outside of 10 au, dust particles must be efficiently transported from the vicinity to distances far away from the star in order to build the large observed dust disks. In such cases, the disk (dust and gas) forms from the inside out. The outward motion of the dust is induced through the radial aerodynamic drag of the radially expanding gas (e.g. Yang & Ciesla 2012). Gas within R_C has a negative radial velocity (towards the star), but the gas close to and beyond R_C viscously spreads outwards. Eventually, the entire gas disk becomes an accretion disk with a negative radial velocity throughout the disk.

The radial motion of the dust depends on its size. The important parameter for dust dynamics is not the particle size but its Stokes number, defined as

$$\text{S}t=\frac{\pi a\rho d}{4{\text{Σ}}_g},$$(1)

where a is the diameter of the dust particle, ρ_d is the particle solid density, and Σ_ɡ is the gas surface density. The radial dust velocity, ${\upsilon }_r^d$, can then be written as

$${\upsilon }_r^d=\frac{2\text{S}t}{1+\text{S}{t}^2}{\upsilon }_t^g+\frac{1}{1+\text{S}{t}^2}{\upsilon }_r^g,$$(2)

where ${\upsilon }_t^g$ and ${\upsilon }_r^g$ are the tangential and radial velocities of the gas relative to a circular Keplerian orbit, respectively. When there is no dust feedback onto the gas, ${\upsilon }_t^g=\eta {\upsilon }_K$ is the difference between the azimuthal gas speed and the Keplerian speed due to the partial pressure support of the gas. The radial velocity of the gas is due to viscosity.

For small dust, when S_t ≪ 1, the radial dust speed is dominated by the radial gas speed (${\upsilon }_r^d\propto {\upsilon }_r^g$, Eq. (2)). Thus, when the dust is small, it initially expands outwards from RC with the gas. Once the dust has grown sufficiently (i.e. St ~ 1), the tangential speed of the gas can become the dominant factor in Eq. (2). Because the gas is sub-Keplerian ${\upsilon }_t^g<0$, the radial dust speed can also become negative once the dust has grown large enough, even if the gas is still in radial expansion. This reflects the fact that dust particles that are large enough feel the headwind of the gas – the dust is moving at Keplerian while the gas is at sub-Keplerian speed. Thus, while the gas can further expand outwards viscously, large dust particles will begin to drift back towards the star.

2.3 Dust growth

Particles grow on a timescale 1/ZΩ, where Z = Σ_ɡ/Σ_ɡ is the local column integrated dust-to-gas ratio, but their growth is limited by the so-called fragmentation barrier (Drążkowska & Alibert 2017). When particles grow, they start to partially decouple from the gas. The turbulence in the disk and the radial drift of particles in the disk then enhance the relative speeds of the dust particles; when this speed becomes greater than the fragmentation velocity, υ_frag, dust particles cannot coagulate further and instead break upon collisions. The largest Stokes number that particles can acquire by coagulation is estimated to be (Drążkowska & Alibert 2017; Birnstiel et al. 2012) the minimum between

$$\text{S}{t}_{\text{frag}}=\frac{0.37{\upsilon }_{\text{frag}}^2}{3\text{Sc\hspace{0.17em}}\alpha c_s^2}$$(3)

and

$$\text{S}{t}_{\text{ddf}}=\frac{0.37{\upsilon }_{\text{frag}}}{2\left|\eta {\upsilon }_K\right|},$$(4)

where a is the gas viscosity parameter, following the assumption that the viscosity $v=\alpha c_s^2/\text{Ω}$ (Shakura & Sunyaev 1973), Sc is the Schmidt number relating viscous angular momentum transfer to turbulent diffusion, and c_s is the local sound speed. Equation (3) comes from the velocity dispersion due to turbulence in the disk, and Eq. (7) comes from the differential radial speed of particles of different Stokes numbers.

The fragmentation velocity υ_frag depends on the material properties. Following the results of laboratory experiments (e.g. Dominik & Tielens 1997; Wada et al. 2007; Blum & Wurm 2008; Teiser & Wurm 2009; Güttler et al. 2010), it is typically assumed that υ_frag = 100 cm s⁻¹ for refractory and silicate particles whereas υ_frag = 1000 cms⁻¹ for icy particles beyond the water snow line. Yet, recent laboratory experiments have shown that ice particles are only ‘sticky’ close to the sublimation temperature and more brittle when the ice is cold (e.g. Musiolik & Wurm 2019). Therefore, we explore an additional fragmentation threshold prescription for icy particles, which is temperature dependent. Similarly, it may be possible that silicate particles become stickier when their temperature is close to sublimation (Pillich et al. 2021) but, awaiting experimental confirmations, we do not yet consider this possibility in our model.

2.4 Planetesimal formation

The currently favoured mechanism for planetesimal formation is through the streaming instability (SI; Youdin & Goodman 2005; Johansen et al. 2007, 2014; Wahlberg Jansson & Johansen 2014, 2017; Simon et al. 2017; Yang et al. 2017; Abod et al. 2019) and subsequent gravitational collapse to form large – the preferred size of 100 km – planetesimals (e.g. Simon et al. 2016; Schäfer et al. 2017; Klahr & Schreiber 2020; Polak & Klahr 2023). The SI is triggered once sufficient dust collects within a certain region of the disk and causes the local dust-to-gas ratio to reach some threshold value (e.g. 0.5; Gole et al. 2020). At that point, clouds of dust particles collapse under their own gravity to form plan-etesimals (e.g. Klahr & Schreiber 2020; Nesvorný et al. 2021; Polak & Klahr 2023).

Previous models that explore the formation of planetesimals within a disk focused on static disks, that is, snapshots of a given disk phase. Such models have been successful in showing that planetesimal formation is particularly favoured in the vicinity of sublimation lines, in particular, the water snow line (e.g. Saito & Sirono 2011; Ida & Guillot 2016; Drążkowska & Alibert 2017; Hyodo et al. 2019, 2021). More recently, these static models were extended to include the temporal evolution of the gas and dust disks and confirm that planetesimal formation at the snow line remains the dominant location for forming a first generation of planetesimals (Drążkowska & Alibert 2017; Drążkowska & Dullemond 2018; Charnoz et al. 2019; Morbidelli et al. 2022). Such evolving disk models capture the expansion phase of the disk and therefore do not rely on a prescribed disk profile, for example the surface density of gas and dust. The addition of the silicate condensation line, in conjunction with a small centrifugal radius, was shown by Morbidelli et al. (2022) to result in planetesimals forming at the silicate line in addition to those forming at the snow line.

Yet, these newer, explicitly time-dependent inside-out formation models are problematic in that they cannot satisfy at least two of our requirements. These disks typically do not result in extended disks (requirement 1) and, by extension, will also struggle to bring CAIs to the outer disk (requirement 3). This shows that a more in-depth investigation is needed, which motivates the present paper.

The reason why the published models fail on requirements 1 and 3 is that the resulting dust disk sizes are merely slightly larger than the location of the water snow line (~ 5 au). This is because particles beyond the snow line rapidly grow and drift back towards the protostar on much shorter timescales due to aerodynamic drag in the tangential direction (e.g. Takeuchi & Lin 2002, 2005).

Thus, the underlying problem is one of the particle sizes and their associated dynamical timescales. Indeed, Eq. (2) tells us that when the dust growth timescale is much shorter than the timescale for particles to be dragged outwards by the gas, dust will be lost into the star efficiently. Therefore, to prevent dust particles from drifting towards the star, we must prevent them from growing to large sizes too fast.

3 Model

We use the previously presented DiskBuild protoplanetary disk model of Morbidelli et al. (2022), which includes dust and gas evolution. Here we summarise the model’s main features and refer the reader to the methods section of Morbidelli et al. (2022) for a detailed model description. We only detail the improvements made for this work.

We typically initiate the model with an empty disk and a protostar with an initial mass of 0.5 M_⊙. This is consistent with a Class-0 protostar. Subsequently, the disk is populated through an infall function describing the amount of mass added to the star-disk system as a function of time and distance to the star. The mass added to the disk is assumed to decay over time as exp(-t/T_infall), where t is time and T_infall is the infall timescale, a free parameter of the model. The time-integrated mass of the infall is scaled to result in a star-disk system with one solar mass.

The green line in Fig. 1 shows an example of the disk mass infall function for T_infall = 100 kyr.

The maximum distance within which material falls into the disk is the centrifugal radius, R_C. As recalled in Sect. 2, the classic recipe for the evolution of R_C over time is derived from the assumption of a rigidly rotating sphere of material (Shu 1977) and is (Hueso & Guillot 2005)

$$R_C\left(t\right)\simeq 53{\left(\frac{\omega }{{10}^{-14}{\text{s}}^{-1}}\right)}^2{\left(\frac{T}{10K}\right)}^{-4}{\left(\frac{M\left(t\right)}{1M_{\odot }}\right)}^3\text{\hspace{0.17em}au,}$$(5)

where ω is the angular speed of the cloud, T is the cloud temperature, and M(t) is the total mass of the star-disk system. For ω = 9 × 10⁻¹⁵ s⁻¹ and T = 15 K, R_C and never exceeds 10 au (orange line in Fig. 2). Foralargerangularspeed of, forexample, ω = 3.1 × 10⁻¹⁴ s⁻¹ the centrifugal radius will grow to 100 au. Therefore, depending on the angular speed of the molecular cloud the centrifugal radius can become very large. As a reference Pignatale et al. (2018) used ω = 1 × 10⁻¹⁴ s⁻¹ and T = 15 K for their study bringing R_C to 10.5 au.

Morbidelli et al. (2022) suggested that the alternative scenario, where R_C remains small throughout the infall process due to magnetic braking of the infalling material, should be appropriate, at least for our Solar System, to aid the formation of planetesimals at two locations within the disk. We thus adopt the prescription of Morbidelli et al. (2022) of

$$R_C\left(t\right)=\frac{0.35}{\sqrt{M_{\star }\left(t\right)}}\text{au,}$$(6)

where M_⋆ is the mass of the protostar in solar masses, M_⊙. We stress that the crucial assumption of Eq. (6) is not its exact form but that R_C remains small, particularly that it remains smaller than the condensation line of silicates and refractories.

There is an ongoing debate over the scale at which this disk forms (e.g. Machida et al. 2014; Masson et al. 2016) and we thus did not constrain ourselves to only exploring scenarios using Eq. (6). Thus, although we mainly present results using that prescription from Morbidelli et al. (2022), we also examine the effects of using the more traditional ‘Shu recipe’ (see the results in Sect. 5.3). In particular, we show results where the R_C grows to 10 and 100 au. The prescription of R_C forms our first main assumption in the model. Material falling closer than 0.05 au (the inner edge of our simulation domain) is assumed to be directly accreted onto the star.

The gas disk evolves under viscous heating and spreading. We use the usual definition of the viscosity v = αH²Ω (or, equivalently $v=\alpha c_s^2/\text{Ω}$), where Ω is the Keplerian frequency and $H=\sqrt{\frac{RT{r}^3}{\mu G{M}_{★}{M}_{\odot }}}$ is the scale height, with R the gas constant, µ the mean molecular weight of the gas, and G the gravitational constant. The scale height is computed self-consistently at each distance, r, of the disk by measuring the temperature, T. The viscosity parameter, α, is a free parameter and varies in time and with radial distance. As discussed in Sect. 2, it is reasonable to assume that α decays over time in a manner proportional to the disk infall function (two examples are shown in Fig. 1). However, the initial value of α – denoted α₀ – is considered a free parameter. A minimum value of a is set at 5 × 10⁻⁵, the order of magnitude of the effective turbulence generated by hydrodynamical mechanisms such as the vertical shear instability (Kumar & Coleman 1993; Urpin & Brandenburg 1998). In addition, at locations in the disk where it is gravitationally unstable or close to instability, the disk develops clumps and waves that also generate an effective viscosity. We take this into account by increasing a in those locations locally (see Eq. (8ff) in Morbidelli et al. 2022).

Of the infalling mass, 1% is considered dust and the rest gas (hydrogen), corresponding to the solar metallicity (Asplund et al. 2009). The dust is further split up into three sub-species: 1) all refractory species with a sublimation temperature above 1400 K, 2) silicates with a sublimation temperature of 1000 K, 3) water/ice with a sublimation temperature of 170 K. In reality, the sublimation temperature for silicates depends on the disk pressure and global chemistry (e.g. the C/O ratio). For instance, Morbidelli et al. (2020) showed that the silicate sublimation temperature could be 1060 K for P = 10⁻⁴ bar and C/O = 1.0. For simplicity, we kept the sublimation temperature of silicates at 1000 K. The species are assumed to have a relative abundance of 0.35/0.35/0.3. When the local disk temperature is above one of these sublimation temperatures, the corresponding dust specie is considered to be in the gaseous form and thus evolves in the same way the overall gas does.

In the part of the disk where a dust specie is in solid form, we track the size of dust particles, or rather its stokes number, St. The model has only one dust size at each radial distance, as in most codes. For dust size distributions that are dominated by the largest size, this is a good approximation and is indeed the result of dust growth models (e.g. Birnstiel et al. 2012; Paruta et al. 2016; Mattsson 2020; Stammler & Birnstiel 2022).

Because of the Eulerian nature of our code, we did not consider only the limit Stokes number given by the fragmentation barriers (3) and (4), where we assume Sc = 0.1 (Morbidelli et al. 2022). We also needed to take into account the fact that particles cannot be so large that they immediately drift out of a given cell. This drift boundary is defined as

$$\text{S}{t}_{\text{drift}}=0.055\frac{{\text{Σ}}_d}{{\text{Σ}}_g}\frac{r\text{Ω}}{\eta {\upsilon }_K},$$(7)

where r is the radial distance to the star. The barriers in Eqs. (4) and (7) are additions compared to the model published in Morbidelli et al. (2022), which only considered Eq. (3). The final particle size is determined from the lowest among St_growth (given by the growth algorithm with timescale ZΩ), St_frag, St_ddf, and St_drift.

We also improved the dust advection treatment in the code. For each cell, we now calculate the flux of particles out of the current cell to the lower and upper neighbouring cell based on the respective dust speed at the edge of the cell. Additionally, we compute the flux of particles from the lower and upper cells to the current cell. Taking into account all four possible loss and gain contributions is important, in particular, at the water snow line, because there the dust size can significantly change from one cell to the next. The particles beyond the snow line may drift towards the star, while those within the snow line may still drift away from the star.

The dust surface density is evolved, taking into account advection and diffusion. The back-reaction from the dust onto the gas is accounted for. At each timestep, the midplane volume density of the dust and gas is calculated. When the ratio of the two exceeds 0.5, we assume that planetesimal formation can occur via the SI in that ring, removing the dust in excess (Gole et al. 2020).

Fig. 1 Mass added to the disk as a function of time shown in units of solar masses per year (in green). In this example, the infall timescale Tinfall = 100 kyr. The yellow lines show the temporal evolution of the viscosity, α, for the two end-member cases, where α0 = 10−1 and α0 = 10−2.

Fig. 2 Centrifugal radius, RC, as a function of time for the two cases in this study. The orange line shows the prescription according to Eq. (5) (ω = 9 × 10−15 s−1 and T = 15 K; Shu 1977), assuming that the angular momentum of infalling material increases rapidly with time. This prescription, inadditionto onewhere RC grows to 100au, ledto the results presented in Sect. 5.3. The green line shows the function of Eq. (6) (Morbidelli et al. 2022), which describes an infall scenario where the infalling material loses angular momentum due to magnetic braking. This scenario was used for most of the cases presented in this work.

4 Model setups and constraints

As described in the introduction, the underlying problem that prevents dust from forming a large disk that extends far beyond the water snow line is that the dust grows too fast. We explore two ways to prevent dust from growing to a size large enough to make it drift towards the star during the expansion phase of the gas disk.

4.1 Expansion speed of the disk

First, a more rapid expansion of the gas disk – which in turn drags the dust particles in the radial direction when the Stokes number is small (Eq. (2)) – can transport dust into more distant regions of the disk before the dust has a chance to grow significantly. Faster expansion of the gas disk should manifest when the gas viscosity (α) is higher or the infall timescale (T_infall) is short. To explore the effect of these two parameters of our model, we varied them.

For the viscosity, we have one free parameter, the initial value of α at the beginning of the simulation, denoted α₀. Once α₀ is set, it decreases as described in Sect. 3 proportional to the mass added to the disk. Because the mass added to the disk decays over time, so will α. We chose to vary α₀ between 0.01 and 0.1 and steps of 0.01. The lower limit is consistent with the nominal case presented in Morbidelli et al. (2022). The upper limit might be considered quite high, but Kuznetsova et al. (2022) showed that for cases where the mass that is added to the disk is a large fraction of the disk mass itself, the disk wide α can reach large values (see their Fig. 8). In particular, when the infalling mass is on the same order as the disk mass, α reaches values of 0.1. Such a mass ratio is reached early in our simulations. Therefore, we believe such a high value of α₀ is plausible for a brief period at the beginning of the simulation. Remember that we let our α decay over time at the same rate as the infalling material decays (Fig. 1).

An increased viscosity has the added benefit of increasing the relative velocities between the dust particles and, therefore, their collision speeds. This results in more fragmentation and, thus, smaller particles, making it easier for the gas to transport the dust to large distances.

Regarding the infall timescale, we tested nine values of T_infall between 15 kyr and 630 kyr. A logarithmic spacing between cases was used. In combination with the ten different α0, we arrive at 90 simulations.

4.2 Fragmentation threshold of the dust

The second way to ensure particles reach greater distances in the disk is more straightforward. In our nominal cases, we follow the assumptions of Morbidelli et al. (2022) and impose a fragmentation threshold of υ_frag = 100 cms⁻¹ for refractory and silicate particles and υ_frag = 1000 cm s⁻¹ for icy particles beyond the water snow line.

However, we also tested a temperature-dependent fragmentation threshold prescription for icy particles:

$${\upsilon }_{\text{frag}}\left(T\right)={\upsilon }_0+{\upsilon }_C\text{Γ}{\left(T\right)}^{\frac{5}{6}},$$(8)

where T is the temperature, υ₀ = 100 cms⁻¹, υ_C = 1600 ms⁻¹ , and

$$\text{Γ}\left(T\right)={\text{Γ}}_C+{\text{Γ}}_{d0}\text{\hspace{0.17em}tanh}\left(\beta \left(T-T_0\right)\right),$$(9)

where Г_C = Г_d0 = 0.25, β = 0.105, and T₀ = 150. These parameters (υ₀, υ_C, Г_C, Г_d0), and T₀) where chosen to match the experimental data presented in Musiolik & Wurm (2019).

Figure 3 shows both the data from Musiolik & Wurm (2019, orange crosses; shifted to account for the different sublimation temperatures between the laboratory and the real disk) and Eq. (8) (light blue line). The fragmentation threshold decreases from 1000 cms⁻¹ to 100 cms⁻¹ between disk temperatures of 170 K and 120 K. The new prescription makes icy particles easier to break in cold regions of the disk. This limits their size and should help transport them to greater distances from the star.

For locations in the disk above the sublimation temperature of170 K(i.e. fordry particles), we retain a fragmentation threshold of 100 cm s⁻¹ whereas for locations with temperature below 170 K we use Eq. (8). We ran two sets of 90 simulations (the variations in α₀ and T_infall) for the nominal fragmentation threshold and the new temperature-dependent fragmentation threshold. The effects from rapidly expanding disks are expected to compound when also applying the new fragmentation threshold.

Fig. 3 Temperature-dependent fragmentation threshold, υfrag, for icy particles according to Eq. (8). The experimental data from Musiolik & Wurm (2019) have been shifted to account for the lower sublimation temperature of water in our model compared to the one in the experiment (~220 K).

4.3 Summary of assumptions

To summarise, there are four main assumptions that we explore in this work. The first is that the centrifugal radius, R_C, either grows to 10 and 100 au, according to Eq. (5) (Shu recipe), or remains small, according to Eq. (6). Our nominal simulations were performed with Eq. (6). Second, the initial disk viscosity, α₀, varies between 0.01 and 0.1. Third, the infall timescale, T_infall, varies between 15 kyr and 630 kyr. Finally, the fragmentation threshold for icy particles is either constant at 1000 cms⁻¹ (the nominal case) or temperature-dependent (according to Eq. (8)).

5 Results

5.1 Temperature-independent fragmentation threshold

First, we present the results from the cases where the nominal fragmentation threshold for dust particles and the small R_C according to Eq. (6) was used. In these cases, particles within the water snow line fragment at 100 cms⁻¹ while those outside at 1000 cms⁻¹.

As discussed in the introduction, the main factor limiting dust transport to large distances is the fast growth and subsequent inward drift of particles once they have crossed the snow line. Already very early on, for example after only 1000 yr, the dust particles just outside the snow line grow to the centimetre scale and effectively stop their outward radial motion. This is shown in panel a1 of Fig. 4, which depicts the results of the case where we have a small viscosity of α₀ = 0.01 and T_infall = 100 kyr (nominal case in Morbidelli et al. 2022). Particles just outside of the water snow line (dashed yellow line) have a size between 0.1 and 1 cm (Fig. 4a₂) and consequently have almost zero radial velocity (Fig. 4a₃). Because the gas continues to spread outwards, the dust and gas disks ‘decouple’, that is, the dust expansion lags the gas expansion. Therefore, even at this very early time, the dust disk is already smaller than the gas disk (fine black dashed line in Fig. 4).

In contrast, when the initial viscosity is much higher, for example α₀ = 0. 1 (Fig. 4b), the dust particles beyond the snow line are roughly an order of magnitude smaller (Fig. 4b₂)and thus retain a positive (outward) motion (Fig. 4b₃).The dust expansion keeps up with the gas expansion, and therefore the two disks retain the same size (Fig. 4b₁).

As expected, disks with larger viscosity expand faster. After 1000 yr of expansion, the gas disk with α₀ = 0.01 has expanded to roughly 4 au (measured where the gas surface density is 1 g cm⁻²). In contrast, the disk with α₀ = 0.1 has reached 10 au and is, therefore, more than double the size of the other (Fig. 4b).

For a given T_infall, the time a disk takes to reach 100 au, denoted T_{100 au}, decreases as the initial viscosity increases (Fig. 5). To measure the size of the disk, we used the location where the gas surface density takes a value of 1 g cm⁻². For the dust, we adopted a value 100 times smaller than that of the gas because of the metallicity of our infalling material being 1% (i.e. a value of 0.01 g cm⁻²). We are aware that this choice is somewhat arbitrary but have found it to be the definition that leads to the easiest and most reliable measure of the disk size, particularly for the dust. Other definitions, for example using the distance containing a certain fraction of the total mass, have proven unstable for the dust.

Figure 5 also shows that there is a transition of the expansion regime. For each value of α₀, the orange star on the corresponding curve indicates the viscous timescale, T_visc, of the disk, to be read on the horizontal axis. T_visc represents the average viscous timescale within 10 au at t = 0 for a disk with an aspect ratio of 6%. When the infall timescale is shorter than the viscous timescale (on the left side of the orange line), the expansion of the disk slows as the infall timescale decreases. In the extreme case where the infall timescale is much shorter than the viscous timescale, the disk’s ability to spread viscously is limited. Thus, the expansion timescale reaches a plateau. This can be clearly seen in the case of the lowest viscosity case.

In contrast, when the infall timescale is larger than the viscous timescale, the expansion of the disk slows with increasing infall timescale. This means the expansion is limited by the amount of material resupplied by the infall. In the most extreme cases T_{100 au} ~ 400 kyr (when α₀ and T_infall are minimal) and T_{100 au} ~ 10 kyr (when α₀ is maximal and T_infall is minimal). We baptise such a rapid expansion, reaching 100 au in just a few tens of thousands of years, the ‘inflationary phase’ of the disk.

Because T_infall in these tests varies by more than one order of magnitude, we might better measure T_{100 au} in units of T_infall. Indeed, the right panel of Fig. 5 shows the expansion time as a fraction of the infall timescale. In this view, we can recognise that for a given α₀, the expansion time as a fraction of T_infall always decreased with increasing T_infall. It is remarkable that if T_infall = T_visc, the value T_{100 au}/T_infall is independent of viscosity (i.e. the orange stars fall on a horizontal line).

Fig. 4 Disk properties for two cases with Tinfall = 100 kyr but different initial viscosities: α0 = 0.01 (a) and α0 = 0.1 (b). Both cases are shown at 1000 yr after disk formation. The top panels (a1 and b1) show the gas surface density (solid purple) and 100 times the dust surface density (solid green) as a function of distance to the star. The vertical dash lines show the centrifugal radius, RC (grey), the condensation lines of refractories (purple) and silicates (green), and the sublimation line of water (yellow). The fine vertical black line shows the distance where the dust surface density has a value of 10−2 g cm−2. The middle panels (a2 and b2) show the dust size, and the bottom panels (a3 and b3) the radial dust speed. A positive radial velocity represents motion away from the star, while a negative velocity is towards the star. Both of the shown simulations assume the nominal fragmentation thresholds of 1 m s−1 and 10 m s−1 for dry and icy particles, respectively.

Fig. 5 Time needed for the gas disk to expand to 100 au, T100 au, as a function of the infall timescale, Tinfall. The size of the disk is measured at the location where the gas surface density is 1 g cm−2. Each curve represents a different value of the initial viscosity, α0. For each value of α0, the orange star on the corresponding curve indicates the viscous timescale, Tvisc, of the disk, to be read on the horizontal axis. Tvisc represents the average viscous timescale within 10 au at t = 0 for a disk with an aspect ratio of 6%. The left panel shows T100 au, and the right panel shows T100 au scaled to Tinfall.

5.1.1 Mass and size of the dust disk

We measured the maximum dust mass a given disk holds 1 au beyond the snow line. To make sure the measurement was not contaminated by the dynamics around the snow line, we chose to exclude the dust mass just outside the snow line. We refer to this part of the disk as the ‘outer disk’. These masses and sizes are illustrated in Fig. 6.

Disks with an initial small viscosity result in small disks that contain little to no dust beyond the snow line (Fig. 6a₁). In these cases, the disks can be as small as 5 au. The most massive disks are formed with the highest viscosity and reach 60 M_⊕ and sizes between 30 and 50 au. For a given viscosity, the infall timescale plays a crucial role in determining the dust mass in the outer disk. The shorter the infall timescale, T_infall, is, the more massive the outer disk is (Fig. 6a₂). Therefore, short T_infall and large α₀ produced the largest and most massive outer disks. These disks thus satisfy our first criteria for good protoplanetary disks of the Solar System.

5.1.2 Planetesimal formation

To address our second criterion for good protoplanetary disks of the Solar System, we evaluate whether planetesimals form and at how many locations in the disk. Figure 7 summarises the mass of planetesimals formed in each of the disks. Because planetesimals typically form at up to two locations in the disk (Fig. 7, right panel), we split the results into ‘rocky’ planetesimals (forming at the silicate condensation line) and ‘icy’ planetesimals forming at or beyond the water snow line.

First, we observe that for most cases with T_infall > 100 kyr, no rocky planetesimals are formed. Second, for rocky planetesimals, there is an optimal viscosity given a T_infall. This is most clearly visible for T_infall = 39 kyr (the third line from the bottom). For this infall timescale, the optimum viscosity to produce rocky planetesimals is α₀ = 0.05. The planetesimal mass decreases for higher and lower values of α₀. When the viscosity is too low, the amount of mass transported to the planetesimal forming region is too small because of the lower radial velocity of the gas, and when the viscosity is too high, the dust cannot settle sufficiently in the midplane to trigger the SI. Third, the mass of icy planetesimals is maximised the larger the viscosity and the shorter the infall timescale. This comes from the fact that those disks are also the most massive beyond the snow line (Fig. 6a). Fourth, a small part of our parameter space (high viscosity and long infall timescales) does not form any planetesimals at any location in the disk. Fifth, the reservoirs of rocky and icy planetesimals have a similar order of magnitude in mass.

5.1.3 CAI transport to the outer disk

For the third criterion for good protoplanetary disks of the Solar System, we track high-temperature condensates. For this purpose, we introduce dust tracers, one for refractory particles that condensate at the refractory line, and a second for refractories that never sublimated. A fraction of the high-temperature condensates will be CAIs, but in our model, we just refer to such particles as potential CAIs because we do not track the full condensation sequence of refractories but rather just treat all refractories as one species of dust. Nevertheless, this lets us determine the locations in the disk that will be enriched or depleted in CAIs.

The ability of the disk to transport CAIs to the outer disk and retain them there depends again on the viscosity of the disk and the infall timescale. In particular, the transport of CAIs is promoted when the centrifugal radius is smaller than the refractory condensation line. If the infall timescale is too long (larger than ~200 kyr for α₀ = 0.05) the disk is rather cold from the beginning, and therefore the refractory condensation line (defined as T = 1400 K) is located inside R_C, and no CAIs are transported to the outer disk (Fig. 8). In contrast, when the infall timescale is short (less than ~100 ky) CAIs are efficiently transported to the outer disk, but then drift back into the inner disk due to the fast evolution of the disk, which transitions to a fully accreting disk within 3–4T_infall. While we show these results for α₀ = 0.05 they are qualitatively the same for other initial viscosities. For larger initial viscosities, the infall timescale where the disk is too cold to create CAIs is shorter (e.g. at T_infall ~ 150 kyr for α₀ = 0.1). Conversely, this transition happens at larger infall timescales when the viscosity is smaller (e.g. at T_infall > 400 kyr for α₀ = 0.01). But in all cases, neither very short nor long T_infall are favoured for the transport of CAIs to the outer disk.

The smaller the initial viscosity is, the larger the fraction of the disk that is populated by CAIs. For example, when α₀ < 0.05 for T_infall = 100 kyr the inner disk gets similarly enriched with CAIs as the outer disk (Fig. 9). When in addition to a low initial viscosity, the infall timescale is also short, then the entire disk is populated by potential CAIs. Such disks would clearly not match the observations. Yet, the larger the initial viscosity, the clearer the divide is between a CAI-enriched outer and CAI-depleted inner disk. The presence of CAIs in outer planetesimals thus suggests a high initial viscosity with the associated rapid expansion phase of the disk. This appears to be consistent with large, kinetic, Si isotopic variations observed in refractory inclusions, which suggest a turbulent environment during condensation (e.g. Marrocchi et al. 2019).

In all of our simulations, we kept the Schmidt number at Sc = 0.1. A higher Schmidt number of, for example, Sc = 1 would aid the transport of CAIs to the outer disk. However, the larger Sc the more the dust will have difficulty settling in the midplane and thus tend to make planetesimal formation more difficult.

Fig. 6 Maximum dust mass 1 au beyond the snow line (SL) as a function of the size of the dust disk, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$, at the time when the disk has reached that maximum mass. The top panels show the dependence on the viscosity, α0, and the bottom row the dependence on the infall timescale, Tinfall. Panels a1 and a2 show the results for the temperature-independent fragmentation threshold, and the panels b1 and b2 show the results for the temperature-dependent fragmentation threshold.

5.2 Temperature-dependent fragmentation threshold

In the case where we impose the temperature-dependent fragmentation threshold beyond the snow line (see Sect. 4 and Fig. 3), we expected that dust fragments more easily and therefore, the outer disk gets populated with more mass. Indeed, all disks now have at least 10 M_⊕ in the outer disk (Fig. 6). Though the disks are, in general, not significantly more massive (10–70 M_⊕ compared to 0–60 M_⊕), the disks with the temperature-dependent fragmentation threshold are much larger (30–80 au instead of 5–50 au). Thus there is, as expected, a general shift to more massive and larger outer disks.

This shift of dust mass from the inner to the outer disk has clear consequences. We now have significantly more icy plan-etesimals than rocky ones (Fig. 10). For some combination of parameters α₀ and T_infall (e.g. 0.07 ≤ α₀ ≤ 0.1 and 40 kyr ≤ T_infall ≤ 100 kyr), a couple of Earth masses of rocky planetes-imals form together with a couple of tens of Earth masses of icy planetesimals. This is in very good agreement with the structure of the Solar System, with massive giant planets’ cores and small terrestrial planets.

Similarly to the temperature-independent fragmentation threshold, there are little to no planetesimals when T_infall > 100 kyr. The delineation is even a bit clearer. Nevertheless, the part of parameter space with two planetesimals rings is roughly equally large irrespective of the fragmentation threshold.

Concerning CAI transport, the overall behaviour is similar to the case with the nominal fragmentation threshold. But, because particles are more easily transported to the outer disk, CAIs also reach much greater distances.

Fig. 7 Mass of planetesimals formed at different locations in the disk. Left panel: Total mass of ‘rocky’ planetesimals as a function of the infall timescale, Tinfall, and the viscosity, α0. Rocky planetesimals are the ones that form around the silicate condensation line. White areas are disks that do not produce any planetesimals, while grey squares indicate disks that produce between zero and two Earth masses of planetesimals. Centre panel: Total mass of ‘icy’ planetesimals for a given disk. Icy planetesimals are the ones formed near the water snow line, typically beyond it. Right panel: Number of locations, i.e. rings, where planetesimals form.

Fig. 8 Ratio between the surface density of high-temperature condensates, Cpot.CAI, and the surface density of all refractory particles, Ctot.ref., for different values of Tinfall and α0 = 0.05. These simulations assume the nominal fragmentation thresholds of 1 m s−1 and 10 m s−1 for dry and icy particles, respectively. The three solid white lines are, from closest to the star to farthest, the sublimation lines of refractories, silicates, and water. The white course dashed line is RC , and the fine dashed line is the dust disk size, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$.

Fig. 9 Ratio between the surface density of high-temperature condensates, Σpot.CAI, and the surface density of all refractory particles, Σtot.ref, for different values of α0 and Tinfall = 100 kyr. These simulations assume the nominal fragmentation thresholds of 1 m s−1 and 10 m s−1 for dry and icy particles, respectively. The three solid white lines are, from closest to the star to farthest, the sublimation lines of refractories, silicates, and water. The white course dashed line is RC, and the fine dashed line is the dust disk size, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$.

5.3 Shu infall

Because our prescription of the infall is somewhat unconventional (i.e. the centrifugal radius R_C ~ 0.35 au; Eq. (6)), we also tested the more common assumption according to Shu (1977). In the ‘Shu case’, the R_C rapidly grows from 1 au to 8 au (Fig. 2, Eq. (5) with ω = 9 × 10⁻¹⁵ s⁻¹ and T = 15 K). This is because the molecular cloud is assumed to be a rigidly rotating body and angular momentum is conserved (i.e. no magnetic braking). Therefore, gas with small angular momentum collapses into the disk first, close to the star. Later, outer shells with larger specific angular momenta fall at greater distances. This behaviour is in contrast to our preferred cases described above, where magnetic braking reduces the angular momentum of the infalling gas to roughly a fixed value independently of the initial angular momentum of the gas in the molecular cloud.

A major consequence of the ‘Shu-type’ infall is connected to the radial gas speed. The disk within Rc is an accretion disk, that is to say, the radial gas velocity, v_r,ɡ, is negative (Fig. 11). Therefore, dust within R_C will also always have a negative radial velocity (v_r,d < 0). Outside of R_C, the disk can spread viscously outwards (${\upsilon }_r^g>0$; Fig. 11), and therefore small dust particles will also have a positive radial motion as long as they do not grow large enough to feel the headwind of the gas and start drifting back towards the star.

We tested two different angular velocities, c , of the molecular cloud. Once with ω = 10⁻¹⁴ s⁻¹ resulting in a maximum R_C of roughly 10 au as shown in Fig. 11 and once with ω = 3.1 × 10⁻¹⁴ s⁻¹ resulting in a maximum R_C of roughly 100 au. The temperature of the molecular clouds is assumed to be 15 K in both cases. We use here the evolution of R_C according to Eq. (3) of Hueso & Guillot (2005). The prescription of the T_infall and α₀ remain the same as above.

In all cases studied, the Shu-type infall has no difficulty producing large and massive disks (Fig. 12). When R_C grows to 10 au, and we use the nominal temperature-independent fragmentation threshold, the disks are between 10 and 100 au and have masses between 2 and 200 M_⊕ (Fig. 12a₁. For the same molecular cloud angular velocity but with the temperature-dependent fragmentation threshold, the disks are overall larger and more massive in particular for the cases with small α₀. The sizes and masses are also confined to 80–150 au and 30–300 M_⊕ (Fig. 12b₁).

When R_C grows to 100 au the disks are even larger and more massive. For the temperature-independent fragmentation threshold, the disks are between 40 and 100 au and have masses between 30 and 600 M_⊕ (Fig. 12a₂). For the temperature-dependent fragmentation threshold, the disk sizes and masses are only weakly dependent on α₀ and T_infall. These disks are between 150 and 400 au and have masses between 300 and 700 M_⊕ (Fig. 12b₂), and therefore very massive and large.

When we prescribe the ‘Shu-infall’ particles in the inner disk (within the water snow line) drift rapidly towards the star (Fig. 11). This does not allow them to pile up at the silicate sublimation line, and therefore no rocky planetesimals are formed in any of the cases (left panels in Fig. 13). Additionally, even at the water snow line, we observe only sparse formation of planetesi-mals (centre panel in Fig. 13). This result is largely independent of which angular velocity of the molecular cloud we used and which fragmentation threshold is applied.

Our results differ from the results found by Drążkowska & Dullemond (2018). We do not find any planetesimal formation during the phase when the snow line moves outwards. This might be caused by the different assumptions of the disk infall prescription. We assume that the mass added to the disk decays over time while a constant function with a sudden cut-off is assumed in Drążkowska & Dullemond (2018). Additionally, we find many fewer planetesimals at the snow line. We believe that Drążkowska & Dullemond (2018) overestimated the amount of water vapour in their disks due to a difference in treatment of the inner disk boundary condition for water vapour from that of hydrogen. This supports planetesimal formation.

Finally, we also studied the transport of CAIs in such disks. As expected no CAIs are able to reach the outer disk, or even the terrestrial planet region (Fig. 14). The example shown in the left panel of Fig. 14 assumes α₀ = 0.05, T_infall = 100 kyr, the temperature-dependent fragmentation threshold, and R_C growing to roughly 10 au but is representative of almost all combinations of α₀ and T_infall. The only exception is for α₀ = 0.01 and T_infall < 25 kyr (right panel of Fig. 14). In this case, some potential CAIs are produced and transported to the outskirts of the disk (at roughly 100 au). For cases where R_C grows to roughly 100 au, the situation is even worse because in none of the cases are there any potential CAIs in the disk. This behaviour is not surprising. The inward motion of the gas prevents any CAIs from being transported to the terrestrial planet region or outer disk.

Our results are broadly consistent with those of Pignatale et al. (2018) in that the fraction of CAIs is largest in the outermost part of the disk (towards the edge of the disk itself). Pignatale et al. (2018) assume a constant function for the infall of material into the disk, whereas we assume a decaying function. Assuming a constant source function results in R_C growing much slower than in our cases. This in turn extends the period during which R_C is smaller than the refractory condensation line. Therefore, CAIs can be produced for longer and transported into more distant regions of the disk. This way the disk generally can be more enhanced with CAIs than in our cases.

Fig. 10 Same as in Fig. 7 but for the temperature-dependent fragmentation threshold.

Fig. 11 Radial velocity of the gas at 100 kyr as a function of the radial distance to the star for the case when Tinfall = 100 kry and α0 = 0.01. The solid line shows the results assuming the infall prescription of Morbidelli et al. (2022), while the dashed line is for the ‘Shu infall’ (Shu 1977). The vertical grey lines denote the position of the centrifugal radius, RC, of the respective case.

Fig. 12 Maximum dust mass 1 au beyond the snow line (SL) as a function of the size of the dust disk, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$, at the time when the disk has reached that maximum mass. Results for the Shu-type infall where RC grows to roughly 10 au are shown in panels a1 and a2, while the ones where RC grows to roughly 100 au are shown in panels b1 and b2. Panels with subscript 1 show the results for the temperature-independent fragmentation threshold, while panels with subscript 2 show the results for the temperature-dependent fragmentation threshold.

6 Discussion and conclusion

Infall of material into protoplanetary disks occurs relatively close to the star – typically a distance much smaller than the observed disk sizes. The disks, therefore, undergo an initial phase of viscous spreading (Lynden-Bell & Pringle 1974; Hueso & Guillot 2005). The dust, on the one hand, is entrained in the outward motion of the gas and, on the other hand, is slowed by the sub-Keplerian motion of the gas (see Eq. (2)), which causes its inward drift. Whether the radial outward entrainment or sub-Keplerian drag dominates the dust motion depends on the particle size.

A key parameter in any protoplanetary disk model is the so-called centrifugal radius, R_C. This is the radius in the disk where the angular momentum is the same as that of the infalling material. If, for example, the pre-stellar cloud rotates as a rigid sphere (Shu 1977), then shells of material closer to the centre collapse first and, having a small specific angular momentum, fall very close to the protostar. Outer shells, with larger angular momenta, will fall at greater distances and in a later stage in disk formation (Shu 1977). In such scenarios, which we refer to as Shu-type infall models, R_C grows with time. Contrary to this, magnetic braking can remove angular momentum from the infalling material. This can cause the material to fall close to the star irrespective of the initial angular momentum of the material.

A disk formation and evolution scenario for the Solar System must satisfy at least the following three requirements:

1. It must develop an extended disk of gas and dust (up to 45 au for dust).

2. In at least two distinct locations in the disk, the dust/gas ratio must be able to rise sufficiently to produce planetesimals and explain the early formation of NC- and CC-iron meteorite parent bodies.

3. Particles that condensed at high temperatures (i.e. CAIs) must be able to reach large heliocentric distances, that is, be transported from the star’s proximity to large distances.

We find that scenarios that employ a Shu-type infall model with an associated large R_C are very successful in achieving requirement 1, as they easily result in large and massive disks. Yet they fail to produce planetesimals at two locations in the disk (requirement 2) and transport CAIs to the outer disk (requirement 3). Therefore, these scenarios are bad candidates for describing the Solar System protoplanetary disk. On the other hand, we show that a disk fed by material with a small R_C can satisfy all three requirements, in particular when the initial viscosity is large and the infall timescale is of the order of or smaller than 100 kyr.

The main results from our nominal disks with a small centrifugal radius, R_C, can be summarised as follows.

1. The larger the initial viscosity, a0, the larger the outer dust disk.

2. The shorter the infall timescale, T_infall, the more massive the outer dust disk.

3. Therefore, an initial inflationary expansion phase is needed to produce large, massive dust disks. The disk can reach a size of 100 au within a few tens of thousands of years.

4. A temperature-dependent fragmentation threshold is more realistic and results in significantly larger and slightly more massive dust disks because particles are more fragile and therefore remain smaller at cold temperatures.

5. No rocky and very few icy planetesimals form when T_infall > 100 kyr.

6. The largest mass of icy planetesimals forms when α₀ > 0.05.

7. There is an optimum a0 that maximises the mass of rocky planetesimals. For example, for T_infall = 39 kyr, it is α₀ = 0.05.

8. The temperature-dependent fragmentation threshold results in more icy than rocky planetesimals (by roughly a factor of 10). This is in contrast to the conventional case, where the two are of the same order. This is a direct consequence of the temperature-dependent fragmentation threshold resulting in more massive outer disks.

Although our disks with a small R_C can satisfy all three of our requirements, there are two additional related requirements that need to be met eventually. Observations show that protoplane-tary disks have a long lifespan, 3–4 million yr (Andrews 2020). All dust in our models (even in the Shu-type infall models) drifts into the star on a timescale of a few hundred thousand years. Therefore, the entire dust disk is lost on that timescale. Not only does this prevent us from explaining long-lived disks, but our disks are also not able to produce a generation of planetesimals late enough to avoid differentiation, because no dust is available at these later times. The retention of a large disk and the production of a population of planetesimals that forms late are two additional requirements to properly describe the protoplanetary disk of the Solar System.

Clearly, our model lacks some additional disk processes that can prevent the loss of dust from the disk. For example, once the disk viscosity is sufficiently small, MHD effects might become dominant and structures (rings and gaps) might appear, impeding dust drift (e.g. Béthune et al. 2016; Riols et al. 2020). This will be the object of future work.

Fig. 13 Same as in Fig. 7 but for the temperature-dependent fragmentation threshold and the Shu-infall model with a maximum RC of roughly 10 au (top panels) and 100 au (lower panels). The results for the temperature-independent fragmentation threshold are qualitatively the same. White areas represent cases without any planetesimal production, while grey areas are cases where the planetesimal mass is between zero and two M⊕.

Fig. 14 Ratio between the surface density of high-temperature condensates, Cpot.CAI, and the surface density of all refractory particles, Ctot.ref.. The case on the left assumes α0 = 0.05 and Tinfall = 100 kyr, and the case on the right assumes α0 = 0.01 and Tinfall = 15 kyr. In both cases, we use the temperature-dependent fragmentation threshold and the RC that grows to roughly 10 au. The three solid white lines are, from closest to the star to farthest, the sublimation lines of refractories, silicates, and water. The dashed white line is RC.

Acknowledgements

We acknowledge the funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement no. 101019380). Additionally, we acknowledge support from programme ANR-20-CE49-0006 (ANR DISKBUILD). We thank Sebastien Charnoz, Yves Marrocchi and Francesco Lovascio for reading the manuscript and providing helpful comments. We thank the anonymous reviewer for their constructive and useful comments that helped us improve the paper.

References

All Figures

	Fig. 1 Mass added to the disk as a function of time shown in units of solar masses per year (in green). In this example, the infall timescale Tinfall = 100 kyr. The yellow lines show the temporal evolution of the viscosity, α, for the two end-member cases, where α0 = 10−1 and α0 = 10−2.
In the text

Fig. 2 Centrifugal radius, RC, as a function of time for the two cases in this study. The orange line shows the prescription according to Eq. (5) (ω = 9 × 10−15 s−1 and T = 15 K; Shu 1977), assuming that the angular momentum of infalling material increases rapidly with time. This prescription, inadditionto onewhere RC grows to 100au, ledto the results presented in Sect. 5.3. The green line shows the function of Eq. (6) (Morbidelli et al. 2022), which describes an infall scenario where the infalling material loses angular momentum due to magnetic braking. This scenario was used for most of the cases presented in this work.

In the text

	Fig. 3 Temperature-dependent fragmentation threshold, υfrag, for icy particles according to Eq. (8). The experimental data from Musiolik & Wurm (2019) have been shifted to account for the lower sublimation temperature of water in our model compared to the one in the experiment (~220 K).
In the text

Fig. 4 Disk properties for two cases with Tinfall = 100 kyr but different initial viscosities: α0 = 0.01 (a) and α0 = 0.1 (b). Both cases are shown at 1000 yr after disk formation. The top panels (a1 and b1) show the gas surface density (solid purple) and 100 times the dust surface density (solid green) as a function of distance to the star. The vertical dash lines show the centrifugal radius, RC (grey), the condensation lines of refractories (purple) and silicates (green), and the sublimation line of water (yellow). The fine vertical black line shows the distance where the dust surface density has a value of 10−2 g cm−2. The middle panels (a2 and b2) show the dust size, and the bottom panels (a3 and b3) the radial dust speed. A positive radial velocity represents motion away from the star, while a negative velocity is towards the star. Both of the shown simulations assume the nominal fragmentation thresholds of 1 m s−1 and 10 m s−1 for dry and icy particles, respectively.

In the text

Fig. 5 Time needed for the gas disk to expand to 100 au, T100 au, as a function of the infall timescale, Tinfall. The size of the disk is measured at the location where the gas surface density is 1 g cm−2. Each curve represents a different value of the initial viscosity, α0. For each value of α0, the orange star on the corresponding curve indicates the viscous timescale, Tvisc, of the disk, to be read on the horizontal axis. Tvisc represents the average viscous timescale within 10 au at t = 0 for a disk with an aspect ratio of 6%. The left panel shows T100 au, and the right panel shows T100 au scaled to Tinfall.

In the text

Fig. 6 Maximum dust mass 1 au beyond the snow line (SL) as a function of the size of the dust disk, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$, at the time when the disk has reached that maximum mass. The top panels show the dependence on the viscosity, α0, and the bottom row the dependence on the infall timescale, Tinfall. Panels a1 and a2 show the results for the temperature-independent fragmentation threshold, and the panels b1 and b2 show the results for the temperature-dependent fragmentation threshold.

In the text

Fig. 7 Mass of planetesimals formed at different locations in the disk. Left panel: Total mass of ‘rocky’ planetesimals as a function of the infall timescale, Tinfall, and the viscosity, α0. Rocky planetesimals are the ones that form around the silicate condensation line. White areas are disks that do not produce any planetesimals, while grey squares indicate disks that produce between zero and two Earth masses of planetesimals. Centre panel: Total mass of ‘icy’ planetesimals for a given disk. Icy planetesimals are the ones formed near the water snow line, typically beyond it. Right panel: Number of locations, i.e. rings, where planetesimals form.

In the text

Fig. 8 Ratio between the surface density of high-temperature condensates, Cpot.CAI, and the surface density of all refractory particles, Ctot.ref., for different values of Tinfall and α0 = 0.05. These simulations assume the nominal fragmentation thresholds of 1 m s−1 and 10 m s−1 for dry and icy particles, respectively. The three solid white lines are, from closest to the star to farthest, the sublimation lines of refractories, silicates, and water. The white course dashed line is RC , and the fine dashed line is the dust disk size, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$.

In the text

Fig. 9 Ratio between the surface density of high-temperature condensates, Σpot.CAI, and the surface density of all refractory particles, Σtot.ref, for different values of α0 and Tinfall = 100 kyr. These simulations assume the nominal fragmentation thresholds of 1 m s−1 and 10 m s−1 for dry and icy particles, respectively. The three solid white lines are, from closest to the star to farthest, the sublimation lines of refractories, silicates, and water. The white course dashed line is RC, and the fine dashed line is the dust disk size, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$.

In the text

	Fig. 10 Same as in Fig. 7 but for the temperature-dependent fragmentation threshold.
In the text

	Fig. 11 Radial velocity of the gas at 100 kyr as a function of the radial distance to the star for the case when Tinfall = 100 kry and α0 = 0.01. The solid line shows the results assuming the infall prescription of Morbidelli et al. (2022), while the dashed line is for the ‘Shu infall’ (Shu 1977). The vertical grey lines denote the position of the centrifugal radius, RC, of the respective case.
In the text

Fig. 12 Maximum dust mass 1 au beyond the snow line (SL) as a function of the size of the dust disk, $R_{0.01\text{\hspace{0.17em}g\hspace{0.17em}\hspace{0.17em}}{\text{cc}}^{-1}}$, at the time when the disk has reached that maximum mass. Results for the Shu-type infall where RC grows to roughly 10 au are shown in panels a1 and a2, while the ones where RC grows to roughly 100 au are shown in panels b1 and b2. Panels with subscript 1 show the results for the temperature-independent fragmentation threshold, while panels with subscript 2 show the results for the temperature-dependent fragmentation threshold.

In the text

Fig. 13 Same as in Fig. 7 but for the temperature-dependent fragmentation threshold and the Shu-infall model with a maximum RC of roughly 10 au (top panels) and 100 au (lower panels). The results for the temperature-independent fragmentation threshold are qualitatively the same. White areas represent cases without any planetesimal production, while grey areas are cases where the planetesimal mass is between zero and two M⊕.

In the text

Fig. 14 Ratio between the surface density of high-temperature condensates, Cpot.CAI, and the surface density of all refractory particles, Ctot.ref.. The case on the left assumes α0 = 0.05 and Tinfall = 100 kyr, and the case on the right assumes α0 = 0.01 and Tinfall = 15 kyr. In both cases, we use the temperature-dependent fragmentation threshold and the RC that grows to roughly 10 au. The three solid white lines are, from closest to the star to farthest, the sublimation lines of refractories, silicates, and water. The dashed white line is RC.

In the text