# Turbulent Dust-trapping Rings as Efficient Sites for Planetesimal Formation

Ziyan Xu Univ Lyon, Univ Lyon 1, ENS de Lyon, CNRS, Centre de Recherche
Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France  Kavli
Institute for Astronomy and Astrophysics, Peking University, Yiheyuan 5,
Haidian Qu, 100871 Beijing, China  Department of Astronomy, Peking University,
Yiheyuan 5, Haidian Qu, 100871 Beijing, China  Xue-Ning Bai Institute for
Advanced Study, Tsinghua University, 100084 Beijing, China  Department of
Astronomy, Tsinghua University, 100084 Beijing, China

###### Abstract

Recent observations of protoplanetary disks (PPDs) in the sub-mm have revealed
the ubiquity of annular substructures, indicative of pebble-sized dust
particles trapped in turbulent ring-like gas pressure bumps. This major
paradigm shift also challenges the leading theory of planetesimal formation
from such pebbles by the streaming instability, which operates in a pressure
gradient and can be suppressed by turbulence. Here we conduct three-
dimensional local shearing-box, non-ideal magnetohydrodynamic (MHD)
simulations of dust trapping in enforced gas pressure bumps including dust
backreaction. Under the moderate level of turbulence generated by the
magnetorotational instability (MRI) with ambipolar diffusion that is suitable
for outer disk conditions, we achieve quasi-steady states of dust trapping
balanced by turbulent diffusion. We find strong dust clumping in all
simulations near the gas pressure maxima, reaching a maximum density well
above the threshold of triggering gravitational collapse to form
planetesimals. A strong pressure bump concentrates dust particles towards bump
center. With a weak pressure bump, dust can also concentrate in secondary
filaments off the bump center due to dust backreaction, but strong clumping
still occurs mainly in the primary ring around the bump center. Our results
reveal dust-trapping rings as robust locations for planetesimal formation in
outer PPDs, while they may possess diverse observational properties.

protoplanetary disks - planet formation - planetesimals - magnetohydrodynamics

††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 INTRODUCTION

Planetesimals are $\sim 1-1000$ km sized bodies in protoplanetary disks (PPDs)
as building blocks of planets. Following the growth of dust grains to mm-cm
sizes in PPDs, planetesimals are believed to form from such grown dust
particles, representing the intermediate stage of planet formation. This stage
is also considered to be the most difficult, as surface forces can no longer
stick particles together upon collisions, nor does self-gravity which
dominates among solids of larger sizes (Chiang & Youdin, 2010). Generally,
planetesimal formation is understood as a two-phase process, where particles
first concentrate into high-density clumps, followed by their gravitational
collapse. What drives the formation of particle clumps thus becomes the key to
understanding planetesimal formation. Regardless of the specific mechanism,
particles of mm-cm sizes strongly interact with gas via aerodynamic drag, and
the gas dynamics in PPDs, especially disk turbulence, must play an important
role.

The leading theory of planetesimal formation is the streaming instability (SI,
Youdin & Goodman, 2005), resulting from the two-way drag force between gas and
dust in the presence of a background pressure gradient in PPDs. Hybrid gas-
dust simulations found that the non-linear evolution of the SI leads to dust
clumping (e.g., Johansen et al., 2009; Bai & Stone, 2010a; Carrera et al.,
2015; Yang et al., 2017), provided the dust abundance $Z$ (ratio of dust to
gas surface density) exceeds some threshold (typically $\gtrsim 0.02$, but see
Li & Youdin, 2021). While highly successful, vast majority of SI simulations
assumed no external turbulence. However, recent studies found that even modest
level of external disk turbulence suppresses the SI (Umurhan et al., 2020;
Chen & Lin, 2020; Gole et al., 2020), challenging the SI paradigm.

Figure 1:  A snapshot of 3D visualization of the dust spatial distribution in
our simulation with strong gas bump (run Z2H) at the time of
$T=1380\Omega_{K}^{-1}$ after inserting dust particles. Dust density is mapped
to brightness, where brighter region indicates higher dust density. Also shown
in the figure are vertically averaged gas density projected to the $x$-$y$
($R$-$\phi$) plane, and vertically and azimuthally averaged radial profiles of
dust and gas density, with amplitudes of the profiles re-scaled for better
illustration.

The outer disk regions ($\gtrsim 20$ AU) are the most accessible to spatially-
resolved disk observations. At sub-mm (where mm-sized dust radiates most
efficiently), the prevalence of ring-like substructures has recently been
firmly established (e.g., Andrews et al., 2018a; Long et al., 2018), widely
believed to be indicative of dust trapping in pressure bumps (Dullemond et
al., 2018)111Unless otherwise noted, “pressure bump” refers to dust trapping
pressure bumps in this paper, i.e., strong enough to overcome the background
pressure gradient and generate local pressure maxima.. Dust trapping in local
pressure maximum is firstly demonstrated by Whipple (1972), and the formation
mechanism of pressure bumps can be planet-disk interaction (see review of
Paardekooper et al., 2022) or non-planet origins such as MHD zonal flows
(e.g., Suriano et al., 2018; Riols et al., 2020; Cui & Bai, 2021). While there
are direct and indirect evidence pointing to weak disk turbulence (Pinte et
al., 2016; Flaherty et al., 2017, 2018; Teague et al., 2018), modeling dust
ring widths suggest modest level of turbulence (Dullemond et al., 2018;
Rosotti et al., 2020). Together, these observational constraints lead to the
emerging picture that most mm-sized dust particles in outer PPDs are most
likely concentrated in weakly-to-modestly turbulent, largely-axisymmetric gas
pressure bumps. The strongly enhanced local dust abundance in such pressure
bumps yields favorable condition for planetesimal formation. However, given
that the SI is unexpected to operate without pressure gradient, and can be
suppressed by external turbulence, it is unclear whether it will further lead
to dust clumping and hence planetesimal formation.

To address the feasibility of planetesimal formation in more realistic
environment of turbulent dust rings, we conduct hybrid gas-particle
simulations in the standard shearing-box framework that mimics a local patch
of the disk. It is anticipated that the magneto-rotational instability (MRI)
is one primary mechanism for driving turbulence in the outer PPDs. Due to very
low level of disk ionization, the MRI is weakened by ambipolar diffusion (AD)
as the dominant non-ideal magnetohydrodynamic (MHD) effect in the outer disk.
Turbulence in our simulations are self-consistently generated by the MRI with
AD. We further implemented a novel forcing prescription to drive a pressure
bump in the center of our simulation box with essentially no impact on disk
turbulence. Our simulations generalize from our earlier work (Xu & Bai, 2022,
hereafter paper I) by imposing and enforcing the pressure bump, allowing us to
achieve a steady-state configuration with dust trapping towards pressure bumps
balanced by dust diffusion in realistic disk turbulence in a self-consistent
and well-controlled manner. They differ from the recent works of (Carrera et
al., 2021, 2022; Carrera & Simon, 2022), who conducted simulations with a
similar procedure to enforce the pressure bump but without external
turbulence, thus cannot achieve such steady-state balance of dust trapping and
diffusion as envisioned in recent disk observations (Dullemond et al., 2018;
Rosotti et al., 2020).

This letter is organized as follows. We first introduce the implementation of
the pressure bump and simulation setup in §2, then present our main results in
§3. In §4, we will further analyze the results, focusing on dust clumping
conditions and ring dynamics, followed by comparison with recent studies and
discussion on observational implications. We summarize and discuss future
prospects in §5.

Table 1: List of simulation runs.

Runs | $Z$ | $N_{x}\times N_{y}\times N_{z}$ | $A_{bump}$ | $H_{d}/H$ | $w_{d}/H$
---|---|---|---|---|---
Z0 | 1e-30 | 512 $\times$ 512 $\times$ 64 | 0.5 | 0.067 | 0.55
Z2 | 0.02 | 512 $\times$ 512 $\times$ 64 | 0.5 | 0.014 | 0.16
Z0H | 1e-30 | 1024 $\times$ 1024 $\times$ 128 | 0.5 | - | -
Z1H | 0.01 | 1024 $\times$ 1024 $\times$ 128 | 0.5 | 0.025 | 0.25
Z2H | 0.02 | 1024 $\times$ 1024 $\times$ 128 | 0.5 | 0.017 | 0.16
Z0w | 1e-30 | 512 $\times$ 512 $\times$ 64 | 0.25 | 0.075 | 0.83
Z1w | 0.01 | 512 $\times$ 512 $\times$ 64 | 0.25 | 0.031 | 0.87
Z2w | 0.02 | 512 $\times$ 512 $\times$ 64 | 0.25 | 0.022 | 0.74

  * •

All simulations have a boxsize of $8H\times 16H\times H$, with
$\beta_{0}=12800$, $Am=2$. Gas bump width $w_{bump}=2H$, with no background
pressure gradient. Particle number is 2 per cell on average, with particle
stopping time $\tau_{s}=0.1$.

  * •

All parameters in this table are initial setups, except that $H_{d}$ and
$w_{d}$ are measured from simulation.

##  2 Method and Simulations

We use the Athena MHD code (Stone et al., 2008) to perform 3D hybrid gas-
particle local shearing-box simulations considering non-ideal MHD with
ambipolar diffusion. The formulation and basic simulation setup is similar to
paper I, which we only describe briefly here. The $x,y,z$ coordinates in a
local shearing-box correspond to radial, azimuthal and vertical directions in
a global disk, as illustrated in Figure 1. We use an isothermal equation of
state with sound speed $c_{s}=1$ and local Keplerian frequency $\Omega_{K}=1$
in code units, which define our length unit of disk scale height
$H=c_{s}/\Omega_{K}=1$. Gas density, velocity and pressure are denoted by
$\rho,{\boldsymbol{v}}$ and $P=\rho c_{s}^{2}$, respectively. Dust particles
are characterized by their dimensionless stopping time $\tau_{s}\equiv
t_{s}\Omega_{0}$, and the height-integrate dust-to-gas mass ratio $Z$ for
background gas surface density. Our simulations focus on the midplane region,
and hence we neglect vertical gravity in the gas (unstratified), but vertical
gravity is included for particles to allow dust settling. We impose a net
vertical magnetic flux characterized by $\beta_{0}=12800$, the ratio of
background gas pressure to magnetic pressure of the net vertical field.
Ambipolar diffusion is characterized by the Ambipolar Elsasser number $Am$,
where we fix $Am=2$ appropriate for typical outer disk conditions.

###  2.1 Implementation of the gas pressure bump

On top of the simulation setup in paper I, we further implement a gas pressure
bump by a prescribed forcing procedure. Let $\rho_{0}=1$ be the constant
background gas density. The forcing will redistribute the density, and we aim
to achieve a Gaussian-like density (and hence pressure) bump in our simulation
box of the following form

$\frac{\rho(x)}{\rho_{0}}=1+A_{\rm bump}\exp\bigg{[}-\frac{x^{2}}{2w_{\rm
bump}^{2}}-0.2\bigg{(}\frac{x^{2}}{2w_{\rm bump}^{2}}\bigg{)}^{4}\bigg{]},\\\
$ (1)

where $A_{\rm bump}$ is the dimensionless amplitude of the pressure bump,
$w_{\rm bump}$ is the effective width of the bump. The first term in the
exponential gives the overall Gaussian profile, while the second term damps
the profile towards the domain edges to satisfy the shearing-periodic radial
boundary condition.

The radial pressure gradient associated with this bump must be balanced by an
azimuthal flow:

$\frac{v^{\prime}_{y,{\rm bump}}(x)}{c_{s}}=\frac{H}{2}\frac{d\ln\rho(x)}{dx}\
,$ (2)

where the prime ′ in $v_{y}$ denotes deviation from the Keplerian velocity
$v_{K}=-(3/2)\Omega_{K}x$. Note that in our simulations, as we focus on
regions near the pressure maxima, we do not impose any background pressure
gradient.

To generate and sustain the bump, we apply an azimuthal forcing in our
simulation, so that the azimuthal velocity relax to the desired profiles
within a timescale of $t_{\rm relax}=60\Omega_{K}^{-1}$:

$\Delta v^{\prime}_{y}(x)=[v_{y,{\rm
bump}}(x)-\overline{v^{\prime}_{y}}(x)][1-\exp(-\Delta t/t_{\rm relax})],$ (3)

where $\overline{v^{\prime}_{y}}(x)$ is azimuthally and vertically averaged.
Our forcing procedure is analogous to that of Carrera et al. (2021), but with
several differences. First, our forcing acts on azimuthal velocity (momentum)
only, thus is mass conserving. It directly mimics the torque exerted on the
disk responsible for making the bump (be it of planet origin, or non-planet
origin as in Cui & Bai, 2021), and the desired density profile is a direct
consequence from this forcing. Second, our forcing is homogeneously applied in
the $y$ and $z$ directions, as opposed to relaxation at individual grid cells
which introduces local friction. Third, our relaxation time is much longer, as
opposed to $\sim 1\Omega_{K}^{-1}$. These facts minimizes the impact of our
forcing procedure to disk turbulence, especially avoiding artificial damping
at small scales.222Our choice of relaxation time is generally consistent with
the typical timescale for magnetic flux evolution (and thus bump evolution)
indicated in Cui & Bai (2021). We further tested $t_{\rm
relax}=30\Omega_{K}^{-1}$ and $120\Omega_{K}^{-1}$ respectively (not shown in
paper), and the result is not sensitive to the change of $t_{\rm relax}$ in
our simulation in this paper.

###  2.2 Simulation Procedures and Parameters

The list of our simulation runs are shown in Table 1. We explore strengths of
the gas pressure bump, presence/absence of dust feedback, different levels of
solid abundance, and simulation resolutions. The simulation setup and
parameters are described as follows.

We first run simulations without dust particles for a time of
$t_{0}=120\Omega_{K}^{-1}$ for the MRI to grow and fully saturate into
turbulence, following our previous work (paper I). The gas pressure bump is
initiated at the beginning of the simulation, with bump width $w_{\rm
bump}=2H$. We explore two cases: strong and weak bumps, with bump amplitude
$A_{\rm bump}=0.5$ and $0.25$ separately, and the weak bump runs are marked by
“w” in the run name.333We have examined that these bump amplitudes are below
the threshold for triggering the Rossby wave instability (RWI) for all
azimuthal wave numbers under global settings (Ono et al., 2016), and we
anticipate the presence of the MRI turbulence also makes our pressure bump
less prone to the RWI. Our choice of bump width generally agrees with
observational constraints (Dullemond et al., 2018; Rosotti et al., 2020), and
both the bump width and amplitude in our simulation are reasonably consistent
with the results of Cui & Bai (2021).

Figure 2:  Time evolution of dust maximum density (upper panel), dust ring
width (middle panel), and dust scale height (lower panel) for all our
simulation runs. Dust maximum density for no-feedback runs (Z0, Z0H, and Z0w)
are scaled to $Z=0.02$.

Dust particles are then inserted into the simulation box at the time of
$t_{0}=120\Omega_{K}^{-1}$ without backreaction to the gas, with random
horizontal positions and a vertical distribution following a Gaussian profile
$\propto\exp(-z^{2}/2H_{\rm d0}^{2})$, with $H_{\rm d0}=0.02H$. We define
$T=T_{\rm sim}-t_{0}$, where $T_{\rm sim}$ is the actual simulation time, and
use $T$ instead of $T_{\rm sim}$ throughout this paper. We use a single
particle species with fixed stopping time $\tau_{s}=0.1$, corresponding to mm-
cm size particles in the midplane regions of the outer disk in standard disk
models (e.g., Chiang & Youdin, 2010). As dust particles passively follow the
gas, settling balanced by turbulent diffusion in the vertical direction makes
them maintain certain scale height (see paper I), while the balance between
trapping toward pressure maxima and turbulent diffusion in the radial
direction allow them to concentrate into a ring with certain ring width
$w_{d}$ (see Figure 2), approaching a steady state after
$T=600\Omega_{K}^{-1}$. Next, we gradually introduce dust feedback to the gas
by increasing dust abundance gradually to the desired value with over a time
interval of $30\Omega_{K}^{-1}$. We consider two global dust abundances
$Z=0.01$ and $Z=0.02$, and these runs are labeled by “$Z$” followed by a
number indicating the $Z$ values. In all simulations, there are on average 2
particles per cell, ensuring sufficient particle statistics for studying dust
dynamics.

We choose simulation box size to be $8H\times 16H\times H$ in $x$, $y$ and $z$
directions respectively. We choose a fiducial resolution of 64/H (in the $x$
and $z$ directions; the resolution in the $y$ direction is half of that). We
further conducted simulations with a higher resolution of 128/H for
comparison, labeled with an “$H$” at the end of the run names.

Figure 3: Snapshots of projected dust density in the $x$-$z$ (upper four
panels) and $x$-$y$ (lower four panels) planes for runs Z0, Z2H, Z0w, Z2w, at
time $T=1380\Omega_{K}^{-1}$, when the secondary filament in Z2w is clearly
seen. Dust density in the $x$-$z$ plane is azimuthally averaged over the whole
domain, while in the $x$-$y$ plane, dust density is vertically averaged around
the midplane ($z=\pm 0.05H$) to highlight dust structures.

##  3 Results

An overview of our simulations is shown in Figure 1, featuring a 3D
visualization of the dust spatial distribution for our run Z2H at time
$T=1380\Omega_{K}^{-1}$. Our forcing procedure generates a wide and steady gas
bump, with additional density fluctuations (i.e., density waves) as a result
of the MRI turbulence (Heinemann & Papaloizou, 2009). Clearly, the dust
concentrates to the center of the gas bump, whose width is much narrower than
the gap bump width. Besides settling and trapping, and the dust distribution
is highly non-smooth and clumpy under background turbulence.

In this section, we present our main simulation results. We measure the
maximum dust density $\rho_{d,{\rm max}}$, normalized by background gas
density $\rho_{0}$, as a proxy for dust clumping. Dust density in no-feedback
runs (Z0 and Z0H) are rescaled to Z=0.02 for better comparison. For reference,
we consider $\rho_{d,{\rm max}}\gtrsim 200\rho_{0}$ as a threshold condition
for triggering gravitational collapse, and hence planetesimal formation, which
is applicable for typical outer disk conditions (see Equation 16 and Section
5.1 in paper I). We further measure the dust ring width $w_{d}$ and dust scale
height $H_{d}$, defined as the rms value of the radial and vertical
coordinates for all dust particles. The time history of $\rho_{d,{\rm
max}}/\rho_{0}$, $w_{d}$ and $H_{d}$ from all simulations is shown in Figure
2. The time-averaged values of $w_{d}$ and $H_{d}$ are presented in Table 1,
averaged within $T=800-1500\Omega_{K}^{-1}$. These are accompanied by Figure
3, where we show the projected dust density in the $x-y$ and $x-z$ planes for
Z0, Z2H, Z02 and Z2w runs at the saturated states. We find that the results
for the strong bump and weak bump cases are distinct, which we discuss
separately below.

###  3.1 The strong bump case: dust clumping in a single thin ring

Top panel of Figure 2 shows that dust density presents modest concentration in
the absence of dust feedback, with $\rho_{d,{\rm max}}/\rho_{0}$ reaching
$\sim 100$, approaching yet barely reach the threshold for triggering
planetesimal formation. This is to due to a combination dust trapping and
intrinsic fluctuations. Reading from the second and third panels of Figure 2,
we find a dust ring width of $w_{d}\sim 0.55H$, and a scale height of
$H_{d}\sim 0.067H$ (comparable to $\sim 0.08H$ found in paper I without
enforced pressure bump). This yields a midplane density enhancement factor of
$\sim 50$ at bump center. Indeed, as shown in Figure 3, the mean midplane
density at bump center for run Z0 (renormalized assuming $Z=0.02$) is about
$\overline{\rho_{d}}\sim 1$. We do not observe obvious dust clumping from the
averaged dust densities along the vertical and azimuthal directions, and thus
attribute an additional factor of $\sim 50-100$ enhancement in $\rho_{d,{\rm
max}}$ from localized density fluctuations. The results are largely
independent of resolution, where $\rho_{d,{\rm max}}$, $w_{d}$ and $H_{d}$ are
consistent with each other between run Z0 and Z0H within fluctuation levels.

Turning on dust feedback leads to a rapid boost in $\rho_{d,{\rm
max}}/\rho_{0}$, well exceeding $10^{3}$ (and hence the Roche density) for all
runs with Z=0.01 (Z1H) and Z=0.02 (Z2 and Z2H), which is therefore expected to
trigger planetesimal formation. There are also visible bright spots seen in
the bottom panel of Figure 3 for run Z2H, suggesting that such clumping is
physical instead of random Poisson fluctuations. There is a trend of higher
resolution leading to higher level of dust clumping, where for runs Z1H and
Z2H, $\rho_{d,{\rm max}}/\rho_{0}\gtrsim 5000$, a factor of $\sim 50$
enhancement compared to the no feedback case. We have further examined the
cumulative probability density distribution (CDF, see Appendix A), which
reveal that dust density is indeed substantially boosted with a significant
fraction (a few tens of percent) residing in over-Roche-density regions.

Somewhat surprisingly, with dust feedback, we observe a dramatic reduction of
the ring width to $w_{d}\approx 0.2H$, $\sim 3$ times narrower than that
without feedback. Moreover, the thickness of the dust layer is also reduced
substantially by a factor of $3-4$. Examining Figure 3 reveal that this
reduction is not homogeneous: there is a strong/spiky dust concentration and
more substantial settling at the bump center (see also the top left panel of
Figure 4 and Figure 3 in paper I), which reduces the measured $w_{d}$ and
$H_{d}$ towards lower values. We note that the reduction of $H_{d}$ was also
reported in paper I, which was mainly attributed to reduction of turbulence
correlation time. However, in the presence of pressure bump, $H_{d}$ is
further reduced by a factor of nearly $\sim 2$. The dramatic reduction of ring
width, as well as further promotion of dust settling, are likely the outcome
of strong dust mass loading in the bump center (see further discussions in
§4.2 and Appendix C). The results are convergent with resolution, as the
values of $w_{d}$ and $H_{d}$ are consistent with each other between runs Z2
and Z2H.

###  3.2 The weak bump case: diverse ring properties

The most prominent feature for the weak bump case compared to the strong bump
case is the dust ring properties. The rest of the results, especially dust
clumping and settling, are similar to the case with a strong bump.

The middple panel of Figure 2 shows that without dust feedback, $w_{d}$ in the
weak bump case is generally higher compared to the strong bump case by a
factor of up to $\sim 2$. This is easily accounted for due to the reduced gas
pressure gradient in the two sides of the gas bump, resulting in weaker dust
trapping. This is also seen in Figure 3, showing a smooth and radially
extended distribution of dust from run Z0w, together with azimuthally extended
stripes akin to the MRI density waves.

When turning on dust feedback, interestingly, the dust ring width $w_{d}$,
measured from the rms radial coordinates, remains similar for all cases.
However, when examining the dust spatial distribution in Figure 3, we find
that instead of concentrating into a single ring at the bump center, dust
dynamics in weak bump is more complicated, potentially forming multiple fine-
scale filaments within the gas bump. In the case of run Z2w, we observe two
dust rings where the primary is close to bump center, and a secondary ring
more than $H$ away from the bump center, in addition to a few minor filaments.
The formation of this more complex ring/filamentary structure is likely caused
by the interplay between dust feedback and gas dynamics within the bump. In
addition, these structures generally stay while slowly evolve over tens to
hundreds of $\Omega_{K}^{-1}$. Details of the formation and evolution of
(multiple) fine-scaled ring structures will be further discussed in §4.2 and
Appendix B.

##  4 Discussion

###  4.1 Clumping conditions

For both the strong and weak bumps, sufficient dust clumping is seen to
trigger planetesimal formation. Here we further discuss the clumping
conditions to gain insight into the clumping mechanism. In doing so, we show
in the top four panels of Figure 4 the radial profiles of vertically and
azimuthally-averaged dust density, gas density and Keplerian-subtracted
azimuthal gas velocity $v^{\prime}_{y}$, at the very midplane of our
simulation ($\pm$ 1 cell). The latter two are useful for quantifying the
pressure gradient, and we further show the time evolution of the radial
profiles of $\overline{v^{\prime}_{y}}$ in four representative runs in the
bottom panels.

#### Dust feedback

Dust feedback is required for dust clumping in our simulations, as no clumping
is seen in our run without feedback (Z0 and Z0w). As seen in the upper panels
of Figure 4, the dust midplane density in the no-feedback cases generally
exhibits a Gaussian-like radial profile, conforming to the standard
expectation of concentration balancing turbulent diffusion (Dullemond et al.,
2018). Including dust feedback, the radial density profiles substantially
deviate from Gaussian, showing density cusps at the location of dust rings,
where dust density is enhanced by a factor of near $100$ relative to the no-
feedback counterparts. The cusps are associated with the locations of dust
clumping, and can be the consequences of dust clumping itself.

#### Solid abundance

Dust clumping occurs in all simulations for mean solid abundance $Z$ (averaged
over the simulation domain) as low as $0.01$. Note that due to dust trapping,
the effective solid abundance (i.e., vertically integrated dust-to-gas mass
ratio) in the bump region is much higher. We may estimate the enhancement
factor to be $\sim 4-8$ in the weak bump and strong bump cases based on the
value of $w_{d}$, corresponding to $Z\gtrsim 0.04$ around bump center in all
our simulations. Therefore, the overall solid abundance $Z$ in the bulk disk
appears less relevant, as dust trapping can easily bring sufficient amount of
dust to the bump region.

#### Local pressure maxima

In the strong bump case (runs Z1H, Z2H), dust clumping primarily occurs at the
central cusp, which is marked in a dashed line in Figure 4. We further see
from the radial profiles of $\rho_{g}$ and $v^{\prime}_{y}$ that this location
exactly coincides with the location of the gas pressure maxima, corresponding
to $v^{\prime}_{y}=0$. In the weak bump case, the primary ring close to the
bump center is also located at the pressure maxima, which in this case is
widened into a plateau with a width of $\sim 0.5H$ as seen from the gas
density profile as well as an extended region with $v^{\prime}_{y}\approx 0$.
The cusp where strong dust clumping occurs is located within the plateau,
again consistent with the requirement of local pressure maxima, or at least a
zero-pressure-gradient region being a necessary condition for strong dust
clumping. These findings are consistent with the results reported in paper I,
where the pressure bumps were not enforced but generated by the MRI zonal
flows (albeit somewhat unrealistic in local unstratified simulations).

#### Weak clumping in the secondary ring

For the weak bump case of Z1w and Z2w, in addition to the primary ring near
the bump center, an off-centered “secondary” ring is identified at
$x\sim-1.2H$. We show in Appendix A that the dust density CDF in the primary
ring is similar to that in the strong bump case, with strong clumping.
However, the dust CDF of the secondary ring is significantly different. The
mean dust density is more than a factor of $\sim 10$ lower, and only a very
small fraction ($\lesssim 3\%$) of particles reside in over-Roche-density
regions (which exhibit as individual spots seen in the bottom right panel of
Figure 3), indicating that clumping is much less efficient. Examining the
density (and hence pressure) profiles in the top right panel of Figure 4, the
location of the secondary ring is clearly not at the pressure maximum.

Overall, for both strong and weak bump cases, dust clumping is tightly related
to the presence of pressure maxima as long as dust backreaction is included.
While secondary rings/filaments can form for the weak bump case, they
generally reside outside of the pressure maxima and do not show
prominent/efficient clumping.

Figure 4: Upper four panels: radial profiles of vertically and azimuthally
averaged dust density ($\rho_{d}$, solid curves in upper two panels), gas
density ($\rho_{g}$, dashed curves in upper two panels), and Keplerian-
subtracted azimuthal gas velocity ($v^{\prime}_{y}$, solid curved in lower two
panels), at the midplane ($\pm 1$ cell from $z=0$), averaged between
$T=800-1500\Omega_{K}^{-1}$. Vertical dashed lines mark the locations of the
dust rings (peak of dust density), with colors corresponding to that used for
radial profiles of each run.
Lower four panels: time evolution of the radial profile of $v^{\prime}_{y}$,
between $T=800-1500\Omega_{K}^{-1}$.

###  4.2 Interplay between dust feedback and gas bump

#### Bump/ring profile

In addition to dust clumping, the interaction between dust and gas also
influences gas bump properties, leading to diversity in dust structures, as we
discuss here.

Dust feedback changes the radial profile of the gas bump. This can be seen in
the yellow region in the time-evolution map of $v^{\prime}_{y}$ at the
midplane in Figure 4 by comparing the upper and lower panels. In the strong
bump case, as well as around the primary ring in the weak bump case, dust
feedback leads to a plateau of Keplerian rotation at the bump center, in
comparison to that in run Z0 and Z0w. This is consistent with the previous
finding in (Huang et al., 2020).

The secondary dust ring, despite not in the pressure maxima, can be sustained
for at least dozens of orbits without drifting towards the bump center.
Examining the profile of $v^{\prime}_{y}$, we find that its value shows a dip
at the secondary dust filament, implying near-Keplerian rotation. This is the
result of feedback from the dusty filaments modifying the rotational profile
of the background gas. It occurs in our weak bump runs likely because the gas
pressure gradient is not strong enough to force the heavily loaded dust
filaments to drift to bump center, resulting in dust trapping.444 One might
argue that the formation of secondary dust rings seen in run Z1w and Z2w be a
consequence of the special timing when we turn on dust feedback due to pre-
existing local density enhancements from the MRI turbulence (e.g., zonal
flows). We thus further perform a simulation where dust feedback is turned on
at a much later time $T_{sim}=1200\Omega_{K}^{-1}$. We still find similar
secondary dust rings/filaments in this case. Such a location shares
similarities with a pressure maxima in that there is almost no radial or
azimuthal drift between gas and dust, and is likely responsible for the weak
clumping discussed in §4.1. These secondary rings also evolve over the
timescale of tens to hundreds of orbits, as we discuss further in Appendix B.

#### Turbulent diffusion

The enhanced dust mass loading in turbulent dust rings also affects the
properties of background MRI turbulence. With feedback included, we find that
the gas rms velocities in both the radial and vertical directions are largely
reduced in the dust ring (see details in Appendix C). Such reductions can
sufficiently explain the results reported in Section 3.1, where the dust rings
in the strong bump case becomes much narrower when dust feedback is included,
and that dust becomes even more settled than that reported in paper I.

###  4.3 Comparison with other recent studies

Recently, Carrera et al. (2021, 2022); Carrera & Simon (2022) conducted a
series of hydrodynamic simulations of the SI, imposing a pressure bump by a
similar (albeit more dissipative) forcing procedure as ours but without
external turbulence. However, the pressure bump in most of their simulations
are not sufficiently strong to overcome the background pressure gradient. As a
result, particles get accumulated (but not trapped) in regions with lower
pressure gradient, making such regions the most favorable for clumping and
planetesimal formation (e.g., Bai & Stone, 2010b). In other words, their
simulations identified the favorable conditions for dust clumping by the
conventional SI in the presence of radial-dependent pressure gradient. Our
simulations, on the other hand, focus on planetesimal formation in dust-
trapping rings, i.e., the pressure maxima. The problem is only well-defined in
the presence of external turbulence, as expected from observations (Dullemond
et al., 2018; Rosotti et al., 2020), otherwise, dust would indefinitely get
concentrated towards ring center over the trapping timescale (this issue may
also apply to the analytic work of Auffinger & Laibe (2018)), forming
planetesimals along the way. Overall, the physical picture in our simulations
is distinctly different from the works of Carrera et al. (2021, 2022), and
likely better reflect outer disk conditions. Together with paper I, we showed
that for the MRI turbulence, the condition for planetesimal formation requires
the presence of pressure maxima, again different from the conventional SI
scenario.

The outer disk may also be subject to turbulence driven by the vertical shear
instability (VSI, Nelson et al., 2013). Using 2D global simulations, Schäfer
et al. (2020) found that the SI can co-exist with the VSI, and the VSI even
facilitates the development of the SI and dust clumping, likely due to the
overdensities induced by the VSI. We note that the overdensities from the VSI
itself is also insufficient to overcome the background pressure gradient. This
is different from the case with the MRI, where we found that dust clumping
requires the presence of pressure maxima, suggesting that dust clumping
conditions likely depends on the nature turbulence. Lehmann & Lin (2022) found
that the presence of an initial pressure bump (but without forcing) further
facilitates dust concentration, though their resolution is insufficient to
capture dust clumping. Moreover, the VSI leads to vortex formation in 3D
(Manger et al., 2020), and hence dust clumping in dust-laden vortices (e.g.,
Raettig et al., 2021), which represents another alternative for planetesimal
formation in turbulent disks. It is worth noting that the VSI and the MRI can
co-exist (Cui & Bai 2022, submitted), and turbulence is dominated by the
MRI/VSI when the ambipolar Elsasser number $Am$ is large/small compared to
unity. It is likely that clumping conditions still depend on the nature of
turbulence (be it due to the VSI or the MRI), although more investigations are
needed to understand the clumping mechanisms, and how the dust dynamics is
affected when both instabilities are present.

###  4.4 Observational implications

#### Dust clumping and fine-scale substructures

There can be multiple observational diagnostics to probe the turbulent dust
rings. First, as such dust rings are clumpy, there can be both optically thick
and optically thin regions within the ring. Observations at sub-mm by ALMA can
hardly resolve such regions, but the clumpiness could be reflected in the
spectral indices, as suggested by Scardoni et al. (2021). Further analysis of
our results with radiative transfer and comparison with ALMA observations may
provide evidence of dust clumping in the observed dust rings.

Second, we have seen the formation of secondary rings/filaments in our weak-
bump simulations. These fine-scale substructures are separated on scales
$\lesssim H$, but are all located within a single gas pressure bump. We note
that generally we do not expect ring-like substructures in the gas component
on scales $\lesssim H$ (Dullemond et al., 2018), whereas observationally,
evidence of such fine-scale ring-like substructures is mounting Jennings et
al. (2022), with hints of multiple rings within previously identified rings.
Our results provide a natural explanation of forming multiple sub-$H$ scaled
substructures within a single gas bump. For instance, they may provide an
alternative explanation to the compact double rings observed in HD 169142 at
57 and 64 AU separately, previously interpreted as the consequence of a
migrating protoplanet (Pérez et al., 2019).

Third, we note that the dust rings from our simulations are not entirely
axisymmetric, with fine-scale substructures along the azimuth. Although our
simulations are local and do not cover the entire azimuth, we notice that
fine-scale, low-contrast asymmetries have been observed in some systems (e.g.,
DM Tau, Hashimoto et al., 2021) and may be common among gapped protoplanetary
disks (van der Marel et al., 2021).

#### Disk turbulence

Dust ring width has recently been used to constrain the level of turbulence
(the $\alpha$ parameter) in disks (Dullemond et al., 2018; Rosotti et al.,
2020). However, as discussed in §4.2, disk turbulence in the presence of a
dust ring can be both anisotropic and non-uniform due to dust feedback.

More specifically, turbulent diffusion in the radial direction in our
simulations is overall stronger than that in the vertical direction, leading
to a large aspect ratio of $w_{d}/H_{d}\sim 10$ for the single dust ring in
run Z2H, which is 5 times of aspect ratio of the gas bump inserted in the
simulation. This effect is even more extreme for the weak bump case, with
$w_{d}/H_{d}\sim 30$ for Z2w, due to the formation of multiple dust filaments.
This flat morphology of dust ring is consistent with recently reported highly
settled disk with relatively extended rings in Oph 163131 (Villenave et al.,
2022), likely suggesting different turbulent levels in radial and vertical
directions.

In addition, dust feedback weakens turbulent diffusion at the dust ring center
(see §4.2). It makes the dust ring narrower in our strong bump case, but the
presence of additional off-centered filaments increases the overall $w_{d}$
(if they are unresolved), which largely cancels out the effect of reduced ring
width. Thus, there is a degeneracy between turbulence level, dust feedback,
and bump amplitude in determining the dust ring width. A single $\alpha$
parameter is not necessarily a good indicator of disk turbulence level, and
special care is needed when interpreting the level of disk turbulence from
observations of dust rings.

#### Prospects for future observations

Overall, our simulations suggest that the dynamics of dust-trapping rings can
be constrained by the presence or absence of fine-scaled substructures.
Simultaneously constraining the gas bump and characterizing the fine-scale
dust substructures are crucial for testify such expectations. High resolution
observations searching for kinematic features (e.g. ExoALMA
555https://www.exoalma.com) will provide essential information on the
properties of the gas pressure bump, while identifying sub-$H$ scaled
substructures requires observations with the longest baseline by ALMA for
nearby disks, and also potentially by future upgrade of ALMA (e.g., Burrill et
al., 2022), ngVLA (e.g., Andrews et al., 2018b) and/or SKA (Ilee et al.,
2020).

##  5 Conclusions and Perspectives

In this paper, we conduct the first controlled experiments of dust dynamics
and clumping in the MRI-turbulent pressure bumps. Our simulations include
realistic level of the MRI turbulence applicable to the ambipolar-diffusion
dominated outer disks, and a carefully-designed forcing scheme to achieve
self-consistent balance between dust trapping and turbulent diffusion. Our
main findings include:

  * •

At solid abundance compatible to solar, we find robust and efficient dust
clumping in radial pressure maxima, with a sizable fraction of dust mass
residing in clumps over the Roche density. Secondary dust filaments may form
outside of the pressure maxima, but shows only weak or no clumping.

  * •

Dust feedback substantially affects ring properties. Strong mass loading in
dust rings reduces turbulent diffusion, making the ring narrower with stronger
dust settling. Weaker pressure bump can result in formation of meta-stable
secondary dust filaments near the bump.

Our work strongly support that dust-trapping rings are robust sites for
planetesimal formation. Our simulations do not include self-gravity, and it is
yet to examine how efficient the dust in dense clumps get converted into
planetesimals. This could result from a competition between self-gravity and
turbulent diffusion (Gerbig et al., 2020), though the internal turbulent
diffusivity is uncertain (but see Klahr & Schreiber (2020, 2021)) and requires
further study.

Our results also have important observational implications. In particular,
depending on the strength of the gas bump, as well as dust abundance, dust
feedback leads to a diversity of dust ring morphologies that may evolve over
time. Conventional estimate of the level of turbulence based on ring width and
dust thickness does not necessarily reflect true level of disk turbulence.
Moreover, our results suggest the presence of fine-scale substructures that
may be resolvable at the highest resolution by ALMA and the future ngVLA.

As a first study, we consider only a single dust species with fixed stopping
time $\tau_{s}=0.1$. The natural next step is to incorporate a dust size
distribution. This will allow us to further examine how particles of different
sizes participate in dust clumping, and to yield more realistic observational
signatures. Moreover, the pressure bumps in our shearing-box simulations rely
on an artificial forcing prescription. Future work should consider more
realistic bump formation scenarios such as magnetic flux concentration and
planet-induced gas bumps in the global context.

## Acknowledgements

We thank Daniel Carrera, Kees Dullemond, Guillaume Laibe and Chao-Chin Yang
for useful discussions, and Greg Herczeg for overseeing the completion of this
work. This work is supported by the National Key R&D Program of China No.
2019YFA0405100, and the China Manned Space Project with NO. CMS-CSST-2021-B09.
ZX Acknowledges the support of the ERC CoG project PODCAST No. 864965 and NSFC
project 11773002. Numerical simulations are conducted on the Orion cluster at
Department of Astronomy, Tsinghua University, and on TianHe-1 (A) at National
Supercomputer Center in Tianjin, China.

## References

  * Andrews et al. (2018a) Andrews, S. M., Terrell, M., Tripathi, A., et al. 2018a, ApJ, 865, 157
  * Andrews et al. (2018b) Andrews, S. M., Wilner, D. J., Macías, E., Carrasco-González, C., & Isella, A. 2018b, in Astronomical Society of the Pacific Conference Series, Vol. 517, Science with a Next Generation Very Large Array, ed. E. Murphy, 137
  * Auffinger & Laibe (2018) Auffinger, J., & Laibe, G. 2018, MNRAS, 473, 796
  * Bai & Stone (2010a) Bai, X.-N., & Stone, J. M. 2010a, ApJ, 722, 1437
  * Bai & Stone (2010b) —. 2010b, ApJ, 722, L220
  * Burrill et al. (2022) Burrill, B. P., Ricci, L., Harter, S. K., Zhang, S., & Zhu, Z. 2022, ApJ, 928, 40
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Carrera & Simon (2022) Carrera, D., & Simon, J. B. 2022, arXiv e-prints, arXiv:2204.14270
  * Carrera et al. (2021) Carrera, D., Simon, J. B., Li, R., Kretke, K. A., & Klahr, H. 2021, AJ, 161, 96
  * Carrera et al. (2022) Carrera, D., Thomas, A. J., Simon, J. B., et al. 2022, ApJ, 927, 52
  * Chen & Lin (2020) Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132
  * 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, 507, 1106
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46
  * 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
  * Gerbig et al. (2020) Gerbig, K., Murray-Clay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91
  * Gole et al. (2020) Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132
  * Hashimoto et al. (2021) Hashimoto, J., Muto, T., Dong, R., et al. 2021, ApJ, 911, 5
  * Heinemann & Papaloizou (2009) Heinemann, T., & Papaloizou, J. C. B. 2009, MNRAS, 397, 64
  * Huang et al. (2020) Huang, P., Li, H., Isella, A., et al. 2020, ApJ, 893, 89
  * Ilee et al. (2020) Ilee, J. D., Hall, C., Walsh, C., et al. 2020, MNRAS, 498, 5116
  * Jennings et al. (2022) Jennings, J., Booth, R. A., Tazzari, M., Clarke, C. J., & Rosotti, G. P. 2022, MNRAS, 509, 2780
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75
  * Klahr & Schreiber (2020) Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54
  * Klahr & Schreiber (2021) —. 2021, ApJ, 911, 9
  * Lehmann & Lin (2022) Lehmann, M., & Lin, M. K. 2022, A&A, 658, A156
  * Li & Youdin (2021) Li, R., & Youdin, A. 2021, arXiv e-prints, arXiv:2105.06042
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Ono et al. (2016) Ono, T., Muto, T., Takeuchi, T., & Nomura, H. 2016, ApJ, 823, 84
  * Paardekooper et al. (2022) Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, arXiv e-prints, arXiv:2203.09595
  * Pérez et al. (2019) Pérez, S., Casassus, S., Baruteau, C., et al. 2019, AJ, 158, 15
  * 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
  * Riols et al. (2020) Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95
  * Rosotti et al. (2020) Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173
  * Scardoni et al. (2021) Scardoni, C. E., Booth, R. A., & Clarke, C. J. 2021, MNRAS, 504, 1495
  * Schäfer et al. (2020) Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  * Suriano et al. (2018) Suriano, S. S., Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2018, MNRAS, 477, 1239
  * Teague et al. (2018) Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133
  * Umurhan et al. (2020) Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4
  * van der Marel et al. (2021) van der Marel, N., Birnstiel, T., Garufi, A., et al. 2021, AJ, 161, 33
  * Villenave et al. (2022) Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Xu & Bai (2022) Xu, Z., & Bai, X.-N. 2022, ApJ, 924, 3
  * 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
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
  * Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81

##  Appendix A Dust clump properties

In order to further analyze the dust properties in dust clumps, we present in
Figure A the cumulative distribution function (CDF) as a function of dust
density for runs Z2H and Z2w. Local dust density $\rho_{d}$ is evenly divided
into 100 bins in logarithmic space between $0.01\rho_{0}$ and 5000 $\rho_{0}$,
and the CDF is obtained by the measuring the probability of a dust particle
residing in a region with local dust density higher than the value in each
bin.

Figure A:  Cumulative distribution function (CDF) of dust density for runs Z2H
and Z2w, as labeled in the titles in each panel. Upper panels show CDF for the
whole simulation domain, whereas the bottom two panels measure the CDFs of run
Z2w for dust particles residing in $x<-0.5H$ and $x>-0.5H$ regions separately,
in order to investigate clump properties in the primary and secondary rings.
The CDFs are measured by averaging within different time intervals in units of
$\Omega_{K}^{-1}$ as marked in different colors (see the legend). The black
curves in each panel shows corresponding runs with no feedback, i.e., Z0 and
Z0w, averaged within $T=720\Omega_{K}^{-1}$ to $1520\Omega_{K}^{-1}$.

Dust density is substantially boosted for both the strong bump and weak bump
cases when feedback is included, compared to the black curves showing the
corresponding no feedback runs (Z0 and Z0w). The CDF is clearly time-
dependent, but after about a hundred local orbits, a significant fraction of
dust particles ($\gtrsim 30\%$) are found to reside in regions that are
$\gtrsim 200$ times background gas density, reaching Roche density. Although
self-gravity is not included in our simulations, such high fraction implies
that planetesimal formation can be highly efficient.

For the weak bump case, we further divide the domain in half at $x=-0.5$ to
separate the primary and secondary rings. We see that strong clumping occurs
mostly in the primary dust ring at the bump center ($x>-0.5H$, see the lower-
right panel), with the CDF very similar to the strong bump case. In the
secondary ring ($x<-0.5H$, lower-left panel), the overall dust density is only
boosted by a small factor $\lesssim 5$ with very clumping. There are only a
couple dust clumps present at later time that appear as outliers in the CDF,
which represent a tiny fraction ($\lesssim 3\%$) of particles. The nature of
these clumps and their potential to form planetesimals deserve future studies.

##  Appendix B Dust ring evolution in long term

In order to examine the lifetime and long-term evolution of dust rings in our
simulation, especially the secondary rings/filaments in the weak bump case, we
continue to run Z1w and Z2w all the way to $T=2880\Omega_{K}^{-1}$
($T_{sim}=3000\Omega_{K}^{-1}$), and show the evolution of dust ring widths
and dust ring profiles in Figure B.

For run Z1w, shown in the upper panel of Figure B, the overall ring width
decreases after $T\sim 1500\Omega_{K}^{-1}$. This can also be seen clearly in
the lower panel of Figure B, showing the evolution of dust radial
distribution, where the secondary ring spreads and later merge into the
primary ring. On the other hand, the secondary ring in run Z2w is more long-
lived, persisting over hundreds of orbits. This is likely due to the higher
inertia in the dust ring which withstands potential disruption by turbulent
diffusion (see also Appendix C) and/or radial drift. It is also in line with
our results in paper I, where dust feedback enhances the MRI zonal flows with
more extended lifetimes. Moreover, we see that additional secondary
rings/filaments may form at off-center positions, and even the primary ring
itself made a small shift in its position over the duration of the long-term
simulation (while still residing in the bump center). The results indicate
that the presence of secondary rings/filaments is likely generic, which
dynamically evolve over the timescale of tens to a few hundreds orbits.

Figure B:  Upper panel: time evolution of dust ring width $w_{d}$ as a
function of time $T=720\Omega_{K}^{-1}$ to $2880\Omega_{K}^{-1}$, featuring
runs Z1w and Z2w at later stages. Lower two panels: time evolution of
azimuthally and vertically averaged radial profiles of dust density at the
midplane ($\pm 1$ cell from $z=0$) for runs Z1w and Z2w, same time interval as
the upper panel.

##  Appendix C Turbulent diffusion in dust rings

Figure C.1 and Figure C.2 show the azimuthal and time averaged gas rms
velocities in the radial and vertical directions, i.e. $\sqrt{\langle
v_{x}^{2}\rangle}$ and $\sqrt{\langle v_{z}^{2}\rangle}$, for the strong and
weak bump cases respectively. They serve as a proxy for estimating the
strength of turbulent diffusion. For run Z0 without dust feedback, gas rms
velocities in both directions are approximately uniform, with $\sqrt{\langle
v_{x}^{2}\rangle}\sim 0.07c_{s}$ and $\sqrt{\langle v_{z}^{2}\rangle}\sim
0.016c_{s}$ (the results from run Z0w are largely the same). With
$w_{d}/w_{\rm bump}=0.28$ and $H_{d}/H=0.067$, our results are overall
consistent with the expectations of $w_{d}/w_{\rm bump}\sim\sqrt{\langle
v_{x}^{2}\rangle/(\Omega_{K}\tau_{s})}$ and $H_{d}/H\sim\sqrt{\langle
v_{z}^{2}\rangle/(\Omega_{K}\tau_{s})}$, modulo correction factors from the
eddy time (Youdin & Lithwick, 2007; Zhu et al., 2015, paper I). They are also
generally consistent with the turbulence parameter $\alpha\sim 10^{-3}$
measured in our simulation.666The MRI turbulence level in our simulations can
be characterized by $\alpha$ parameter (Shakura & Sunyaev, 1973), and is
obtained by summing over the Maxwell and Reynold stresses normalized by
thermal pressure. Our simulation without dust feedback (Z0) gives time-
averaged $\alpha=1.08\times 10^{-3}$ , consistent with previous studies (Zhu
et al., 2015; Xu et al., 2017, paper I).

With feedback included in run Z2H, we find that the radial gas rms velocity is
largely reduced in the dust ring, which is clearly visible in Figure C.1. We
obtain $\sqrt{\langle v_{x}^{2}\rangle}\sim 0.07c_{s}$ by averaging within the
whole domain, which is consistent to the value measured in Z0. In the the dust
ring, however, we find $\sqrt{\langle v_{x,{\rm ring}}^{2}\rangle}\sim
0.02c_{s}$, about 3 times lower than that outside of the ring. This reduction
in $\sqrt{\langle v_{x}^{2}\rangle}$ generally matches the reduction in ring
width reported earlier. In the vertical direction, we note that in paper I, we
found that dust feedback leads to a strong reduction in the turbulent eddy
time which promotes settling. In this work, we find that in addition, with
more enhanced dust mass loading, the vertical rms velocity $\sqrt{\langle
v_{z}^{2}\rangle}$ is also reduced. Although the reduction is less dramatic
than that in the radial direction, it is sufficient to account for the
additional settling found in Figure 2. The reduction in both $\sqrt{\langle
v_{x}^{2}\rangle}$ and $\sqrt{\langle v_{z}^{2}\rangle}$ is also seen in run
Z2w (see right panels of Figure C.2).

Figure C.1: Azimuthally averaged spatial distribution of dust density (upper
panels in each plot), gas rms velocity on radial direction (middles panels)
and vertical direction (lower panels), time-averaged within
$T=800-1500\Omega_{K}^{-1}$, for runs Z0 (left panels) and Z2H (right panels).
Note that this figure is not to the real aspect ratio.

Figure C.2: Same as Figure C.1, but for runs Z0w and Z2w.


# Increasing Supply Chain Resiliency Through Equilibrium Pricing and
Stipulating Transportation Quota Regulation

Mostafa Pazoki Edwards School of Business, University of Saskatchewan,
Saskatoon, SK, Canada. Email<EMAIL_ADDRESS>Hamed Samarghandi
_Corresponding Author_ ; Edwards School of Business, University of
Saskatchewan, Saskatoon, SK, Canada. Email<EMAIL_ADDRESS>Mehdi
Behroozi Department of Mechanical and Industrial Engineering, Northeastern
University, Boston, MA, USA. Email<EMAIL_ADDRESS>
###### Abstract

Supply chain disruption can occur for a variety of reasons, including natural
disasters or market dynamics. If the disruption is profound and with dire
consequences for the economy, the regulators may decide to intervene to
minimize the impact for the betterment of the society. This paper investigates
the minimum quota regulation on transportation amounts, stipulated by the
government in a market where transportation capacity is below total production
and profitability levels differ significantly among different products. In
North America, an interesting example can happen in rail transportation
market, where the rail capacity is used for a variety of products and
commodities such as oil and grains. This research assumes that there is a
shipping company with limited capacity which will ship a group of products
with heterogeneous transportation and production costs and prices.
Mathematical problems for the market players as well as the government are
presented, solutions are proposed, and implemented in a framed Canadian case
study. Subsequently, the conditions that justify government intervention are
identified, and an algorithm to obtain the optimum minimum quota is presented.

Keywords: Rail Transportation Pricing Optimization; Game Theory; Minimum Quota
Regulation; Resilient Commodities Supply Chain

##  1 Introduction

Many societal requirements and special circumstances may necessitate
government’s limited intervention in free markets. The most recent and major
example of such benevolent interventions is the unprecedented policies of the
governments around the world in the wake of the COVID-19 Pandemic. The field
of transportation and supply chain can be heavily disrupted by crises of this
type and can benefit from policymakers’ intervention. Supply chain disruptions
can occur due to several other forms of events such as a natural disaster in a
major port, an obstruction in the Suez Canal, chronic shortages in products
that are either not economical to produce or their technology is not
available, strikes impacting major suppliers of certain products, and
sanctions or blockades against a major producer country of certain products or
commodities. Such disruptions could happen in the supply chain of almost all
products regardless of how basic or advanced they might be.

It is also possible for these disruptions to happen due to the free markets
internal dynamics. These may occur when a purely profit-driven supply chain
values certain products more than others while the disfavored products,
although non-profitable for the supply chain and transportation companies, are
vital for the society. Furthermore, time sensitivity of such situations or
perishability of the products may not allow us to wait for the market to
regulate itself. Imposing a minimum shipment quota for certain products may be
necessary for the betterment of society or resiliency reasons. Two interesting
cases happened in Canada’s (CanadaCase) and North Dakota’s (USNYT; USNBCNews)
rail transportation systems. Here, we focus on the Canadian case.

In 2013, the Canadian railroad capacity was strained. Grain harvested, which
needed to be shipped to the consumption markets by rail, reached an all-time
high. Soaring global oil prices and limited pipeline capacity caused the
producers to consider shipping oil to the international markets by rail.
Meanwhile, an abnormally cold winter stalled the trains. A carry-over of 18
million tonnes of crops was expected in 2014, six million tonnes higher than
the average of 12 million. Moreover, oil producers’ profit margins were high;
hence, they were able to allocate most of the rail transportation capacity to
oil and oil-derivatives. This situation urged the Federal Government to ensure
that carry-over volumes of grain will be shipped to the markets.

In May 2014, after a series of consultation with various interest groups and
stakeholders including farmers, oil producers, various other commodity
producers such as chemical and lumber producers, transportation companies,
public representatives, etc., the Federal Government passed “Canada’s Fair
Rail for Grain Farmers Act” (Canada2014Act) to reduce backlogs at rail
terminals and grain elevators, and to ensure that farmers can transport their
grain to the markets during the peak price season. The act states that
Canadian National Railway111https://www.cn.ca/en (CN) and Canadian Pacific
Railway222https://www.cpr.ca/en (CP) must each ship 500,000 tonnes of grain
per week. Motivated by the described event, this article examines the minimum
quota regulation stipulated on transportation companies, where a group of
suppliers are competing over the transportation capacity.

This paper assumes that there are several suppliers producing distinct
products and commodities. These commodities must be shipped from the point of
production to the point of sale or consumption using the same mode of
transportation, which, in our case, is rail. Since the number of available
trains and shipping containers are limited, the suppliers compete with each
other to gain access to transportation. If the profits for one of the products
is much higher than the others, the producers will offer dominating
transportation prices to the transportation provider, essentially cutting the
producers of the other products off the market.

This was the case in 2013, when oil producers offered much higher prices for
rail transportation compared to grain producers, which left the farmers in a
weak position. The grain backlog and unavailability of transportation costed
farmers between two and four billion dollars. In case one or more of these
left-off products are of strategic importance, government will be willing to
intervene in favor of the suppliers if they cannot suggest a competitive
transportation price to the shipping firms. As a result of the described
situation, the Federal Government mandated a quota (5,500 carloads per week at
the time) for grain transportation and a penalty for non-compliance through “
_Canada’s Fair Rail for Farmers Act_.”

Once the regulation was in place, all the stakeholders expressed
disappointment. According to grain producers, neither the assigned quota was
enough to remove the backlog, nor the penalty was enough to convince the
transporters to meet the quota. The farmers believed that the oil producers
were able to cover the penalty for the transportation companies. On the other
hand, the oil producers accused the government of intervening in the free
market of transportation and tilting the ground against the profitable
commodity.

Transportation capacity, either defined for people’s movement in a
transportation network or moving commodities, depends almost completely on
node capacity and route capacity. In general, transportation capacity is the
_maximum throughput_ of a system performed within the permitted constraints
(Kasikitwiwat2005). Therefore, for the rail transportation, capacity could
mean the minimum of total car capacity (appropriate for handling a specific
product), the rail road capacity, and temporary stocking capacity in
warehouses, which also includes loading/unloading capacity. In this paper,
rail transportation capacity is measured by the number of carloads of products
shipped during a certain amount of time, which In North America, is reported
by the transporters on a weekly basis, and is available from various
sources333For example, see the weekly reports of the Association of American
Railroads (AAR), accessed on September 5, 2022 at https://www.aar.org/data-
center/rail-traffic-data. When the demand for shipping a product is more than
the addressed transportation capacity for a limited amount of time, the
suppliers need to stock the excess products and pay the required expenses.
However, if the total supply and demand for the products are constantly more
than transportation capacity, as was the case for the described Canadian
episode, outside (non-market) intervention is required to guarantee some
transportation share for the suppliers of all the various products and
commodities requiring transportation. Having this in mind, we seek to address
the following research questions regarding the public policy and
transportation regulation:

######  Q 1.

The government wants to know when it should intervene in the market in favor
of the producers of one of the commodities.


# Limit spectral measures of matrix distributions of metric triples

A. M. Vershik , F. V. Petrov St Petersburg Department of V. A. Steklov
Mathematical Institute of Russian Academy of Sciences, St Petersburg State
University, Institute for Information Transmission Problems of Russian Academy
of Sciences (Kharkevich Institute); supported by RSF grant project no.
21-11-00152St Petersburg Department of V. A. Steklov Mathematical Institute of
Russian Academy of Sciences, St Petersburg State University; supported by RSF
grant project no. 21-11-00152

###### Abstract

A notion of the limit spectral measure of a metric triple (i.e., a metric
measure space) is defined. If the metric is square integrable, then the limit
spectral measure is deterministic and coinsides with the spectrum of the
integral operator in $L^{2}(\mu)$ with kernel $\rho$. We construct an example
in which there is no deterministic spectral measure.

##  1 Definitions

Consider a _metric triple_ $\tau=(X,\mu,\rho)$ (here triple is for space-
measure-metric), where $(X,\mu)$ is standard continuous probability measure
space (Lebesgue – Rokhlin space), and $(X,\rho)$ is a complete separable
metric (Polish) space, and the measure $\mu$ is a Borel and non-degenerate
(=with a full support) measure. For the metric triple $\tau=(X,\mu,\rho)$
consider the so-called _matrix distribution_ [2], i.e. a measure on infinite
matrices of distances

$M_{\tau}=(\rho(x_{i},x_{j}))_{i,j=1}^{\infty},$

where $\\{x_{i}\\}_{i=1,\dots}$ is an infinite sequence of independent random
variables distributed by the measure $\mu$. This measure has important
properties and, first of all, it is invariant and ergodic with respect to the
infinite symmetric group acting on matrices by simultaneous substitutions of
rows and columns. Before discussing other properties of matrix distributions,
let us formulate a theorem on the classification of $mm$-spaces.

The exact terminology has not yet been worked out: Gromov’s notion [3]
“$mm$-space” — that is, a space with a measure and a metric, has not been
defined quite precisely. The notion of an “admissible triple” $(X,\mu,\rho)$
which we introduced earlier in [11] has exactly defined meaning, but the name
does not seem expressive enough. We propose the term ”metric triple”; note
that the word “metric” applies to both measures and metrics (see also the
survey [12]).

An isomorphism between two metric triples $\tau_{1}(X_{1},\mu_{1},\rho_{1})$
and $\tau_{2}=(X_{2},\mu_{2},\rho_{2})$ is an invertible measure-preserving
isometry:

$T:X_{1}\to X_{2},T_{*}(\mu_{1})=\mu_{2},\rho_{2}(Tx,Ty)=\rho_{1}(x,y).$

Denote the matrix distributions of these triples by $M_{\tau_{i}},i=1,2$.

The measures $\mu_{1},\mu_{2}$, if not stated otherwise, are assumed to be
continuous and non-degenerated (that is, positive on all open sets).

######  Theorem 1 (Gromov– Vershik).

Two metric triples $\tau_{1}$ and $\tau_{2}$ are isomorphic if and only if
their matrix distributions, as measures on matrices, coincide:
$M_{\tau_{1}}=M_{\tau_{2}}$.

Concerning this theorem, see [3, 2], and the survey [12].

Thus, every metric triple by this theorem is uniquely determib ned by its
matrix distribution upto isomorphism.

For each metric triple $\tau=(X,\mu,\rho)$, we may try to associate the
integral operator $I_{\tau}$ for which the metric is the kernel:

$(I_{\tau}f)(y)=\int\rho(x,y)f(x)d\mu(x).$

If the metric $\rho$ is square summable with repect to the measure $\mu$:
$\int_{X\times X}\rho(x,y)^{2}d\mu(x)d\mu(y)<\infty,$ then this operator is
well defined in the space $L^{2}_{\mu}(X)$, it is self-adjoint and is a
Hilbert – Schmidt operator. Thus, it is compact and has a purely point
spectrum, which is a square summable sequence of real numbers. But if the
metric is only summable, then the natural domain of $I_{\tau}$ is rather
intricate, and dealing with spectrum requires caution.

For example, under the assumption of integrability, in particular in the case
of compactness or finiteness of the diameter of the metric space, this
operator is known to be defined in the space of measurable bounded functions
$L^{\infty}$. It seems to be difficult to express in terms of metric the
possibility of considering this operator as a (at least unbounded closed)
operator in $L^{2}_{\mu}(X)$ and say something about its spectrum.

But our approach of studying metric triples concerns not this potentially
definable operator, but matrix distributions of the space.

For a matrix $M_{n}$, an upper minor of order $n$ of a random matrix
$M_{\tau}$, consider the corresponding spectral measure on the real line,
which we think about as the sum of unit delta-measures in the eigenvalues of
$M_{n}$ (multiplicity respected). On teh set $\rm{Meas}$ of measures on the
line consider the topology defined by the functionals $\mu\to\int fd\mu$,
where $f$ is a continuous function on the line with a bounded support not
containing 0.

######  Definition 1.

The spectral measure of the metric triple is the limit in the space
$\rm{Meas}({\mathbb{R}})$, if it exists, of a sequence of measures of the form
$d\mu_{n}(x\cdot a_{n})$, where $a_{n}$ is some normalizing sequence of
positive numbers.

Our goal is to find out the conditions on the metric that guarantee the
existence of the spectral measure , its computation and, most importantly, to
investigate whichinformation about the properties of the triple does spectral
measure retain. At present these questions have not been studied yet, if we do
not have in mind the case of finiteness of the integral of the square of the
metric. As far as we know, even the answer to the question posed earlier by
the first author as to whether the $mm$-space is determined by a spectral
measure, is open. Most likely, the general answer is negative, but it is not
clear whether it happens to be positive in any sufficiently rich class of
metrics. It is interesting to relate metric invariants (e.g., entropy) and
spectral invariants of metric triples.

In numerous works the problem similar to the above one has been studied, but
not for metrics and $mm$-spaces, but for symmetric kernels of general form.
The main consequence for our case is that for a metric integrable with the
square the spectral measure (to be normalized by $a_{n}\equiv n$) coincides
with the purely point spectrum of the operator $I_{\tau}$ (see above), that
is, the sum of the delta-measures in its eigenvalues, cf. [1][Theorem 3.1].

##  2 Main question and example

Our main example shows that there may be no deterministic limit spectral
measure. The authors are not aware of similar examples in the theory of
U-statistics, but we also consider them, apparently for the first time, in
application to matrix distributions.

The main question is: _for which metric triples the limit spectral measure (as
a measure on $\rm{Meas}$) is a non-degenerate measure?_

We give a simple example of the situation, not calculating the spectral
measure completely, but showing only nondeterminacy of some functional on the
limit measure.

Namely, we consider the behavior of the trace of the square of the minor of
oder $n$ and show that it already has a nondegenerate distribution in the
limit.

This distribution is related to the behavior of $U$-statistics of the form

$\sum_{i<j}h(x_{i},x_{j}),$

where $x_{i}$ are i.i.d., but the first moment of a non-negative random
variable $h(x_{1},x_{2})$ is infinite. We have been unable to find such a
general statement in the extensive literature on $U$-statistics (see [9, 5] et
al.)

######  Theorem 2.

Consider as a metric triple a real line with the usual Euclidean metric and as
a measure — the standard Cauchy distribution $\pi^{-1}\frac{dx}{1+x^{2}}$.
Then the normlaized by $n^{3}$ traces of the squares of the minors of the
matrix distribution converge to some non-degenerate $\alpha$-stable
distribution with $\alpha=1/2$.

###### Proof.

The trace of the square of our random $n\times n$-matrix $M_{n}$ equals to the
sum of squares of its entries, i.e.,

$\frac{1}{2}\rm trM_{n}^{2}=\sum_{1\leqslant i<j\leqslant
n}|x_{i}-x_{j}|^{2}=n\sum x_{i}^{2}-(\sum x_{i})^{2},$

where $(x_{1},\ldots,x_{n})$ is Cauchy i.i.d. sample. The random variable
$\frac{1}{n}\sum x_{i}$ again has the standard Cauchy distribution, thus
$(\sum x_{i})^{2}$ is of order $n^{2}$, in particular, $n^{-3}(\sum
x_{i})^{2}$ converges to 0 in probability. On the other hand, $x_{i}^{2}$ is a
non-negative random variable with tails of order

${\rm prob}(x_{1}^{2}>t)={\rm
prob}(|x_{1}|>\sqrt{t})\sim\frac{1}{\pi\sqrt{t}},t\to\infty.$

It is known [10] that in this case $n^{-2}(x_{1}^{2}+\ldots+x_{n}^{2})$
converges in distribution to the appropriate $1/2$-stable law. Thus so does
$\frac{1}{2n^{3}}\rm{T\\!r}M_{n}^{2}$. ∎

It would be interesting to describe the distribution not only of the sum of
squares of eigenvalues, but also of arbitrary powers and the eigenvalues
themselves.

##  3 Formulation of the general problem and relation to random matrices

In [4], at the suggestion of the first author, preliminary calculations of
spectra of matrix distributions were done for the most natural triples:
$n$-dimensional spheres with Euclidean metric normalized by a Lebesgue
measure. The results of the calculations indicated a significant difference in
the spectral measures already for spheres of small dimensions and not very
large orders of minors. Unfortunately, the authors do not know whether these
experiments were continued. Of course, in these cases due to compactness the
limit spectrum is deterministic and should coincide with the spectrum of the
integral operator $I$, (which was not clear then) and which apparently has
been known for a long time.

The fundamental question is whether there are metric triples (with a non-
integrable metric) whose limit spectrum is similar to the limit spectra of
random matrix theory (i.e. semicircular and similar laws.) To clarify this
difficult question we should first compare the theory of matrix distributions
with the theory of random matrices.

Matrix distributions as measures on the space of infinite matrices are
$S_{\infty}$-invariant and ergodic with respect to the infinite symmetric
group. It makes sense to compare the collection of such measures with a
similar collection of $U_{\infty}$-invariant ergodic measures on of the space
of infinite Hermitian matrices, which correspond to Wishart statistics, cf.
[7, 6, 8]. The common feature of these and other measures is that they are
parameterized by similar sets of frequencies, which differs from measures on
matrix spaces with independent elements (GOE, GUE, white noise.)

## References

  * [1] V. Koltchinskii, E. Giné. Random matrix approximation of spectra of integral operators. _Bernoulli_ 6:1, 113–167, 2000.
  * [2] A. M. Vershik. Random metric spaces and universality. Russian Math. Surveys, 59:2 (2004), 259–295
  * [3] M. Gromov. Metric Structures for Riemannian and Non-Riemannian Spaces. Birkhauser, Boston, 1999.
  * [4] E. Bogomolny, O. Bohigas, C Schmit. Spectral properties of distance matrices. _J. Phys. A_ , 36:12, 3595–3616, 2003.
  * [5] V. Bentkus, F. Götze. Optimal bounds in non-gaussian limit theorems for $U$-statistics. _Ann. Prob._ 27:1, 454–521, 1999.
  * [6] G. Olshanski, A. Vershik. Ergodic unitarily invariant measures on the space of infinite Hermitian matrices, in _Contemporary Mathematical Physics, Amer. Math. Soc. Transl. Ser. 2, vol. 175, Amer. Math. Soc., Providence, RI_ , 137–175, 1996.
  * [7] D. Pickrell. Mackey analysis of infinite classical motion groups. _Pacific J. Math._ 150:1, 139–166, 1991.
  * [8] A. M. Vershik, F. V. Petrov. Central measures of continuous graded graphs: the case of distinct frequencies. _Europ. J. Math._ 8:2, 481–493, 2022.
  * [9] V. S. Koroljuk and Yu. V. Borovskich. Theory of U-statistics. Springer, 1994.
  * [10] I. A. Ibragimov, Yu. V. Linnik. Independent and Stationary Sequences of Random Variables. Wolters-Noordhoff., 1971.
  * [11] A. M. Vershik, P. B. Zatitskiy, F. V. Petrov. Geometry and dynamics of admissible metrics in measure spaces. _Central European Journal of Mathematics_ , 11:3, 379–400, 2013.
  * [12] A. M. Vershik, G. A. Veprev, P. B. Zatitskiy. Dynamics of metrics in the measure spaces and scaling entropy (in print, Russian Math Surveys).