# Dust Settling and Clumping in MRI Turbulent Outer Protoplanetary Disks

Ziyan Xu Kavli Institute for Astronomy and Astrophysics, Peking University,
100871 Beijing, China  Department of Astronomy, Peking University, 100871
Beijing, China  Xue-Ning Bai Institute for Advanced Study, Tsinghua
University, 100084 Beijing, China  Department of Astronomy, Tsinghua
University, 100084 Beijing, China

###### Abstract

Planetesimal formation is a crucial yet poorly understood process in planet
formation. It is widely believed that planetesimal formation is the outcome of
dust clumping by the streaming instability (SI). However, recent analytical
and numerical studies have shown that the SI can be damped or suppressed by
external turbulence, and at least the outer regions of protoplanetary disks
are likely weakly turbulent due to magneto-rotational instability (MRI). We
conduct high-resolution local shearing-box simulations of hybrid particle-gas
magnetohydrodynamics (MHD), incorporating ambipolar diffusion as the dominant
non-ideal MHD effect, applicable to outer disk regions. We first show that
dust backreaction enhances dust settling towards the midplane by reducing
turbulence correlation time. Under modest level of MRI turbulence, we find
that dust clumping is in fact easier than the conventional SI case, in the
sense that the threshold of solid abundance for clumping is lower. The key to
dust clumping includes dust backreaction and the presence of local pressure
maxima, which in our work is formed by the MRI zonal flows overcoming
background pressure gradient. Overall, our results support planetesimal
formation in the MRI-turbulent outer protoplanetary disks, especially in ring-
like substructures.

protoplanetary disks - planet formation - planetesimals - magnetohydrodynamics

##  1 INTRODUCTION

The process of planet formation encompasses a growth of more than 13 orders of
magnitude in size. At the lower end, micron-sized dust particles stick
together and coagulate into mm to cm sized solids. At the higher end, km-sized
planetesimals grow into planetary cores, and eventually into terrestrial or
giant planets of up to $10^{5}$ km in size. The intermediate stage involves
the formation of km-sized planetesimals out of from mm-cm sized particles, and
it is perhaps the least understood problem in planet formation (e.g., see
review by Chiang & Youdin (2010)).

###  1.1 Overview of planetesimal formation theory

The obstacle in planetesimal formation arises primarily because of the
existence of growth barriers. First, direct growth of dust beyond $\sim$mm-cm
size has been found to be difficult, and encounters fragmentation and bouncing
barriers (e.g., Blum & Wurm, 2008; Güttler et al., 2010; Zsom et al., 2010).
While various scenarios have been proposed to mitigate the problem (e.g.,
Windmark et al., 2012; Okuzumi et al., 2012; Kataoka et al., 2013), they
generally rely on not-well-known dust surface properties, and more recent
models can yield very different conclusions (e.g., Okuzumi & Tazaki, 2019).
Second, even dust can grow larger, they suffer from the so-called “meter-size
barrier” (Weidenschilling, 1977), or more generally the drift barrier (e.g.,
Birnstiel et al., 2012), so that dust must grow sufficiently quickly to bypass
this barrier before they are lost through radial drift.

It is thus generally believed that planetesimal formation must proceed through
some collective effect, where the responsible physical mechanism first gathers
dust particles into dense clumps, which then collapse under their own self-
gravity. Under this framework, the key to understanding planetesimal formation
is to identify and characterize mechanisms that can strongly concentrate
particles, typically of mm-cm in size as the end product of grain growth.

Particles with a broad range of sizes interact with gaseous protoplanetary
disks (PPDs) via gas drag, characterized by the stopping time $t_{\rm stop}$.
It is conventional to cast the stopping time in dimensionless form as
$\tau_{s}\equiv\Omega t_{\rm stop}$, where $\Omega$ is the local orbital
frequency of the disk. Particles are considered strongly coupled to the gas
for $\tau_{s}\ll 1$, and the coupling is marginal for $\tau_{s}\sim 1$. Such
interactions lead to the aforementioned radial drift, where an outward radial
gas pressure gradient typically present in PPDs is felt by the gas but not
dust. The gas thus rotates at sub-Keplerian velocity, while particles feel a
head wind and hence drift radially inward. More generically, particles always
drift towards higher pressure, and the effect is the most prominent for
marginally coupled particles with $\tau_{s}\sim 1$. Particles with mm-cm sizes
are typically marginally coupled with $\tau_{s}\sim 0.1-1$ in the outer PPDs
(a few tens of AU), and more strongly coupled in the inner disks ($\sim
10^{-3}$ in the inner few AU). Therefore, mutual interaction between dust and
gas is the key piece of physics in any theory of planetesimal formation.

Multiple planetesimal formation mechanisms have been proposed. Early studies
focused on gravitational instability (GI) of a thin dust layer (Safronov &
Zvjagina, 1969; Goldreich & Ward, 1973). However, it requires extreme level of
settling, and the dust layer itself can be subject to the Kelvin-Helmholtz
instability (KHI, Weidenschilling, 1980). Triggering the GI requires
unrealistically high dust-to-gas ratio (Sekiya, 1998; Youdin & Shu, 2002),
although it can be mitigated to less extreme conditions (Chiang, 2008;
Barranco, 2009; Lee et al., 2010a, b). For strongly-coupled dust,
concentration may be achieved through secular GI (e.g., Youdin, 2011;
Takahashi & Inutsuka, 2014; Tominaga et al., 2018) and two-component viscous
GI (Tominaga et al., 2019), resulting from gas drag and self-gravity over
secular timescales, with the latter further involving the action of gas
viscosity. They generally require very low level of external turbulence
(viscosity) to operate. On the other hand, turbulence contains eddies of
different sizes which themselves may selectively concentrate particles whose
stopping time is comparable to eddy turnover time (e.g., Cuzzi et al., 2001;
Pan et al., 2011; Hopkins, 2016; Hartlep & Cuzzi, 2020). In addition, specific
hydrodynamic instabilities can also result in strong vortices that act as dust
traps (e.g., Raettig et al., 2021).

Over the past decade, perhaps the most popular and promising scenario for
planetesimal formation is the streaming instability (SI, Goodman & Pindor,
2000; Youdin & Goodman, 2005). The key ingredient behind the SI is the
incorporation of backreaction from dust to gas, leading to a two-fluid type
instability. It arises at the cost of the free energy from relative drift
between dust and gas (which ultimately arises from the pressure gradient in
the disk) rather than orbital shear. It is commonly interpreted as dust pile
up causing traffic jams (Youdin & Johansen, 2007; Johansen & Youdin, 2007),
more recent investigations have revealed the richness of the instability in
its linear behaviors (Jacquet et al., 2011; Lin & Youdin, 2017; Jaupart &
Laibe, 2020; Pan & Yu, 2020; Pan, 2020; Lin, 2021), and can be considered as a
special case of the resonant drag instability (Squire & Hopkins, 2018, 2020).

Non-linear simulations particle-gas dynamics in (vertically-stratified) PPDs
have revealed that upon incorporating dust backreaction, the SI leads to the
formation of dense dust clumps with sufficiently high density to drive
gravitational collapse (Johansen et al., 2009b, 2014; Simon et al., 2016;
Schäfer et al., 2017; Abod et al., 2019). The strength of particle
concentration depends on physical conditions in PPDs, particularly particle
stopping time, the solid abundance, and the background gas pressure gradient
in the disk (Bai & Stone, 2010a, c; Carrera et al., 2015; Yang et al., 2017).
With these simulation results, the SI scenario has been widely applied to
model planetesimal formation in global disk context (e.g. Drażkowska &
Dullemond, 2014; Schoonenberg & Ormel, 2017).

###  1.2 The issue with disk turbulence

The aforementioned studies of the SI typically consider dust dynamics without
external turbulence. However, recent studies have shown that the linear growth
of the SI can be damped or fully suppressed by even moderate level (Shakura-
Sunyaev $\alpha\gtrsim 10^{-4}$, Shakura & Sunyaev, 1973) of disk viscosity
(Umurhan et al., 2020; Chen & Lin, 2020), a result that has been qualitatively
confirmed from non-linear SI simulations with external driven turbulence (Gole
et al., 2020). This is a serious problem that endangers the applicability of
the SI to planetesimal formation.

Observational inference of disk turbulence have been inconclusive, where there
is evidence for both strong and weak turbulence. Direct measurements of
turbulent line broadening have yielded either strong (Flaherty et al., 2020)
or non-detectable (Flaherty et al., 2017, 2018; Teague et al., 2018) level of
turbulence. Level of turbulence can also be inferred from measuring dust
distributions. In particular, dust layer thickness is determined by a balance
of dust settling from vertical gravity and turbulent diffusion. In the case of
the HL Tau disk (ALMA Partnership et al., 2015), the sharpness of dust rings
indicates a geometrically thin dusty disk corresponding to weak turbulence
with $\alpha\sim$ a few times $10^{-4}$ (Pinte et al., 2016). Recent disk
observations have also revealed rich substructures dominated by dust rings
(e.g. Andrews et al., 2018; Long et al., 2018), where the width of dust rings
result from dust trapping balancing turbulent diffusion. Modeling and
measuring dust and gas ring width have led to a conclusion of relatively
strong turbulence with $\alpha/\tau_{s}\gtrsim 10^{-2}$, at least in disk
outer regions (Dullemond et al., 2018; Rosotti et al., 2020). Despite some
controversies, the inferred turbulence levels are typically weak, but still
well above $\alpha\sim 10^{-4}$, which would have suppressed the SI.

An important caveat behind the conundrum above is that turbulence may not be
simply approximated as a viscosity, nor mimicked by artificial turbulence
driving. Depending on the source of turbulence driving, it may share
substantially different properties. For instance, it was found that the SI can
coexist with the vertical shear instability (VSI, Nelson et al., 2013) which,
despite yielding a turbulence level of $\alpha\sim 10^{-3}$, can
constructively act with the SI to promote dust clumping (Schäfer et al.,
2020). Therefore, it is of crucial importance to study the interplay of the SI
with more realistic turbulence.

###  1.3 This work

In this work, we focus on the outer regions of PPDs, which are accessible to
spatially-resolved observations. The most widely invoked mechanism for
generating turbulence is the magnetorotational instability (MRI, Balbus &
Hawley (1991)), whereas long as gas and magnetic fields are well coupled, or
in the ideal magnetohydrodynamics (MHD) regime, it generates vigorous
turbulence that efficiently transports angular momentum and drives disk
accretion. However, PPDs are weakly ionized, making gas and magnetic fields
poorly coupled, resulting in three non-ideal MHD effects, namely, Ohmic
resistivity, the Hall effect, and ambipolar diffusion (AD) (e.g. Wardle, 2007;
Bai, 2011a; Armitage, 2015; Lesur, 2020). The three effects dominate in
progressively lower density regions, where the outer disk is largely dominated
by AD. Simulations have revealed that the MRI is damped by AD (Bai & Stone,
2011), but not fully suppressed given the expected level of ionization in the
outer disk (Simon et al., 2013b, a; Bai, 2015), with $\alpha$ on the order of
$10^{-3}$. It is our goal to study the interplay of the SI with the MRI in
this regime.

We note that there are several similar studies in the literature considering
particle concentration in MRI turbulence. Majority of the studies considered
the MRI in ideal MHD, with low resolution (typically $\lesssim 16$ cells per
scale height), well below typical resolutions needed to properly resolve the
SI (e.g. Balsara et al., 2009; Tilley et al., 2010). Exceptions include
Johansen et al. (2007, 2011), who conducted high-resolution local MRI
simulations in ideal MHD and found particle concentration in local pressure
bumps (known as zonal flows, Johansen et al. (2009a)), leading to planetesimal
formation. Such zonal flows can be long-lived and likely play a significant
role in concentrating particles (Dittrich et al., 2013). Yang et al. (2018)
conducted large-box (but relatively low resolution) simulations with Ohmic
resistivity in the midplane region (and hence modest turbulent damping), and
also reported particle clumping despite that particles do not strongly settle
to the midplane region. However, we anticipate the properties of the MRI
turbulence to be very different in outer PPDs in the presence of AD.

Here, we perform very high resolution 3D simulations of particle-gas mixture
in the local shearing-box framework. Both hydrodynamic and non-ideal MHD with
AD are considered. Our simulations have both reasonable box size to
accommodate the MRI turbulence and sufficient resolution to reasonably resolve
the SI, if present, which allow us to study the interplay between realistic
MRI turbulence and dust dynamics in the outer PPDs. We are particularly
interested in particle settling and diffusion in vertical direction, as well
as particle clumping. By comparison between runs with and without magnetic
field, as well as with and without dust backreaction, we aim to decipher the
potential mechanisms that control dust diffusion, settling and concentration,
and whether the SI still operates to form planetesimals.

This paper is organized as follows. In §2, we describe basic formulations and
simulation setup. In §3, we overview the simulation results. In §4, §5, and
§6, we analyze the simulation results, focusing on particle vertical settling
and clumping respectively. The results are further analyzed and discussed in
§7 to address the physics behind particle settling behaviour and particle
clumping, as well as implications in realistic disks. We summarize and discuss
future prospects in §8.

Figure 1: A snapshot of 3D visualization of the dust spatial distribution in
our non-ideal MHD shearing-box simulation (run ADZ2) at the time of
$T=120\Omega_{K}^{-1}$ after inserting dust particles, with projections in
three orthogonal directions. Dust density is mapped to brightness, where
brighter region indicates higher dust density.

##  2 Method and Simulations

###  2.1 Formulation and simulation framework

We use the Athena MHD code (Stone et al., 2008), a higher-order Godunov code
with constrained transport to preserve the divergence-free property for
magnetic fields (Gardiner & Stone, 2005, 2008) to perform 3D local MHD
simulations in the shearing sheet approximation (Goldreich & Lynden-Bell,
1965). We use a local reference frame located at a fiducial radius corotating
with Keplerian frequency $\Omega_{K}$. The dynamical equations are written in
Cartesian coordinates, with $\hat{\boldsymbol{x}}$, $\hat{\boldsymbol{y}}$ and
$\hat{\boldsymbol{z}}$ denoting unit vectors pointing to the radial, azimuthal
and vertical directions, and ${\boldsymbol{\Omega}}_{K}$ is along the
$\hat{\boldsymbol{z}}$ direction. The gas density and velocity are denoted by
$\rho_{g}$ and ${\boldsymbol{u}}$ separately in this non-inertial frame. The
HLLD Riemann solver (Miyoshi & Kusano, 2005) with third order reconstruction,
combined with the CTU integrator are used in all MHD simulations. Dust is
implemented as Lagrangian particles coupled with gas via aerodynamic drag, and
the backreaction from the particles to gas is included, as implemented in (Bai
& Stone, 2010b). The drag force between gas and dust particles is
characterized by particle stopping time $t_{\rm stop}$. For a particle $i$
with velocity of ${\boldsymbol{v}_{i}}$, the drag force on the particle is
$({\boldsymbol{u}}-{\boldsymbol{v}_{i}})/t_{\rm stop}$ per particle unit mass.
With magnetic field denoted by ${\boldsymbol{B}}$, the MHD equations are
written as:

$\frac{\partial\rho_{g}}{\partial t}+\nabla\cdot(\rho_{g}{\boldsymbol{u}})=0,$
(1) $\begin{split}&\frac{\partial\rho_{g}{\boldsymbol{u}}}{\partial
t}+\nabla\cdot(\rho_{g}{\boldsymbol{u}}^{T}{\boldsymbol{u}}+{\sf T})\\\
&=\rho_{g}[2{\boldsymbol{u}}\times{\boldsymbol{\Omega}_{K}}+3\Omega_{K}^{2}x{\hat{\boldsymbol{x}}}+\epsilon\frac{{\bar{\boldsymbol{v}}}-{\boldsymbol{u}}}{t_{\rm
stop}}],\end{split}$ (2) $\frac{\partial{\boldsymbol{B}}}{\partial
t}=\nabla\times\left[{\boldsymbol{u}}\times{\boldsymbol{B}}+\frac{({\boldsymbol{J}}\times{\boldsymbol{B}})\times{\boldsymbol{B}}}{\gamma\rho_{i}\rho_{g}}\right],$
(3)

where ${\sf T}$ denotes the total stress tensor:

${\sf T}=(P+B^{2}/2){\sf I}-{\boldsymbol{B}}^{T}{\boldsymbol{B}},$ (4)

with ${\sf I}$ being the identity tensor and $P$ is the gas pressure. We
assume an isothermal equation of state for the gas so that
$P=\rho_{g}c_{s}^{2}$, where $c_{s}$ is the isothermal sound speed. The disk
scale height is given by $H_{g}=c_{s}/\Omega_{K}$. The units for magnetic
field is chosen such that magnetic permeability $\mu=1$, and hence magnetic
pressure is given by $B^{2}/2$, and current density is given by
${\boldsymbol{J}}=\nabla\times{\boldsymbol{B}}$.

The first and second terms on the right side of equation (LABEL:eq.MHD2)
represent Coriolis force and radial tidal gravity, separately. The third term
on the right side of (LABEL:eq.MHD2) is the backreaction from dust particles
to gas where $\epsilon$ and ${\bar{\boldsymbol{v}}}$ are the local particle
mass density and mean particle velocity, obtained by interpolation. The disk
vertical gravity in the gas is ignored in our simulations and hence the gas is
unstratified (but dust is stratified, see later). This is justified as we are
mainly interested in regions within one gas scale height about the midplane.

In the induction Equation (3), the last term corresponds to ambipolar
diffusion (AD), where $\gamma$ is the coefficient of momentum exchange in ion-
neutral collisions and $\rho_{i}$ denotes the ion density. AD is paraterized
by Elsässer number:

$Am=\frac{\gamma\rho_{i}}{\Omega_{K}},$ (5)

with $Am\to\infty$ corresponds to ideal MHD regime, where the magnetic field
is frozen to the gas. The gas dynamics in the disk is significantly affected
by AD with $Am\lesssim 10$ (Bai & Stone, 2011). Bai (2011a, b) found that
$Am\sim 1$ is widely applicable at the outer PPDs ($\gtrsim$ 30 AU) where AD
is the dominant non-ideal MHD effect.

Without stratification, initial gas density $\rho_{0}$ is uniform. We thus set
$c_{s}=\Omega_{K}=H_{g}=\rho_{0}=1$ in code units. Periodic boundary
conditions are used in both vertical and azimuthal directions, while in the
radial direction, we use shearing periodic boundary conditions (Hawley et al.,
1995). An orbital advection scheme is used which improves the accuracy of the
simulations (Stone & Gardiner, 2010), in which we subtract Keplerian shear
${\boldsymbol{u}}_{K}\equiv-(3/2)\Omega_{K}x\hat{\boldsymbol{y}}$ and only
evolve the deviation
${\boldsymbol{u}}^{\prime}\equiv{\boldsymbol{u}}-{\boldsymbol{u}}_{K}$.

Fiducially, we impose a net vertical magnetic field $B_{0}$ in our MHD
simulations. The strength of $B_{0}$ is characterized plasma $\beta$,

$\beta_{0}\equiv\frac{\rho_{0}c_{s}^{2}}{B_{0}^{2}/2},$ (6)

corresponding to the ratio of gas pressure $\rho_{0}c_{s}^{2}$ to the magnetic
pressure of the net vertical field. The net vertical magnetic field is needed
for the MRI turbulence to be sustained under strong AD in outer disk regions
(Bai & Stone, 2011; Simon et al., 2013b), and $\beta_{0}\sim 10^{4}$ is
required for achieving accretion rates consistent with observations (Simon et
al., 2013a; Bai, 2015).

The equation of motion for particle $i$ can then be written as

$\frac{d{\boldsymbol{v}}_{i}}{dt}=-2\eta
v_{K}\Omega_{K}{\hat{\boldsymbol{x}}}+2{\boldsymbol{v}}_{i}\times{\boldsymbol{\Omega}_{K}}+3\Omega_{K}^{2}x_{i}{\hat{\boldsymbol{x}}}-\Omega_{K}^{2}z_{i}{\hat{\boldsymbol{z}}}-\frac{{{\boldsymbol{v}}_{i}}-{\boldsymbol{u}}}{t_{\rm
stop}}.$ (7)

In the above, the first term on the right side $-2\eta
v_{K}\Omega_{K}{\hat{\boldsymbol{x}}}$ represents a force pointing inward in
order to mimic the outward gas radial pressure gradient (see Bai & Stone
(2010b)), which is equivalent to but more convenient than directly
implementing the pressure gradient to the gas in Equation (LABEL:eq.MHD2). The
parameter $\Delta v_{K}\equiv\eta v_{K}$ represents the deviation from
Keplerian rotation due to radial pressure gradient, which is characterized by
$\Pi\equiv\eta v_{K}/c_{s}$ (Bai & Stone, 2010c). In the absence of gas drag,
$\Delta v_{K}$ gives the relative velocity between gas and dust, i.e., the
headwind. With our treatment, gas would stay at Keplerian while dust would
travel at super-Keplerian speeds, maintaining the relative drift between the
two components.

Table 1: List of simulation runs.

Run name | $\Pi$ | $Z$ | $N_{x}\times N_{y}\times N_{z}$ | $\beta$ | $Am$
---|---|---|---|---|---
hydZ2 | 0.1 | 0.02 | 480 $\times$ 480 $\times$ 480 | $\infty$ | -
hydZ4 | 0.1 | 0.04 | 480 $\times$ 480 $\times$ 480 | $\infty$ | -
ADZ0 | 0.1 | $10^{-5}$ | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ05 | 0.1 | 0.005 | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ1 | 0.1 | 0.01 | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ2 | 0.1 | 0.02 | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ4 | 0.1 | 0.04 | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ1H | 0.1 | 0.01 | 384 $\times$ 384 $\times$ 384 | 12800 | 2
ADZ2H | 0.1 | 0.02 | 384 $\times$ 384 $\times$ 384 | 12800 | 2
ADZ2P0 | 0 | 0.02 | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ2P005 | 0.05 | 0.02 | 256 $\times$ 256 $\times$ 256 | 12800 | 2
ADZ2P02 | 0.2 | 0.02 | 256 $\times$ 256 $\times$ 256 | 12800 | 2

  * •

All simulations have a particle number of 1 per cell on average, with particle
stopping time $\tau_{s}=0.1$.

  * •

The simulation box size is $H\times 2H\times H$ for MHD runs, and $0.6H\times
1.2H\times 0.6H$ for hydrodynamic runs.

The second and third terms on the right side of equation (7) again represent
Coriolis force and radial tidal gravity, separately. Unlike for the gas,
vertical gravity for dust particles is included in the simulation given by
$-\Omega_{K}^{2}z_{i}{\hat{\boldsymbol{z}}}$, thus allowing particles to
settle towards the midplane. The last term represents the drag force on
individual particles, where gas velocities are interpolated to individual
particle positions. It is convenient to use a dimensionless stopping time
defined as $\tau_{s}\equiv\Omega_{K}t_{\rm stop}$, where particles with
$\tau_{s}\ll 1$ are strongly coupled to the gas while particles with
$\tau_{s}\gg 1$ are loosely coupled to the gas.

We use the second-order semi-implicit solver to integrate particles, with the
triangular shaped cloud (TSC) scheme for interpolation (for dust integrator)
and deposits (for dust backreaction), as described in Bai & Stone (2010b).

###  2.2 Simulation setup and parameters

The list of our simulation runs is shown in Table 1. Fiducially, we perform
non-ideal MHD simulations (labeled by “AD”) with particle feedback included,
and an additional MHD run with sufficiently low solid abundance so that the
dust feedback is negligible. We refer to the latter simulation as “no
feedback” (ADZ0), and rescale the particle density as appropriate when
comparing with other runs with dust feedback. We also perform pure
hydrodynamic simulations (labeled by “hyd”) with particle feedback (i.e. the
standard SI simulations) for reference. In addition, we explore different
levels of solid abundance, pressure gradient, turbulence levels and simulation
resolutions. The simulation setup and parameters are described as follows.

As a first study, our simulations consider a single particle species with
fixed stopping time $\tau_{s}=0.1$ for all our simulations. At outer disk
regions, gas drag is in the Epstein regime with $t_{\rm
stop}=\rho_{s}a/\rho_{g}c_{s}$ (Epstein, 1924), where $\rho_{s}$ and $a$ are
the solid density and physical size of dust particles. Fixing $\tau_{s}$ is
justified because gas density undergoes very little variation in unstratified
MRI simulation box with AD (Bai & Stone, 2011). In standard disk models,
$\tau_{s}=0.1$ corresponds to mm-cm size particles in the midplane regions of
the outer disk (e.g., Chiang & Youdin, 2010), which are also the most active
in driving the SI (e.g., Bai & Stone, 2010a).

We consider a range of particle abundances, measured by the height-integrated
dust-to-gas mass ratio $Z$, ranging from $Z=0.005$ to $Z=0.04$. We also
consider a run with $Z=10^{-5}$ where there is essentially no particle
backreaction. These runs are labeled by “$Z$” followed by a number indicating
the $Z$ values. Note that in interpreting $Z$, we assume the gas surface
density is given by $\Sigma_{g}=\sqrt{2\pi}\rho_{0}H_{g}$ as if the gas is
vertically stratified. The height-integrated dust mean surface density is thus
$\Sigma_{d}=Z\Sigma_{g}$. This particle abundance is then evenly divided to
individual particle “mass density” in the simulations.

For gas radial pressure gradient, we choose $\Pi=0.1$ by default, which is
generally applicable at outer disk ($\gtrsim$ 30 AU) 111For standard minimum-
mass solar nebular (MMSN, Weidenschilling, 1977) model with
$\Sigma(R)=1700(R/{\rm AU})^{-3/2}\ {\rm g}\ {\rm cm}^{-2}$ and
$T(R)=280(R/{\rm AU})^{-1/2}{\rm K}$, we have (e.g., Xu et al., 2017)
$\frac{\Delta v_{K}}{c_{s}}\approx 0.127\left(\frac{R}{30{\rm
AU}}\right)^{1/4}\ .$ (8) The observationally inferred disk surface
distribution in the main body of the disk is shallower ($\Sigma\sim R^{-1}$
instead of $R^{-3/2}$, e.g., Andrews et al. (2009, 2010), which lowers $\Delta
v_{K}/c_{s}$ to $\sim 0.1$ at 30 AU.. We further consider simulations with
different values of $\Pi=$ 0, 0.05, 0.1, 0.2 for comparison. These runs are
labeled by “$P$” followed by a number indicating the $\Pi$ values.

In our MRI simulations, we fix $Am=2$ for our simulations. This makes
turbulence slightly stronger than that with $Am=1$ so that we can better
resolve the dust particle layer, while still maintaining a realistic
turbulence level. We also fix the net vertical field with
$\beta_{0}=1.28\times 10^{4}$. This yields a most unstable MRI wavelength to
be about $\sim 0.12H$ (Bai & Stone (2011), see their Equation 22).

We choose simulation box size to be $H\times 2H\times H$ in $x$, $y$ and $z$
dimensions respectively, with a fiducial grid size of $256^{3}$. Thus, we can
well resolve the MRI with $\sim 30$ cells per most unstable wavelengths. Note
that the resolution in the $y-$ axis is half of that in other dimensions due
to the anisotropic nature of the MRI turbulence (e.g., Hawley et al., 2011).
Such resolution is also sufficient to resolve the SI given $\Pi=0.1$, with
$\sim 26$ cells in a characteristic length of $\Pi H$ for the SI. To our
knowledge, this resolution is also higher (if not much higher) than all
previous MRI simulations with particles (closest to ours is Johansen et al.
(2011), where their high-resolution run resolves the length of $\Pi H$ with
$\sim 20$ cells). We also consider a higher resolution run with a grid size of
$384^{3}$ for comparison. For reference, we also conduct two pure hydrodynamic
simulations, with smaller simulation box size of $0.6H\times 1.2H\times 0.6H$
and higher resolution with $480^{3}$ grid points. The main reason for using
even higher resolution is to properly resolve the thinner dust layer in this
case.

For MHD simulations, we first run simulations without particles for a time of
$t_{0}=120\Omega_{K}^{-1}$ with initially seeded small-amplitude random
velocity perturbations in order for the MRI to grow and fully saturate into
turbulence. Particles are then inserted into the simulation box with random
positions but with vertical distribution following a Gaussian profile
$\propto\exp(-z^{2}/2H_{\rm p0}^{2})$, with $H_{\rm p0}=0.02H_{g}$. The
particle velocity is initialized according to Nakagawa-Sekiya-Hayashi (NSH)
equilibrium (Nakagawa et al., 1986). For pure hydrodynamic simulations, we use
the same procedure except that particles are injected at the beginning. For
the remaining of this paper, we use the time $T$ after particle injection to
denote simulation time. In all simulations, there are on average one particle
per cell. As particles settle towards the midplane, allowing us to build up
sufficient particle statistics in the simulations to study the dynamics and
clumping in the particle layer.

Figure 2:  Time evolution of particle scale height (upper panel) and particle
maximum density (lower panel) for simulation runs with $\Pi=0.1$ and fiducial
resolution, namely, runs ADZ0, ADZ05, ADZ1, ADZ2, ADZ4, hydZ2, and hydZ4. “No
FB” stands for no feedback. The particle maximum density for run ADZ0 (shown
with black lines) is scaled to $Z=0.02$. The horizontal axis represents time
after inserting particles.

##  3 Overview of simulations

Before presenting more detailed analysis and discussions of our simulation
results, an overview of our simulations is shown in Figure 1, with a 3D
visualization of the dust spatial distribution in our fiducial non-ideal MHD
simulation (ADZ2) at the time of $T=120\Omega_{K}^{-1}$ after particle
injection. The figure highlights that dust settles to the midplane,
maintaining a layer with scale height $H_{p}$ due to a balance between
vertical gravity and turbulent diffusion. The dust layer is also radially
segregated, which is due to the formation of zonal flows, as will be discussed
later.

#### Dust Settling

The upper panel of Figure 2 shows the evolution of the normalized particle
scale height $H_{p}/H$ over time, for all the runs with pressure gradient of
$\Pi=0.1$, where $H_{p}$ is measured using the rms value of the vertical
coordinate for all the particles.

For the MRI turbulent run without particle feedback, the particle layer is
puffed up from the initial setup and saturates within the time of $\sim
30\Omega_{K}^{-1}$. The resulting $H_{p}\sim 0.08H$ is comparable with that
from the non-ideal MHD simulation with AD for $\tau_{s}=0.1$ in Zhu et al.
(2015) and Xu et al. (2017). With particle feedback included in the
simulations, the particle scale height in the steady state is lower than that
without feedback, ranging from $H_{p}\sim 0.03H$ (resolved by $\sim$16 cells)
to $H_{p}\sim 0.06H$. This reduction of dust layer thickness indicates that
particle feedback enhances dust settling in MRI turbulent disks.

For comparison, our pure hydrodynamic runs saturate at the time of $\sim
150\Omega_{K}^{-1}$, with $H_{p}\sim 0.03H$ and $H_{p}\sim 0.02H$ for hydZ2
and hydZ4, separately. The $H_{p}$ in the no-feedback (pure MRI) run is $\sim
3$ times higher than in the hydrodynamic (pure SI) runs, while with feedback
included, $H_{p}$ in AD case is only slightly/modestly higher than in the
hydrodynamic case. We will further discuss the vertical transport of dust
particles and the role of dust feedback in dust settling in §4.

#### Dust Concentration

In addition to the formation of dust filament seen in Figure 1, we further
show in the lower panel of Figure 2 the particle maximum density in the
simulations as a function of time.

With dust feedback included, particle maximum density is significantly
enhanced in our MHD simulations, especially when the solid abundance reaches
$Z\gtrsim 0.02$ (ADZ2 and ADZ4), with local particle density reaching up to
1000 times background gas density at the midplane, sufficient to trigger
planetesimal formation if self-gravity is included. For comparison, particle
clumping with similar particle maximum density is seen at $Z=0.04$ for the
pure hydrodynamic case, but not for the $Z=0.02$ case. For AD runs with lower
solid abundances (ADZ1 and ADZ05), particle concentration is much weaker and
the maximum particle local density is rarely exceed a hundred times of the
midplane gas density.

We will further analyze particle concentration in §5, and discuss the role of
pressure variations induced by zonal flows on particle concentration.

#### Level of turbulence

Both the dust layer thickness and the concentration of dust particles can be
largely affected by the turbulence level in the disk. The level of MRI
turbulence in our simulation can be characterized by $\alpha$ parameter for
disk angular momentum transport (Shakura & Sunyaev, 1973), which is obtained
by evaluating the time- and volume-averaged sum of the Maxwell stress and
Reynolds stress, normalized by thermal pressure:

$\alpha\equiv\frac{-B_{x}B_{y}+\rho v_{x}v_{y}}{\rho c_{s}^{2}}.$ (9)

Our non-ideal MHD simulation without feedback (ADZ0) gives $\alpha=8.3\times
10^{-4}$, which is consistent with previous studies (Zhu et al., 2015; Xu et
al., 2017). On the other hand, more direct characterization of level of
turbulence arises from simply computing the rms velocities for each velocity
components. We find from the ADZ0 run that $\langle
v_{x}^{2}/c_{s}^{2}\rangle=1.2\times 10^{-3}$, $\langle
v_{y}^{2}/c_{s}^{2}\rangle=1.3\times 10^{-3}$, and $\langle
v_{z}^{2}/c_{s}^{2}\rangle=6.2\times 10^{-4}$, comparable with the results
from Zhu et al. (2015). Simulations with dust feedback included generally give
such rms velocities that are close to that from run ADZ0.

##  4 Particle vertical transport

One important finding from our simulations is the reduction of particle scale
height when considering particle feedback. More quantitatively, here we list
the $H_{p}$ measured in our simulations, time-averaged from
$T=60\Omega_{K}^{-1}$ to the end of each run, being 0.082$H$, 0.055$H$,
0.051$H$, 0.037$H$ and 0.034$H$ for runs ADZ0 (no feedback), ADZ05, ADZ1, ADZ2
and ADZ4, respectively. In this section, we diagnose in more detail about the
cause of this reduction, from two separate perspectives, discussed in the
following two subsections.

As a preliminary, we expect particle scale height $H_{p}$ to be related to
both the turbulence properties and the particle stopping time $\tau_{s}$, and
can be written as (Youdin & Lithwick, 2007)

$H_{p}=\sqrt{\frac{\langle
v_{z}\rangle^{2}\tau_{e}}{\Omega^{2}\tau_{s}}}\cdot\sqrt{\frac{\tau_{s}+\tau_{e}}{\tau_{s}+\tau_{e}+\tau_{s}\tau_{e}^{2}}},$
(10)

where the dimensionless eddy time $\tau_{e}\equiv\Omega_{K}t_{\rm eddy}$ is
defined for the convenience of characterizing the turbulence eddy time $t_{\rm
eddy}$ (formally defined in §4.2). The turbulence properties control particle
scale height with both the gas vertical rms velocity $\langle v_{z}\rangle$
and the eddy time $t_{\rm eddy}$. The gas vertical diffusion coefficient can
be written as

$D_{g,z}=\langle v_{z}\rangle^{2}t_{\rm eddy}.$ (11)

For dust particles with Stokes number $St\equiv\tau_{s}/\tau_{e}\ll 1$ and
$\tau_{s}\tau_{e}\ll 1$, the second term of Equation (10) is approximately
unity, and Equation (10) is reduced to

$H_{p}\sim\sqrt{\frac{D_{g,z}}{\Omega\tau_{s}}}=\sqrt{\frac{\alpha_{z}}{\tau_{s}}}H,$
(12)

where $\alpha_{z}\equiv D_{g,z}/(c_{s}H)$ is a dimensionless parameter defined
for characterizing vertical turbulent diffusion.

###  4.1 Vertical profiles

If the dust vertical profile is Gaussian, it is uniquely characterized by
$H_{p}$. However, this is not necessarily the case. In Figure 3, we show the
vertical profiles of particle density (upper panel) and the gas vertical rms
velocity (lower panel), in order to further study particle vertical transport.

For MHD simulation without particle feedback (ADZ0), the particle density
profile generally agrees with the Gaussian profile, especially close to the
disk midplane. For run ADZ05, Where particle feedback is included but no
particle clumping is seen ($\max(\rho_{p})/\rho_{0}$ always lower than 100),
the dust density profile is overall Gaussian (except for more extended wings),
but becomes narrower than that for run ADZ0.

Despite this reduction in particle scale height, the lower panel of Figure 3
shows that the gas vertical rms velocity $\langle v^{2}_{z}\rangle$ in non-
ideal MHD simulations with and without dust feedback is highly consistent
throughout the disk, with $\langle v^{2}_{z}\rangle/c_{s}^{2}\sim 0.0006$.
This indicates that dust feedback is not controlling the particle layer
thickness in MRI turbulent disks by affecting the $\langle v^{2}_{z}\rangle$.
Recalling Equation (10), we thus need to further examine effect of the
turbulence eddy time (see next subsection).

For run ADZ2 (and similarly ADZ4 but not shown), there is further reduction in
$H_{p}$, but more interestingly, the vertical density profile of particles is
clearly non-Gaussian, as seen in Figure 2. A cusp is present at the midplane
in the density vertical profile, and besides the cusp, the density profile is
roughly exponential rather than Gaussian. Note that this density profile is
measured over time $T=100-200\Omega_{K}^{-1}$ after particle injection, where
strong dust clumping is present, and such deviations are likely the direct
consequence of dust clumping. However, dust clumping and additional settling
does not react back to the gas strong enough to affect the vertical rms
velocity of the gas, and we see that $\langle v^{2}_{z}\rangle$ remains
vertically uniform throughout the disk.

The results obtained above are in contrast with our pure hydrodynamic
simulations. In run hydZ2, the vertical dust diffusion is generated by the SI
turbulence itself. We see that the vertical rms velocity peaks at the
midplane, which is natural as most particles that drive the SI reside in the
midplane region, with $\langle v^{2}_{z}\rangle/c_{s}^{2}\sim 0.0003$,
slightly lower than that of the MHD runs. It then falls off away from the
midplane 222Away from the midplane, the increase of $\langle v^{2}_{z}\rangle$
is likely attributed to the numerical artifact from periodic vertical boundary
conditions, as discussed in Li et al. (2018), but as very few particles reside
at larger vertical heights, its impact on particle dynamics is minimum.. As a
result, the particle vertical profile becomes more extended than the Gaussian
profile at the midplane, but rapidly drops away from the midplane. This result
is similar to those from previous SI simulations (e.g. Bai & Stone, 2010a).
Overall, comparison of our hydrodynamic and MRI runs shows that the presence
of external turbulence has a strong effect on dust dynamics.

Figure 3:  Vertical profiles of particle density (upper panel) and gas
vertical rms velocity (lower panel) for runs hydZ2, ADZ0, ADZ05, and ADZ2.
Solid lines are time and horizontally averaged vertical profiles in the
saturate states of the runs. Dashed lines in the upper panel are the Gaussian
profiles assuming time-averaged particle scale heights from each run.

###  4.2 The effect of $t_{\rm eddy}$

The fact that dust settles more strongly when turning on feedback, while
$\langle v^{2}_{z}\rangle$ remains largely unchanged raises further question
to the origin of the reduced particle vertical diffusion. In addition to the
mean square gas vertical velocity $\langle v_{z}^{2}\rangle$, particle
diffusion in turbulent gas disk is also determined by the correlation time of
turbulent fluctuations, defined as

$t_{\rm eddy}\equiv\int_{0}^{\infty}R_{zz}(\tau)/R_{zz}(0)d\tau,$ (13)

where $R_{zz}(\tau)\equiv\langle v_{z}(\tau)v_{z}(0)\rangle$ is the auto
correlation function for $v_{z}$ at a timescale of $\tau$. Although
conventional ideal MHD simulations suggest $t_{\rm eddy}\sim\Omega_{K}^{-1}$
(Fromang & Papaloizou, 2006; Carballido et al., 2011), Zhu et al. (2015) found
that the correlation time in the vertical direction is increased in non-ideal
MHD simulations with AD.

In order to investigate the effect of particle feedback on the turbulent
correlation time, we calculate the auto-correlation functions for $v_{z}$ in
our simulations. We first output our simulation data with time intervals down
to 0.1$\Omega_{K}^{-1}$, and shift the data along the azimuthal direction with
a distance of 1.5$\Omega_{K}xt$ to correct for Keplerian shear. The auto-
correlation function $R_{zz}(\Delta t)$ is then derived by averaging
$v_{z}(t)v_{z}(t+\Delta t)$ over time and space. The results are shown in
Figure 4.

For pure hydrodynamic simulations, the integration of the auto correlation
function gives $t_{\rm eddy}\approx 1.5\Omega_{K}^{-1}$, which is consistent
with the conventional expectation of $t_{\rm eddy}\sim\Omega_{K}^{-1}$. For
non-ideal MHD simulation without particle feedback, the auto correlation
function extends to longer timescales and the correlation time is increased to
$t_{\rm eddy}\approx 5.0\Omega_{K}^{-1}$ due to the ambipolar diffusion,
overall consistent with Zhu et al. (2015). When particle is included in the
non-ideal MHD simulation, we find that the auto correlation function drops
rapidly at $\Delta t>2\Omega_{K}^{-1}$, and the resulting $t_{\rm eddy}\approx
1.5\Omega_{K}^{-1}$, close to the measured eddy time in the pure hydrodynamic
simulation. With the measured $t_{\rm eddy}$ and $\langle v_{z}\rangle^{2}$,
the vertical diffusion coefficients can be calculated with Equation (11). Our
calculation gives $D_{g,z}\sim 0.003$ for ADZ0, and $D_{g,z}\sim 9\times
10^{-4}$ for ADZ2. For pure hydrodynamic case, turbulent diffusion is
inhomogeneous and is strongest at midplane with $\langle
v^{2}_{z}\rangle/c_{s}^{2}\sim 3\times 10^{-4}$, corresponding to $D_{g,z}\sim
4.5\times 10^{-4}$. Following equation (12), our inferred $D_{g,z}$ would give
$H_{p}\approx 0.17H$ for ADZ0, and $H_{p}\approx 0.095H$ for ADZ2. Both of
these are higher than $H_{p}$ measured in the simulation by a factor of $\sim
2$, which may reflect the complex nature of realistic turbulence, but the
values are consistent with the fact that the particle layer thickness is
reduced by $\sim$ half in ADZ2 compared to ADZ0.

We conclude that in MRI turbulent disks under AD, particle feedback likely
enhances particle settling by reducing turbulence correlation time, instead of
controlling the gas rms velocity. We emphasize that our results are applicable
only in the MRI turbulence with AD, where $t_{\rm eddy}$ is increased compared
to ideal MHD case (Zhu et al., 2015). We have also tested (but not presented)
simulations with ideal MHD, and such reduction in eddy time is not observed
when dust feedback is included. This comparison further indicates that dust
feedback reduces the turbulence correlation time by canceling out the
enhancement of correlation time due to AD.

Figure 4:  Auto-correlation functions of gas vertical velocity fluctuations
for runs ADZ0, ADZ2, and hydZ2 as a function of time scale. Without dust
feedback, the eddy time for run ADZ0 is $t_{\rm eddy}\sim 5\Omega_{K}^{-1}$.
With dust feedback in ADZ2, the eddy time is reduced and is comparable to the
pure hydrodynamic case ($t_{\rm eddy}\sim 1.5\Omega_{K}^{-1}$).

##  5 Particle Concentration

###  5.1 Particle clumping

While our simulations do not include self-gravity, we expect that
gravitational collapse can be triggered to form planetesimals when dust
particles are sufficiently concentrated in high-density clumps. The standard
criterion for gravitational collapse is that self-gravity of the dust clump
must overcome the stellar tidal shear, i.e., dust density in the clump must
exceed the Roche density:

${\rho_{R}}=\frac{9M_{\star}}{4\pi a^{3}},$ (14)

where $M_{\star}$ denotes the mass of the central star, and $a$ is the
distance from the central star. Assuming a standard MMSN disk, the gas density
at the midplane is given by (Hayashi, 1981)

$\rho_{0}=1.4\times 10^{-9}(a/AU)^{-2.75}{\rm g~{}cm}^{-3}.$ (15)

Combining Equations (14) and (15), the corresponding local dust-to-gas ratio
for a solar-mass star is

$\frac{\rho_{R}}{\rho_{0}}=3.03\times 10^{2}(a/AU)^{-0.25}{\rm g~{}cm}^{-3},$
(16)

At 1 AU, ${\rho_{R}}/{\rho_{0}}$ is approximately 300, and
${\rho_{R}}/{\rho_{0}}\sim 130$ at outer disk around 30 AU.

In the lower panel of Figure 2, the evolution of particle maximum density
indicates strong clumping for solid abundances of $Z=0.02$ and $0.04$ in the
presence of MRI turbulence and dust feedback, as ${\rm max}(\rho_{p})$ reaches
$10^{3}$ larger than the gas density at midplane, well exceeding the estimated
Roche density. For lower solid abundances of Z=0.01 and Z=0.005, particle
concentration is much weaker and the maximum solid density is generally lower
than the Roche density. Also note that in run ADZ2, the particle maximum
density temporarily drops after $T\sim 200\Omega_{K}^{-1}$, indicating that
the clump is disrupted. This may be attributed to the fact that dust self-
gravity is not included in our simulations, and that $Z=0.02$ is very close to
the threshold solid abundance for dust clumping (see next subsection for more
clues). We will further discuss the survival of dust clumps and collapse
criterion in §7.3.

It is useful to compare the results of particle clumping in our MHD
simulations with those in pure hydrodynamic simulations. In our pure
hydrodynamic simulation, particle clumping is seen in for $Z=0.04$, but not
seen for $Z=0.02$. Surprisingly, this critical $Z$ for clumping is higher than
that in our MHD simulations. In general, dust abundance $Z$ above solar is
needed to trigger clumping from the SI (Johansen et al., 2009b; Bai & Stone,
2010a). Assuming a background gas pressure gradient of $\Pi=0.05$, the
clumping threshold for particles with $\tau_{s}\approx 0.1$ is about $Z\gtrsim
0.02$ (Carrera et al., 2015, Yang et al., 2017, but see Li & Youdin, 2021 who
found a lower $Z$ threshold with improved numerics). However, our simulation
adopts a higher pressure gradient of $\Pi=0.1$ that is likely more applicable
to the outer disk, and our result of dust clumping threshold at $Z=0.04$ is
generally consistent with the study of pressure gradient versus critical solid
abundance (Bai & Stone, 2010c; Sekiya & Onishi, 2018).

Overall, our simulations show that under more realistic MRI turbulence at the
outer disk, particle clumping generally requires lower solid abundance (i.e.,
easier!) compared to conventional SI studies in pure hydrodynamic case, and
the trend that clumping favors higher Z still holds in our MHD simulations.

Figure 5:  Time evolution of the azimuthally and vertically averaged radial
profiles of gas density (left), Keplerian-subtracted gas azimuthal velocity
(middle), and dust density (right) at the midplane ($\pm 1$ cell from $z=0$),
for the simulation run ADZ0 (upper panels) and ADZ2 (lower panels). Dust
densities for ADZ0 run is scaled to Z=0.02 for easy comparison. Note that
$T=0$ stands for the time of particle injection, and our first particle output
is at $T=10\Omega_{K}^{-1}$.

###  5.2 Clumping vs. pressure variations

Apparently, our finding on particle clumping challenges results from recent
studies (e.g. Umurhan et al., 2020; Chen & Lin, 2020; Gole et al., 2020) that
turbulence damp the SI and thus suppress dust clumping. To further diagnose
the reason why MRI turbulence promotes clumping, and the clumping process in
general, we present in Figure 5 the time evolution of azimuthally and
vertically averaged radial profiles of gas density, gas azimuthal velocity,
and dust density for run ADZ0 and ADZ2. The formation and evolution of the
dust clumps in ADZ2 are seen in the lower-right panel at $T=80\sim
180\Omega_{K}^{-1}$ and $T=330\sim 380\Omega_{K}^{-1}$, corresponding to the
periods of dust maximum density reaching $\sim 1000$ in Figure 2.

One important feature in the MRI turbulence is the spontaneous formation of
long-lived (tens of orbits) zonal flows, corresponding to large-scale, quasi-
axisymmetric density (and hence pressure) variations (Fromang & Nelson, 2005;
Johansen et al., 2006, 2009a). These large-scale gas pressure variations can
be seen in the left panels of Figure 5. The pressure variations change the
pressure gradient profile, alter particle radial drift and potentially trap
particles at the local pressure maxima. Note that because we implement
background pressure gradient by a radial force, the actual local pressure
maximum, if present, will appear slightly to the negative-$x$ (left) side of
the density maximum.

A more straightforward way to examine the correspondence of particle clumping
and pressure variations is by analyzing the gas azimuthal velocity. A positive
pressure gradient leads to super-Keplerian rotation, corresponding to
$v_{y}>0$ in our simulation with Keplerian shear subtracted. A negative
pressure gradient corresponds to $v_{y}<0$. Hence, the local pressure maximum
is marked by $v_{y}=0$ regions, with $v_{y}>0$ on the left side and $v_{y}<0$
on the right. The gas azimuthal velocity $v_{y}$ versus time is shown in the
middle panels of Figure 5 333Our formulation mimics the outward pressure
gradient in the gas by an inward radial force on the particles (see §2.1),
making gas azimuthal velocity $\eta v_{K}=\Pi c_{s}$ larger than the actual
case. In Figure 5, we have subtracted this component to show $v_{y}$ in the
actual situation..

In ADZ0 where no clumping is present, we see that pressure variation from the
zonal flow is not strong enough to overcome the background pressure gradient,
and gas azimuthal velocity $v_{y}$ is always negative (indicating sub-
Keplerian rotation) throughout the simulation box. With dust feedback included
in ADZ2, on the other hand, pressure variation is enhanced. It successfully
overcomes the background gradient to generate a pressure maximum ($v_{y}=0$
region) to the left of the density peak, and particle clumping is seen in
regions with $v_{y}=0$. Although it is less clear why this is the case, our
results indicate that dust feedback can promote clumping by enhancing the gas
pressure bump in the zonal flow. Moreover, we can infer that the disruption of
the clump around $T=200\Omega_{K}^{-1}$ is associated with the weakening of
the zonal flow, where a pressure maxima is no longer sustained.

Overall, we have identified that the formation of pressure maxima and dust
feedback are the two essential elements for dust clumping, and that dust
feedback appears to enhance zonal flows and promote the formation of pressure
maxima. Dust clumping and planetesimal formation thus follow a two-stage
process with particle concentration in pressure maxima followed by clumping
due dust feedback.

What is yet to be clarified is the interpretation on the role of dust
feedback. It is commonly taken into granted that incorporation of dust
feedback leads to clumping by the SI. However, dust feedback does not
necessarily lead to clumping (solid abundance must exceed some threshold), and
on the other hand, the SI is not supposed to operate at pressure maxima (its
free energy comes from the pressure gradient). While small pressure gradient
promotes dust clumping for pure SI (Bai & Stone, 2010c), we do not observe
dust clumping when pressure variation from the zonal flow is close but not yet
to overcome background pressure gradient (e.g., around time $T\sim
200\Omega_{K}^{-1}$ in run ADZ2 where $\rho_{\rm max}$ plunges, also see §7.1
and Figure 9). Alternatively, one may also argue that our results are in line
with the notion that the SI is suppressed by modest-to-strong level of
turbulence (Gole et al., 2020). Yet, the requirement of dust feedback for
clumping remains to be addressed.

Similar results were obtained in ideal MHD simulations of Johansen et al.
(2011), where they attribute the clumping to dust trapping in pressure “bumps”
generated by the zonal flow with the SI being a secondary effect.
Nevertheless, zonal flows in their ideal MHD simulations appear more
intermittent with lower amplitudes than ours, and they did not address whether
pressure variations in their simulations actually overcome background pressure
gradient to become bumps. Our analysis from more realistic setting with AD-
dominated MRI turbulence in the outer PPDs allows us to clearly identify the
occurrence of clumping with pressure maxima, and we will present further study
in §7.1 which will strengthen our conclusions.

###  5.3 Properties of particle clumps

In the lower panel of Figure 2, we identify clumping events by the maximum
local dust density. In order to further analyze the properties in the particle
clumps, we show in Figure 6 the distribution function of dust density. The
local dust density $\rho_{p}$ is uniformly divided into 1000 bins in
logarithmic space between $0.01\rho_{0}$ and $1000\rho_{0}$, and the
distribution function is measured by counting the number of dust particles
located in a region with local $\rho_{p}$ in each bin.

We take the first clumping period in ADZ2 ($T=0-210\Omega_{K}^{-1}$) for
analysis, and measure the distribution function every $40\Omega_{K}^{-1}$
(except that the first measurement is taken at $T=10\Omega_{K}^{-1}$ to avoid
numerical artifacts at the beginning of inserting particles). Each
distribution function is time-averaged within 10$\Omega_{K}^{-1}$, in order to
reasonably reduce the noise level. For comparison, we also present a
distribution function of run ADZ0 measured at $T=11-20\Omega_{K}^{-1}$, same
time period as the first function measured in ADZ2. Due to the absence of dust
clumping, the distribution function for ADZ0 is generally consistent
throughout simulation time.

The whole process of clump formation and evolution in ADZ2 is well illustrated
by the $\rho_{p}$ distribution function. The run starts with a smooth Poisson-
like distribution prior to clumping ($T=11-20\Omega_{K}^{-1}$), with width
reflecting the density fluctuations in turbulence. The distribution function
peaks at $\rho_{p}\sim\rho_{0}$, consistent with the expected dust density at
midplane that we initially set up in the simulation, $\rho_{p0}\sim
Z\rho_{0}H/H_{p0}\sim\rho_{0}$, where $H_{p0}=0.02$ is the initial particle
scale height set in the simulation. The distribution maximizes at
$\rho_{p}\sim 10\rho_{0}$, consistent with the $max(\rho_{p})/\rho_{0}$ value
shown in Figure 2.

The profile then flattens and shifts towards higher density as the clumping
starts ($T=81-90\Omega_{K}^{-1}$). During the time of steady, maximum clumping
(e.g. $T=121-130\Omega_{K}^{-1}$), the density distribution peaks at
$\rho_{p}\sim 100\rho_{0}$, and maximizes at $\rho_{p}\sim 1000\rho_{0}$,
again consistent with Figure 2. After the period of clumping, the particle
density distribution shifts back after $T=161-170\Omega_{K}^{-1}$, with the
distribution profile similar to that before clumping. This profile before and
after clumping is also consistent with the density distribution profile in
ADZ0, indicating that dust feedback itself has little effect on particle
density distribution properties, and the profile/shape of distribution
function is mainly controlled by the presence of strong clumping events.

Overall, during the clumping event, the bulk of the dust density distribution
shifts towards higher values, deviates from Poisson profile, and peaks at
close to the high-density end of the distribution. This further indicates that
in addition to forming the densest clumps of $\rho_{p}\sim 1000\rho_{0}$, a
large fraction of particles also reside in relatively dense clumping regions
($\sim 100\rho_{0}$).

Figure 6:  Distribution function of dust density in runs ADZ0 (no feedback)
and ADZ2 ($Z=0.02$). For run ADZ2, the distribution function is measured at
different simulation times $T$ in the unit of $\Omega_{K}^{-1}$, and marked in
different colors (see legend). Each distribution function is time-averaged
within $10\Omega_{K}^{-1}$. For run ADZ0, the distribution function is
measured and time-averaged between $T=11-20\Omega_{K}^{-1}$.

##  6 Radial transport of dust particles

Figure 7:  Left/middle: probability distribution for particle radial
displacements in run ADZ0 (no feedback)/run ADZ2 (Z=0.02). The probability
distributions are measured at different simulation times $T$ in the unit of
$\Omega_{K}^{-1}$, and marked in different colors (see legend). Upper right:
time evolution of the variance of the particle radial displacement. Lower
right: time evolution of the mean particle radial displacement. The black
dashed lines in the right panels show the linear fit to the results from run
ADZ0.

In this section, we analyze the radial transport of dust particles, including
the role of dust clumping on dust transport. In doing so, we track the radial
trajectory $x_{i}(t)$ of a sample of particles in the simulation. From the
time of particles injection $t_{0}$, we trace the radial displacement of
particles $x_{i}(t_{0}+\Delta t)-x_{i}(t_{0})$, and reset particle positions
when they cross radial boundaries in order to maintain continuous particle
trajectories.

The left and middle panels of Figure 7 show the distribution of particle
radial displacement over time of $T=180\Omega_{K}^{-1}$ after injecting
particles for ADZ0 and ADZ2 runs respectively, corresponding to the first
round of particle clumping in run ADZ2 (see the lower panel of Figure 2).
Without dust feedback, the bulk of particles drift radially inward, and
diffuse to exhibit broader radial width over time. In the following
subsections, we analyze particle radial drift and diffusion separately.

###  6.1 Radial drift

In laminar disks, for a sub-Keplerian speed $\Delta v_{K}$, dust particles
suffer from radial drift. For particles with dimensionless stopping time of
$\tau_{s}$, and ignore particle backreaction, the radial and azimuthal drift
velocity relative to the local Keplerian velocity is given by (Nakagawa et
al., 1986)

$v_{r}=-2\frac{\tau_{s}}{\tau_{s}^{2}+1}\Delta v_{K}\ ,\quad
v_{\phi}=-\frac{1}{\tau_{s}^{2}+1}\Delta v_{K}\ .$ (17)

For our choice of $\Delta v_{K}=0.1c_{s}$, dust particles are expected to
drift inward with $v_{r}=-0.02c_{s}$.

We measure the mean particle radial displacement $\langle\Delta x\rangle$ as a
function of time for ADZ0 and ADZ2 runs, shown in the lower right panel of
Figure 7. For ADZ0, the particle displacement increases almost linearly with
time, and a linear regression fit estimates the averaged radial drift velocity
of $v_{r}=-0.017c_{s}$, very close to the theoretical expectation of
$-0.02c_{s}$. We also see that the drift speed also slightly deviate from the
linear fit, likely owing to weak zonal flows that slightly alter the radial
pressure profile.

For run ADZ2, when dust feedback is considered, the $\langle\Delta x\rangle-T$
curve does not closely follow the theoretical expectation anymore. Around
$T\sim 50-150\Omega_{K}^{-1}$, corresponding to the time period with strongest
clumping and particle trapping in zonal flows in the simulation, the radial
drift is strongly suppressed. Before and after this period, the slope of the
$\langle\Delta x\rangle-T$ curve remains similar to that of ADZ0.

This reduction in dust radial drift in run ADZ2 is also seen in the middle
panel of Figure 7. Comparing to the smooth shift of particle distributions in
run ADZ0, the particle radial distribution in run ADZ2 shows substantial
overlap around $T=60-150\Omega_{K}^{-1}$, indicating the bulk of the particles
are trapped and the radial drift is halted.

###  6.2 Radial diffusion

Under the assumption of random walk, the dust radial diffusion coefficient
$D_{x}$ can be measured by

$D_{x}=\frac{1}{2}\frac{d\sigma_{x}^{2}}{dt},$ (18)

where $\sigma_{x}$ is the rms of particle radial displacement distribution. We
measure the $\sigma_{x}^{2}$ as a function of time for runs ADZ0 and ADZ2,
presented in the upper right panel of Figure 7.

For ADZ0, the $\sigma_{x}^{2}-T$ curve closely follows the linear relation,
and by fitting the slope, we obtain $D_{x}\approx 1.2\times 10^{-3}$. Assuming
the gas diffusion coefficient on the radial direction $D_{g,x}$ is largely the
same (given $\tau_{s}\ll 1$), the value is lower than $D_{g,z}=0.003$ that we
measured for ADZ0. This is overall consistent with the result of Zhu et al.
(2015), where $D_{g,x}$ is generally lower (by a factor of $\sim 2$) than
$D_{g,z}$ in the MRI turbulence with AD. But we also note with the reduction
of turbulent eddy time in run AD2 (§4.2), the value of $D_{g,z}$ and $D_{g,x}$
becomes comparable.

When dust feedback is included, with the presence of dust clumping, particles
no longer follow normal diffusion. The $\sigma_{x}^{2}-T$ curve thus depends
on the starting time of measuring the displacement, and is no longer well
characterized by a diffusion coefficient. In run ADZ2, the $\sigma_{x}^{2}-T$
curve is largely consistent with that of ADZ0 before $T\sim
50\Omega_{K}^{-1}$, prior to strong particle clumping. During the clumping
period ($T=50-150\Omega_{K}^{-1}$), the $\sigma_{x}^{2}-T$ curve is flatter,
mainly due to the bulk of particles trapped in the clump and pressure maxima
in the zonal flow, similar to our finding in particle radial drift in §6.1. At
the end of clumping, $\sigma_{x}^{2}$ is strongly enhanced. As zonal flow
weakens, many particles leak from the left side of the pressure bump, re-enter
the box from the right side of the bump, and immediately experience a stronger
pressure gradient to gain a faster radial drift velocity than the background,
which enlarges $\sigma_{x}^{2}$. This can also be seen in the middle panel of
Figure 7, where particles with largest radial displacement in run ADZ02 start
to catch up with those in run ADZ0 at the final time of
$T=180\Omega_{K}^{-1}$. This enhancement of particle drift also contributes to
the dust radial drift velocity that is slightly faster than the ADZ0 case
before and after dust clumping period of $T\sim 50-150\Omega_{K}^{-1}$ (see
lower right panel of Figure 7).

Overall, our analysis of dust radial transport indicates that pressure maxima
from zonal flows, even being temporal, can efficiently trap particles and
overcome turbulent diffusion.

##  7 Discussion

###  7.1 The effect of background pressure gradient

Figure 8:  Time evolution of particle maximum density since particle
injection, similar to the lower panel of Figure 2, but for additional
simulation runs. Upper panel: simulations for $Z=0.02$, compared with the
higher resolution run (ADZ2H). Middle panel: simulations for $Z=0.01$,
compared with higher resolution run (ADZ1H). Lower panel: simulations with
different background pressure gradients. “No FB” stands for no feedback.

In §5, we show that the particle clumping mechanism in our simulation is
different from the conventional SI studies. Particle clumping under the
conventional SI scenario is expected to be triggered in the presence of a weak
background pressure gradient, while dust clumping under MRI turbulence occurs
in the presence of a local pressure maximum, where the pressure variation
induced by zonal flow overcomes the background pressure gradient. In order to
further examine this clumping mechanism, we perform runs with different
background pressure gradients for comparison.

The lower panel of Figure 8 shows the time evolution of the local maximum
particle density in simulations with pressure gradients of $\Pi=0$, $0.05$,
$0.1$, and $0.2$, separately. No particle clumping is seen under $\Pi=0.2$,
the strongest background pressure gradient that we test, and clumping is
stronger towards lower pressure gradient. Especially, the strongest particle
clumping is seen in our simulation in the absence of a background pressure
gradient. We note that the overall trend is similar to that with the SI: the
critical $Z$ required for clumping increases monotonically with the background
pressure gradient $\Pi$ (Bai & Stone, 2010c). Nevertheless, the SI does not
directly operate in pressure bumps as its free energy arises from background
pressure gradient. Thus, particle clumping seen in our $\Pi=0$ run further
highlights the different clumping mechanism in our case under MRI turbulence.

In §5.2, we illustrated from Figure 5 that the clumping occurs at the local
gas pressure maxima, or regions where gas rotation is Keplerian (i.e.
$v_{y}=0$). Similar to Figure 5, the upper and middle panels of Figure 9
further show the relation between dust clumping and the local pressure maxima,
for our runs with different background pressure gradients.

For run ADZ2P0 without background pressure gradient, a pressure variation is
formed at $T=0-30\Omega_{K}^{-1}$ at $x\sim-0.25H$ with an amplitude of $\sim
0.03\rho_{0}$, strong enough to generate a pressure maximum (see the
corresponding $v_{y}=0$ region in the middle panel). Although this pressure
variation evolves and disappears at later time, it triggers particle clumping,
and the clump stays relatively steadily even after the pressure maxima
disappeared at later time. Similarly, a second weak particle concentration is
seen at $x\sim 0.1H-0.2H$, with corresponding gas pressure maximum and
Keplerian rotation ($v_{y}=0$) region seen in the simulation box. A third
$v_{y}=0$ region is seen at $x\sim 0.4H$. However, this region does not
correspond to local pressure maximum, but instead a pressure minimum, as the
gas rotates at a sub-Keplerian velocity at left side and super-Keplerian at
the right. Thus, this $v_{y}=0$ region does not lead to particle clumping.

In our simulation run with higher pressure gradient $\Pi=0.2$, a pressure bump
that is largely similar to that in the ADZ2P0 run is formed, also at
$T=0-30\Omega_{K}^{-1}$, $x\sim-0.25H$, with the same amplitude of $\sim
0.03\rho_{0}$. However, no robust particle clumping is seen in the right
panel, simply because the pressure bump does not overcome the stronger
background pressure gradient, and $v_{y}$ is overall sub-Keplerian throughout
the simulation. The only region with $v_{y}$ relatively close to Keplerian
velocity at later stage of the simulation. It leads to a weak particle
concentration, as is seen in a small particle density enhancement shown in the
right middle panel of Figure 9.

This exploration of pressure gradients again illustrates that the key reason
for dust clumping is the presence of the pressure maxima. The change of
background pressure gradient affects particle clumping by altering the overall
radial profile of gas pressure and hence the presence of the pressure maxima.

Figure 9:  Similar to Figure 5, but for simulation runs of ADZ2P0 (top),
ADZ2P02 (middle), and ADZ2H (bottom).

###  7.2 The effect of resolution

Our fiducial resolution of 256/H fairly well resolves the particle layer
($H_{p}\sim 0.03H$ for ADZ2) with $\sim$ 16 cells, as well as the width of the
particle clumping region ($\lesssim 0.1H$, see e.g. Figure 5) with $\sim 25$
cells. In our simulations with higher resolution of 384/H, dust clumping
properties can be slightly different. The upper panel of Figure 8 shows that
for Z=0.02, higher resolution may lead to weaker clumping, where particles
reach a maximum density similar to the fiducial resolution, but the clumping
lasts for a shorter time. For Z=0.01 (middle panel), particle clumping is not
seen in the simulation with resolution of 384/H, whereas in fiducial
resolution clumping occurs briefly.

The effect of resolution on particle clumping is further shown in the lower
panel of Figure 9, the space-time plot for run ADZ2H. Similar to previously
demonstrated (see 5.2 and 7.1), dust clumping is coincident both in time and
space with the presence of local pressure maxima ($v_{y}=0$ regions). Compared
to run ADZ2, zonal flow in the high resolution ADZ2H run is overall weaker,
and the resulting amplitude of the pressure variation is on average lower. As
a result, the different outcomes in dust clumping we observe in runs with
different resolutions likely primarily reflect the non-convergence of zonal
flow properties with resolution. We note that in reality, zonal flows
properties are tied to disk evolution at global scale, and are not necessarily
well characterized in local simulations (see further discussion in §7.4).
Therefore, we do not pursue further discussion about this non-convergence, but
mainly highlight the connection between dust clumping and presence of a
pressure maxima.

###  7.3 Further collapsing criterion and the effect of turbulence diffusion

Without external turbulent diffusion, SI simulations with self-gravity show
that clumping is followed by gravitational collapse into planetesimals (e.g.,
Simon et al., 2016; Abod et al., 2019; Li et al., 2019). Our simulations lack
self-gravity. However, with external turbulence, the gravitational collapse of
dust clumps may also be suppressed by turbulent diffusion.

Gerbig et al. (2020) investigated the collapsing criterion considering the
competition between self-gravity and turbulent diffusion by both KHI and SI,
in addition to the conventional criterion of self-gravity overcoming the
stellar tidal shear (i.e., the Roche criterion). In order to overcome the
tidal shear, the dust clump radius $r_{\rm clump}$ must not exceed the Hill
length scale, given by their equation (13)

$r_{\rm clump}<\frac{2\pi}{9}\frac{G\Sigma_{d,c}}{\Omega_{K}^{2}}\equiv L_{\rm
Hill}\ ,$ (19)

where $\Sigma_{d,c}$ is the surface density of the dust clump.

In order to overcome the turbulent diffusion, the clump radius $r_{\rm clump}$
must be large enough so that the timescale for particles to be diffused over
$r_{\rm clump}$ is longer than the gravitational collapsing timescale.
Assuming turbulent diffusion coefficient is characterized by $\alpha$, their
equation (9) gives

$r_{\rm clump}>\frac{\alpha}{\tau_{s}}\frac{c_{s}^{2}}{2\pi
G\Sigma_{d,c}}\equiv L_{\rm diff}\ .$ (20)

Combining equations (19) and (20) leads to

$\frac{3}{2}\sqrt{\frac{\alpha}{\tau_{s}}}<\frac{\pi
G\Sigma_{d,c}}{c_{s}\Omega_{K}}\ .$ (21)

By adopting Toomre $Q$ parameter (Toomre, 1964)

$Q\equiv\frac{c_{s}\Omega_{K}}{\pi G\Sigma_{g}}\ ,$ (22)

and the local metallicity enhancement

$\epsilon\equiv\frac{\Sigma_{d,c}}{\Sigma_{d}}=\frac{\Sigma_{d,c}}{Z\Sigma_{g}},$
(23)

their equation (70) gives a collapsing criterion of

$Q<Q_{crit}=\frac{2\epsilon
Z}{3}\sqrt{\frac{\tau_{s}}{\alpha}}=\frac{2\epsilon Z}{3}\frac{H}{H_{p}}\ .$
(24)

In our simulations, the maximum particle density in runs with particle
clumping (e.g. ADZ2) is on the order of $10^{3}$, around 100 times higher than
the runs without clumping (e.g. ADZ05, see the lower panel of Figure 2).
Adopting this enhancement in particle maximum density of $\epsilon=100$,
together with $Z=0.02$, and $H/H_{p}\sim 20$, analytical prediction by Gerbig
et al. (2020) gives a critical collapsing value of $Q_{crit}\sim 27$,
comparable to the Toomre $Q\sim 30$ in MMSN (Weidenschilling 1977, Hayashi
1981).444The surface density slope in MMSN is steep, with $\Sigma\sim
R^{-3/2}$, reaching $\Sigma\approx 10$g cm-2 at 30AU. More realistic disk
surface density profiles are shallower, e.g., Andrews et al., 2009, leading to
higher $\Sigma$ at 30 AU and hence smaller $Q$ values. Considering the fact
that the diffusion $\alpha$ in equation (24) is expected to be measured within
the dust clump, which could largely deviate from the global $\alpha$ estimated
from $H_{p}/H$, our estimation here gives an lower limit of $Q_{crit}$.

Under SI-induced turbulence, Klahr & Schreiber (2020, 2021) tested the
collapsing criterion using simulation with box size comparable to the size of
a dust clump, measuring particle diffusivity in the clumping region. They
suggested that the level of turbulence within the SI-clumps sets the size of
planetesimals. However, as clumping occurs in pressure bumps in our
simulations without much free energy to drive the SI, source of turbulence in
our (unresolved) clumps is likely dominated by MRI. Though we are not in a
position to infer the exact turbulence level at this moment, our study may
open up a new avenue for understanding the size distribution of planetesimals.

###  7.4 Towards planetesimal formation in realistic disks (and rings)

In §1.2, we raised the conundrum where external turbulence is disfavored by
conventional understanding of planetesimal formation by the SI, yet
observations indicate modest level of turbulence by modeling the width of ALMA
rings. Our simulations represent the first step towards resolving this
conundrum: we have found strong evidence that particle clumping occurs in
pressure maxima formed by the MRI zonal flows. To some extent, our results are
consistent with the suppression of SI by external turbulence, as we do not
find particle clumping occuring in regions with finite pressure gradient. It
remains to identify the precise mechanism for particle clumping in pressure
maxima, where the SI is not expected to operate, but at least we observe that
dust feedback does play a role in enhancing zonal flows and promote clumping
(see also analysis by Auffinger & Laibe, 2018).

Our results highlight the importance of investigating planetesimal formation
in turbulent pressure bumps. The recent work of Carrera et al. (2021), who
conducted hydrodynamic simulations of the SI, is towards this direction.
However, most of their imposed pressure variations were insufficient to
overcome background pressure gradient, and planetesimal formation in their
simulations mostly occurs by the SI in regions with low pressure gradient,
instead of at the point of pressure maximum (if present). This is no longer
the case in our simulations with external MRI turbulence, where clumping
occurs just at the pressure maxima. We also note that in the presence of a
real pressure bump but without external turbulence, while dust can drive some
turbulence and form planetesimals by the SI as they drift towards the pressure
maximum, as also found in Carrera et al. (2021), such processes are transient
and remaining dust particles would eventually settle and concentrate
indefinitely towards the midplane of the bump center as free energy from
radial drift is exhausted. Such a state is inconsistent with the finite width
of ALMA rings. Therefore, we anticipate that in reality, planetesimal
formation likely occurs in turbulent pressure bumps. While we do not impose
pressure bumps on our own, the MRI zonal flows provides a natural mechanism to
create the pressure variations, potentially overcoming the background pressure
gradient to form pressure bumps, serving as an initial step towards more
realistic understandings.

One major caveat in our simulations is the use of local shearing-box framework
without vertical stratification, making the zonal flow properties subject to
large uncertainties. Recent global simulations of the MRI with AD confirmed
the ubiquity of zonal flows in more realistic disk setting (Cui & Bai, 2021).
The zonal flows are of modest amplitudes (but higher than ours) with pressure
bumps separated by a few disk scale heights, likely compatible with some ALMA
ring systems that show finer substructures (e.g. Jennings et al., 2021).
Therefore, zonal flows from in our simulations are likely not far from more
realistic ones, while future investigations of dust dynamics in more realistic
pressure bumps are needed to gain further understandings.

##  8 Conclusions and Perspectives

In this paper, we conduct high-resolution local shearing-box simulations of
hybrid particle-gas MHD to study dust dynamics in weakly turbulent PPDs. We
target disk outer regions with turbulence driven by the MRI, incorporating AD
as the dominant non-ideal MHD effect, which provides a relatively low but
realistic turbulent environment. Dust backreaction is included, which is
believed to be the key ingredient for dust concentration and planetesimal
formation in conventional SI theory without external turbulence. One of the
main motivations is to resolve the recently raised conundrum that while the SI
is widely believed to be the most promising mechanism for planetesimal
formation, it can be damped or suppressed by external turbulence, as recent
theoretical and numerical studies have shown (Umurhan et al., 2020; Chen &
Lin, 2020; Gole et al., 2020). Such turbulence is likely present in PPDs from
both observational (Dullemond et al., 2018; Rosotti et al., 2020) and
theoretical (e.g., Simon et al., 2013a; Cui & Bai, 2021) perspectives, which
fundamentally challenges our understandings of planetesimal formation.

This conundrum is resolved based on our simulation results, which are
summarized below:

  * 1.

Dust feedback promotes dust settling towards the midplane in MRI turbulent
disks with AD. This enhancement of settling is due to dust feedback reducing
the turbulence correlation time.

  * 2.

Dust clumping under AD-dominated MRI turbulence is generally easier compared
to the conventional SI case under pure hydrodynamics. Dust clumping favors
higher solid abundance $Z$ in the disk, similar to the case of pure SI, but
the threshold $Z$ is lower in AD-dominated MRI turbulence.

  * 3.

The necessary conditions for dust clumping under MRI turbulence include dust
feedback and the presence of gas pressure maxima.

  * 4.

In the AD-dominated MRI turbulence, pressure maxima is achieved by the MRI
zonal flow that is strong enough to overcome background pressure gradient.
Dust feedback enhances the amplitudes of the zonal flow and promotes clumping.

Our results suggest that dust clumping and planetesimal formation is a two-
stage process, where dust concentration by pressure maxima followed by
clumping due to dust feedback. Nevertheless, how dust feedback leads to
clumping is debatable. While it is commonly assumed to be due to the SI, the
SI is not expected to operate at the pressure maxima, nor under modest-to-
strong level of turbulence. The exact role of dust feedback, and potentially
the SI, on particle clumping in the MRI turbulence, remains open for future
investigations.

As a first study towards dust dynamics under more realistic disk environments,
our work is subject to several limitations. First, the properties of the zonal
flows are not necessarily well captured in local unstratified shearing-box
simulations (e.g., Bai & Stone, 2014). Fortunately, recent full 3D global
simulations confirm the existence of such zonal flows that do form pressure
bumps (Cui & Bai, 2021), with bump separation of a few $H$. While
incorporating dust to such global simulations with similar resolution can be
extremely computationally costly, this can be mitigated by conducting local
simulations and drive a pressure bump by forcing, to be shown in our immediate
follow up paper. Second, our simulations have neglected self-gravity and hence
cannot directly follow the collapse of dust clumps towards planetesimals.
Introducing self-gravity is another natural next step, which would break the
scale-free nature of our simulations, and we defer for future studies.
Finally, we consider single particle species with fixed stopping time in this
work. Further explorations should include a dust size distribution, where dust
of different sizes likely act and feedback differently in turbulent MRI disks
with/without pressure bumps. Including a dust size distribution will also
allow us to compare and contrast with results of multi-species streaming
instability without external turbulence (e.g., Krapp et al., 2019; Schaffer et
al., 2021; Zhu & Yang, 2021).

We thank Chao-Chin Yang, Min-Kai Lin, and Guillaume Laibe for constructive
conversations and discussions, and Greg Herczeg for overseeing the completion
of this work. We thank the referee, Daniel Carrera, for a detailed and
constructive report and helpful discussions. This project is supported by the
National Key R&D Program of China (No.2019YFA0405100). ZX acknowledges the
support of NSFC project 11773002. Numerical simulations are conducted on
TianHe-1 (A) at National Supercomputer Center in Tianjin, China, and on the
Orion cluster at Department of Astronomy, Tsinghua University.

## References

  * Abod et al. (2019) Abod, C. P., Simon, J. B., Li, R., et al. 2019, ApJ, 883, 192
  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3
  * Andrews et al. (2018) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018, ApJ, 865, 157
  * Andrews et al. (2009) Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502
  * Andrews et al. (2010) —. 2010, ApJ, 723, 1241
  * Armitage (2015) Armitage, P. J. 2015, arXiv e-prints, arXiv:1509.06382
  * Auffinger & Laibe (2018) Auffinger, J., & Laibe, G. 2018, MNRAS, 473, 796
  * Bai (2011a) Bai, X.-N. 2011a, ApJ, 739, 50
  * Bai (2011b) —. 2011b, ApJ, 739, 51
  * Bai (2015) —. 2015, ApJ, 798, 84
  * Bai & Stone (2010a) Bai, X.-N., & Stone, J. M. 2010a, ApJ, 722, 1437
  * Bai & Stone (2010b) —. 2010b, ApJS, 190, 297
  * Bai & Stone (2010c) —. 2010c, ApJ, 722, L220
  * Bai & Stone (2011) —. 2011, ApJ, 736, 144
  * Bai & Stone (2014) —. 2014, ApJ, 796, 31
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Balsara et al. (2009) Balsara, D. S., Tilley, D. A., Rettig, T., & Brittain, S. D. 2009, MNRAS, 397, 24
  * Barranco (2009) Barranco, J. A. 2009, ApJ, 691, 907
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  * Blum & Wurm (2008) Blum, J., & Wurm, G. 2008, ARA&A, 46, 21
  * Carballido et al. (2011) Carballido, A., Bai, X.-N., & Cuzzi, J. N. 2011, MNRAS, 415, 93
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Carrera et al. (2021) Carrera, D., Simon, J. B., Li, R., Kretke, K. A., & Klahr, H. 2021, AJ, 161, 96
  * Chen & Lin (2020) Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132
  * Chiang (2008) Chiang, E. 2008, ApJ, 675, 1549
  * Chiang & Youdin (2010) Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493
  * Cui & Bai (2021) Cui, C., & Bai, X.-N. 2021, MNRAS, arXiv:2106.10167
  * Cuzzi et al. (2001) Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496
  * Dittrich et al. (2013) Dittrich, K., Klahr, H., & Johansen, A. 2013, ApJ, 763, 117
  * Drażkowska & Dullemond (2014) Drażkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46
  * Epstein (1924) Epstein, P. S. 1924, Physical Review, 23, 710
  * Flaherty et al. (2020) Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109
  * Flaherty et al. (2018) Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117
  * Flaherty et al. (2017) Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150
  * Fromang & Nelson (2005) Fromang, S., & Nelson, R. P. 2005, MNRAS, 364, L81
  * Fromang & Papaloizou (2006) Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751
  * Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509
  * Gardiner & Stone (2008) —. 2008, Journal of Computational Physics, 227, 4123
  * Gerbig et al. (2020) Gerbig, K., Murray-Clay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125
  * Goldreich & Ward (1973) Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132
  * Goodman & Pindor (2000) Goodman, J., & Pindor, B. 2000, Icarus, 148, 537
  * Güttler et al. (2010) Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56
  * Hartlep & Cuzzi (2020) Hartlep, T., & Cuzzi, J. N. 2020, ApJ, 892, 120
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  * Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  * Hopkins (2016) Hopkins, P. F. 2016, MNRAS, 455, 89
  * Jacquet et al. (2011) Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591
  * Jaupart & Laibe (2020) Jaupart, E., & Laibe, G. 2020, MNRAS, 492, 4591
  * Jennings et al. (2021) Jennings, J., Booth, R. A., Tazzari, M., Clarke, C. J., & Rosotti, G. P. 2021, arXiv e-prints, arXiv:2103.02392
  * 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
  * Johansen et al. (2006) Johansen, A., Klahr, H., & Henning, T. 2006, ApJ, 636, 1121
  * Johansen et al. (2011) —. 2011, A&A, 529, A62
  * Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. 2007, ApJ, 662, 627
  * Johansen et al. (2009a) Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269
  * Johansen et al. (2009b) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009b, ApJ, 704, L75
  * Kataoka et al. (2013) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4
  * Klahr & Schreiber (2020) Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54
  * Klahr & Schreiber (2021) —. 2021, ApJ, 911, 9
  * Krapp et al. (2019) Krapp, L., Benítez-Llambay, P., Gressel, O., & Pessah, M. E. 2019, ApJ, 878, L30
  * Lee et al. (2010a) Lee, A. T., Chiang, E., Asay-Davis, X., & Barranco, J. 2010a, ApJ, 718, 1367
  * Lee et al. (2010b) —. 2010b, ApJ, 725, 1938
  * Lesur (2020) Lesur, G. 2020, arXiv e-prints, arXiv:2007.15967
  * Li & Youdin (2021) Li, R., & Youdin, A. 2021, arXiv e-prints, arXiv:2105.06042
  * Li et al. (2018) Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14
  * Li et al. (2019) —. 2019, ApJ, 885, 69
  * Lin (2021) Lin, M.-K. 2021, ApJ, 907, 64
  * Lin & Youdin (2017) Lin, M.-K., & Youdin, A. N. 2017, ApJ, 849, 129
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17
  * Miyoshi & Kusano (2005) Miyoshi, T., & Kusano, K. 2005, in AGU Fall Meeting Abstracts, Vol. 2005, SM51B–1295
  * Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Okuzumi et al. (2012) Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
  * Okuzumi & Tazaki (2019) Okuzumi, S., & Tazaki, R. 2019, ApJ, 878, 132
  * Pan (2020) Pan, L. 2020, ApJ, 898, 8
  * Pan et al. (2011) Pan, L., Padoan, P., Scalo, J., Kritsuk, A. G., & Norman, M. L. 2011, ApJ, 740, 6
  * Pan & Yu (2020) Pan, L., & Yu, C. 2020, ApJ, 898, 7
  * Pinte et al. (2016) Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25
  * Raettig et al. (2021) Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92
  * Rosotti et al. (2020) Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173
  * Safronov & Zvjagina (1969) Safronov, V. S., & Zvjagina, E. V. 1969, Icarus, 10, 109
  * Schäfer et al. (2020) Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190
  * Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69
  * Schaffer et al. (2021) Schaffer, N., Johansen, A., & Lambrechts, M. 2021, arXiv e-prints, arXiv:2106.15302
  * Schoonenberg & Ormel (2017) Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21
  * Sekiya (1998) Sekiya, M. 1998, Icarus, 133, 298
  * Sekiya & Onishi (2018) Sekiya, M., & Onishi, I. K. 2018, ApJ, 860, 140
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55
  * Simon et al. (2013a) Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ, 775, 73
  * Simon et al. (2013b) Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ, 764, 66
  * Squire & Hopkins (2018) Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011
  * Squire & Hopkins (2020) —. 2020, MNRAS, 498, 1239
  * Stone & Gardiner (2010) Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  * Takahashi & Inutsuka (2014) Takahashi, S. Z., & Inutsuka, S.-i. 2014, ApJ, 794, 55
  * Teague et al. (2018) Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133
  * Tilley et al. (2010) Tilley, D. A., Balsara, D. S., Brittain, S. D., & Rettig, T. 2010, MNRAS, 403, 211
  * Tominaga et al. (2018) Tominaga, R. T., Inutsuka, S.-i., & Takahashi, S. Z. 2018, PASJ, 70, 3
  * Tominaga et al. (2019) Tominaga, R. T., Takahashi, S. Z., & Inutsuka, S.-i. 2019, ApJ, 881, 53
  * Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217
  * Umurhan et al. (2020) Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4
  * Wardle (2007) Wardle, M. 2007, Ap&SS, 311, 35
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Weidenschilling (1980) —. 1980, Icarus, 44, 172
  * Windmark et al. (2012) Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73
  * Xu et al. (2017) Xu, Z., Bai, X.-N., & Murray-Clay, R. A. 2017, ApJ, 847, 52
  * Yang et al. (2017) Yang, C. C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27
  * Youdin & Johansen (2007) Youdin, A., & Johansen, A. 2007, ApJ, 662, 613
  * Youdin (2011) Youdin, A. N. 2011, ApJ, 731, 99
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  * Youdin & Shu (2002) Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494
  * Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81
  * Zhu & Yang (2021) Zhu, Z., & Yang, C.-C. 2021, MNRAS, 501, 467
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57