Dust Settling and Clumping in MRI Turbulent Outer Protoplanetary Disks

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.


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)).

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 ∼mm-cm size has been found to be difficult, and encounters fragmentation and Corresponding author: Ziyan Xu, Xue-Ning Bai ziyanx@pku.edu.cn, xbai@tsinghua.edu.cn 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 stop . It is conventional to cast the stopping time in dimensionless form as τ s ≡ Ωt stop , where Ω is the local orbital frequency of the disk. Particles are considered strongly coupled to the gas for τ s 1, and the coupling is marginal for τ s ∼ 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 τ s ∼ 1. Particles with mm-cm sizes are typically marginally coupled with τ s ∼ 0.1 − 1 in the outer PPDs (a few tens of AU), and more strongly coupled in the inner disks (∼ 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 , 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. 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(Johansen et al. , 2014Simon 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).

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 α 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(Flaherty et al. , 2018Teague 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 α ∼ 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 α/τ s 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 α ∼ 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 α ∼ 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.

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 α 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 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. ( , 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.

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 Ω K . The dynamical equations are written in Cartesian coordinates, withx, y andẑ denoting unit vectors pointing to the radial, azimuthal and vertical directions, and Ω K is along the 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Ω −1 K after inserting dust particles, with projections in three orthogonal directions. Dust density is mapped to brightness, where brighter region indicates higher dust density. z direction. The gas density and velocity are denoted by ρ g and 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 stop . For a particle i with velocity of v i , the drag force on the particle is (u − v i )/t stop per particle unit mass. With magnetic field denoted by B, the MHD equations are written as: where T denotes the total stress tensor: with I being the identity tensor and P is the gas pressure. We assume an isothermal equation of state for the gas so that P = ρ g c 2 s , where c s is the isothermal sound speed. The disk scale height is given by H g = c s /Ω K . The units for magnetic field is chosen such that magnetic permeability µ = 1, and hence magnetic pressure is given by B 2 /2, and current density is given by J = ∇ × B.
The first and second terms on the right side of equation (2) represent Coriolis force and radial tidal gravity, separately. The third term on the right side of (2) is the backreaction from dust particles to gas where andv 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 γ is the coefficient of momentum exchange in ion-neutral collisions and ρ i denotes the ion density. AD is paraterized by Elsässer number: with Am → ∞ 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 10 (Bai & Stone 2011). Bai (2011a,b) found that Am ∼ 1 is widely applicable at the outer PPDs ( 30 AU) where AD is the dominant non-ideal MHD effect. Without stratification, initial gas density ρ 0 is uniform. We thus set c s = Ω K = H g = ρ 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 u K ≡ −(3/2)Ω K xŷ and only evolve the deviation u ≡ u − u K .
Fiducially, we impose a net vertical magnetic field B 0 in our MHD simulations. The strength of B 0 is characterized plasma β, corresponding to the ratio of gas pressure ρ 0 c 2 s 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 β 0 ∼ 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 (7) In the above, the first term on the right side −2ηv K Ω Kx 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 (2). The parameter ∆v K ≡ ηv K represents the deviation from Keplerian rotation due to radial pressure gradient, which is characterized by Π ≡ ηv K /c s (Bai & Stone 2010c). In the absence of gas drag, ∆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.
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 −Ω 2 K z iẑ , 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 τ s ≡ Ω K t stop , where particles with τ s 1 are strongly coupled to the gas while particles with τ s 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).

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 τ s = 0.1 for all our simulations. At outer disk regions, gas drag is in the Epstein regime with t stop = ρ s a/ρ g c s (Epstein 1924), where ρ s and a are the solid density and physical size of dust particles. Fixing τ s is justified because gas density undergoes very little variation in unstratified MRI simulation box with AD (Bai & Stone 2011). In standard disk models, τ 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 Σ g = √ 2πρ 0 H g as if the gas is vertically stratified. The height-integrated dust mean surface density is thus Σ d = ZΣ g . This particle abundance is then evenly divided to individual particle "mass density" in the simulations.
For gas radial pressure gradient, we choose Π = 0.1 by default, which is generally applicable at outer disk ( 30 AU) 1 . We further consider simulations with different values of Π = 0, 0.05, 0.1, 0.2 for comparison. These runs are labeled by "P " followed by a number indicating the Π 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 β 0 = 1.28 × 10 4 . This yields a most unstable MRI wavelength to be about ∼ 0.12H (Bai & Stone (2011), see their Equation 22).
We choose simulation box size to be H × 2H × H in x, y and z dimensions respectively, with a fiducial grid size of 256 3 . Thus, we can well resolve the MRI with ∼ 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 Π = 0.1, with ∼ 26 cells in a characteristic length of Π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 ΠH with ∼ 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 × 1.2H × 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Ω −1 K 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 ∝ exp(−z 2 /2H 2 p0 ), with H p0 = 0.02H g . The particle velocity is initialized according to Nakagawa-Sekiya-Hayashi (NSH) equilibrium (Nakagawa et al. 1986). For pure hydrodynamic The observationally inferred disk surface distribution in the main body of the disk is shallower (Σ ∼ R −1 instead of R −3/2 , e.g., Andrews et al. (2009Andrews et al. ( , 2010, which lowers ∆v K /cs to ∼ 0.1 at 30 AU.
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.

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Ω −1 K 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 Π = 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 ∼ 30Ω −1 K . The resulting H p ∼ 0.08H is comparable with that from the non-ideal MHD simulation with AD for τ 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 ∼ 0.03H (resolved by ∼16 cells) to H p ∼ 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 ∼ 150Ω −1 K , with H p ∼ 0.03H and H p ∼ 0.02H for hydZ2 and hydZ4, separately. The H p in the no-feedback (pure MRI) run is ∼ 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. Figure 2. Time evolution of particle scale height (upper panel) and particle maximum density (lower panel) for simulation runs with Π = 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.
With dust feedback included, particle maximum density is significantly enhanced in our MHD simulations, especially when the solid abundance reaches Z 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 af-fected by the turbulence level in the disk. The level of MRI turbulence in our simulation can be characterized by α 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: Our non-ideal MHD simulation without feedback (ADZ0) gives α = 8.3 × 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 v 2 x /c 2 s = 1.2 × 10 −3 , v 2 y /c 2 s = 1.3 × 10 −3 , and v 2 z /c 2 s = 6.2 × 10 −4 , comparable with the results from Zhu et al. (2015). Simulations with dust feedback

Xu & Bai
included generally give such rms velocities that are close to that from run ADZ0.

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Ω −1 K to the end of each run, being 0.082H, 0.055H, 0.051H, 0.037H and 0.034H 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 τ s , and can be written as (Youdin & Lithwick 2007) where the dimensionless eddy time τ e ≡ Ω K t eddy is defined for the convenience of characterizing the turbulence eddy time t eddy (formally defined in §4.2). The turbulence properties control particle scale height with both the gas vertical rms velocity v z and the eddy time t eddy . The gas vertical diffusion coefficient can be written as For dust particles with Stokes number St ≡ τ s /τ e 1 and τ s τ e 1, the second term of Equation (10) is approximately unity, and Equation (10) is reduced to where α z ≡ D g,z /(c s H) is a dimensionless parameter defined for characterizing vertical turbulent diffusion.

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(ρ p )/ρ 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 v 2 z in non-ideal MHD simulations with and without dust feedback is highly consistent throughout the disk, with v 2 z /c 2 s ∼ 0.0006. This indicates that dust feedback is not controlling the particle layer thickness in MRI turbulent disks by affecting the v 2 z . 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Ω −1 K 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 v 2 z 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 v 2 z /c 2 s ∼ 0.0003, slightly lower than that of the MHD runs. It then falls off away from the midplane 2 . 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.

The effect of t eddy
The fact that dust settles more strongly when turning on feedback, while v 2 z remains largely unchanged raises further question to the origin of the reduced particle vertical diffusion. In addition to the mean square gas vertical velocity v 2 z , particle diffusion in turbulent gas disk is also determined by the correlation time of turbulent fluctuations, defined as where R zz (τ ) ≡ v z (τ )v z (0) is the auto correlation function for v z at a timescale of τ . Although conventional ideal MHD simulations suggest t eddy ∼ Ω −1 K (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 autocorrelation functions for v z in our simulations. We first output our simulation data with time intervals down to 0.1Ω −1 K , and shift the data along the azimuthal direction with a distance of 1.5Ω K xt to correct for Keplerian shear. The auto-correlation function R zz (∆t) is then derived by averaging v z (t)v z (t + ∆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 eddy ≈ 1.5Ω −1 K , which is consistent with the conventional expectation of t eddy ∼ Ω −1 K . For non-ideal MHD simulation without particle feedback, the auto correlation function extends to longer timescales and the correlation time is increased to t eddy ≈ 5.0Ω −1 K 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 ∆t > 2Ω −1 K , and the resulting t eddy ≈ 1.5Ω −1 K , close to the measured eddy time in the pure hydrodynamic simulation. With the measured t eddy and v z 2 , the vertical diffusion coefficients can be calculated with Equation (11). Our calculation gives D g,z ∼ 0.003 for ADZ0, and D g,z ∼ 9 × 10 −4 for ADZ2. For pure hydrodynamic case, turbulent diffusion is inhomogeneous and is strongest at midplane with v 2 z /c 2 s ∼ 3 × 10 −4 , corresponding to D g,z ∼ 4.5 × 10 −4 . Following equation (12), our inferred D g,z would give H p ≈ 0.17H for ADZ0, and H p ≈ 0.095H for ADZ2. Both of these are higher than H p measured in the simulation by a factor of ∼ 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 ∼ 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 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.

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: 10 1 10 0 10 1 10 2 Timescale t / 1  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 eddy ∼ 5Ω −1 K . With dust feedback in ADZ2, the eddy time is reduced and is comparable to the pure hydrodynamic case (t eddy ∼ 1.5Ω −1 K ).
where M 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) ρ 0 = 1.4 × 10 −9 (a/AU ) −2.75 g cm −3 .
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 max(ρ 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 ∼ 200Ω −1 K , 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 hydro-dynamic 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 Π = 0.05, the clumping threshold for particles with τ s ≈ 0.1 is about Z 0.02 (Carrera et al. 2015, but see Li & Youdin 2021 who found a lower Z threshold with improved numerics). However, our simulation adopts a higher pressure gradient of Π = 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.

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 ∼ 180Ω −1 K and T = 330 ∼ 380Ω −1 K , corresponding to the periods of dust maximum density reaching ∼ 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, quasiaxisymmetric density (and hence pressure) variations (Fromang & Nelson 2005;Johansen et al. 2006Johansen et al. , 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. 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 (±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Ω −1 K .
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 3 . 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 (indi-3 Our formulation mimics the outward pressure gradient in the gas by an inward radial force on the particles (see §2.1), making gas azimuthal velocity ηv K = Πcs larger than the actual case. In Figure 5, we have subtracted this component to show vy in the actual situation.

Xu & Bai
cating 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Ω −1 K 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 ∼ 200Ω −1 K in run ADZ2 where ρ 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 max-ima, and we will present further study in §7.1 which will strengthen our conclusions.

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 ρ p is uniformly divided into 1000 bins in logarithmic space between 0.01ρ 0 and 1000ρ 0 , and the distribution function is measured by counting the number of dust particles located in a region with local ρ p in each bin.
We take the first clumping period in ADZ2 (T = 0 − 210Ω −1 K ) for analysis, and measure the distribution function every 40Ω −1 K (except that the first measurement is taken at T = 10Ω −1 K to avoid numerical artifacts at the beginning of inserting particles). Each distribution function is time-averaged within 10Ω −1 K , in order to reasonably reduce the noise level. For comparison, we also present a distribution function of run ADZ0 measured at T = 11 − 20Ω −1 K , 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 ρ p distribution function. The run starts with a smooth Poisson-like distribution prior to clumping (T = 11 − 20Ω −1 K ), with width reflecting the density fluctuations in turbulence. The distribution function peaks at ρ p ∼ ρ 0 , consistent with the expected dust density at midplane that we initially set up in the simulation, ρ p0 ∼ Zρ 0 H/H p0 ∼ ρ 0 , where H p0 = 0.02 is the initial particle scale height set in the simulation. The distribution maximizes at ρ p ∼ 10ρ 0 , consistent with the max(ρ p )/ρ 0 value shown in Figure  2.
The profile then flattens and shifts towards higher density as the clumping starts (T = 81 − 90Ω −1 K ). During the time of steady, maximum clumping (e.g. T = 121 − 130Ω −1 K ), the density distribution peaks at ρ p ∼ 100ρ 0 , and maximizes at ρ p ∼ 1000ρ 0 , again consistent with Figure 2. After the period of clumping, the particle density distribution shifts back after T = 161 − 170Ω −1 K , 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.  No FB  T=11 -20  T=41 -50  T=81 -90  T=121 -130  T=161 -170 T=201 -210 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 Ω −1 K , and marked in different colors (see legend). Each distribution function is time-averaged within 10Ω −1 K . For run ADZ0, the distribution function is measured and time-averaged between T = 11 − 20Ω −1 K .
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 ρ p ∼ 1000ρ 0 , a large fraction of particles also reside in relatively dense clumping regions (∼ 100ρ 0 ).

RADIAL TRANSPORT OF DUST PARTICLES
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 + ∆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Ω −1 K 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.

Radial drift
In laminar disks, for a sub-Keplerian speed ∆v K , dust particles suffer from radial drift. For particles with dimensionless stopping time of τ s , and ignore particle backreaction, the radial and azimuthal drift velocity relative to the local Keplerian velocity is given by (Nakagawa et al. 1986 For our choice of ∆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 ∆x 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 ∆x − T curve does not closely follow the theoretical expectation anymore. Around T ∼ 50 − 150Ω −1 K , 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 ∆x − 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Ω −1 K , indicating the bulk of the particles are trapped and the radial drift is halted.

Radial diffusion
Under the assumption of random walk, the dust radial diffusion coefficient D x can be measured by where σ x is the rms of particle radial displacement distribution. We measure the σ 2 x as a function of time for runs ADZ0 and ADZ2, presented in the upper right panel of Figure 7.
For ADZ0, the σ 2 x − T curve closely follows the linear relation, and by fitting the slope, we obtain D x ≈ 1.2 × 10 −3 . Assuming the gas diffusion coefficient on the radial direction D g,x is largely the same (given τ s 1), the value is lower than D g,z = 0.003 that we measured for ADZ0. This is overall consistent with the result of  T=30  T=60  T=90  T=120  T=150  T=180   0  25  50  75  100  125  150  175  x /H 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 Ω −1 K , 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. Zhu et al. (2015), where D g,x is generally lower (by a factor of ∼ 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 σ 2 x −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 σ 2 x − T curve is largely consistent with that of ADZ0 before T ∼ 50Ω −1 K , prior to strong particle clumping. During the clumping period (T = 50 − 150Ω −1 K ), the σ 2 x − 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, σ 2 x 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 σ 2 x . 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Ω −1 K . 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 ∼ 50 − 150Ω −1 K (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. 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 Π = 0, 0.05, 0.1, and 0.2, separately. No particle clumping is seen under Π = 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 Π (Bai & Stone 2010c 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. free energy arises from background pressure gradient. Thus, particle clumping seen in our Π = 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Ω −1 K at x ∼ −0.25H with an amplitude of ∼ 0.03ρ 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 ∼ 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 ∼ 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 Π = 0.2, a pressure bump that is largely similar to that in the ADZ2P0 run is formed, also at T = 0 − 30Ω −1 K , x ∼ −0.25H, with the same amplitude of ∼ 0.03ρ 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.

The effect of resolution
Our fiducial resolution of 256/H fairly well resolves the particle layer (H p ∼ 0.03H for ADZ2) with ∼ 16 cells, as well as the width of the particle clumping region ( 0.1H, see e.g. Figure 5) with ∼ 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.

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 clump must not exceed the Hill length scale, given by their equation (13) r clump < 2π 9 where Σ d,c is the surface density of the dust clump. In order to overcome the turbulent diffusion, the clump radius r clump must be large enough so that the timescale for particles to be diffused over r clump is longer than the gravitational collapsing timescale. Assuming turbulent diffusion coefficient is characterized by α, their equation (9) gives Combining equations (19) and (20) leads to By adopting Toomre Q parameter (Toomre 1964) and the local metallicity enhancement their equation (70) gives a collapsing criterion of 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 = 100, together with Z = 0.02, and H/H p ∼ 20, analytical prediction by Gerbig et al. (2020) gives a critical collapsing value of Q crit ∼ 27, comparable to the Toomre Q ∼ 30 in MMSN (Weidenschilling 1977, Hayashi 1981. 4 Considering the fact that the diffusion α in equation (24) is expected to be measured within the dust clump, which could largely deviate from the global α 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.

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.

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 nonideal 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 to-wards 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).