# An inflationary disk phase to explain extended protoplanetary dust disks

Raphael Marschall  CNRS, Observatoire de la Côte d’Azur, Laboratoire J.-L.
Lagrange, CS 34229, 06304 Nice Cedex 4, France  Alessandro Morbidelli  CNRS,
Observatoire de la Côte d’Azur, Laboratoire J.-L. Lagrange, CS 34229, 06304
Nice Cedex 4, France

(June 30, 2022)

###### Abstract

Context: Understanding planetesimal formation is an essential first step to
understanding planet formation. The distribution of these first solid bodies
will drive the locations where planetary embryos can grow, eventually leading
to fully-fledged planets.
Aim: We seek to understand the parameter space of possible protoplanetary disk
formation and evolution models of our Solar System. A good protoplanetary disk
scenario for the Solar System must meet at least the following three criteria:
1) It must produce an extended gas and dust disk (e.g., 45 au for the dust);
2) within the disk, the local dust-to-gas ratio in at least two distinct
locations must sufficiently increase to explain the early formation of the
parent bodies of non-carbonaceous and carbonaceous iron meteorite; and 3) dust
particles, which have condensed at high temperatures (i.e., calcium–aluminium-
rich inclusion, CAIs), must be transported to the outer disk. Though able to
satisfy a combination of these three criteria, current protoplanetary disk
models have not been successful in recreating all three features. We aim to
find scenarios that satisfy all three criteria.
Methods: In this study, we use a 1D disk model that tracks the evolution of
the gas and dust disks. Planetesimals are formed within the disk at locations
where the streaming instability can be triggered. We explore a large parameter
space to study the effect of the disk viscosity, the timescale of infall of
material into the disk, the distance within which material is deposited into
the disk, and the fragmentation threshold of dust particles.
Results: We find that scenarios with a large initial disk viscosity
($\alpha>0.05$), relatively short infall timescale
($T_{\texttt{infall}}<100-200$ kyr), and a small centrifugal radius
($R_{C}\sim 0.4$ au; the distance within which material falls into the disk)
result in disks that satisfy all three criteria for a good protoplanetary disk
of the Solar System. The large initial viscosity and short infall timescale
result in a rapid initial expansion of the disk, which we dub the inflationary
phase of the disk. Furthermore, a temperature-dependent fragmentation
threshold, which mimics that cold icy particles break more easily, results in
larger and more massive disks. This in turn, results in more “icy” than
“rocky” planetesimals. Such scenarios are also better in line with our Solar
System, which has small terrestrial planets and massive giant planet cores.
Finally, we find that scenarios with large $R_{C}$ cannot transport CAIs to
the outer disk and do not produce planetesimals at two locations within the
disk.

##  1 Introduction

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

Observations of protoplanetary dust disks show two distinct properties: they
are large and long-lasting. Their sizes range from $10-500$ au with typical
sizes $\sim 30$ au (Tripathi et al., 2017; Andrews et al., 2018; Hendler et
al., 2020), and 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
these were 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 (Nesvorný et al.,
2017; Nesvorný, 2018). Furthermore, 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 have likely formed late (e.g., Nakashima et al.,
2015; Nimmo et al., 2018; Neumann et al., 2018), i.e., several million years
after the formation of the first solids, the so-called calcium–aluminium-rich
inclusion (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 $\mathrm{CO_{2}}$ and CO (e.g.,
Eberhardt et al., 1987; Morse et al., 2015; Gasc et al., 2017). The presence
of these highly volatile species also in very large comets ($\sim 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 of they formed late, at a time when most short-lived
radioactive elements (e.g. 26Al) 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-3; 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 both samples from differentiated and un-differentiated parent
bodies. The latter formed significantly later – up to 5 million years 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 will focus in this work on the first generation of planetesimal,
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 its large size.

Focusing on the first generation of planetesimals, the differentiated parent
bodies of iron meteorites, we find that these can be divided into two
isotopically distinct groups akin to carbonaceous chondrites (CC) and non-
carbonaceous chondrites (NC) (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 have formed simultaneously, they must
form at distinctly different locations in the disk that can have a different
disk composition. Therefore, our second requirement for a good model of the
Solar System disk is that it produces planetesimals at two distinct locations
in the disk.

Finally, the oldest Solar System solids, CAIs, are thought to have formed as
high temperature condensates very close (few tenths of an au) to the proto-Sun
(Scott & Krot, 2003). The age of CAIs sets what is usually considered time
zero of Solar System formation (see review by Chaussidon & Liu, 2015). Their
age is $4,567.30\pm 0.16$ million years 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 4,568.7 Myr (Piralla et al., 2023; Desch et al.,
2022). The duration of CAI formation appears to be very short, from $\sim 100$
kyr (Connelly et al., 2012) to just $\sim 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 formed the furthest away
are more enriched with them. This implies that these high-temperature
condensates have been transported efficiently to the outer disk, so that the
latter 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 present already 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 properties:

  1. 1.

it must develop an extended disk of gas and dust (up to 45 au for the dust);

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

  3. 3.

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

In this work, we try to build such a scenario. In section 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
(Sec. 3) before discussing the model setup (Sec. 4). In particular, we will
describe four assumptions’ influence on satisfying the Solar System
constraints. These are i) the centrifugal radius, $R_{C}$; ii) the initial
viscosity of the disk, $\alpha_{0}$; iii) the infall timescale of material
onto the disk, $T_{\text{infall}}$; and iv) the effect of a temperature-
dependent fragmentation threshold for icy particles. Our results are presented
in Sec. 5.

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

##  2 Key processes in disk evolution and planetesimal formation

As anticipated in the introduction, we start discussing key processes in the
formation and evolution of the disk and planetesimal accretion, focusing on
the unknowns we will parametrise and test in our models.

### 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, e.g., 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 proto-star. 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 pre-stellar 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 proto-star, sketched in Fig.
19). The formation of such small disks is observed in some
magnetohydrodynamics (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 will 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^{5}$ y, with large uncertainties (Larson, 1969; Vaytet
et al., 2018; Wurster et al., 2021, 2022), so that 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
don’t 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 suffer strong Raynold stresses that act as an effective viscosity
(Kuznetsova et al., 2022). Thus, it seems legitimate to assume that a disk
which is still accreting mass has a viscosity proportional to the mass infall
rate, but the proportionality factor is poorly constrained, and therefore we
will consider different values in our study.

### Motion of dust particles within the disk

For disks forming with a small $R_{C}$ where, e.g., 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\Sigma_{g}}\quad,$ (1)

where $a$ is the diameter of the dust particle, $\rho_{d}$ is the particle
solid density, and $\Sigma_{g}$ is the gas surface density. The radial dust
velocity, $v_{r}^{d}$, can then be written as

$v_{r}^{d}=\frac{2{\text{S}t}}{1+{\text{S}t}^{2}}v_{t}^{g}+\frac{1}{1+{\text{S}t}^{2}}v_{r}^{g}\quad,$
(2)

where $v_{t}^{g}$ and $v_{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, $v_{t}^{g}=\eta v_{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 ${\text{S}t}\ll 1$, the radial dust speed is dominated by
the radial gas speed ($v_{r}^{d}\propto v_{r}^{g}$, Eq. 2). Thus, when the
dust is small, it initially expands outwards from $R_{C}$ with the gas. Once
the dust has grown sufficiently (i.e., ${\text{S}t}\sim 1$), the tangential
speed of the gas can become the dominant factor in Eq. 2. Because the gas is
sub-keplerian $v_{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.

### Dust growth

Particles grow on a timescale $1/Z\Omega$, where $Z=\Sigma_{d}/\Sigma_{g}$ 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 among dust particles and when the latter is larger than the
fragmentation velocity $v_{\text{frag}}$, dust particles cannot coagulate
further but rather break upon collisions.

The largest Stokes number that particles can acquire by coagulation is
estimated to be (Drążkowska & Alibert, 2017) the minimum between:

${\text{S}t}_{\text{frag}}=\frac{0.37v_{\text{frag}}^{2}}{3\texttt{Sc}\leavevmode\nobreak\
\alpha c_{s}^{2}}\quad,$ (3)

and

${\text{S}t}_{\text{ddf}}=\frac{0.37v_{\text{frag}}}{2|\eta v_{K}|}\quad,$ (4)

where $\alpha$ is the gas viscosity parameter, following the assumption that
the viscosity $\nu=\alpha c_{s}^{2}/\Omega$ (Shakura & Sunyaev, 1973), Sc is
the Schmidt number relating viscous angular momentum transfer to turbulent
diffusion, and $c_{s}$ is the local sound speed. Eq. (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 $v_{\text{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 $v_{\text{frag}}=100$ cm/s
for refractory and silicate particles whereas $v_{\text{frag}}=1,000$ cm/s for
icy particles beyond the water snowline. 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 will explore an additional fragmentation threshold
prescription for icy particles, which is temperature dependent. Similarly, it
may be possible that silicate particles become more sticky when their
temperature is close to sublimation (Pillich et al., 2021) but, awaiting
experimental confirmations, we don’t yet consider this possibility in our
model.

### 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 planetesimals (e.g., Klahr &
Schreiber, 2020; Nesvorný et al., 2021; Polak & Klahr, 2023).

Previous models exploring the formation of planetesimals within a disk have
focused on static disks, i.e., 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
snowline (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 snowline 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, e.g., 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
exhibit the problem that they cannot satisfy at least two of our requirements.
These disks typically don’t 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 snowline ($\sim 5$ au). This is because particles beyond the snowline
rapidly grow and drift back towards the proto-star on much shorter time scales
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, equation 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 proto-star with an
initial mass of $0.5M_{\odot}$. 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_{\text{infall}})$, where $t$ is time and $T_{\text{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 Figure 1 shows an example of the disk mass infall function
for $T_{\text{infall}}=100$ kyr.

Figure 1: The mass added to the disk as a function of time is shown in green
and units of solar masses per year. In this example, the infall timescale,
$T_{\text{infall}}=100$ kyr. The yellow lines show the temporal evolution of
the viscosity, $\alpha$, for the two end-member cases where
$\alpha_{0}=10^{-1}$ and $\alpha_{0}=10^{-2}$. Figure 2: The centrifugal
radius, $R_{C}$, is shown as a function of time for the two cases in this
study. The orange line shows the prescription according to Eq. 5
($\omega=9\times 10^{-15}$ s-1 and $T=15$ K; Shu, 1977), assuming that the
angular momentum of infalling material increases rapidly with time. This
prescription, in addition to one where $R_{C}$ grows to 100 au, has been used
for results in Sec. 5.3. The green line shows the function of Eq. 6
(Morbidelli et al., 2022) describing an infall scenario where the infalling
material loses angular momentum due to magnetic braking. The latter was used
for most cases presented in this work.

The maximum distance within which material falls into the disk is the
centrifugal radius, $R_{C}$. As recalled in section 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}(t)\simeq
53\left(\frac{\omega}{10^{-14}\text{s}^{-1}}\right)^{2}\left(\frac{T}{10\text{K}}\right)^{-4}\left(\frac{M(t)}{1M_{\odot}}\right)^{3}\text{au}\quad,$
(5)

where $\omega$ 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
$\omega=9\times 10^{-15}$ s-1 and $T=15$ K, $R_{C}$ and never exceeds $10$ au
(orange line in Fig. 2). For a larger angular speed of, e.g.,
$\omega=3.1\times 10^{-14}$ s-1 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 $\omega=1\times 10^{-14}$ s-1 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}(t)=\frac{0.35}{\sqrt{M_{\star}(t)}}\,\text{au}\quad,$ (6)

where $M_{\star}$ is the mass of the proto-star in solar masses, $M_{\odot}$.
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 don’t constrain
ourselves to only exploring scenarios using Eq. 6. Thus, although we mainly
present results using that prescription from Morbidelli et al. (2022) we will
also examine the effects of using the more traditional “Shu recipe” (see
results in Sec. 5.3). In particular, we will show results where the $R_{C}$
grows to 10 and 100 au, respectively. 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 $\nu=\alpha H^{2}\Omega$ (or, equivalently
$\nu=\alpha c_{s}^{2}/\Omega)$, where $\Omega$ is the keplerian frequency and
$H=\sqrt{\frac{RTr^{3}}{\mu GM_{\star}M_{\odot}}}$ is the scale height, with
$R$ the gas constant, $\mu$ 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, $\alpha$, is a free parameter and varies in time and with radial
distance. As discussed in Sect. 2, it is reasonable to assume that $\alpha$
decays over time in a manner proportional to the disk infall function (two
examples are shown in Fig. 1). However, the initial value of $\alpha$ –
denoted $\alpha_{0}$ – is considered a free parameter. A minimum value of
$\alpha$ is set at $5\times 10^{-5}$, 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 $\alpha$ 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 $1,400$ K, 2) silicates with a sublimation
temperature of $1,000$ 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 1,060
K for $P=10^{-4}$ bar and C/O=1.0. For simplicity, we have kept the
sublimation temperature of silicates at $1,000$ 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, ${\text{S}t}$. 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 don’t just consider the limit
Stokes number given by the fragmentation barriers (3) and (4), where we assume
$\texttt{Sc}=0.1$ (Morbidelli et al., 2022), but also need to consider that
particles cannot be as large that they immediately drift out of a given cell.
This drift boundary is defined as

${\text{S}t}_{\text{drift}}=0.055\frac{\Sigma_{d}}{\Sigma_{g}}\frac{r\Omega}{\eta
v_{K}}\quad,$ (7)

where $r$ is the radial distance to the star. The barriers (4) and (7) are
additions to the model compared to the one published in Morbidelli et al.
(2022), which only considered (3).

The final particle size is determined through the minimum among
${\text{S}t}_{\text{growth}}$, given by the growth algorithm with timescale
$Z\Omega$, and ${\text{S}t}_{\text{frag}}$, ${\text{S}t}_{\text{ddf}}$ and
${\text{S}t}_{\text{drift}}$.

We have 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/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/upper cell to the current cell. Taking into account all four possible
loss/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 streaming instability in that ring, removing the dust in
excess (Gole et al., 2020).

##  4 Model setups and constraints

As described in the introduction, the underlying problem that prevents dust
from forming a large disk which extends far beyond the water snowline is that
it grows too fast. We will 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 ($\alpha$) is higher or the infall timescale
($T_{\text{infall}}$) is short. To explore the effect of these two parameters
of our model, we have varied them.

For the viscosity, we have one free parameter, the initial value of $\alpha$
at the beginning of the simulation, denoted $\alpha_{0}$. Once $\alpha_{0}$ is
set, it decreases as described in Sec. 3 proportional to the mass added to the
disk. Because the mass added to the disk decays over time, so will $\alpha$.
We have chosen to vary $\alpha_{0}$ 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 $\alpha$
can reach large values (see their Fig. 8). In particular, when the infalling
mass is on the same order as the disk mass, $\alpha$ reaches values of 0.1.
Such a mass ratio is reached early in our simulations. Therefore, we believe
such a high value of $\alpha_{0}$ is plausible for a brief period at the
beginning of the simulation. Remember that we let our $\alpha$ 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 have tested nine values of
$T_{\text{infall}}$ between $15$ kyr and $630$ kyr. A logarithmic spacing
between cases was used. In combination with the ten different $\alpha_{0}$, we
arrive at 90 simulations.

###  4.2 Fragmentation threshold of the dust

The second way to ensure particles reach larger 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 $v_{\text{frag}}=100$
cm/s for refractory and silicate particles and $v_{\text{frag}}=1,000$ cm/s
for icy particles beyond the water snowline.

Figure 3: The temperature-dependent fragmentation threshold,
$v_{\text{frag}}$, for icy particles according to Eq. 8

is shown. The experimental data from Musiolik & Wurm (2019) have been shifted
to lower temperatures to account for the lower sublimation temperature of
water in our model compared to the one in the experiment ($\sim 220$ K).

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

$v_{\text{frag}}(T)=v_{0}+v_{C}\Gamma(T)^{\frac{5}{6}}\quad,$ (8)

where $T$ is the temperature, $v_{0}=100$ cm/s, $v_{C}=1,600$ m/s, and

$\Gamma(T)=\Gamma_{C}+\Gamma_{d0}\tanh(\beta(T-T_{0}))\quad,$ (9)

where $\Gamma_{C}=\Gamma_{d0}=0.25$, $\beta=0.105$, and $T_{0}=150$. These
parameters ($v_{0}$, $v_{C}$, $\Gamma_{C}$, $\Gamma_{d0}$, and $T_{0}$) 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 $1,000$ cm/s to $100$ cm/s 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 larger distances from the star.

For locations in the disk above the sublimation temperature of $170$ K, i.e.,
for dry particles, we retain a fragmentation threshold of $100$ cm/s whereas
for locations with temperature below $170$ K we use Eq. 8. We have run two
sets of 90 simulations (the variations in $\alpha_{0}$ and
$T_{\text{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.

###  4.3 Summary of assumptions

To summarise, there are four main assumptions that we will explore in this
work:

  1. 1.

The centrifugal radius, $R_{C}$: either according to Eq. 5 (‘Shu recipe’)
growing to 10 and 100 au respectively or Eq. 6 where $R_{C}$ remains small.
Our nominal simulations are performed with Eq. 6.

  2. 2.

Variation of the initial disk viscosity, $\alpha_{0}$, between 0.01 and 0.1.

  3. 3.

Variation of the infall timescale $T_{\text{infall}}$ between $15$ kyr and
$630$ kyr.

  4. 4.

The fragmentation threshold for icy particles: either constant at 1,000 cm/s
(nominal case) or temperature-dependent according to Eq. 8.

##  5 Results

###  5.1 Temperature independent fragmentation threshold

Figure 4: The disk properties for two case with $T_{\text{infall}}=100$ kyr
but different initial viscosity (a with $\alpha_{0}=0.01$ and b with
$\alpha_{0}=0.1$) are shown. Both cases are shown at 1,000 years 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,
$R_{C}$ (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/cm2. 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 and 10 m/s for dry and icy particles, respectively.

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 snowline fragment at $100$
cm/s while those outside at $1,000$ cm/s.

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 snowline. Already very early on, e.g., after only
1,000 years, the dust particles just outside the snowline grow to the
centimetre scale and effectively stop their outward radial motion. This is
shown in panel a1 of Figure 4, which depicts the results of the case where we
have a small viscosity of $\alpha_{0}=0.01$ and $T_{\text{infall}}=100$ kyr
(nominal case in Morbidelli et al., 2022). Particles just outside of the water
snowline (dashed yellow line) have a size between 0.1 and 1 cm (Fig. 4a2) and
consequently have almost zero radial velocity (Fig. 4a3). Because the gas
continues to spread outwards, the dust and gas disks “decouple”, i.e., the
dust expansion lags the one of the gas. 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, e.g., $\alpha_{0}=0.1$
(Fig. 4b), the dust particles beyond the snowline are roughly an order of
magnitude smaller (Fig. 4b2) and thus retain a positive/outward motion (Fig.
4b3). The dust expansion keeps up with the gas expansion, and therefore the
two disks retain the same size (Fig. 4b1).

As expected, disks with larger viscosity expand faster. After 1,000 years of
expansion, the gas disk with $\alpha_{0}=0.01$ has expanded to roughly 4 au
(measured where the gas surface density is 1 g/cm2). In contrast, the disk
with $\alpha_{0}=0.1$ has reached 10 au and is, therefore, more than double
the size of the other (Fig. 4b).

Figure 5: The panes show the time needed for the gas disk to expand to 100 au,
$T_{\text{100 au}}$, as a function of the infall timescale,
$T_{\text{infall}}$. The size of the disk is measured at the location where
the gas surface density is 1 g/cm2. Each curve of different colour represents
a different value of the initial viscosity, $\alpha_{0}$. For each value of
$\alpha_{0}$ the orange star on the corresponding curve indicates the viscous
timescale, $T_{\text{visc}}$, of the disk, to be read on the horizontal axis.
$T_{\text{visc}}$ represents the average viscous timescale within 10 au at
$t=0$ for a disk with an aspect ratio of $6\%$. The left panel shows
$T_{\text{100 au}}$ while the right panel shows $T_{\text{100 au}}$ scaled to
$T_{\text{infall}}$

For a given $T_{\text{infall}}$, the time a disk takes to reach 100 au,
denoted $T_{\text{100 au}}$, decreases as the initial viscosity increases
(Fig. 5). To measure the size of the disk, we have used the location where the
gas surface density takes a value of 1 g/cm2. For the dust, we have adopted a
value 100 times smaller than the gas because of the metallicity of our
infalling material being $1\%$, i.e., a value of 0.01 g/cm2. 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, e.g., 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 $\alpha_{0}$, the orange star on the corresponding curve
indicates the viscous timescale, $T_{\text{visc}}$, of the disk, to be read on
the horizontal axis. $T_{\text{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_{\text{100 au}}\sim 400$ kyr (when $\alpha_{0}$
and $T_{\text{infall}}$ are minimal) and $T_{\text{100 au}}\sim 10$ kyr (when
$\alpha_{0}$ is maximal and $T_{\text{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_{\text{infall}}$ in these tests varies by more than one order of
magnitude, we might better measure $T_{\text{100 au}}$ in units of
$T_{\text{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 $\alpha_{0}$, the expansion time as a fraction of
$T_{\text{infall}}$ always decreased with increasing $T_{\text{infall}}$. It
is remarkable that if $T_{\text{infall}}=T_{\text{visc}}$, the value
$T_{\text{100 au}}/T_{\text{infall}}$ is independent of viscosity (i.e. the
orange stars fall on a horizontal line).

Figure 6: The maximum dust mass 1 au beyond the snowline (SL) is shown as a
function of the size of the dust disk, $R_{\text{0.01 g/cc}}$, at the time
when the disk has reached that maximum mass. The top panels show the
dependency on the viscosity, $\alpha_{0}$, while the bottom row the dependency
on the infall timescale, $T_{\text{infall}}$. Panels a shows the results for
the temperature-independent fragmentation threshold, while the b panels show
the results for the temperature-dependent fragmentation threshold.

####  5.1.1 Mass and size of the dust disk

We have measured the maximum dust mass a given disk holds 1 au beyond the
snowline. To make sure the measurement was not contaminated by the dynamics
around the snowline, we chose to exclude the dust mass just outside the
snowline. We will 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 snowline (Fig. 6a1). 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_{\text{infall}}$,
is, the more massive the outer disk is (Fig. 6a2). Therefore, short
$T_{\text{infall}}$ and large $\alpha_{0}$ produced the largest and most
massive outer disks. These disks thus satisfy our first criteria for good
protoplanetary disks of the Solar System.

Figure 7: The panels show the mass of planetesimals formed at different
locations in the disk. The left panel shows the total mass of “rocky”
planetesimals as a function of the infall timescale, $T_{\text{infall}}$, and
the viscosity, $\alpha_{0}$. Rocky planetesimals are the ones that form around
the silicate condensation line. White areas are disks that don’t produce any
planetesimals, while grey squares indicate disks that produce between zero and
two Earth masses of planetesimals. The centre panel show the total mass of
“icy” planetesimals for a given disk. Icy planetesimals are the ones formed
around the water snowline, typically outside of it. The right panel shows the
number of locations, i.e., rings, where planetesimals form.

####  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 have split the results into “rocky” planetesimals
(forming at the silicate condensation line) and “icy” planetesimals forming
at/outside of the water snowline.

First, we observe that for most cases with $T_{\text{infall}}>100$ kyr, no
“rocky” planetesimals are formed. Second, for “rocky” planetesimals, there is
an optimal viscosity given a $T_{\text{infall}}$. This is most clearly visible
for $T_{\text{infall}}=39$ kys (the third line from the bottom). For this
infall timescale, the optimum viscosity to produce “rocky” planetesimals is
$\alpha_{0}=0.05$. The planetesimal mass decreases for higher and lower values
of $\alpha_{0}$. 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 snowline (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 will 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.

Figure 8: The panels show the ratio between the surface density of high-
temperature condensates, $\Sigma_{\text{pot.CAI}}$, and the surface density of
all refractory particles, $\Sigma_{\text{tot.ref.}}$ for different values of
$T_{\text{infall}}$ and $\alpha_{0}=0.05$. These simulations assumed the
nominal fragmentation thresholds of 1 m/s and 10 m/s 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 $R_{C}$, and the fine dashed line is the dust
disk size, $R_{\text{0.01 g/cc}}$.

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 time scale. 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 $\sim 200$ ky for $\alpha_{0}=0.05$) the disk is rather cold
from the beginning, and therefore the refractory condensation line (defined as
$T=1,400$ 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 $\sim 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_{\text{infall}}$.
While we show these results for $\alpha_{0}=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_{\text{infall}}\sim 150$ kyr for $\alpha_{0}=0.1$). Conversely, this
transition happens at larger infall timescales when the viscosity is smaller
(e.g., at $T_{\text{infall}}>400$ kyr for $\alpha_{0}=0.01$). But in all
cases, neither very short nor long $T_{\text{infall}}$ are favoured for the
transport of CAIs to the outer disk.

Figure 9: The panels show the ratio between the surface density of high-
temperature condensates, $\Sigma_{\text{pot.CAI}}$, and the surface density of
all refractory particles, $\Sigma_{\text{tot.ref.}}$ for different values of
$\alpha_{0}$ and $T_{\text{infall}}=100$ kyr. These simulations assumed the
nominal fragmentation thresholds of 1 m/s and 10 m/s 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 $R_{C}$, and the fine dashed line is the dust
disk size, $R_{\text{0.01 g/cc}}$.

The smaller the initial viscosity is, the larger the fraction of the disk that
is populated by CAIs. For example, when $\alpha_{0}<0.05$ for
$T_{\text{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 have kept the Schmidt number at
$\texttt{Sc}=0.1$. A higher Schmidt number of, e.g., $\texttt{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.

###  5.2 Temperature dependent fragmentation threshold

In the case where we impose the temperature-dependent fragmentation threshold
beyond the snowline (see Sec. 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.

Figure 10: The same panels are shown as in Fig. 7 but for the temperature-
dependent fragmentation threshold.

This shift of dust mass from the inner to the outer disk has clear
consequences. We now have significantly more “icy” planetesimals than “rocky”
ones (Fig. 10). For some combination of parameters $\alpha_{0}$ and
$T_{\text{infall}}$ (e.g., $0.07\leq\alpha_{0}\leq 0.1$ and $40\text{ kyr}\leq
T_{\text{infall}}\leq 100$ kyr), a couple of Earth masses of “rocky”
planetesimals 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_{\text{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 larger distances.

###  5.3 Shu infall

Figure 11: The radial velocity of the gas at $100$ kyr is shown as a function
of the radial distance to the star for the case when $T_{\text{infall}}=100$
kry and $\alpha_{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, $R_{C}$, of the respective case.

Because our prescription of the infall is somewhat unconventional, i.e., the
centrifugal radius, $R_{C}\sim 0.35$ au (Eq 6), we have 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 $\omega=9\times 10^{-15}$
s-1 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 larger 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 $R_{C}$ is an accretion disk, i.e., the radial gas
velocity, $v_{r,g}$, 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 ($v_{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 have tested two different angular velocities, $\omega$, of the molecular
cloud. Once with $\omega=10^{-14}$ s-1 resulting in a maximum $R_{C}$ of
roughly $10$ au as shown in Fig. 11 and once with $\omega=3.1\times 10^{-14}$
s-1 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_{\text{infall}}$ and $\alpha_{0}$ remain the same as
above.

Figure 12: The maximum dust mass 1 au beyond the snowline (SL) is shown as a
function of the size of the dust disk, $R_{\text{0.01 g/cc}}$, at the time
when the disk has reached that maximum mass. Results for the “Shu-type” infall
where $R_{C}$ grows to roughly $10$ au are shown in panels a, while the ones
where $R_{C}$ grows to roughly $100$ au are shown in panels b. 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 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. 12a1). 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 $\alpha_{0}$. The sizes and masses are
also confined to 80-150 au and 30-300 M⊕ (Fig. 12b1).

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. 12a2). For the
temperature-dependent fragmentation threshold, the disk sizes and masses are
only weakly dependent on $\alpha_{0}$ and $T_{\text{infall}}$. These disks are
between 150 and 400 au and have masses between 300 and 700 M⊕ (Fig. 12b2), and
therefore very massive and large.

Figure 13: The same panels are shown as in Fig. 7 but for the temperature-
dependent fragmentation threshold and the “Shu-infall” model with a maximum
$R_{C}$ of roughly 10 au (top panels) and 100 au (lower panels) respectively.
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⊕.

When we prescribe the “Shu-infall” particles in the inner disk (within the
water snowline) 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 snowline, we observe only sparse formation of
planetesimals (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 much 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 to that of hydrogen. This supports
planetesimal formation.

Figure 14: Show the ratio between the surface density of high-temperature
condensates, $\Sigma_{\text{pot.CAI}}$, and the surface density of all
refractory particles, $\Sigma_{\text{tot.ref.}}$. The case on the left assumes
$\alpha_{0}=0.05$, $T_{\text{infall}}=100$ kyr, and the case on the right
assumes $\alpha_{0}=0.01$, $T_{\text{infall}}=15$ kyr. In both cases, we’ve
used the temperature-dependent fragmentation threshold and the $R_{C}$ 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 white dashed line is $R_{C}$.

Finally, we have 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
$\alpha_{0}=0.05$, $T_{\text{infall}}=100$ kyr, the temperature-dependent
fragmentation threshold, and $R_{C}$ growing to roughly 10 au but is
representative of almost all combinations of $\alpha_{0}$ and
$T_{\text{infall}}$. The only exception is for $\alpha_{0}=0.01$ and
$T_{\text{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.

##  6 Discussion and conclusion

Infall of material into protoplanetary disks occurs more or less close to the
star – typically much less 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 down 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 e.g., 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 proto-star. Outer shells, with larger angular momentum, will fall
at larger distances and in a later stage in disk formation (Shu, 1977). In
such scenarios, $R_{C}$ grows with time and we refer to them as “Shu-type”
infall models. Contrary to this, magnetic braking can remove angular momentum
from the infalling material. This can cause the material to fall close to the
star irrespectively of the initial angular momentum of the material.

In the introduction we have described that a disk formation and evolution
scenario for the Solar System must satisfy at least the following three
requirements:

  1. 1.

it must develop an extended disk of gas and dust (up 45 au for dust);

  2. 2.

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

  3. 3.

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

We found that scenarios using 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 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, the infall timescale is of
order 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. 1.

The larger is the initial viscosity, $\alpha_{0}$, the larger is the outer
dust disk.

  2. 2.

The shorter is the infall timescale, $T_{\text{infall}}$, the more massive is
the outer dust disk.

  3. 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. 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. 5.

No “rocky” and very few “icy” planetesimals form when $T_{\text{infall}}>100$
kyr.

  6. 6.

The largest mass of “icy” planetesimals forms when $\alpha_{0}>0.05$.

  7. 7.

There is an optimum $\alpha_{0}$ that maximises the mass of “rocky”
planetesimals. For example, for $T_{\text{infall}}=39$ kyr it is
$\alpha_{0}=0.05$.

  8. 8.

The temperature-dependent fragmentation threshold results in more “icy” than
“rocky” planetesimals (roughly by a factor of 10) than in 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 the three requirements we
had put forth at the beginning, there are two additional related requirements
that will need to be met eventually but cannot at this point. Observations
show that protoplanetary disks are long-lived, i.e., $3-4$ million years
(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 for a good 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 magneto-hydrodynamic (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.

## Acknowledgments

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

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192, doi: 10.3847/1538-4357/ab40a3
  * Amelin et al. (2010) Amelin, Y., Kaltenbach, A., Iizuka, T., et al. 2010, Earth and Planetary Science Letters, 300, 343, doi: 10.1016/j.epsl.2010.10.015
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483, doi: 10.1146/annurev-astro-031220-010302
  * Andrews et al. (2018) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157, doi: 10.3847/1538-4357/aadd9f
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Barenfeld et al. (2017) Barenfeld, S. A., Carpenter, J. M., Sargent, A. I., Isella, A., & Ricci, L. 2017, ApJ, 851, 85, doi: 10.3847/1538-4357/aa989d
  * Berthier et al. (2020) Berthier, J., Descamps, P., Vachier, F., et al. 2020, Icarus, 352, 113990, doi: 10.1016/j.icarus.2020.113990
  * Béthune et al. (2016) Béthune, W., Lesur, G., & Ferreira, J. 2016, A&A, 589, A87, doi: 10.1051/0004-6361/201527874
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148, doi: 10.1051/0004-6361/201118136
  * Blum & Wurm (2008) Blum, J., & Wurm, G. 2008, ARA&A, 46, 21, doi: 10.1146/annurev.astro.46.060407.145152
  * Bouvier & Wadhwa (2010) Bouvier, A., & Wadhwa, M. 2010, Nature Geoscience, 3, 637, doi: 10.1038/ngeo941
  * Brownlee et al. (2006) Brownlee, D., Tsou, P., Aléon, J., et al. 2006, Science, 314, 1711, doi: 10.1126/science.1135840
  * Capria et al. (2000) Capria, M. T., Coradini, A., De Sanctis, M. C., & Orosei, R. 2000, A&A, 357, 359
  * Chambers (2001) Chambers, J. E. 2001, Icarus, 152, 205, doi: 10.1006/icar.2001.6639
  * Charnoz et al. (2019) Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019, A&A, 627, A50, doi: 10.1051/0004-6361/201833216
  * Chaussidon & Liu (2015) Chaussidon, M., & Liu, M.-C. 2015, Washington DC American Geophysical Union Geophysical Monograph Series, 212, 1, doi: 10.1002/9781118860359.ch1
  * Connelly et al. (2012) Connelly, J. N., Bizzarro, M., Krot, A. N., et al. 2012, Science, 338, 651, doi: 10.1126/science.1226919
  * Desch et al. (2022) Desch, S. J., Dunlap, D. R., Dunham, E. T., Williams, C. D., & Mane, P. 2022, arXiv e-prints, arXiv:2212.00390, doi: 10.48550/arXiv.2212.00390
  * Dominik & Tielens (1997) Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647, doi: 10.1086/303996
  * Drążkowska & Alibert (2017) Drążkowska, J., & Alibert, Y. 2017, A&A, 608, A92, doi: 10.1051/0004-6361/201731491
  * Drążkowska & Dullemond (2018) Drążkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62, doi: 10.1051/0004-6361/201732221
  * Eberhardt et al. (1987) Eberhardt, P., Krankowsky, D., Schulte, W., et al. 1987, A&A, 187, 481
  * Gasc et al. (2017) Gasc, S., Altwegg, K., Balsiger, H., et al. 2017, MNRAS, 469, S108, doi: 10.1093/mnras/stx1412
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132, doi: 10.3847/1538-4357/abc334
  * Groussin et al. (2019) Groussin, O., Attree, N., Brouet, Y., et al. 2019, Space Science Reviews, 215, 29, doi: 10.1007/s11214-019-0594-x
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56, doi: 10.1051/0004-6361/200912852
  * Hendler et al. (2020) Hendler, N., Pascucci, I., Pinilla, P., et al. 2020, ApJ, 895, 126, doi: 10.3847/1538-4357/ab70ba
  * Hueso & Guillot (2005) Hueso, R., & Guillot, T. 2005, A&A, 442, 703, doi: 10.1051/0004-6361:20041905
  * Hyodo et al. (2021) Hyodo, R., Guillot, T., Ida, S., Okuzumi, S., & Youdin, A. N. 2021, A&A, 646, A14, doi: 10.1051/0004-6361/202039894
  * Hyodo et al. (2019) Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90, doi: 10.1051/0004-6361/201935935
  * Ida & Guillot (2016) Ida, S., & Guillot, T. 2016, A&A, 596, L3, doi: 10.1051/0004-6361/201629680
  * Jacobsen et al. (2008) Jacobsen, B., Yin, Q.-z., Moynier, F., et al. 2008, Earth and Planetary Science Letters, 272, 353, doi: 10.1016/j.epsl.2008.05.003
  * Jacquet et al. (2011) Jacquet, E., Fromang, S., & Gounelle, M. 2011, A&A, 526, L8, doi: 10.1051/0004-6361/201016118
  * Johansen et al. (2014) Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 547, doi: 10.2458/azu_uapress_9780816531240-ch024
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022, doi: 10.1038/nature06086
  * Kelley et al. (2022) Kelley, M. S. P., Kokotanekova, R., Holt, C. E., et al. 2022, ApJ, 933, L44, doi: 10.3847/2041-8213/ac7bec
  * Klahr & Schreiber (2020) Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54, doi: 10.3847/1538-4357/abac58
  * Kruijer et al. (2017) Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proceedings of the National Academy of Science, 114, 6712, doi: 10.1073/pnas.1704461114
  * Kumar & Coleman (1993) Kumar, S., & Coleman, C. S. 1993, MNRAS, 260, 323, doi: 10.1093/mnras/260.2.323
  * Kuznetsova et al. (2022) Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, ApJ, 928, 92, doi: 10.3847/1538-4357/ac54a8
  * Larson (1969) Larson, R. B. 1969, MNRAS, 145, 271, doi: 10.1093/mnras/145.3.271
  * Lee et al. (2021) Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021, A&A, 648, A101, doi: 10.1051/0004-6361/202038105
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Machida & Basu (2019) Machida, M. N., & Basu, S. 2019, ApJ, 876, 149, doi: 10.3847/1538-4357/ab18a7
  * Machida et al. (2014) Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2014, MNRAS, 438, 2278, doi: 10.1093/mnras/stt2343
  * Machida & Matsumoto (2011) Machida, M. N., & Matsumoto, T. 2011, MNRAS, 413, 2767, doi: 10.1111/j.1365-2966.2011.18349.x
  * Marrocchi et al. (2019) Marrocchi, Y., Villeneuve, J., Jacquet, E., Piralla, M., & Chaussidon, M. 2019, Proceedings of the National Academy of Science, 116, 23461, doi: 10.1073/pnas.1912479116
  * Masson et al. (2016) Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32, doi: 10.1051/0004-6361/201526371
  * Mattsson (2020) Mattsson, L. 2020, MNRAS, 499, 6035, doi: 10.1093/mnras/staa3114
  * Morbidelli et al. (2022) Morbidelli, A., Baillié, K., Batygin, K., et al. 2022, Nature Astronomy, 6, 72, doi: 10.1038/s41550-021-01517-7
  * Morbidelli et al. (2020) Morbidelli, A., Libourel, G., Palme, H., Jacobson, S. A., & Rubie, D. C. 2020, Earth and Planetary Science Letters, 538, 116220, doi: 10.1016/j.epsl.2020.116220
  * Morse et al. (2015) Morse, A., Mousis, O., Sheridan, S., et al. 2015, A&A, 583, A42, doi: 10.1051/0004-6361/201526624
  * Musiolik & Wurm (2019) Musiolik, G., & Wurm, G. 2019, ApJ, 873, 58, doi: 10.3847/1538-4357/ab0428
  * Nakashima et al. (2015) Nakashima, D., Ushikubo, T., Kita, N. T., et al. 2015, Earth and Planetary Science Letters, 410, 54, doi: 10.1016/j.epsl.2014.11.020
  * Nesvorný (2018) Nesvorný, D. 2018, ARA&A, 56, 137, doi: 10.1146/annurev-astro-081817-052028
  * Nesvorný et al. (2021) Nesvorný, D., Li, R., Simon, J. B., et al. 2021, The Planetary Science Journal, 2, 27, doi: 10.3847/PSJ/abd858
  * Nesvorný et al. (2017) Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, ApJ, 845, 27, doi: 10.3847/1538-4357/aa7cf6
  * Nesvorný et al. (2022) Nesvorný, D., Vokrouhlický, D., & Fraser, W. C. 2022, AJ, 163, 137, doi: 10.3847/1538-3881/ac4bc9
  * Neumann et al. (2018) Neumann, W., Kruijer, T. S., Breuer, D., & Kleine, T. 2018, Journal of Geophysical Research (Planets), 123, 421, doi: 10.1002/2017JE005411
  * Nimmo et al. (2018) Nimmo, F., Kretke, K., Ida, S., Matsumura, S., & Kleine, T. 2018, Space Sci. Rev., 214, 101, doi: 10.1007/s11214-018-0533-2
  * Paruta et al. (2016) Paruta, P., Hendrix, T., & Keppens, R. 2016, Astronomy and Computing, 16, 155, doi: 10.1016/j.ascom.2016.05.002
  * Pignatale et al. (2018) Pignatale, F. C., Charnoz, S., Chaussidon, M., & Jacquet, E. 2018, ApJ, 867, L23, doi: 10.3847/2041-8213/aaeb22
  * Pillich et al. (2021) Pillich, C., Bogdan, T., Landers, J., Wurm, G., & Wende, H. 2021, A&A, 652, A106, doi: 10.1051/0004-6361/202140601
  * Pineda et al. (2020) Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nature Astronomy, 4, 1158, doi: 10.1038/s41550-020-1150-z
  * Piralla et al. (2023) Piralla, M., Villeneuve, J., Schnuriger, N., Bekaert, D. V., & Marrocchi, Y. 2023, Icarus, 394, 115427, doi: 10.1016/j.icarus.2023.115427
  * Polak & Klahr (2023) Polak, B., & Klahr, H. 2023, ApJ, 943, 125, doi: 10.3847/1538-4357/aca58f
  * Preusker et al. (2017) Preusker, F., Scholten, F., Matz, K. D., et al. 2017, A&A, 607, L1, doi: 10.1051/0004-6361/201731798
  * Riols et al. (2020) Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95, doi: 10.1051/0004-6361/201937418
  * Rodenkirch & Dullemond (2022) Rodenkirch, P. J., & Dullemond, C. P. 2022, A&A, 659, A42, doi: 10.1051/0004-6361/202142571
  * Ruíz-Rodríguez et al. (2018) Ruíz-Rodríguez, D., Cieza, L. A., Williams, J. P., et al. 2018, MNRAS, 478, 3674, doi: 10.1093/mnras/sty1351
  * Saito & Sirono (2011) Saito, E., & Sirono, S.-i. 2011, ApJ, 728, 20, doi: 10.1088/0004-637X/728/1/20
  * Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69, doi: 10.1051/0004-6361/201629561
  * Scott & Krot (2003) Scott, E. R. D., & Krot, A. N. 2003, Treatise on Geochemistry, 1, 711, doi: 10.1016/B0-08-043751-6/01145-2
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488, doi: 10.1086/155274
  * Shu et al. (2001) Shu, F. H., Shang, H., Gounelle, M., Glassgold, A. E., & Lee, T. 2001, ApJ, 548, 1029, doi: 10.1086/319018
  * Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55, doi: 10.3847/0004-637X/822/1/55
  * Simon et al. (2017) Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12, doi: 10.3847/2041-8213/aa8c79
  * Spencer et al. (2020) Spencer, J. R., Stern, S. A., Moore, J. M., et al. 2020, Science, 367, aay3999, doi: 10.1126/science.aay3999
  * Spitzer et al. (2021) Spitzer, F., Burkhardt, C., Nimmo, F., & Kleine, T. 2021, Earth and Planetary Science Letters, 576, 117211, doi: 10.1016/j.epsl.2021.117211
  * Stammler & Birnstiel (2022) Stammler, S. M., & Birnstiel, T. 2022, ApJ, 935, 35, doi: 10.3847/1538-4357/ac7d58
  * Takeuchi & Lin (2002) Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344, doi: 10.1086/344437
  * Takeuchi & Lin (2005) —. 2005, ApJ, 623, 482, doi: 10.1086/428378
  * Teiser & Wurm (2009) Teiser, J., & Wurm, G. 2009, MNRAS, 393, 1584, doi: 10.1111/j.1365-2966.2008.14289.x
  * Tobin et al. (2010) Tobin, J. J., Hartmann, L., & Loinard, L. 2010, ApJ, 722, L12, doi: 10.1088/2041-8205/722/1/L12
  * Tripathi et al. (2017) Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44, doi: 10.3847/1538-4357/aa7c62
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399, doi: 10.1046/j.1365-8711.1998.01118.x
  * Vaytet et al. (2018) Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5, doi: 10.1051/0004-6361/201732075
  * Wada et al. (2007) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320, doi: 10.1086/514332
  * Wahlberg Jansson & Johansen (2014) Wahlberg Jansson, K., & Johansen, A. 2014, A&A, 570, A47, doi: 10.1051/0004-6361/201424369
  * Wahlberg Jansson & Johansen (2017) —. 2017, MNRAS, 469, S149, doi: 10.1093/mnras/stx1470
  * Walsh et al. (2011) Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, Nature, 475, 206, doi: 10.1038/nature10201
  * Warren (2011) Warren, P. H. 2011, Geochim. Cosmochim. Acta, 75, 6912, doi: 10.1016/j.gca.2011.09.011
  * Wurster et al. (2021) Wurster, J., Bate, M. R., & Bonnell, I. A. 2021, MNRAS, 507, 2354, doi: 10.1093/mnras/stab2296
  * Wurster et al. (2022) Wurster, J., Bate, M. R., Price, D. J., & Bonnell, I. A. 2022, MNRAS, 511, 746, doi: 10.1093/mnras/stac123
  * Yang et al. (2017) Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80, doi: 10.1051/0004-6361/201630106
  * Yang & Ciesla (2012) Yang, L., & Ciesla, F. J. 2012, , 47, 99, doi: 10.1111/j.1945-5100.2011.01315.x
  * Yen et al. (2019) Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, ApJ, 880, 69, doi: 10.3847/1538-4357/ab29f8
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459, doi: 10.1086/426895
  * Zolensky et al. (2006) Zolensky, M. E., Zega, T. J., Yano, H., et al. 2006, Science, 314, 1735, doi: 10.1126/science.1135842