# Anisotropic Infall and Substructure formation in Embedded Disks

Aleksandra Kuznetsova NHFP Sagan Postdoctoral Fellow  American Museum of
Natural History, 200 Central Park West, New York, New York, 10024, USA  Jaehan
Bae Department of Astronomy, University of Florida, 211 Bryant Space Science
Center, Gainesville, FL 32611-2055 USA  Lee Hartmann Department of Astronomy,
University of Michigan, 1085 S. University Ave., Ann Arbor, MI 48109
Mordecai-Mark Mac Low American Museum of Natural History, 200 Central Park
West, New York, New York, 10024, USA

###### Abstract

The filamentary nature of accretion streams found around embedded sources
suggest that protostellar disks experience heterogenous infall from the star-
forming environment, consistent with the accretion behavior onto star-forming
cores in top-down star-cluster formation simulations. This may produce disk
substructures in the form of rings, gaps, and spirals continuing to be
identified by high-resolution imaging surveys in both embedded Class 0/I and
later Class II sources. We present a parameter study of anisotropic infall,
informed by the properties of accretion flows onto protostellar cores in
numerical simulations, and varying the relative specific angular momentum of
incoming flows as well as their flow geometry. Our results show that
anisotropic infall perturbs the disk and readily launches the Rossby wave
instability (RWI). It forms vortices at the inner and outer edge of the infall
zone where material is deposited. These vortices drive spiral waves and
angular momentum transport, with some models able to drive stresses on the
order of $\alpha\sim 10^{-2}$. The resulting azimuthal shear forms robust
pressure bumps that act as barriers to radial drift of dust grains, as
demonstrated by post-processing calculations of drift-dominated dust
evolution. We discuss how a self-consistent model of anisotropic infall can
account for the formation of millimeter rings in the outer disk as well as
producing compact dust disks, consistent with observations of embedded
sources.

hydrodynamics — instabilities — planetary systems: formation — planetary
systems: protoplanetary disks

††journal: ApJ

##  1 Introduction

High-resolution imaging studies of protoplanetary disks (canonically Class II
sources) in the continuum have found that many of the brightest sources are
not well described by a smooth emission profile. For instance, the objects
mapped in the Disk Substructures at High Angular Resolution Project (DSHARP)
survey are rife with substructure in the form of rings, gaps, and spirals
(Andrews et al., 2018). Several studies have linked the location of ringed
structures with that of potential protoplanets and the indirect kinematic
effects of such planets (e.g. Teague et al., 2018; Pinte et al., 2019),
motivating the search for an assumed population of embedded planets (Bae et
al., 2018; Disk Dynamics Collaboration et al., 2020). While direct detection
of protoplanet candidates is still in its infancy, cases like the discovery
and characterization of embedded protoplanet candidates in PDS 70 (Keppler et
al., 2018; Benisty et al., 2021), lend credence to a causal relationship
between planets and disk substructure.

Studies of planet-disk interaction have initially assumed that gaps were
predominantly opened by relatively high (Jupiter) mass planets (Malik et al.,
2015), but recent work has shown that planets of even lower (Super-Earth)
masses can not only efficiently dynamically perturb disks to create rings and
gaps, but that a single planet can create the appearance of multiple
concentric rings within the disk (Dong et al., 2017; Bae et al., 2017).
However, dust evolution models have shown that ringed structures observed in
the millimeter continuum can be produced by any mechanism that results in the
trapping of grains in pressure bumps – regions in the disk where the pressure
gradient is perturbed locally, eventually arresting the inward radial drift of
dust (Pinilla et al., 2012). The concentration of dust in pressure bumps may
trigger planetesimal formation; for example, through the streaming instability
(Carrera et al., 2021), such that the coincidence of putative embedded
protoplanets and disk substructures could also be explained by rings serving
as sites of planet formation.

Recently, the VLA/ALMA Nascent Disk and Multiplicity (VANDAM) survey targeting
earlier phase (Class 0/I) protostellar disks, still embedded within their
infalling envelopes, has also identified sources with continuum substructure
(Sheehan et al., 2020). In some younger disks, like in the case of HL Tau
which has prominent concentric ring structures (ALMA Partnership et al.,
2015), molecular line emission observations have associated flows of dense
infalling gas with the presence of rings and gaps (Yen et al., 2019). Direct
observations and characterization of infalling streams of filamentary gas
originating outside of the dense protostellar core (Pineda et al., 2020) have
helped _streamers_ newly enter the vocabulary, although non-isotropy in the
envelopes of early phase protostellar objects had been observed for some time
before (Tobin et al., 2010).

Given that embedded objects are regularly experiencing infall of nascent cloud
material and that such infall is likely to be filamentary and anisotropic,
infall-induced pressure perturbations could be a reliable mechanism for
seeding substructure at early disk phases. Numerical simulations have already
shown that the infall of material onto existing disks has important dynamical
consequences. Studies have found that infalling material landing on disks
could effectively launch the Rossby wave instability (RWI), forming dust-
trapping vortices and spiral arms (Bae et al., 2015), aid in disk
fragmentation by inducing gravitational instability (GI) (Kratter et al.,
2008), drive high rates of angular momentum transport (Lesur et al., 2015),
and potentially power accretion outbursts at early times (Vorobyov & Basu,
2010; Bae et al., 2014).

Numerical simulations of star cluster formation that follow star formation
starting from the nascent molecular cloud have found that anisotropic infall
is the norm rather than the exception (e.g. Smith et al., 2011; Kuffmeier et
al., 2017; Kuznetsova et al., 2019)—a natural consequence of the inherent
heterogeneity within filamentary star-forming environments. In Kuznetsova et
al. (2020, K20 hereafter) we were able to statistically characterize the
accretion behavior for all centers of gravitational collapse in our
simulations by measuring the flux of mass and angular momentum into patches of
cells surrounding sink particles.

We found that the accretion of cloud material onto these protostellar systems
was not continuous, instead proceeding in several discrete episodes between
periods of quiescence. For any given core, accretion events, often associated
with spatially distinct reservoirs, have uncorrelated specific angular
momenta, variable in direction and magnitude.

The results of these and similar numerical studies call into question the
traditional analytic descriptions of disk formation and evolution through
self-similar uniformly rotating collapse (i.e. Terebey et al., 1984).

Analytic descriptions of protostellar infall typically model the deposition of
mass and momentum onto disks by assuming parabolic trajectories of fluid
parcels from an isotropic, uniformly rotating, spherical envelope (see:
Ulrich, 1976; Cassen & Moosman, 1981).

In the self-similar framework, the specific angular momentum of infalling
material is monotonically increasing with time as outer layers of the cloud
fall in, gradually building the disk outward. In this paper, we present a
modified infall model, informed by the results of K20, that takes into account
stochastic variations in the angular momentum of material accreted into
envelopes, as well as the filamentary and episodic nature of envelope infall.
This anisotropic model is used here in a parameter study across infall
properties, investigating the relative efficiency of transport and structure
formation. We demonstrate how the evolution of the infall-induced RWI affects
structure formation as a function of the infall properties, particularly
examining the effects on disk transport (§3.2) and dust trapping (§3.4). We
show that our model of anisotropic infall results in robust structure
formation and predicts a variety of ring and gap architectures, including a
significant likelihood for compact dust disks.

##  2 Model Description

###  2.1 FARGO3D

The simulations in this paper are performed on a two-dimensional cylindrical
grid ($r$-$\phi$) with the _FARGO3D_ (GPU) code (Benítez-Llambay & Masset,
2016). _FARGO3D_ solves the basic hydrodynamic equations with user-defined
source terms, $\dot{\Sigma}_{\rm in}$ and $\mathbf{F}_{\rm in}$, to represent
the injection of infalling mass and momentum:

$\frac{\partial\Sigma}{\partial
t}+\nabla\cdot(\Sigma\mathbf{u})=\dot{\Sigma}_{\rm in},$ (1)
$\Sigma\left(\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}\right)=-\nabla
P-\Sigma\nabla\Phi+\nabla\cdot\Pi+\mathbf{F}_{\rm in},$ (2)

where $P=\Sigma c_{\rm s}^{2}$ is the isothermal pressure, $\Phi$ is the
central star’s gravitational potential, and $\Pi$ is the viscous stress
tensor.

###  2.2 Disk Setup

The initial disk for every run has a $\Sigma\propto{R}^{-1}$ power-law surface
density profile, with an exponential cutoff at the nominal disk radius,
$R_{\rm d}=60\ \mathrm{au}$, described by

$\Sigma(R)=\Sigma_{0}{\left(\frac{R}{R_{0}}\right)}^{-1}e^{-R/R_{\rm d}},$ (3)

where the inner radius, $R_{0}=1\ \mathrm{au}$, is set to $\Sigma_{0}=1700\
\mathrm{g\ {cm}^{-2}}$ for an initial disk mass of $M_{d}=0.07M_{\odot}$.

The radial domain starts at $R_{0}$ and extends out to $R_{\rm max}=260\
\mathrm{au}$. As such, we do not directly study accretion onto the central
star. At the inner boundary, we use the wave-killing boundary prescription (De
Val-Borro et al., 2006) to avoid wave reflection, while we allow outflow from
the outer boundary. We adopt a temperature scaling consistent with heating by
stellar irradiation such that $T\propto{R}^{-1/2}$, corresponding to a disk
aspect ratio $H/R=0.046\ (R/R_{0})^{1/4}$, which sets the sound speed $c_{\rm
s}$ at each location in the radial grid The isothermal equation of state
neglects local heating through shocks or other internal heating processes.

The hydrodynamic simulations are performed on a cylindrical grid of dimension
$N_{r}=2048,N_{\theta}=1024$. For the purposes of capturing the development of
local pressure gradients in a global simulation, we use a logarithmic radial
grid corresponding to $\Delta R/R=0.0027$. In numerical studies of RWI
vortices in protoplanetary disks, the radial half width of vortex structures
depends on the local scale height and varies from $\Delta r_{\rm v}=0.5$–$2H$
(Ono et al., 2018). With the chosen $\Delta R/R$, there are at minimum 16
cells across the smallest vortices.

We use the $\alpha$-disk prescription in _FARGO3D_ to set the viscosity
$\nu=\alpha c_{\rm s}H$. To ensure that any viscous transport occurs as a
result of infall-induced mechanical stresses, the background level of $\alpha$
is set to $\alpha=10^{-7}$.

###  2.3 Modeling Infall

We calculate the mass and momentum deposition of infalling material by
adapting the analytic models in Ulrich (1976, UCM hereafter); Cassen & Moosman
(1981, UCM hereafter), which calculate the flux of material landing on the
disk based on parabolic trajectories of fluid parcels. The fluid parcels are
assumed to inherit the specific angular momentum of the infalling cloud
$j_{c}$ which sets the maximum landing radius for material, the centrifugal
radius

$R_{\rm c}=\frac{j_{\rm c}^{2}}{GM_{*}}.$ (4)

In the UCM model, fluid parcels travel along streamlines in $(r,\theta)$
starting from a range of incidence angles $\theta_{0}=[0,\pi]$ until landing
on the disk surface where $\theta=\pi$. The landing radius of a specific
streamline with an incidence angle $\theta_{0}$ scales as $R=R_{\rm
c}\sin^{2}\theta_{0}$. Fully spherical infall will deposit material from
$R=[0,R_{\rm c}]$. The instantaneous mass infall rate of the UCM models is
then:

$\dot{\Sigma}_{\mathrm{UCM}}(R,t)=\frac{\dot{M}_{\rm in}}{4\pi RR_{\rm
c}}\left(1-\frac{R}{R_{\rm c}}\right)^{-1/2},$ (5)

where $\dot{M}_{\rm in}$ is the total mass infall rate of the cloud.

We modify the UCM models to allow for filamentary infall by introducing the
parameter $R_{\rm in}$ to denote the inner boundary of the infall zone, which
has width $\Delta R=R_{\rm c}-R_{\rm in}$. For the purpose of our parameter
study, the modified equations are scaled such that the total integrated mass
infall rate is always set to the total cloud infall rate $\dot{M}_{\rm
in}=1.0\times 10^{-6}\ \mathrm{M_{\odot}\ {yr}^{-1}}$, regardless of how
filamentary the accretion is. The modified infall rate

$\begin{split}\dot{\Sigma}_{\rm
in}(R,t)&=\dot{\Sigma}_{\mathrm{UCM}}(R,t)\left(1-\frac{R_{\rm in}}{R_{\rm
c}}\right)^{-1/2}\\\ &=\dot{\Sigma}_{\mathrm{UCM}}(R,t)\left(\frac{\Delta
R}{R_{\rm c}}\right)^{-1/2},\\\ \end{split}$ (6)

is the surface density source term in equation 1. The relative width of the
annulus of deposited material is set by the fraction of solid angle subtended
by the flow streamlines, $f$, such that

$f^{2}=\frac{\Delta R}{R_{\rm c}}.$ (7)

Momentum conservation when material with velocity $\mathbf{u}_{\rm in}$
impacts disk gas, with flow velocity $\mathbf{u}$ allows us to derive the
force imparted by the infalling material (the shear source term in Equation 2)

$\mathbf{F}_{\rm in}(R,t)=\dot{\Sigma}_{\rm in}(\mathbf{u}_{\rm
in}(R,t)-\mathbf{u}(R,t)).$ (8)

The components of the radial and azimuthal flow velocities are given by the
UCM projectile velocities:

$\begin{split}u_{R,{\rm in}}(R,t)&=-\left(\frac{GM_{*}}{R}\right)^{1/2},\\\
u_{\phi,{\rm in}}(R,t)&=\left(\frac{GM_{*}}{R_{\rm c}}\right)^{1/2}.\\\
\end{split}$ (9)

For each episode of infall, set to a duration of $t_{\rm ep}=1.5\times 10^{4}\
\mathrm{yr}$, the simulation is evolved for an additional period of $2.0\times
10^{4}\ \mathrm{yr}$ without infall, in order to investigate the longevity of
infall-produced structures during quiescence. The disk gains $\sim 0.015\
\mathrm{M_{\odot}}$ of material during the course of infall, about $20\%$ of
its initial mass.

The models used in this parameter study are detailed in Table 1. Their
properties are set by the inner and outer boundaries of the infall zone,
$R_{\rm in}$ and $R_{\rm c}$, and the relative change in the specific angular
momentum of the envelope,

$\delta j=(j-j_{0})/j_{0}.$ (10)

The parameters were chosen based on the results of K20 to comprise a mix of
filamentary and spherical infall scenarios and represent the probable range of
values of $\delta j$.

We assume that the initial disk surface density profile (Equation 3) with
radius $R_{\rm d}$ was set by a prior episode of infall with specific angular
momentum $j_{0}=(GM_{*}R_{\rm d})^{1/2}$. Thus, in our model, the centrifugal
radius of infall, $R_{\rm c}$, parameterizes the relative change in specific
angular momentum

$\delta j=\left(\frac{R_{\rm c}}{R_{\rm d}}\right)^{1/2}-1.$ (11)

The results of K20 suggest that for a typical envelope accretion event of
similar duration to $t_{\rm ep}$, the probable range of relative angular
momentum change is $|\delta j|=0.005-0.3$. The most commonly occurring
episodic accretion events have $\delta j=\pm 0.1$; these are denoted as ave.
in Table 1. The K20 runs show that 60% of accretion events are primarily 1- or
2-dimensional, occurring in streams or sheets. In this paper, we classify runs
using the width parameter $f^{2}$ defined in Equation (7), with runs having
$f^{2}\lesssim 0.4$ described as stream-like, those with as sheet-like, and
those with $f^{2}>0.75$ as cloud-like. We choose runs from K20 that cover the
relevant parameter space of anisotropic infall conditions, as shown in Figure
1, so as to ensure each type of geometry is represented for every choice of
$\delta j$.

Figure 1:  The parameter space of anistropic infall explored across the
variation of $\delta j$ in the infalling envelope and the filling factor of
infall $f^{2}=\Delta R/R_{\rm c}$, where the positions of runs in the
parameter space used in this paper are indicated by black edged circles. Runs
corresponding to the more typical accretion events in K20 are shown as bolder
circles. For reference, the fiducial run _rd_sheet_ is marked with a central
dot.

The infall properties, $\delta j$ and $f$, can be expressed as a function of
the characteristic radii of the infall and initial disk radius: $R_{\rm
c}/R_{\rm d}$ and $\Delta R/R_{\rm c}$. In this paper, we discuss results
across the parameter space in terms of whether infall occurs in the inner disk
$R_{\rm c}<R_{\rm d}$ or outer disk $R_{\rm c}>R_{\rm d}$ and if accretion is
stream-, sheet-, or cloud-like.

Model |  $R_{\rm in}$ [au] |  $R_{\rm c}$ [au] | $\delta j$ | $f$
---|---|---|---|---
$\delta j$ = 0
rd_stream | 45 | 60 | 0 | stream
rd_sheet | 30 | 60 | 0 | sheet
rd_cloud | 15 | 60 | 0 | cloud
$-\delta j$
45_stream | 30 | 45 | ave. | stream
45_sheet | 15 | 45 | ave. | sheet
45_cloud | 5 | 45 | ave. | cloud
30_stream | 20 | 30 | high | stream
30_sheet | 15 | 30 | high | sheet
30_cloud | 5 | 30 | high | cloud
$+\delta j$
75_stream | 60 | 75 | ave. | stream
75_sheet | 45 | 75 | ave. | sheet
75_sheet | 30 | 75 | ave. | sheet
75_cloud | 5 | 75 | ave. | cloud
105_stream | 75 | 105 | high | stream
105_sheet | 60 | 105 | high | sheet
105_cloud | 30 | 105 | high | cloud
Table 1: Parameter Study Models

###  2.4 Post-processing Dust Transport

One of our goals is to determine whether infall episodes generate ring-like
substructures in the dust density by perturbing the radial pressure gradient
of the disk. In order to do this, we post-process each run at the end of the
infall episode with a 1D dust transport algorithm.

We parameterize the pressure gradient with

$\eta=\frac{1}{2}\left(\frac{c_{\rm s}}{v_{\rm k}}\right)^{2}\frac{d\ln
P}{d\ln R},$ (12)

where $\eta$ can be shown to be the fractional deviation of the azimuthal
velocity from the Keplerian velocity $v_{\rm k}$,

$\eta=(u_{\phi}-v_{\rm k})/v_{\rm k}$ (13)

As defined, for a typical unperturbed disk pressure gradient, $\eta$ is
negative,so the gas in the disk orbits at a sub-Keplerian speed, supported by
outward gas pressure.

Dust particles, however, do not have internal pressure, and their motions are
dictated by how well-coupled they are to the gas flow. The dust-gas coupling
can be described by the Stokes number $\mathrm{St}$, the dimensionless ratio
of the drag stopping time to the local orbital time. Particles with
$\mathrm{St}\ll 1$ move with the gas, while those with $\mathrm{St}\gg 1$
orbit closer to Keplerian speeds. For physical conditions within a typical
protoplanetary disk in vertical hydrostatic equilibrium, the Stokes number of
particles in the midplane as a function of the local gas surface density, ,
and the grain size $a$ is (Birnstiel et al., 2016)

$\mathrm{St}=\frac{a\rho_{\rm s}}{\Sigma_{\rm g}}\frac{\pi}{2},$ (14)

where $\rho_{\rm s}$ is the solid material density of the grains. Following
Weidenschilling (1977), the radial drift velocity of dust particles in an
accretion disk depends primarily on $\eta$ and $\mathrm{St}$ as

$u_{\rm d}=\frac{\eta v_{\rm k}}{\mathrm{St}+1/\mathrm{St}}.$ (15)

For simplicity, we adopt an advection dominated model of dust transport, such
that the local change in the dust surface density, $\Sigma_{\rm d}$, driven by
the drift velocity and gradients in the surface density and drift velocity, is

$\begin{split}\frac{\partial\Sigma_{\rm d}}{\partial
t}&=-\frac{1}{R}\frac{\partial}{\partial R}\left(R\Sigma_{\rm d}u_{\rm
d}\right)\\\ &=-\Sigma_{\rm d}\frac{u_{\rm d}}{R}-\Sigma_{\rm
d}\frac{\partial{u_{\rm d}}}{\partial R}-\frac{\partial\Sigma_{\rm
d}}{\partial R}u_{\rm d}.\end{split}$ (16)

In this paper, we calculate the dust transport during the post-infall period
due to radial drift for

particle sizes in sixteen logarithmically spaced bins

in the range $\ \mathrm{[\mu m,cm]}$, at a solid density of $\rho_{\rm s}=1.2\
\mathrm{g\ cm^{-3}}$.

For the disk setup used here, all particle sizes have $\mathrm{St}<1$ across
the disk. We assume an initial ISM dust-to-gas ratio $\epsilon=\Sigma_{\rm
d}/\Sigma_{\rm g}=0.01$ at $t=t_{\rm ep}$, using the azimuthal average of
simulation outputs as the input, and numerically integrate the radial dust
transport equation (16) with a first order 1D conservative upwind scheme,
calculating $u_{\rm d}$ at each timestep from the corresponding simulation
output. These calculations are presented in §3.4 as a conservative estimate on
the relative range and distribution of dust substructures possible as a result
of this infall model. Dedicated studies of dust evolution for cases of
interest will be done in future work.

##  3 Infall-induced Effects

###  3.1 Infall and the Rossby Wave Instability

Figure 2: Growth of the maximum non-axisymmetric surface density perturbations
$\delta\Sigma/\overline{\Sigma}$ scaled with the orbital timescale at the
inner radius $\tau_{\rm in}$ shows the most unstable modes launch the RWI and
how they proceed through the linear phase and onto non-linear growth. The
color scale of the lines matches the parameters in Figure 1, scaling with
$R_{\rm c}/R_{\rm d}$. Similar $\Delta R/R_{\rm c}$ have similar launch times
across $R_{\rm c}/R_{\rm d}$.

The modeled streams add mass and angular momentum within the prescribed infall
zone $[R_{\rm in},R_{\rm c}]$, with a radial dependence described by Equations
(6) and (8). While the magnitude and extent of the effect on the disk across
the infall zone will vary based on the parameters of the infall, at both
$R_{\rm in}$ and $R_{\rm c}$, where there is a transition between the
unperturbed disk profile and the altered region of the infall zone, the disk
will develop local bumps in the surface density and pressure gradients. These
perturbations create prime conditions for the development of the RWI.

The traditional stability criterion of the RWI derived from linear analysis
(Lovelace et al., 1999), requires that the disk vortensity profile

$q(R)=\frac{d/dR(R^{4}\Omega^{2})}{2R^{3}\Sigma\Omega}$ (17)

have at least one extremum. Practically, depending on the initial disk
structure, this means that local perturbations in the surface density and
rotational velocity profiles both contribute to the development of the RWI.
The non-linear growth and saturation of the RWI has been studied extensively
with local box simulations mapping the growth rates of unstable regimes across
parameter spaces of the shape, width, and amplitude of local perturbations
(i.e. Ono et al., 2018). In this study, we focus on investigating the long-
term effects of the RWI on the disk structure and dynamics in the context of
our prescribed physical scenario, rather than characterizing the properties of
the RWI itself.

In Figure 2, we show the development of the phases of the RWI in terms of the
growth of the fastest growing mode by plotting the maximum non-axisymmetric
perturbation to the mean surface density $\delta\Sigma/\overline{\Sigma}$ over
time, where $\delta\Sigma=\Sigma-\overline{\Sigma}$ and $\overline{\Sigma}$
denotes the azimuthal average of the disk surface density. The launch of the
RWI occurs when the magnitude of non-axisymmetric perturbations grows past the
initial noise level such that $(\delta\Sigma/\overline{\Sigma})_{\rm max}\sim
10^{-4}$.

During the linear regime, the growth rate—the slope of the curves in Figure
2—is roughly constant. The instability saturates and becomes non-
linear,resulting in the perturbations no longer growing as they approach the
local disk average. Given that the global model has two regions where the RWI
can develop and in some cases these regions have separations small enough for
interaction, some non-linear behavior such as mode-mode coupling can be
expected during what would traditionally be called the linear phase. For
example, runs where the infall zone lies outside the disk radius $R_{\rm d}$
show fluctuations in the growth rate. In fact, we can not expect completely
analogous behavior to dedicated RWI studies, given that during infall,
perturbations to the surface density and azimuthal velocity are consistently
reinforced.

Figure 3: The growth of the disk surface density in the central 100 au of the
disk during each epoch: the linear phase of the RWI, the non-linear phase of
the RWI, the end of the infall episode, and post-infall. (Left:) The fiducial
run _rd_sheet_ where $\delta j=0$ such that $R_{\rm c}=R_{\rm d}$ and infall
is sheet-like. (Center:) Comparison of inner and outer disk infall runs,
_45_sheet_ and _75_sheet_ , for which $R_{\rm c}<R_{\rm d}$ and $R_{\rm
c}>R_{\rm d}$ from top to bottom. (Right:) Comparison of runs _rd_cloud_ and
_rd_stream_ in which the annulus width $\Delta R$ is widened or narrowed
compared to the fiducial run. All panels have colorbars scaled to that of the
fiducial run for each respective epoch. The boundaries of the infall zone are
annotated with white dotted circles.

While the vortex formation characteristic of the RWI proceeds more rapidly at
$R_{\rm in}$, the growth of the RWI is not simply dependent on the orbital
timescales at $R_{\rm in}$, given the presence of two zones in the disk that
are liable to be unstable to the RWI. We find that the average growth rate
$\gamma$ in terms of the orbital timescale at the outer radius, $\tau_{\rm
c}=2\pi\Omega^{-1}_{\rm k}(R_{\rm c})$, varies across filling factor as
$\gamma\sim 0.6,1,$ and $3\ \tau_{\rm c}^{-1}$ for stream-, sheet-, and cloud-
like runs respectively. Likewise, for each filling factor we find the minimum
saturation time for the RWI, $t_{\rm RWI}=8,5,$ and $2\ \tau_{\rm c}$.

We use the nominal scaling of $t_{\rm RWI}$ and our chosen infall episode
length ($t_{\rm ep}=1.5\times 10^{4}\ \mathrm{yrs}$) to designate four main
epochs of interest in our analysis and their representative times:

  * •

the linear phase of the RWI: $t\lesssim t_{\rm RWI}$

  * •

the non-linear phase: $t\sim 1.5\ t_{\rm RWI}$

  * •

the end of the infall episode: $t=t_{\rm ep}$

  * •

post-infall: $t\lesssim 3.5\times 10^{4}\ \mathrm{yrs}\simeq 2.33\ t_{\rm ep}$

In Figure 3, we show the evolution of perturbations to the 2D surface density
$\Delta\Sigma/\Sigma_{0}$ over the course of the four epochs. The leftmost
panels show the results from the fiducial run _rd_sheet_ , a sheet-like infall
with $R_{\rm in}=30\ \mathrm{au}$ and $R_{\rm c}=R_{\rm d}=60\ \mathrm{au}$.
During the linear phase, perturbations are prominent within the infall zone,
with vortices forming at $R_{\rm in}$. As the RWI becomes non-linear, vortices
at the inner edge of the infall zone merge and spiral waves are driven outside
of $R_{\rm c}$. At the end of infall and post-infall, $R_{\rm in}$ hosts a
single vortex—its spiral wake visible in contrast to the gaseous bump in the
infall zone. The other panels in Figure 3 demonstrate how this behavior varies
relative to the fiducial run if infalling material has less or more angular
momentum (center) or is more or less filamentary (right). The most dramatic
departures in the morphology of structures are achieved by narrowing or
widening the ring. Narrow annuli of deposition, in addition to leading to a
higher average surface density in the infall zone, create large-scale non-
axisymmetric structures.

Figure 4: Evolution of infall zone vortices and spiral features at each epoch:
the linear phase of the RWI, the non-linear phase of the RWI, the end of the
infall episode, and post-infall for (left) the fiducial run _rd_sheet_ such
that $R_{\rm c}=R_{\rm d}$ and infall is sheet-like; (center) runs _45_sheet_
and _75_sheet_ for inner and outer disk infall respectively; and (right) runs
_rd_cloud_ and _rd_stream_ in which the annulus width $\Delta R$ is widened or
narrowed compared to the fiducial run. Each panel is in $\phi$-$r$
coordinates, 60 au wide, and centered on the infall zone. Vortex centers have
normalized, density-weighted vorticity $\delta\omega\sim 1$, while spiral
features are visible as alternating stripes of $\delta\omega\pm 0.1$.

RWI vortices are visible in Figure 3 as localized surface density
perturbations at the edges of the infall zones. We use the vorticity,
$\omega=\nabla\times\bf{u}$, to define a normalized density-weighted vorticity
factor,

$\delta\omega=\Sigma(\omega-\omega_{0})/(\Sigma_{0}\omega_{0})$ (18)

plotted in Figure 4 at the four main epochs for the same set of selected runs
as in Figure 3 in $\phi$-$r$ space. With this normalization, the centers of
vortices have $\delta\omega\sim 1$ and alternating regions of
$\delta\omega\sim\pm 0.1$ highlight the peaks and troughs of spiral waves.

During the linear phase, individual vortices form at the edges of the infall
zones, each vortex appearing to be the foot point for a spiral wave. With
vortex formation occurring at both the inner and outer edges of the infall
zone, spiral density waves are launched from both $R_{\rm in}$ and $R_{\rm
c}$, resulting in interference within the infall zone itself. As the
instability grows, smaller denser vortices merge into larger, more diffuse
structures, with fewer spiral arms. Merged vortices are wider at outer radii,
their radial extent scaling with the scale height H at vortex center.

While most vortices at the inner edge of the infall zone merge into a single
vortex within the infall episode period, runs with the same $R_{\rm c}$ do not
necessarily share the same timescales for vortex merging. The average number
of total vortices across each run depends on the orbital timescale of the
midpoint between $R_{\rm in}$ and $R_{\rm c}$, suggesting that interactions
between modes in narrower rings of infall can affect the evolution of vortex
structures. Even so, interference between spiral waves within the infall zone
can still be seen for the widest separations.

Through the mechanism of the RWI, annular infall has implications for two
major disk processes: inducing turbulence and thereby transport of disk
materials (§3.2) and the creation of pressure bumps liable to trap dust
particles by arresting radial drift (§3.4). We explore the longevity of
structures across our parameter space after infall in §3.3.

###  3.2 Turbulence and Transport

During infall, the bulk of mass transport is driven directly by infall as
shown by the average radial mass flux within the disk:

$\dot{M}(R)=2\pi R\overline{\Sigma u_{R}}$ (19)

for our comparative set of models (Figure 5). Inward mass flux—where
$\dot{M}(R)<0$ —is at most 20–30% of the total mass infall rate $\dot{M}_{\rm
in}=10^{-6}M_{\odot}\ \ \mathrm{yr}^{-1}$.

Exterior to $R_{\rm c}$, mass flows outward, contributing to spreading of the
gas ring accumulating within the infall zone. Interior to $R_{\rm in}$, in the
innermost regions of the disk, the radial mass flux is dominated by wave
transport, particularly for models with inner disk infall such that $R_{\rm
c}<R_{\rm d}$. We expect that spiral waves may transport angular momentum
outward, thus allowing further mass accretion to occur. However, during the
infall stage mass flux in these inner regions is at most 10% of $\dot{M}_{\rm
in}$. At this stage, the time average of the mass flux within the outer extent
of the infall zone is between $\dot{M}(<R_{\rm c})\sim
10^{-8}$–$10^{-7}M_{\odot}\ \mathrm{yr}^{-1}$ for all disks, with models of
stream-like infall experiencing the most fluctuation within that range.

Figure 5: The average radial mass flux in units of total mass infall rate,
$\dot{M}_{\rm in}=10^{-6}M_{\odot}\ \mathrm{yr^{-1}}$, across the inner 100 au
at the end of infall—time averaged across 0.3 kyr—where negative values
indicate inward mass transport. (Left:) the fiducial run, sheet-like accretion
with $R_{\rm c}=R_{\rm d}$. (Center:) sheet-like runs with inner disk $R_{\rm
c}<R_{\rm d}$ and outer disk $R_{\rm c}>R_{\rm d}$ infall. (Right:) cloud- and
stream- like infall with $R_{\rm c}=R_{\rm d}$.

In addition to the injection of momentum in the radial direction, the impact
of infalling material onto the disk generates azimuthal shear (see Eq. 9).
Shear stresses, also called Reynolds stresses, can drive turbulence. We
calculate the turbulent viscosity parameter $\alpha$ due to Reynolds stresses
at each radius:

${}\alpha_{\rm Rey}(R)=\frac{\int\Sigma\delta u_{R}\delta u_{\phi}\ R\
d\phi}{\int\Sigma c_{\rm s}^{2}\ R\ d\phi}$ (20)

where $\delta u_{R}$ and $\delta u_{\phi}$ are the residual velocities from
the azimuthal average, defined by $\delta u=u-\overline{u}$. The disk-wide
$\alpha_{\rm Rey}$ is computed as the radial average of Equation (20). The
total $\alpha$ – a sum of the induced Reynolds $\alpha$ and the background
viscosity parameter: $\alpha=\alpha_{\rm Rey}+10^{-7}$ – is an average
estimate for the turbulent stresses acting over time at a particular radius in
the disk; the shear stresses produced by spiral waves will vary over azimuth.

In the prescribed infall model, azimuthal shear is introduced across the
entire infall zone, with the largest velocity shear initially at the inner
edge of the infall zone $R_{\rm in}$. Once the RWI develops, spiral arms
driven by vortices contribute to the total mechanical shear stress. For
stream-like infall with narrow infall zones, interference between spiral waves
drives turbulence within the infall zone (Figure 6), raising the overall level
of $\alpha$.

In the case of inner disk infall where $R_{\rm c}/R_{\rm d}<1$, $\alpha$
decreases outside of the infall zone, in direct contrast with outer disk
infall ( $R_{\rm c}/R_{\rm d}>1$), where $\alpha$ grows outward, peaking much
further out in the disk, well outside of the infall zone. The fiducial case
where $R_{\rm c}=R_{\rm d}$ produces the most constant distribution of
$\alpha$ across the disk.

Figure 6: Measured turbulence due to infall across three zones of the disk,
$R\leq R_{\rm in}$, $R_{\rm in}<R\leq R_{\rm c}$, and $R>R_{\rm c}$, shown
here with each zone scaled to the same apparent width for illustrative
purposes. Panels from top to bottom show $\alpha$ at $t=t_{\rm ep}$ for inner
disk ($R_{\rm c}<R_{\rm d}$), fiducial ($R_{\rm c}=R_{\rm d}$), and outer disk
($R_{\rm c}>R_{\rm d}$) infall, respectively.

Stream and sheet-like infall produces, on average, larger $\alpha$ values than
more spherical deposition, but the location of the centrifugal radius $R_{\rm
c}/R_{\rm d}$ plays the most important role in determining the disk-wide
$\alpha$. Outer disk infall drives the largest amount of turbulence across the
largest fractions of the disk, with values reaching $\alpha\sim
10^{-3}$–$10^{-2}$, around two orders of magnitude greater than the typical
values with inner disk infall.

During infall, the average disk-wide $\alpha$ grows on similar timescales to
the RWI, with the initial growth scaling as $\tau_{\rm in}$. As the
instability saturates, the magnitude of the disk-wide $\alpha$ approaches a
constant value with time (Figure 7). While more annular concentrated
deposition develops robust $\alpha$ on slower timescales, this stream-like
infall is on the whole just as effective at generating turbulence, likely in
part due to the elevated turbulence generated within the infall zone for those
cases.

Figure 7:  The disk-wide measure of the turbulent $\alpha$ parameter, over the
course of simulation runtime. Circles mark $t_{\rm RWI}$ for each run. Figure
8: The disk-wide $\alpha$ (circles) and maximum $\alpha$ (triangles) at
$t=t_{\rm ep}$ for each model, plotted against the model’s mass loading
parameter $\log\langle\Sigma_{in}\rangle-\log\Sigma_{0}(R_{\rm c})$.
Corresponding values from the standard run in Bae et al. (2015) are plotted as
starred markers. The black and gray dotted lines have slopes of 2 and 1,
respectively. Points are colored by their outer infall radius $R_{\rm c}$
following Figure 1, .

The strength of the RWI driven turbulence depends on the size of the initial
perturbations produced by infall mass loading. Previous studies find the
amount of angular momentum transport by spiral waves launched from vortices is
directly correlated to their size (Paardekooper et al., 2010), which is
influenced by both the scale height at vortex center and the relative
amplitude of the initial surface density perturbation (Ono et al., 2018), both
of which scale with the centrifugal radius of infall in this work. To measure
the initial perturbations produced by infall, we define a quantity for the
relative mass loading

$\langle\Sigma_{\rm in}\rangle/\Sigma_{0}(R_{\rm c})$, the ratio of the
average surface density perturbation per dynamical time

$\langle\Sigma_{\rm in}\rangle=\frac{\dot{M}_{\rm in}}{4\pi(R_{\rm
c}^{2}-R_{\rm in}^{2})}\tau_{\rm c}$ (21)

to the initial disk surface density at the centrifugal radius
$\Sigma_{0}(R_{\rm c})$. Plotted in log-log space in Figure 8, the disk-wide
$\alpha$ at the end of infall has a tight relationship with the relative mass
loading, where the black dotted reference line has the form:

$\log\alpha=2\ [\log\langle\Sigma_{\rm in}\rangle-\log\Sigma_{0}(R_{\rm
c})]-0.5.$ (22)

Compared to the disk-wide $\alpha$ at the same simulation time, the
relationship to mass loading of the maximum $\alpha$ achieved in the disk has
more scatter; the gray reference line for those points has a slope of unity,
for comparison. We calculate the mass-loading ratio for the standard infall
model in Bae et al. (2015) and find that their measured $\alpha\sim
10^{-3}$–$10^{-2}$, is also consistent with this scaling.

Overall, we find that in these models the magnitude of measured turbulence
driven by infall is heavily dependent on whether it occurs primarily in the
inner or outer disk, since the maximum achievable $\alpha$ is directly
correlated with the outermost radius of the infall zone $R_{\rm c}$. Outer
disk infall drives the highest magnitude turbulence across the largest
fraction of the disk area.

The relationship in Figure 8 is

thus consistent with the scaling with $R_{\rm c}$ expected from Paardekooper
et al. (2010) and Ono et al. (2018).

###  3.3 Post-Infall Conditions

Each model was evolved for an additional $2\times 10^{4}\ \mathrm{yr}$ after
$t_{\rm ep}$ in order to investigate the longevity of flow structures produced
by infall. While some structures show evidence of spreading, particularly for
cases of outer disk infall, the radial profile of the surface density does not
evolve significantly in the post-infall phase (see Figure 9). However, as
shown in Figure 4, the post-infall phase marks the start of vortex migration,
where individual vortices at $R_{\rm in}$ or $R_{\rm c}$ can be seen to have
drifted inward from their initial positions. This basic migration behavior is
consistent with other numerical simulations of Rossby vortices (Paardekooper
et al., 2010), and has been attributed to a net exchange of angular momentum
between the disk and vortex transported by spiral density waves due to a
difference between the angular momentum flux at the inner and outer boundary
of the vortex.

Figure 9: The surface density (top), pressure gradient parameterized as $\eta$
(center), and measured $\alpha$ parameter (bottom) for the fiducial run shown
at the end of infall and post-infall, $t=1.5\times 10^{4}\ \mathrm{yrs}$ and
$t=3.5\times 10^{4}\ \mathrm{yrs}$, in dark blue and light blue dashed lines,
respectively. The initial profiles for $\Sigma_{0}$ and $\eta_{0}$ are shown
in gray.

We estimate the inward migration rate of vortices by tracking the change in
the radial location of the innermost peak in vorticity at $R_{\rm in}$ and
$R_{\rm c}$, calculating the average migration speed over the post-infall
period. As in previous studies, there is a correlation between migration rates
and the scale height at vortex center $H_{\rm v}$. At $R_{\rm in}$, inward
vortex migration rates

$\begin{split}\tau_{\rm v}&\sim 2\times 10^{-8}\ \mathrm{au\
yr^{-1}}\left(\frac{H_{\rm v}}{0.05\ \mathrm{au}}\right),\\\ &\sim 2\times
10^{-8}\ \mathrm{au\ yr^{-1}}\left(\frac{R_{\rm in}}{1\
\mathrm{au}}\right)^{5/2}\end{split}$ (23)

have the same scaling with $H_{\rm v}$ to that of isolated vortices in
Paardekooper et al. (2010). At $R_{\rm c}$, however, we find that inward
migration rates are generally slower than for inner infall zone vortices at
comparable radii, but have a much steeper relationship with scale height,

$\begin{split}\tau_{\rm v}&\sim 3\times 10^{-11}\ \mathrm{au\
yr^{-1}}\left(\frac{H_{\rm v}}{0.05\ \mathrm{au}}\right)^{3}\\\ &\sim 3\times
10^{-11}\ \mathrm{au\ yr^{-1}}\left(\frac{R_{\rm c}}{1\
\mathrm{au}}\right)^{15/4}.\end{split}$ (24)

During infall, the disk surface density is consistently perturbed, leading to
the formation of multiple generations of RWI vortices, particularly for
vortices at outer radii (Figure 4). In our case, angular momentum exchange is
complicated by interactions between multiple generations of vortices, such
that inward migration speeds of the first generation vortices are not direct
proxies for angular momentum transport.

In Paardekooper et al. (2010), pressure bumps are argued to impede vortex
migration, but our post-infall calculations are more consistent with Ono et
al. (2018), who find that pressure bumps induced self-consistently by the RWI
do not prevent vortex migration. In this study, we find that RWI induced
pressure perturbations tend to migrate along with their associated vortex.

However, the stalling effect of external pressure and surface density bumps on
migration discussed in Paardekooper et al. (2010) could become important for
the long-term evolution of vortices produced at $R_{\rm c}$ and migrating
through the infall zone.

We do note that these comparisons to local box simulations are largely
qualitative, as dedicated RWI studies induce the RWI with constant amplitude
density perturbations, have constant numerical viscosity, and do not include
the effects of additional velocity shear due to the impact of infalling flow
streams; these additional effects reinforce the vortensity profile of the disk
and locally alter the turbulent $\alpha$.

Based on the range of migration speeds $\sim 10^{-5}$–$10^{-3}\ \mathrm{au\
yr^{-1}}$ , the fastest migration timescales would be on the order of several
hundred kiloyears, with more typical estimates on the order of several
megayears.

During infall, we can expect the turbulent $\alpha$, in part, to be locally
driven by the momentum deposition of infalling material, particularly for
stream-like infall. Post-infall, without additional turbulent forcing, we find
that the average disk $\alpha$ is sustained at (or surpasses) its end of
episode value for infall parameters $R_{\rm c}>R_{\rm d}$ (Figure 7).

The increase in $\alpha$ for outer disk infall is associated with infall
models that host multiple generations of vortices, producing additional shear
stresses.

However, $\alpha$ does decay post-infall for episodes of inner disk infall.
The decay times for the value of the disk-wide $\alpha$ to drop by a factor of
two, range from $\tau_{\alpha}=1$–$4\times 10^{4}\ \mathrm{yr}$. This decrease
for inner disk infall models largely occurs at radii $R>R_{\rm c}$ (Figure 9).

###  3.4 Effects on Dust Evolution

After infall, the dynamical landscape of the disk has been significantly
altered. Following the model in Birnstiel et al. (2016), we look at the
relevant timescales for basic dust evolution. The speed of dust transport by
turbulent diffusion is characterized by the diffusion timescale

$\tau_{\mathrm{diff}}=L^{2}\frac{(1+\mathrm{St}^{2})}{\nu},$ (25)

where the characteristic length scale is $L$ and the viscosity is set by the
total measured $\alpha$. Transport due to the inward radial drift of
particles, induced by the pressure gradient (see Equation 12) is described by
the radial drift timescale

$\tau_{\mathrm{drift}}=L\frac{(\mathrm{St}+\mathrm{St}^{-1})}{\eta v_{\rm
k}}.$ (26)

We also compare the collisional growth timescale due to particle sticking

$\tau_{\mathrm{growth}}=\frac{a\rho_{\rm s}}{\rho_{\rm d}\Delta v},$ (27)

where we assume relative collisional velocities due to turbulence

$\Delta v=c_{\rm
s}\left(\frac{3\alpha}{\mathrm{St}+\mathrm{St}^{-1}}\right)^{1/2}$ (28)

from Ormel & Cuzzi (2007) and use the dust volume density with a vertical
distribution set by settling (Birnstiel et al., 2016)

$\rho_{\rm
d}=\frac{\Sigma_{g}\epsilon}{H\sqrt{2\pi}}\left(\frac{\alpha}{\alpha+\mathrm{St}}\right)^{1/2}.$
(29)

In Figure 10, we compare these timescales for dust evolution at the start of
the simulation and at $t=t_{\rm ep}$. We assume an initial dust-to-gas ratio
$\epsilon=0.01$ equal to the ISM value for three sizes of dust: 1 $\mu$m, 100
$\mu$m, and 1 mm. Using $L=0.1R$ as the operative length scale, we compare the
times over which one might expect changes in the radial dust distribution of
at least $10\%$.

Figure 10:  Diffusion (magenta), radial drift (orange), and collisional growth
(navy) timescales at $t=0$ (left) and $t=t_{\rm ep}$ (right) for sheet-like
infall with $R_{\rm c}=R_{\rm d}$ for 1 $\mu$m (top), 100 $\mu$m (middle), and
1 mm (bottom) dust. Small particles initially show long time-scales for drift
and diffusion, but after the infall episode, disk perturbations and infall-
induced turbulence shorten $\tau_{\rm diff}$, shorten $\tau_{\rm growth}$ in
the outer disk for micron-sized particles, and lengthen $\tau_{\rm drift}$
within the pressure bump region. Dotted lines show (vertical) the locations of
$R_{\rm in}$ and $R_{\rm c}$ and (horizontal) the characteristic dynamical
timescale.

Comparing to the timescales over which we evolve the dynamics, on the order of
$\sim 10^{4}\ \mathrm{yr}$, micron-sized grains would initially show limited
growth at disk radii above $60\ \mathrm{au}$. However, with infall, the growth
timescale is shortened for small dust out to nearly $120\ \mathrm{au}$ for the
fiducial run shown here. Larger grains at outer radii have fairly rapid drift
timescales, so inward drift can push grains into growth dominated regimes.
Without a means to arrest radial drift, this leads to the standard picture of
drift limited growth. In the case of infall, however, rapid inward drift from
the outer disk could lead to pileup within pressure bumps, potentially leading
to accelerated growth there. The diffusion timescales also become relatively
more important within the inner disk for small grains, although still on
longer timescales than the post-infall period we simulate here. While Figure
10 shows approximate timescales, it suggests that the post-infall disk
environment could be more conducive to planet formation than the typical
steady-state disk.

Moving beyond estimates using the drift velocity, drift transport depends on
multiple components. In Equation 16, perturbations to the dust surface density
due to particle drift depend on (1) the local drift velocity (2) the
differential drift velocity, because of a traffic jam-like effect, and (3) the
gradient in the surface density that dictates the bulk transport of material
along the flow. The main pressure bump corresponds to where $u_{\rm d}>0$, but
most of the contribution to local drift arises from gradients in $u_{\rm d}$,
which peak near both $R_{\rm in}$ and $R_{\rm c}$, rather than just at the
inward edge of the infall zone.

Note that in this infall model, perturbations to the pressure gradient occur
not only due to perturbations in the gas surface density, but also as a direct
consequence of the velocity shear induced by infalling material (Eq. 8).Thus,
substructure in the gas surface density will not exactly match the amplitude
and location of substructure in the dust surface density—a potential
discriminating factor between rings induced by infall and those directly
induced by planets, further discussed in § 4.

###  3.5 Rings and gaps

We calculate the amplitude of relative perturbations in the dust density
(§2.4) as the fractional change in dust surface density due to drift relative
to an initially smooth surface density profile

$A=\frac{\Sigma}{\Sigma_{0}}-1$ (30)

at $t=t_{ep}+15\ \mathrm{kyr}$. The results presented here are a conservative
estimate of the _relative_ locations and amplitudes of potential dust
substructures across our parameter space of simulations and possible dust
sizes, without any model assumptions about the growth or fragmentation
behavior of dust populations.

Figure 11:  (Top:) Relative amplitude of dust surface-density perturbation $A$
(Eq. 30) as a result of local dust concentration in pressure bumps for
particle sizes in 16 bins. (Bottom:) Relative locations of rings of
millimeter-sized dust (orange lines) compared to the gas surface density
(teal), normalized to the same amplitude. In each column, the middle plot
(dashed line) shows the sheet-like infall of the top panel with the stream-
like and cloud-like models with the same $R_{\rm c}$ shown for comparison.
(Left:) The fiducial run, sheet-like accretion with $R_{\rm c}=R_{\rm d}=60\
\mathrm{au}$. (Center:) Inner disk, with $R_{\rm c}=45\ \mathrm{au}$. (Right:)
outer disk, with $R_{\rm c}=75\ \mathrm{au}$. The maximum grain size due to
fragmentation $a_{\rm frag}$ is shown with a dashed white line. $R_{\rm in}$
for each model and $R_{\rm c}$ are denoted by gray dotted lines.

In Figure 11, we show the amplitude of perturbations versus radius for
selected runs post-infall to demonstrate the relative locations of rings and
gaps in the dust and gas for a range of dust sizes. For reference, we plot the
size corresponding to the fragmentation limit due to turbulent collisions
(Birnstiel et al., 2016)

$a_{\rm frag}\simeq 0.08\frac{\Sigma_{\rm
g}}{\rho_{s}\alpha}\left(\frac{v_{\rm frag}}{c_{\rm s}}\right)^{2}$ (31)

as a dashed white line, where $v_{\rm frag}$ is the fragmentation velocity,
for which we use the threshold velocity $10\ \mathrm{m~{}s^{-1}}$ for
predominantly icy grains in the outer disk (Gundlach & Blum, 2015). The
turbulent fragmentation size limit extends down to millimeter sizes for the
case of outer disk infall; the estimates of the perturbation amplitudes in
those cases should be taken as an upper limit. The outer disk infall case,
where $\alpha$ is constant outside of the infall zone (Figure 6), also
illustrates a physically self-consistent scenario in which the fragmentation
barrier is tracing out the same region in $a$-$R$ space as the radial drift
line.

As expected, the largest dust sizes experience the most overall inward drift.
In most cases, this would lead to fairly compact dust distributions, if not
for the robust pressure bumps that serve as drift barriers. Infall-induced
pressure variations are strongest at $R_{\rm in}$ compared to the pressure
bumps at $R_{\rm c}$, with the exception for cases of stream-like infall.
Large grains concentrate in rings near both $R_{\rm in}$ and $R_{\rm c}$,
migrating inward with time, such that most rings are several astronomical
units inward of their respective infall zone edges. In the dust profiles in
the bottom panel of Figure 11, the outermost peak at $\sim 100\ \mathrm{au}$
corresponds to the location of the dust line, demarcating the progression of
radial drift for millimeter sized particles over the post-infall period.

The width of dust ring structures is correlated with their location in the
disk. Consistent with dedicated studies of the RWI (Ono et al., 2018), we find
the width of structures correlates with the disk scale height at peak center
where the peak full width at half maximum is roughly half of the scale height
at its center. Relative amplitudes of the dust surface density perturbations
have opposite trends depending on whether structures are in the outer or inner
disk, where the trend for $A$ increases for $R<R_{\rm d}$ and decreases for
$R>R_{\rm d}$. Relative to the gas, peaks in the dust are narrow and
concentrated near the edges of the infall zone, whereas the gas surface
density is dominated by a broad ring.

##  4 Discussion

###  4.1 Hallmarks of Heterogeneous Infall

The traditional analytic framework that informs our understanding of disk
formation and evolution (e.g. Terebey et al., 1984) does not produce
structured disks. Such infall of uniformly rotating spherical envelopes
produces smooth disks that grow outward with time such that $R_{\rm c}\propto
t^{3}$. Previous studies of 1D dust evolution in such viscously evolving disks
during the initial self-similar collapse stage find the smoothly varying
infall of the traditional model to have little impact in the overall outcome
of dust population growth (Birnstiel et al., 2010).

Numerical simulations that capture the onset of disk instabilities have
established that the infall of envelope material can have dynamical
consequences relevant for disk transport and planet formation processes (Lesur
et al., 2015; Bae et al., 2015). However, Bae et al. (2015) show that even
with the ability to incite the RWI, monotonically increasing $R_{\rm c}$ to
mimic self-similar disk formation can still stall the instability when the
migration rate outpaces the growth of the instability.

Without an understanding of how model assumptions fit into a larger disk
evolution framework it is difficult to contextualize what a typical model
system looks like and determine a range of probable effects. In this work we
focus on a parameter exploration of a self-consistent model for the infall of
star-forming material in order to constrain the hallmarks of these model
systems.

The growing evidence for heterogeneous accretion onto cores both in
observations (Tobin et al., 2010; Yen et al., 2019; Pineda et al., 2020) and
numerical simulations (Smith et al., 2011; Kuffmeier et al., 2017; Kuznetsova
et al., 2020) suggests that we should not merely view anisotropic infall as a
potential mechanism for generating substructure, but rather acknowledge that
infall is likely to be anisotropic in some fashion and that this will
naturally lead to the formation of disk substructure, playing an important
role in disk evolution and planet formation.

Figure 12: The relative location of millimeter-sized dust rings and kinematic
deviations from the initial azimuthal gas velocity profile at a time of
$1.5\times 10^{4}\ \mathrm{yrs}$ post-infall, with the original locations of
$R_{\rm in}$ and $R_{\rm c}$ denoted by dotted annuli (left), are such that
the minimum of the dust gap, the location of which is marked in both panels by
a vertical gray dotted line, near the inner ring corresponds with the peak in
super-Keplerian velocity as expected for localized vortex structures (center),
but the apparent gap in front of the outer ring is centered where deviations
are zero, with the pattern more akin to that produced by spiral wakes
associated with a planetary companion (right). Bottom center and right panels
show sample radial profiles of the velocity deviation (light blue), with the
navy bold curve corresponding to the profile at the vortex azimuth $\phi_{\rm
v}$

As a caveat, the range of parameters used here is based on the results from
numerical simulations that, due to limitations in resolution, are biased
toward higher mass cores; thus, our choice of initial star and disk model is
kept constant, but consistent with the kinds of systems whose substructures
have been well-characterized in the radio (e.g. Andrews et al., 2018; Sheehan
et al., 2020) which also tend to be higher mass and relatively bright systems.

From our study, we have highlighted some of the main features of a system that
is undergoing or has undergone an infall episode: a broad gaseous bump across
the infall zone with a prominent drift barrier ($\eta>0$) at $R_{\rm in}$ due
to the presence of an inner vortex and a weaker barrier due to a vortex at
$R_{\rm c}$. Such vortices launch spiral waves, which can entrain dust in
narrower structures (Bae et al., 2015), and have been shown to produce
contrasts similar to spiral waves driven by low mass planets (up to tens of
Earth masses) (Huang et al., 2019).

Post-infall, differential radial drift contributes to the pile-up of mm-size
dust at the edges of the infall zone. For infall in the outer disk, this
results in the formation of rings outside of the nominal disk radius, $R_{\rm
d}$. In Figure 12, we show an example comparison of the locations of
millimeter-sized dust rings and gaps with velocity deviations in the gas. This
comparison demonstrates that vortices produce kinematic deviations like those
expected for planets (Disk Dynamics Collaboration et al., 2020), where spiral
wakes produce strong deviations on the order of $5\%$. Vortices rotate
opposite to the flows perturbed by planetary companions and thus their
kinematic signatures are expected to be easily differentiated from planetary
perturbers. While this is true for the vortices near the inner ring, the
velocity deviations near the outer ring at $R_{\rm c}$ are centered on the gap
between the two rings, potentially producing the appearance of deviations akin
to those from a planetary perturber.

The predicted appearance of dust substructures due to the presence of a gap-
opening planet has been well-described (Dong & Fung, 2017), and analysis of
the substructures in later phase sources from the DSHARP survey have yielded
an inferred population of embedded planets (Zhang et al., 2018) with planet
mass, $\alpha$, and $H/R$ used as free parameters. However, without accurate
tracers of both the dust and the gas, the presence of continuum rings alone is
not unique to the disk-planet interaction scenario. For instance, studies of
the molecular gas of the substructured source GM Auriga have found the
potential for the gravitational instability (Schwarz et al., 2021) and complex
structure in the surrounding gas pointing to late-stage infall (Huang et al.,
2021).

The one-dimensional dust evolution calculation presented here neglects non-
axisymmetry, but RWI vortices have been shown to be effective for trapping
grains of at least sub-mm sizes in non-axisymmetric structures (Meheut et al.,
2012; Barge et al., 2017). Rings and non-axisymmetric vortex structures are
likely to be most prominent for outer disk infall, given that the radial
extent of vortices is dependent on the scale height at their location. This
trend is consistent with the prevalence of non-axisymmetric features observed
largely on the outskirts of disks (van der Marel et al., 2020). Planetary
companions have been invoked to explain the presence of gaps and rings,
however generating the large scale structures seen in van der Marel et al.
(2020) with a planetary perturber requires sub-stellar mass companions
($M>15M_{\rm J}$).

###  4.2 Anisotropic infall makes small disks?

One could expect that as dust ring architectures are relatively common for the
annular infall modeled in this study, they should be prevalent among embedded
sources. However, the VANDAM survey of Class 0/I sources in Orion is dominated
by compact unstructured sources (Tobin et al., 2020); whereas substructures in
the sample appear exclusively as large-scale exterior rings (Sheehan et al.,
2020). The typical effective dust radii of these systems is less than $100\
\mathrm{au}$, with most single sources having inferred disk sizes less than
even $50\ \mathrm{au}$, smaller than that expected from a rotating collapse
scenario.

Figure 13:  surface density profiles (gray) for inner vs outer disk infall,
smoothed over $10\ \mathrm{au}$ (left). The effective disk size $R_{\rm eff}$
denotes where the mass enclosed is $68\%$ of the total mass. Dashed lines
denote surface density power-law profiles, where outer disk infall is best fit
by the initial $-1$ profile (black dashed line) and inner disk infall is more
consistent with a steeper $-3/2$ power-law (maroon dashed line). The mean
effective disk size is 55 au, systematically smaller than $R_{\rm d}$, with a
weak dependence on the $R_{\rm c}$ of the infall zone (right).

Looking at the evolution of millimeter grain surface density due to radial
drift (Figure 11), both outer and inner disk infall show similar relative
enhancement of the dust density at the pressure bumps. However, outer disk
infall, which occurs outside of the disk critical radius $R_{\rm d}$, deposits
material in an initially less dense part of the disk over a larger surface
area. A relative enhancement of a factor of two is still a fairly low absolute
surface density when it occurs in the outer disk, but a more significant
concentration of mass in the inner disk.

In Figure 13, we compare dust disk surface density profiles, smoothed over
$10\ \mathrm{au}$, for inner and outer disk infall. As an analogue to radii
inferred from continuum observations, we denote the location of the effective
radius $R_{\rm eff}$ as where the enclosed mass of the smoothed profile
contain $68\%$ of the total mass, such that the unperturbed initial surface
density profile has $R_{\rm eff}=R_{\rm d}$. Inner disk infall produces
systematically smaller disks as dust drift creates pile-ups near $R_{\rm c}$.
In addition, surface density profiles are more consistent with a steeper
power-law $\Sigma\propto R^{-3/2}$. Outer disk infall is less likely to affect
the overall fit of the surface density profile, tracing out the original disk
radius $R_{\rm d}$, as perturbations to the surface density profile then occur
well outside of $R_{\rm eff}$. Anisotropic infall, where inner and outer disk
infall is just as likely, accounting for radial drift, produces mean dust disk
sizes $R_{\rm eff}\sim 55\ \mathrm{au}$; consistent with an observed
population of protostellar sources dominated by compact dust disks.
Substructured disks in the VANDAM survey tend to be larger than the mean disk
sizes in the sample, given that their radii are influenced by the location of
the outer dust ring. As outer disk infall closer to $R_{\rm d}$ is less likely
to be resolved as a ring-like structure, large scale exterior substructures
like those in (Sheehan et al., 2020) should in fact be rarer among the
population of observed early phase objects.

###  4.3 Prospects for Disk Evolution

This study evolves the disk post-infall for several tens of kiloyears,
relatively short compared to the megayear timescales over which we expect
protoplanetary disks to evolve. The results of this study provide a starting
point for further investigation into the long-term consequences of anistropic
infall on disk evolution.

Notably, our results show a variety of radial profiles of the $\alpha$
parameter across the disk, increasing outside the infall zone for outer disk
infall and decreasing for inner disk infall (Figure 6). In the steady-state
$\alpha$-disk model , which assumes a constant disk-wide $\alpha$ to
parameterize viscosity, the viscous timescale $t_{\nu}\sim R^{2}/(\alpha
c_{\rm s}H)$ depends on radius as $t_{\nu}\sim R$ for a disk with a
temperature structure set by stellar heating $T\propto R^{1/2}$. Given the
radial dependence in $\alpha$ that we find for outer disk infall, one could
expect the same viscous timescale at all radii, compared to inner disk infall
which would only have comparable transport times interior to the infall zone.

Bae et al. (2014) has shown that isotropic infall could power accretion
outbursts on short timescales during disk formation. Accelerated transport
across the whole disk due to anisotropic infall in the outer disk has
consequences for the protostellar accretion and outburst behavior of early
phase sources. Characterizing transport and outburst behavior due to
anisotropic infall will require treatment of heating and cooling and the
inclusion of disk self-gravity to capture fragmentation.

Within the parameter space and initial disk parameters explored here, we do
not expect disk self-gravity to play an important role. At the end of infall,
the Toomre Q parameter is well above the marginal stability criterion. In our
study, by the time Q reaches a value at which self-gravity could affect the
stability of high mode number vortices $Q<R/H$ (Lovelace & Hohlfeld, 2013),
the instability has already saturated. However, the chosen $\dot{M}_{\rm in}$
is marginal, and higher infall rates will be in the regime at which self-
gravity will need to be considered for both the RWI and gravitational
instability.

Bae et al. (2015) show that quenching of the infall-induced RWI and
dissipation of vortex structures can occur in the self-gravitating regime.
Population level predictions for a full range of evolutionary outcomes, such
as those originally put forth in Kratter et al. (2008), including the effects
of anisotropic infall will require including disk self-gravity and expanding
the parameter space.

Even as anisotropic infall creates regions of enhanced angular momentum
transport, the results demonstrate that it also sets up barriers to dust
transport by drift, robustly trapping most particle sizes at $R_{\rm in}$.
Over longer timescales, the pressure bump at $R_{\rm in}$ could effectively
segregate dust populations in the disk. Multiple episodes of infall from
different accretion reservoirs may lead to compositional heterogeneity within
the disk, a proposed mechanism for accounting for the isotopic variations in
solar system chondrites (Piani et al., 2021). Investigating how anisotropic
infall may affect the long-term compositional evolution of protoplanet bodies
will require more sophisticated modeling of the dust evolution, accounting for
growth and fragmentation over longer timescales.

##  5 Summary

Growing evidence from both observations and numerical simulations suggests
that disks are fed by anisotropic streams of infalling material, counter to
the self-similar theoretical approach for disk formation and evolution. We
performed a parameter study of anistropic infall onto disks, with model
parameters informed by the properties of protostellar core accretion events
from large-scale star cluster formation simulations (Kuznetsova et al., 2020).
The parameter space subtends accretion flow geometry from the nearly spherical
to the stream-like and landing zones within the inner disk and in the outer
disk. Results show that anisotropic infall significantly perturbs the disk,
generating perturbations strong enough to incite the Rossby wave instability
within a few dynamical times at the outermost radius of infall.

Infalling material generates shear stresses, with the highest measured
viscosity parameter $\alpha\sim 10^{-2}$ occurring in cases of outer disk
infall. How locally generated mechanical turbulent stresses affect long-term
disk transport and evolution and potentially induce outbursts is a subject of
further study.

Infall-induced azimuthal shear is robust near the edges of the infall zone,
allowing entrapment of the largest dust particles, with the potential for
accelerated planet formation at early times. We post-process simulations,
calculating the one-dimensional dust evolution due to radial drift. We find
that ring formation due to differential drift can occur rapidly in the post-
infall phase, with extended substructures likely for infall occurring in the
outer disk. Infall in the inner disk, equally likely in our theoretical
framework, accumulates dust in the infall zone, resulting in a compact dust
disk, consistent with measurements of small sub-millimeter dust radii in
embedded sources.

The anisotropic infall model presented here results in the ready formation and
reinforcement of pressure bumps in the disk during the earlier phases of
protostellar disk evolution. The presence of multiple drift barriers likely
has consequences on the long term evolution of the disk dynamics and on the
bulk composition and sizes of planetary bodies segregated by drift barriers.

Support for A.K. was provided by NASA through the NASA Hubble Fellowship grant
#HST-HF2-51463.001-A awarded by the Space Telescope Science Institute, which
is operated by the Association of Universities for Research in Astronomy,
Incorporated, under NASA contract NAS5-26555. M.-M.M.L. was partially
supported by NSF grant AST18-15461. This work used the HPC _Mendel_ cluster at
the American Museum of Natural History developed with National Science
Foundation (NSF) Campus Cyberinfrastructure support through Award#1925590 and
the Extreme Science and Engineering Discovery Environment (XSEDE) _Bridges_ at
the Pittsburgh Supercomputing Center under allocation AST190025, which is
supported by NSF grant number ACI-1548562.

##  Appendix A Effect of Numerical Resolution

The simulations in this study were performed with a resolution of $N_{\rm
R}=2048,N_{\phi}=1024$, on a logarithmic radial grid corresponding to a
constant $\Delta R/R\sim 0.0027$. Here, we compare radially logarithmic grids
of varying resolution, $N_{\rm R}=256,512,1024,2048$ with
$N_{\phi}=128,256,512,1024$, respectively, for the fiducial infall model used
above: $R_{\rm in}=30\ \mathrm{au},R_{\rm c}=60\ \mathrm{au}$. When comparing
the development of the RWI over time using the maximum relative non-
axisymmetric perturbation strength, as in Figure 2, we find that the
resolution affects the launch time of the instability, as well as the
perturbation strength during the non-linear phase (Figure 14).

Figure 14: The evolution of the maximum non-axisymmetric surface density
perturbation over the simulation time, as in Figure 2, for the numerical
resolutions tested, with the fiducial run ($N_{R}=2048$) shown as a bold black
line. The launch time of the instability is delayed for lower resolutions. The
inset region marked by the gray box in the left panel is magnified on the
right and shown on a linear scale.

Notably, for the lowest resolution shown for comparison $N_{\rm
R}=512,N_{\phi}=256$, the perturbation strength decays after infall ends. At
that resolution, the RWI at $R_{\rm c}$ is quenched, producing no long-lived
vortices. This effect can be seen in the evolution of the viscosity parameter
over time for the lowest resolution case, wherein the disk-wide $\alpha$ does
not see the same secondary growth period after $4\ \mathrm{kyr}$ corresponding
to the growth of the RWI at $R_{\rm c}$. Compared to the overall trend of the
decay of the disk-wide $\alpha$ of the fiducial run, shown as a smoothed
trend-line in Figure 15, the deviations introduced by the time variation of
the fiducial disk-wide $\alpha$ is usually on the level of 10–25%. The
magnitude of errors introduced at lower resolutions can be twice that of the
noise level in the fiducial case.

Figure 15: The time evolution of the disk-wide $\alpha$ parameter at different
numerical resolutions, with the fiducial run ($N_{R}=2048$) shown as a bold
black line. The inset region marked by the gray box in the left panel is
magnified on the right and shown on a linear scale. (inset) A smoothed trend
line corresponding to the fiducial run is overplotted as a dashed white line,
with contours for $10\%,25\%,$ and $50\%$ deviation from the trend shown in
shades of violet.

Across the metrics introduced in Figure 9, we show in Figure 16 that quenching
of the RWI by low numerical resolution at the outer radius results in a
relatively lower surface density (and pressure) bump near $R_{\rm c}$. At
higher resolutions, the relative perturbation in the pressure gradient
compared to its initial value $\eta-\eta_{0}$ is not strongly affected by
resolution.

However, the value of the large scale pressure gradient represented by $\eta$
does depend on how well the pressure gradient is resolved. Though $\eta$ sets
the value of the dust drift velocity $u_{\rm d}$ and associated drift
timescales, the formation of dust gaps and rings depends much more on the
differential drift than it does the absolute value of the drift velocity. The
primary effect of under-resolving the pressure gradient is realized at smaller
radii, where the direction of the pressure gradient is reversed, resulting in
artificially super-Keplerian flow in the inner disk where $\eta>0$.

Figure 16:  (left) Metrics of disk dynamics shown in Figure 9 across different
numerical resolutions selected during the infall period, $t<t_{\rm ep}$, from
top to bottom: surface density, pressure gradient parameter $\eta$, and
measured $\alpha$ parameters. (right) Deviations in same metrics, from top to
bottom: the perturbation to the surface density
$\delta\Sigma=\Sigma/\Sigma_{0}-1$, perturbation to the pressure gradient
$\eta-\eta_{0}$, and perturbation relative to the disk-wide $\alpha$ parameter
measured for the fiducial run at the selected time, $\alpha_{\rm
fid}=3.86\times 10^{-4}$.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership, T., Brogan, C., Pérez, L., et al. 2015, The Astrophysical Journal Letters, 808, L3
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, The Astrophysical Journal Letters, 869, L41
  * Bae et al. (2015) Bae, J., Hartmann, L., & Zhu, Z. 2015, The Astrophysical Journal, 805, 15
  * Bae et al. (2014) Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, The Astrophysical Journal, 795, 61
  * Bae et al. (2018) Bae, J., Pinilla, P., & Birnstiel, T. 2018, ApJ, 864, L26, doi: 10.3847/2041-8213/aadd51
  * Bae et al. (2017) Bae, J., Zhu, Z., & Hartmann, L. 2017, The Astrophysical Journal, 850, 201
  * Barge et al. (2017) Barge, P., Ricci, L., Carilli, C. L., & Previn-Ratnasingam, R. 2017, Astronomy & Astrophysics, 605, A122
  * Benisty et al. (2021) Benisty, M., Bae, J., Facchini, S., et al. 2021, The Astrophysical Journal Letters, 916, L2
  * Benítez-Llambay & Masset (2016) Benítez-Llambay, P., & Masset, F. S. 2016, The Astrophysical Journal Supplement Series, 223, 11
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C., & Brauer, F. 2010, Astronomy & Astrophysics, 513, A79
  * Birnstiel et al. (2016) Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Science Reviews, 205, 41
  * Carrera et al. (2021) Carrera, D., Thomas, A., Simon, J. B., et al. 2021, arXiv preprint arXiv:2108.08315
  * Cassen & Moosman (1981) Cassen, P. M., & Moosman, A. 1981, Icarus, 48, 353
  * De Val-Borro et al. (2006) De Val-Borro, M., Edgar, R., Artymowicz, P., et al. 2006, Monthly Notices of the Royal Astronomical Society, 370, 529
  * Disk Dynamics Collaboration et al. (2020) Disk Dynamics Collaboration, A., Armitage, P. J., Bae, J., et al. 2020, arXiv preprint arXiv:2009.04345
  * Dong & Fung (2017) Dong, R., & Fung, J. 2017, The Astrophysical Journal, 835, 146
  * Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, The Astrophysical Journal, 843, 127
  * Gundlach & Blum (2015) Gundlach, B., & Blum, J. 2015, The Astrophysical Journal, 798, 34
  * Huang et al. (2021) Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, arXiv e-prints, arXiv:2109.06224.  https://arxiv.org/abs/2109.06224
  * Huang et al. (2019) Huang, P., Dong, R., Li, H., Li, S., & Ji, J. 2019, The Astrophysical Journal Letters, 883, L39
  * Keppler et al. (2018) Keppler, M., Benisty, M., Muller, A., et al. 2018, Astronomy and Astrophysics, 617
  * Kratter et al. (2008) Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, The Astrophysical Journal, 681, 375
  * Kuffmeier et al. (2017) Kuffmeier, M., Haugboelle, T., & Nordlund, A. 2017, The Astrophysical Journal, 846, 7
  * Kuznetsova et al. (2019) Kuznetsova, A., Hartmann, L., & Heitsch, F. 2019, The Astrophysical Journal, 876, 33
  * Kuznetsova et al. (2020) —. 2020, The Astrophysical Journal, 893, 73
  * Lesur et al. (2015) Lesur, G., Hennebelle, P., & Fromang, S. 2015, Astronomy & Astrophysics, 582, L9
  * Lovelace & Hohlfeld (2013) Lovelace, R., & Hohlfeld, R. 2013, Monthly Notices of the Royal Astronomical Society, 429, 529
  * Lovelace et al. (1999) Lovelace, R., Li, H., Colgate, S., & Nelson, A. 1999, The Astrophysical Journal, 513, 805
  * Malik et al. (2015) Malik, M., Meru, F., Mayer, L., & Meyer, M. 2015, The Astrophysical Journal, 802, 56
  * Meheut et al. (2012) Meheut, H., Meliani, Z., Varniére, P., & Benz, W. 2012, Astronomy and Astrophysics, 545
  * Ono et al. (2018) Ono, T., Muto, T., Tomida, K., & Zhu, Z. 2018, The Astrophysical Journal, 864, 70
  * Ormel & Cuzzi (2007) Ormel, C., & Cuzzi, J. 2007, Astronomy & Astrophysics, 466, 413
  * Paardekooper et al. (2010) Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. 2010, The Astrophysical Journal, 725, 146
  * Piani et al. (2021) Piani, L., Marrocchi, Y., Vacher, L. G., Yurimoto, H., & Bizzarro, M. 2021, Earth and Planetary Science Letters, 567, 117008
  * Pineda et al. (2020) Pineda, J. E., Segura-Cox, D. M., Caselli, P., et al. 2020, Nature Astronomy, 4, 1158
  * Pinilla et al. (2012) Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, Astronomy & Astrophysics, 538, A114
  * Pinte et al. (2019) Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nature Astronomy, 3, 1109
  * Schwarz et al. (2021) Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, arXiv preprint arXiv:2109.06228
  * Sheehan et al. (2020) Sheehan, P. D., Tobin, J. J., Federman, S., Megeath, S. T., & Looney, L. W. 2020, The Astrophysical Journal, 902, 141
  * Smith et al. (2011) Smith, R. J., Glover, S. C. O., Bonnell, I. A., Clark, P. C., & Klessen, R. S. 2011, Monthly Notices of the Royal Astronomical Society, 411, 1354
  * Teague et al. (2018) Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, The Astrophysical Journal Letters, 860, L12
  * Terebey et al. (1984) Terebey, S., Shu, F. H.-S., & Cassen, P. M. 1984, The Astrophysical Journal, 286, 529
  * Tobin et al. (2010) Tobin, J. J., Hartmann, L., Chiang, H.-F., & Looney, L. W. 2010, The Astrophysical Journal, 712, 1010
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, The Astrophysical Journal, 890, 130
  * Ulrich (1976) Ulrich, R. K. 1976, The Astrophysical Journal, 210, 377
  * van der Marel et al. (2020) van der Marel, N., Birnstiel, T., Garufi, A., et al. 2020, The Astronomical Journal, 161, 33
  * Vorobyov & Basu (2010) Vorobyov, E. I., & Basu, S. 2010, The Astrophysical Journal, 719, 1896
  * Weidenschilling (1977) Weidenschilling, S. 1977, Monthly Notices of the Royal Astronomical Society, 180, 57
  * Yen et al. (2019) Yen, H.-W., Gu, P.-G., Hirano, N., et al. 2019, The Astrophysical Journal, 880, 69
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, The Astrophysical Journal Letters, 869, L47