11institutetext: ISAS/JAXA, Sagamihara, Kanagawa, Japan (11email:
<EMAIL_ADDRESS>22institutetext: Earth-Life Science Institute, Tokyo Institute
of Technology, Meguro-ku, Tokyo 152-8550, Japan 33institutetext: Laboratoire
J.-L. Lagrange, Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
F-06304 Nice, France

# A ”no-drift” runaway pile-up of pebbles in protoplanetary disks in which
midplane turbulence increases with radius

Ryuki Hyodo 11 Shigeru Ida 22 Tristan Guillot 33

(DRAFT: )

###### Abstract

Context. A notable challenge of planet formation is to find a path to directly
form planetesimals from small particles.

Aims. We aim to understand how drifting pebbles pile up in a protoplanetary
disk with a nonuniform turbulence structure.

Methods. We consider a disk structure in which the midplane turbulence
viscosity increases with the radius in protoplanetary disks, such as in the
outer region of a dead zone. We perform 1D diffusion-advection simulations of
pebbles that include back-reaction (the inertia) to the radial drift and the
vertical and radial diffusions of pebbles for a given pebble-to-gas mass flux.

Results. We report a new mechanism, the ”no-drift” runaway pile-up, that leads
to a runaway accumulation of pebbles in disks, thus favoring the formation of
planetesimals by streaming and/or gravitational instabilities. This occurs
when pebbles drifting in from the outer disk and entering a dead zone
experience a decrease in vertical turbulence. The scale height of the pebble
subdisk then decreases, and, for small enough values of the turbulence in the
dead zone and high values of the pebble-to-gas flux ratio, the back-reaction
of pebbles on gas leads to a significant decrease in their drift velocity and
thus their progressive accumulation. This occurs when the ratio of the flux of
pebbles to that of the gas is large enough that the effect dominates over any
Kelvin-Helmholtz shear instability. This process is independent of the
existence of a pressure bump.

###### Key Words.:

Planets and satellites: formation, Planet-disk interactions, Accretion,
accretion disks

Figure 1: Schematic illustration of pebble drift and its pile-up within a
protoplanetary disk with a dead zone. The disk gas accretion is characterized
by $\alpha$-parameter $a_{\rm acc}$, while the midplane diffusivity, being a
dead zone ($\alpha_{\rm dead}\ll\alpha_{\rm acc}$) in the inner region, is
characterized by $\alpha_{\rm mid}$. During the inward drift of pebbles, the
pebble scale height $H_{\rm p}$ decreases as $H_{\rm p}\propto\alpha_{\rm
mid}^{1/2}$ until a KH instability prevents it from becoming smaller. A
smaller $H_{\rm p}$ leads to an elevated local midplane concentration of
pebbles within a thinner midplane layer. The elevated midplane pebble-to-gas
ratio causes the back-reaction to be more effective in reducing the radial
drift velocity of pebbles. Such a physical interplay with a sufficiently large
pebble-to-gas mass flux results in a progressive accumulation of pebbles in a
runaway fashion (i.e., the ”no-drift” runaway pile-up).

##  1 Introduction

Forming planetesimals directly from small particles in an evolving
protoplanetary disk is a major challenge in planet formation due to the
”growth barrier” (Blum & Wurm 2000; Zsom et al. 2010) and the ”drift barrier”
(Whipple 1972; Weidenschilling 1977). Streaming instability (SI) may be a
promising mechanism for forming planetesimals directly from pebbles, although
it requires special conditions: the Stokes number $\tau_{\rm s}\gtrsim 0.01$
and a local elevated solid-to-gas ratio ($Z\equiv\rho_{\rm p}/\rho_{\rm
g}\gtrsim 1$, where $\rho_{\rm p}$ and $\rho_{\rm g}$ are the spatial
densities of pebbles and gas, respectively; Youdin & Shu 2004; Youdin &
Goodman 2005; Carrera et al. 2015). How, when, and whether such conditions are
met in an evolving protoplanetary disk are still a matter of intense debate.

Different mechanisms have been proposed for planetesimal formation via the
local enhancement of materials. They include the diffusive redistribution and
recondensation of water vapor outside the snow line (e.g., Stevenson & Lunine
1988; Ciesla & Cuzzi 2006; Ros & Johansen 2013), and detailed numerical
simulations have been performed (e.g., Dra̧żkowska & Alibert 2017;
Schoonenberg & Ormel 2017; Hyodo et al. 2019; Gárate et al. 2020; Hyodo et al.
2021). Other mechanisms include pile-ups at pressure maxima. These pressure
maxima can be created by the planet’s gravity (e.g., Dipierro & Laibe 2017;
Kanagawa et al. 2018), by a sharp change in the local magneto-rotational
instability- (MRI-) driven turbulence structure at an evaporation front (e.g.,
Kretke & Lin 2007; Brauer et al. 2008), by a change in the gas accretion
velocity at the inner (e.g., Chatterjee & Tan 2014; Ueda et al. 2019; Charnoz
et al. 2019) and outer (e.g., Pinilla et al. 2016) edges of a dead zone,
and/or by the change in the gas profile due to gas-pebble friction (back-
reaction) combined with pebble growth (Gonzalez et al. 2017). However, a
positive pressure gradient may not be likely because the non-
magnetohydrodynamic (non-MHD) effects erase a clear transition between dead
and active zones (e.g., Mori et al. 2017). The back-reaction of solids onto
the gas would smooth out the bump (e.g., Taki et al. 2016; Kanagawa et al.
2018). Further studies are required to assess these likelihoods in an evolving
disk.

Over the last few decades, MHD simulations have shown that a region where
ionization is too low for the MRI to operate (i.e., a ”dead zone”; Gammie
(1996)) might ubiquitously exist in the inner part of the disk midplane, and
that only surface layers are magnetically active, supporting accretion (Fig.
1; Gressel et al. 2015; Simon et al. 2015; Bai & Stone 2013; Bai et al. 2016;
Mori et al. 2017). The transition between active and dead zones might be
smooth (i.e., the $\alpha$-parameter in the disk midplane $\alpha_{\rm mid}$
increases as a function of the distance to the star $r$), and $\alpha_{\rm
dead}=10^{-5}-10^{-3}$ is reported for the $\alpha$-parameter within a dead
zone (e.g., Gressel et al. 2015; Simon et al. 2015; Mori et al. 2019). We note
that the sound waves propagated from MRI-active surface layers (e.g., Okuzumi
& Hirose 2011; Yang et al. 2018) and the vertical shear instability (VSI;
e.g., Stoll & Kley 2016; Flock et al. 2020) could induce vertical mixing with
an equivalent $\alpha$ for vertical diffusivity up to $\sim 10^{-3}$, even
within a dead zone. These results imply nonuniform disk turbulence structures.

In this letter, we consider a toy model of a protoplanetary disk with a dead
zone located at the inner region of the disk. We consider different
$\alpha$-parameters for the gas accretion and for the pebble motion. Including
the effects of the back-reaction of pebbles onto gas, which slows the drift
velocity of pebbles, we study how drifting pebbles within a dead zone pile up.
We demonstrate that a runaway pile-up of pebbles occurs for a sufficiently
large pebble-to-gas mass flux $F_{\rm p/g}$ when pebbles reach a critical low-
level of turbulence; we call this the ”no-drift” (ND) runaway pile-up (the
summary of this mechanism is shown by a schematic illustration in Fig. 1). It
is worth mentioning that the newly reported ND mechanism does not require the
snow line or the pressure maxima.

In Section 2, we describe the disk models and settings of our 1D simulations
that solve the diffusion-advection of drifting pebbles. In Section 3, we show
the results of our 1D simulations as well as analytical arguments. In Section
4, we summarize this letter.

##  2 Models and settings

###  2.1 Gas structure

Here, we adopt a toy model where the gas accretion toward the central star is
characterized by a non-dimensional $\alpha$-parameter $\alpha_{\rm acc}$ (Fig.
1). Using the prescription of the classical $\alpha$-accretion disk model
(Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974), the surface density of
the gas is given as

$\Sigma_{\rm g}=\frac{\dot{M}_{\rm g}}{3\pi\nu_{\rm acc}}=\frac{\dot{M}_{\rm
g}}{3\pi\alpha_{\rm acc}c_{\rm s}^{2}\Omega_{\rm K}^{-1}},$ (1)

where $\dot{M}_{\rm g}$ and $\nu_{\rm acc}=\alpha_{\rm acc}c_{\rm
s}^{2}\Omega_{\rm K}^{-1}$ are the gas mass accretion rate and the effective
viscosity, respectively ($c_{\rm s}$ is the gas sound velocity and
$\Omega_{\rm K}$ is the Keplerian orbital frequency).

The gas rotates at sub-Keplerian speed. The degree of the deviation of the gas
rotation frequency from that of Keplerian $\eta$ is given by

$\eta\equiv\frac{\Omega_{\rm K}-\Omega}{\Omega_{\rm
K}}=-\frac{1}{2}\frac{\partial\ln P_{\rm g}}{\partial\ln r}\left(\frac{H_{\rm
g}}{r}\right)^{2}=C_{\rm\eta}\left(\frac{H_{\rm g}}{r}\right)^{2},$ (2)

where $\Omega$, $P_{\rm g}$, and $H_{\rm g}$ are the gas orbital frequency,
the gas pressure, and the gas scale height, respectively, and $C_{\rm\eta}$ is
defined as

$C_{\rm\eta}\equiv-\frac{1}{2}\frac{\partial\ln P_{\rm g}}{\partial\ln r},$
(3)

which depends on the temperature profile (e.g., $C_{\rm\eta}=11/8$ for
$T\propto r^{-1/2}$).

Using the $\alpha$-disk prescription, the gas accretion velocity $v_{\rm g}$
is written as

$\displaystyle v_{\rm g}$ $\displaystyle=-\frac{3\nu_{\rm
acc}}{2r}=-\frac{3\alpha_{\rm acc}H_{\rm g}^{2}\Omega_{\rm
K}}{2r}=-\frac{3\alpha_{\rm acc}}{2}\left(\frac{H_{\rm g}}{r}\right)^{2}v_{\rm
K},$ (4) $\displaystyle\simeq-\frac{3}{2}\alpha_{\rm acc}\eta v_{\rm
K}\left(-\frac{1}{2}\frac{d\ln P_{\rm g}}{d\ln
r}\right)^{-1}=-\frac{3}{2}\alpha_{\rm acc}\eta v_{\rm K}C_{\rm\eta}^{-1},$
(5)

where a negative sign indicates accretion toward the central star.

###  2.2 Pebbles in the disk midplane

The radial drift velocity of pebbles, including the effects of gas-solid
friction $-$ drift back-reaction (hereafter Drift-BKR) $-$ is given as (Ida &
Guillot 2016; Schoonenberg & Ormel 2017; Hyodo et al. 2019)

$\displaystyle v_{\rm p}=-\frac{\Lambda}{1+\Lambda^{2}\tau_{\rm
s}^{2}}\left(2\tau_{\rm s}\Lambda\eta v_{\rm K}-v_{\rm g}\right),$ (6)

where $\tau_{\rm s}$ is the Stokes number of pebbles. $\Lambda\equiv\rho_{\rm
g}/(\rho_{\rm g}+\rho_{\rm p})=1/(1+Z)$ characterizes the strength of the
back-reaction due to the pile-up of pebbles, where $Z\equiv\rho_{\rm
p}/\rho_{\rm g}$ is the midplane pebble-to-gas density ratio ($\rho_{\rm
g}=\Sigma_{\rm g}/\sqrt{2\pi}H_{\rm g}$ and $\rho_{\rm p}=\Sigma_{\rm
p}/\sqrt{2\pi}H_{\rm p}$, where $\Sigma_{\rm p}$ and $H_{\rm p}$ are the
surface density and scale height of pebbles, respectively).

The disk midplane could be too weakly ionized for the MRI to operate (the dead
zone with extremely weak turbulence; Gammie (1996)). Thus, the disk midplane
could have very different turbulence from that of gas accretion, which
controls the disk surface density (i.e., $\alpha_{\rm acc}$). In this letter,
we use an effective viscous parameter $\alpha_{\rm mid}$ to characterize the
radial and vertical diffusion processes of pebbles near the disk midplane
(Fig. 1).

The scale height of pebbles characterizes the degree of pebble concentration
in the disk midplane ($\rho_{\rm p}\propto H^{-1}_{\rm p}$). In the steady
state, the scale height of pebbles is regulated by the vertical turbulent
stirring (Dubrulle et al. 1995; Youdin & Lithwick 2007; Okuzumi et al. 2012;
Hyodo et al. 2019) as

$\displaystyle H_{\rm p,tur}=\left(1+\frac{\tau_{\rm s}}{\alpha_{\rm
mid}\left(1+Z\right)^{-K}}\right)^{-1/2}H_{\rm g},$ (7)

where a coefficient $K$ characterizes the strength of the back-reaction onto
the diffusivity (Hyodo et al. 2019; Ida et al. 2021) and $K=1$ is used (the
choice of the $K$ value does not alter the conclusion – see Appendix A for the
$K=0$ case).

For a small $\alpha_{\rm mid}$ (i.e., for a small $H_{\rm p,tur}$), a vertical
shear Kelvin-Helmholtz (KH) instability prevents a further decrease in the
pebble scale height. This minimum scale height, $H_{\rm p,KH}$, is (Hyodo et
al. 2021)

$\displaystyle H_{\rm p,KH}\simeq
Ri^{1/2}\frac{Z^{1/2}}{\left(1+Z\right)^{3/2}}C_{\rm\eta}\left(\frac{H_{\rm
g}}{r}\right)H_{\rm g}=Ri^{1/2}\frac{Z^{1/2}}{\left(1+Z\right)^{3/2}}\eta r,$
(8)

which is valid for $Z\lesssim 1$, and $Z\gg 1$ indicates a gravitational
collapse. Thus, the scale height of pebbles $H_{\rm p}$ is given as

$\displaystyle H_{\rm p}=\max\left\\{H_{\rm p,tur},H_{\rm p,KH}\right\\}.$ (9)

Figure 2: $\alpha_{\rm mid}/\alpha_{\rm acc}$ (top panels) and $Z=\rho_{\rm
p}/\rho_{\rm g}$ (bottom panels) as a function of the distance to a star in
the 1D numerical simulations. Included are the cases where $\alpha_{\rm
acc}=10^{-2}$ and $\tau_{\rm p}=0.1$. From left to right, the panels show the
cases where $(F_{\rm p/g},\alpha_{\rm dead},r^{*},\Delta r_{\rm tra})$ =
$(0.11,10^{-4},9{\rm\,au},5{\rm\,au})$,
$(0.17,10^{-4},9{\rm\,au},5{\rm\,au})$,
$(0.47,10^{-4},9{\rm\,au},5{\rm\,au})$, and
$(0.17,10^{-4},5{\rm\,au},0.5{\rm\,au})$. The gray lines in the top panels
represent the analytical critical $\alpha_{\rm mid}/\alpha_{\rm acc}$ below
which the ND runaway pile-up is expected to occur for a given $F_{\rm p/g}$
(Eq. (23) with $F_{\rm p/g}>F_{\rm p/g,crit2}$). The gray lines in the bottom
panels show the critical $Z$ (i.e., $Z_{\rm crit}=1$) above which the ND
runaway pile-up occurs (Section 3.2). The blue, green, and red lines in the
bottom panels are those at $t=3\times 10^{4}$ yrs, $t=1\times 10^{5}$ yrs, and
$t=3\times 10^{5}$ yrs from the beginning of the calculations, respectively.
Panel (a) shows the system that reaches a steady-state, while the other cases
(panels (b), (c), and (d)) show the ND runaway pile-ups.

###  2.3 Numerical settings

We performed 1D diffusion-advection simulations that included the back-
reaction (the inertia) to radial drift of pebbles that slows the pebble drift
velocity as pile-up proceeds (see details in Hyodo et al. 2019, 2021). The
drift velocity of pebbles is given by Eq. (6), while the diffusivity is given
as $D_{\rm p}=\alpha_{\rm mid}c_{\rm s}^{2}\Omega_{\rm
K}^{-1}\Lambda^{K}/(1+{\rm\tau_{\rm s}^{2}})$. We included back-reaction onto
the diffusivity of pebbles (with $K=1$). We assumed that pebbles drift outside
the snow line, and the sublimation of pebbles was neglected (i.e., the dead
zone exists beyond the snow line). The Stokes number of pebbles was set to be
constant: $\tau_{\rm s}=0.1$ (Okuzumi et al. 2016; Ida & Guillot 2016). The
disk midplane temperature was set to $T(r)=150{\,\rm K}\times(r/3{\,\rm
au})^{-\beta}$ ($\beta=1/2$). The surface density of the gas (molecular weight
of $\mu_{\rm g}=2.34$) is described by $\Sigma_{\rm g}=\dot{M}_{\rm
g}/3\pi\nu_{\rm g}$, where $\nu_{\rm g}=\alpha_{\rm acc}c_{\rm
s}^{2}\Omega_{\rm K}^{-1}$. We used $\alpha_{\rm acc}=10^{-2}$ and
$\dot{M}_{\rm g}=10^{-8}$ $M_{\odot}$/year. These led to $C_{\rm\eta}=11/8$
for $\Sigma_{\rm g}\propto r^{-1}$ and $T\propto r^{-1/2}$.

We considered a dead zone in the inner region of the disk midplane, and we
used a non-dimensional turbulence parameter in the midplane $\alpha_{\rm
mid}$, which differs from the one that characterizes gas accretion,
$\alpha_{\rm acc}$ (Eq. (1)). The midplane $\alpha_{\rm mid}$ is modeled as

$\alpha_{\rm mid}(r)=\alpha_{\rm dead}+\left(\frac{\alpha_{\rm
acc}-\alpha_{\rm
dead}}{2}\right)\left[\mathrm{erf}\left(3+\frac{6\left(r-r^{*}\right)}{\Delta
r_{\rm tra}}\right)+1\right],$ (10)

where: $\alpha_{\rm acc}$ and $\alpha_{\rm dead}$ are those outside and inside
a dead zone; $r^{*}$ is the innermost radial distance, where $\alpha_{\rm
mid}=\alpha_{\rm acc}$; and $\Delta r_{\rm tra}$ is the radial width of the
transition from $\alpha_{\rm mid}=\alpha_{\rm acc}$ to $\alpha_{\rm
mid}=\alpha_{\rm dead}$. As shown below, the ND mode occurs for an arbitrary
choice of the dead zone structure, that is, irrespective of a sharp or a
smooth change between the active and dead zones (i.e., an arbitrary choice of
$\alpha_{\rm dead}$, $\alpha_{\rm acc}$, $r^{*}$, and $\Delta r_{\rm tra}$),
as long as it satisfies that $\alpha_{\rm mid}$ is smaller than a threshold
value (Section 3.2.2).

At the beginning of the 1D simulations, we set the pebble-to-gas mass flux
$F_{\rm p/g}$ at the outer boundary ($r_{\rm out}=15$ au), and $F_{\rm p/g}$
at $r_{\rm out}$ was fixed throughout the simulations. We note that Elbakyan
et al. (2020) performed a time-dependent simulation of disk formation to find
a high $F_{\rm p/g}$ variation between $\mathcal{O}(10^{-4})$ to
$\mathcal{O}(1)$ in a single disk evolution (see also Ida et al. 2021).

##  3 Results

###  3.1 Numerical results

Figure 2 shows the results of our 1D simulations for different combinations of
$F_{\rm p/g}$ and the midplane turbulence structures (i.e., different choices
of $\alpha_{\rm dead}$, $r^{*}$, and $\Delta r_{\rm tra}$ in Eq. (10)) for the
case of $\alpha_{\rm acc}=10^{-2}$. The top panels show given turbulence
structures (i.e., $\alpha_{\rm mid}(r)/\alpha_{\rm acc}$) and the bottom
panels show the resultant midplane pebble-to-gas ratio in their spatial
densities (i.e., $\rho_{\rm p}/\rho_{\rm g}$). The gray lines in the top
panels are the analytically derived critical $\alpha_{\rm mid}/\alpha_{\rm
acc}$ below which a runaway pile-up of pebbles is expected to occur for a
given $F_{\rm p/g}$ (Eq. (23) with $\tau_{\rm s}=0.1$).

For a small $F_{\rm p/g}$ (e.g., $F_{\rm p/g}=0.11$; panel (a) in Fig. 2), the
system reaches a steady state. For this given $F_{\rm p/g}$, analytical
arguments (the gray line) predict the requirement of a smaller $\alpha_{\rm
mid}/\alpha_{\rm acc}$ for a runaway pile-up of pebbles to occur (i.e., the
black line in Fig. 2 needs to be below the gray line).

For a larger $F_{\rm p/g}$ (e.g., $F_{\rm p/g}=0.17$; panel (b) in Fig. 2)
that has the same $\alpha_{\rm mid}/\alpha_{\rm acc}$ structure as that of
panel (a), pebbles progressively pile up over time. Continuous pile-up occurs
as the drift velocity of pebbles progressively decreases due to strong Drift-
BKR and as the pile-up efficiently proceeds. For the $F_{\rm p/g}$ given here,
the gray line overlaps with the black line (i.e., the analytical requirement
for the runaway pile-up is met) and the analytical arguments (Section 3.2)
show a very good accordance with 1D numerical simulations. For an even larger
$F_{\rm p/g}$ with the same $\alpha_{\rm mid}(r)/\alpha_{\rm acc}$ as those in
panels (a) and (b), a runaway pile-up of pebbles occurs more quickly and
significantly (panel (c) in Fig. 2).

Comparing panels (b) and (d) in Fig. 2 (the same $F_{\rm p/g}$, $\alpha_{\rm
dead}$, and $\alpha_{\rm acc}$ but different shapes and locations of the
transition between active and dead zones), 1D simulations demonstrate that a
runaway pile-up of pebbles occurs irrespective of the shape and location of
the dead zone as long as $F_{\rm p/g}$ is sufficiently large. This shows the
robustness of this new mechanism of runaway pebble pile-up.

Numerical results show that when the change in the midplane turbulence on $r$
is smooth (i.e., $\alpha_{\rm mid}$), the pile-up of pebbles efficiently
propagates toward the outer region in the disk (compare panels (b) and (d) in
Fig. 2). This indicates that, when the ND runaway pile-up occurs at a given
radial distance, a progressive pile-up of pebbles could occur even at a
greater radial distance where the ND conditions are not met, leading to a
global pile-up of pebbles for the SI to operate.

When $\alpha_{\rm acc}=10^{-3}$, 1D simulations only show a runaway pile-up of
pebbles for a very large value of $F_{\rm p/g}$ ($F_{\rm p/g}\gtrsim 0.8$).
This is because the surface density of the gas increases (Eq. (1)) for a
smaller $\alpha_{\rm acc}$, and $\rho_{\rm p}/\rho_{\rm g}$ correspondingly
decreases.

###  3.2 Analytical arguments

Figure 3: Midplane pebble-to-gas ratio $Z\equiv\rho_{\rm p}/\rho_{\rm g}$ in
the $\alpha_{\rm mid}/\alpha_{\rm acc}-F_{\rm p/g}$ space. Included are the
cases where $r=5$ au, $C_{\rm\eta}=11/8$ for $T\propto r^{-1/2}$, $\tau_{\rm
s}=0.1$, and $Ri=0.5$. The left and right panels are the cases where
$\alpha_{\rm acc}=10^{-2}$ and $\alpha_{\rm acc}=10^{-3}$, respectively. The
color contours are obtained by directly solving Eq. (11). The red regions
indicate the no-drift runaway pile-up. The diagonal dashed black line (Eq.
(16)) and the horizontal dashed black line (Eq. (20)) are the analytically
derived critical $F_{\rm p/g,crit}$ above which the ND runaway pile-up takes
place. The critical $Z$ above which the ND runaway pile-up occurs (i.e.,
$Z_{\rm crit}=1$) is shown by the white lines. The parameter map is very
weakly dependent on the radial distance to the star (see Eqs. (16) and (20)),
and the ND mode can occur irrespective of the shape and position of a dead
zone.

####  3.2.1 Midplane concentration of pebbles

Below, we discuss the concentration of pebbles in the disk midplane
considering Drift-BKR. Our arguments are made using the analytical equation
presented below (Eq. (11)). We solve this equation using two distinct
approaches: a direct numerical approach and an analytical approach with some
approximation. In both ways, we identify a parameter space ($F_{\rm
p/g}-\alpha_{\rm acc}-\alpha_{\rm mid}$ space) where the ND runaway pile-up of
drifting pebbles could take place.

The concentration of pebbles at the midplane is written as

$Z\equiv\frac{\rho_{\rm p}}{\rho_{\rm g}}=\frac{\Sigma_{\rm p}}{\Sigma_{\rm
g}}h_{\rm p/g}^{-1}=\frac{v_{\rm g}}{v_{\rm p}}h_{\rm p/g}^{-1}F_{\rm p/g},$
(11)

where $h_{\rm p/g}\equiv H_{\rm p}/H_{\rm g}$ and $v_{\rm p}$ are functions of
$\Lambda(Z)$. We directly find a solution (i.e., solve for $Z$) for Eq. (11)
for a given $F_{\rm p/g}$ and $\alpha_{\rm mid}$. The color contours in Fig. 3
show directly obtained $\rho_{\rm p}/\rho_{\rm g}$ and the red region is where
no steady-state solution is found (i.e., the radial drift of pebbles is halted
due to Drift-BKR). The dashed lines in Fig. 3 are analytically derived
critical $F_{\rm p/g,crit}$, above which no steady-state solution is found for
$\rho_{\rm p}/\rho_{\rm g}$ (i.e., the ND runaway pile-up). The analytical
estimations reproduce the direct solutions well.

The $F_{\rm p/g,crit}$ can be divided into two different regimes. In the first
regime, $H_{\rm p,tur}>H_{\rm p,KH,max}$, where $H_{\rm p,KH,max}$ is the
maximum scale height of pebbles regulated by a KH instability. In this case,
using Eqs. (5) and (6), the vertically averaged metallicity of pebbles is
given as

$\displaystyle Z_{\rm\Sigma}$ $\displaystyle\equiv\frac{\Sigma_{\rm
p}}{\Sigma_{\rm g}}=F_{\rm p/g}\times\frac{v_{\rm g}}{v_{\rm p}},$ (12)
$\displaystyle\simeq F_{\rm p/g}\times\frac{-\frac{3}{2}\alpha_{\rm acc}\eta
v_{\rm K}C_{\rm\eta}^{-1}}{-\left(\frac{\Lambda}{1+\Lambda^{2}\tau_{\rm
s}^{2}}\right)\left(2\Lambda\tau_{\rm s}+\frac{3\alpha_{\rm
acc}}{2C_{\rm\eta}}\right)\eta v_{\rm K}},$ (13)
$\displaystyle\simeq\frac{3\alpha_{\rm acc}F_{\rm p/g}}{4\tau_{\rm
s}C_{\rm\eta}\Lambda^{2}}=\frac{3\alpha_{\rm acc}F_{\rm p/g}}{4\tau_{\rm
s}C_{\rm\eta}}\left(1+h_{\rm p/g}^{-1}Z_{\rm\Sigma}\right)^{2},$ (14)

where an approximation is made for $\alpha_{\rm acc}\ll\tau_{\rm s}\ll 1$,
$Z\ll 1$, $\Lambda\simeq 1$. Solving this equation gives

$Z_{\rm\Sigma}=\frac{\left(1-2ab\right)\pm\sqrt{1-4ab}}{2ab^{2}},$ (15)

where $a=\frac{3\alpha_{\rm acc}F_{\rm p/g}}{4\tau_{\rm s}C_{\rm\eta}}$ and
$b=h_{\rm p/g}^{-1}$. The $Z_{\rm\Sigma}$ has a real solution when $1-4ab>0$
and a critical $F_{\rm p/g,crit1}$ for the first regime is given as

$\displaystyle F_{\rm p/g,crit1}$
$\displaystyle=\frac{C_{\rm\eta}}{3}\frac{\tau_{\rm s}}{\alpha_{\rm
acc}}h_{\rm p/g}\simeq\frac{\left(\alpha_{\rm mid}\tau_{\rm
s}\right)^{1/2}C_{\rm\eta}}{3\alpha_{\rm acc}},$ (16) $\displaystyle\simeq
0.15\times\left(\frac{\alpha_{\rm
acc}}{10^{-2}}\right)^{-1}\left(\frac{\alpha_{\rm
mid}}{10^{-4}}\right)^{1/2}\left(\frac{\tau_{\rm
s}}{0.1}\right)^{1/2}\left(\frac{C_{\rm\eta}}{11/8}\right)$ (17)
$\displaystyle\hskip 110.00017pt\left[{\rm for\,}H_{\rm p,tur}>H_{\rm
p,KH,max}\right],$

where $h_{\rm p/g}\simeq\left(\tau_{\rm s}/\left(\alpha_{\rm
mid}(1+Z)^{-K}\right)\right)^{-1/2}\simeq\left(\tau_{\rm s}/\alpha_{\rm
mid}\right)^{-1/2}$ for $Z\ll 1$ and $\alpha_{\rm mid}\ll 1+\tau_{\rm s}$ (Eq.
(9)). The above criteria for $F_{\rm p/g,crit1}$ are independent of $r$, and
thus the ND mechanism occurs irrespective of the position and shape of a dead
zone as long as it satisfies its criteria (i.e., Eq. (16)). The diagonal
dashed black line in Fig. 3 shows Eq. (16) and is in accordance with the
direct solutions of Eq. (11) (color contours in Fig. 3). The critical $Z$
above which the ND occurs, $Z_{\rm crit}$, is obtained by inserting Eq. (16)
into Eq. (11), and $Z_{\rm crit}\sim 1$ for the first regime. Analytical
$Z_{\rm crit}=1$ agrees with the direct solution (color contours in Fig. 3).
We note that because $F_{\rm p/g,crit1}$ is proportional to $C_{\rm\eta}$, it
increases with a steeper $\Sigma_{\rm g}$ gradient at an outer boundary of a
dead zone. However, the ND mode does not need a sharp boundary. If the
boundary is very gradual, the change in $C_{\rm\eta}$ is negligible.

The second regime appears when a KH instability plays a role, which
corresponds to $H_{\rm p,tur}<H_{\rm p,KH,max}$. Thus, in this case, the scale
height of pebbles is described by $H_{\rm p,KH}$, which is independent of
$\alpha_{\rm mid}$. As $\rho_{\rm p}/\rho_{\rm g}$ increases from zero,
$H_{\rm p,KH}$ initially increases and reaches its maximum $H_{\rm p,KH,max}$
at $Z=\rho_{\rm p}/\rho_{\rm g}=1/2$. As $\rho_{\rm p}/\rho_{\rm g}$ further
increases, $H_{\rm p,KH}$ decreases111We note that $H_{\rm p,KH}$ is a mean
square root of the vertical height of particles but not the truncated height
of the distribution.. By smoothly connecting with the critical value of
$Z_{\rm crit}=1$ in the first regime (i.e., $Z_{\rm crit}=1$ as well in the
second regime, which agrees with the direct solution; see the color contours
in Fig. 3), the critical $F_{\rm p/g}$ in the second regime $F_{\rm
p/g,crit2}$ adopts $H_{\rm p,KH}^{Z=1}$ ($H_{\rm p,KH}$ with $Z=1$), and
$Z_{\rm crit}$ is given as

$\displaystyle Z_{\rm crit}$ $\displaystyle=\frac{\rho_{\rm p}}{\rho_{\rm
g}}\mid_{H_{\rm p,KH}^{Z=1}}=\frac{\Sigma_{\rm p}}{\Sigma_{\rm g}}\left(h_{\rm
p/g,KH}^{Z=1}\right)^{-1},$ (18) $\displaystyle\simeq\frac{3\alpha_{\rm
acc}F_{\rm p/g,crit2}}{4\tau_{\rm s}C_{\rm\eta}(\Lambda_{\rm
KH}^{Z=1})^{2}}\left(h_{\rm p/g,KH}^{Z=1}\right)^{-1},$ (19)

where $h_{\rm p/g,KH}^{Z=1}\equiv H_{\rm p,KH}^{Z=1}/H_{\rm g}$ and an
approximation is made for $\alpha_{\rm acc}\ll\tau_{\rm s}\ll 1$ and
$\Lambda\simeq 1$. The $\Lambda_{\rm KH}^{Z=1}$ is $\Lambda$ using $H_{\rm
p,KH}^{Z=1}$. Here, $Z=1$ and $\Lambda_{\rm KH}^{Z=1}=1/2$. Thus, $F_{\rm
p/g,crit2}$ with $C_{\rm\eta}=11/8$ is given as

$\displaystyle F_{\rm p/g,crit2}$
$\displaystyle=\frac{C_{\rm\eta}}{3}\left(\frac{\tau_{\rm s}}{\alpha_{\rm
acc}}\right)h_{\rm p/g,KH}^{Z=1}=\frac{11}{24}\left(\frac{\tau_{\rm
s}}{\alpha_{\rm acc}}\right)h_{\rm p/g,KH}^{Z=1},$ (20) $\displaystyle\simeq
0.06\times\left(\frac{\alpha_{\rm
acc}}{10^{-2}}\right)^{-1}\left(\frac{\tau_{\rm
s}}{0.1}\right)\left(\frac{Ri}{0.5}\right)^{1/2}\left(\frac{H_{\rm
g}/r}{0.04}\right),$ (21) $\displaystyle\hskip 110.00017pt\left[{\rm
for\,}H_{\rm p,tur}<H_{\rm p,KH,max}\right]$

which is independent of $\alpha_{\rm mid}$. We note that the above criterion
for $F_{\rm p/g,crit2}$ has a very weak dependence on $r$ through $H_{\rm
g}/r$. Equation (20) is plotted by a horizontal dashed black line in Fig. 3,
and it shows a good consistency with the direct solution of Eq. (11).
Including a KH instability (i.e., the second regime and Eq. (20)) prevents the
scale height of pebbles from becoming increasingly smaller as $\alpha_{\rm
mid}$ decreases, reducing the parameter range of the ND runaway pile-up222The
vertical shear streaming instability (VSSI; Lin 2020) might affect the pebble
scale height when the midplane diffusivity is very small. If the VSSI
dominated over a KH instability for the scale height of pebbles, the boundary
of the ND regime would be regulated by the VSSI..

####  3.2.2 Criteria for the ND runaway pile-up

Summarizing the analytical arguments above, the scale height of pebbles
decreases with decreasing $\alpha_{\rm mid}$ until a KH instability plays a
role. As pebbles drift inward from the active to the dead zones in the disk
midplane, $\alpha_{\rm mid}$ becomes smaller with decreasing distance to the
star $r$, while keeping a constant $F_{\rm p/g}$. This implies that the
evolutionary path in the $\alpha_{\rm mid}/\alpha_{\rm acc}-F_{\rm p/g}$ space
(Fig. 3) is a horizontal shift from right to left for a given fixed $F_{\rm
p/g}$ as pebbles drift inward. Therefore, for the ND mode to take place, the
following two conditions need to be met.

First, $F_{\rm p/g}$ needs to be sufficiently large to satisfy $F_{\rm
p/g}\geq F_{\rm p/g,crit2}$. As $F_{\rm p/g,crit2}\propto 1/\alpha_{\rm acc}$
(Eq. (20)), the ND mode is limited to relatively large $F_{\rm p/g}$ values
for the $\alpha_{\rm acc}=10^{-3}$ case ($F_{\rm p/g}\gtrsim 0.8$ at $5$ au),
while only moderate $F_{\rm p/g}\gtrsim 0.08$ is required in the case of
$\alpha_{\rm acc}=10^{-2}$.

Second, a sufficiently small $\alpha_{\rm mid}/\alpha_{\rm acc}$ is required
in order to have a sufficiently small scale height of pebbles. Rewriting Eq.
(16) gives a critical $\alpha_{\rm mid}/\alpha_{\rm acc}$ below which the ND
occurs for a given $F_{\rm p/g}$ as

$\displaystyle\frac{\alpha_{\rm mid,crit}}{\alpha_{\rm acc}}$
$\displaystyle\equiv\left(\frac{3F_{\rm
p/g}}{C_{\rm\eta}}\right)^{2}\alpha_{\rm acc}\tau_{\rm s}^{-1},$ (22)
$\displaystyle\simeq 4.76\times 10^{-3}\left(\frac{F_{\rm
p/g}}{0.1}\right)^{2}\left(\frac{C_{\rm\eta}}{11/8}\right)^{-2}\left(\frac{\alpha_{\rm
acc}}{10^{-2}}\right)\left(\frac{\tau_{\rm s}}{0.1}\right)^{-1},$ (23)

which is independent of $r$. These two criteria are in very good accordance
with the 1D simulations (Section 3.1).

##  4 Conclusions

In this letter, we studied how drifting pebbles pile up in a protoplanetary
disk with a nonuniform turbulence structure, including the back-reaction of
pebbles onto the gas that slows the drift velocity of pebbles as pile-up
proceeds. We considered that gas accretion is regulated by an
$\alpha$-parameter that is distinct from that of midplane turbulence. In the
disk midplane, the turbulence strength can decrease with decreasing radial
distance (e.g., the case where the inner region of the disk is an MRI-inactive
dead-zone and the outer region is MRI-active). Thus, drifting pebbles are
further concentrated in the midplane as the scale height decreases.

We demonstrated a new mechanism for a runaway pile-up of pebbles to occur for
a moderate $F_{\rm p/g}$ when pebbles reach a critical level of low turbulence
in the dead zone via a continuous slowing of their drift velocity due to the
back-reaction (the ”no-drift” (ND) runaway pile-up). The ND runaway pile-up
occurs irrespective of the shape of the dead zone. Our results imply that SI
could also occur as a natural consequence of pebble drift in a protoplanetary
disk with a dead zone.

The minimum scale height of pebbles in a dead zone is a critical parameter
that characterizes the ND mode. Further studies that include detailed physical
processes – such as a KH instability, VSI, vertical shear streaming
instability, and the velocity fluctuation by the sound waves propagated from
MRI-active surface layers – are required.

###### Acknowledgements.

We thank Satoshi Okuzumi and Shoji Mori for helpful comments and fruitful
discussion. R.H. acknowledges the financial support of JSPS Grants-in-Aid
(JP17J01269, 18K13600). R.H. also acknowledges JAXA’s International Top Young
program. S.I. acknowledges the financial support (JSPS Kakenhi 15H02065, MEXT
Kakenhi 18H05438). T.G. was partially supported by a JSPS Long Term Fellowship
at the University of Tokyo.

## References

  * Bai & Stone (2013) Bai, X.-N. & Stone, J. M. 2013, ApJ, 769, 76
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152
  * Blum & Wurm (2000) Blum, J. & Wurm, G. 2000, Icarus, 143, 138
  * Brauer et al. (2008) Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1
  * Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  * Charnoz et al. (2019) Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019, A&A, 627, A50
  * Chatterjee & Tan (2014) Chatterjee, S. & Tan, J. C. 2014, ApJ, 780, 53
  * Ciesla & Cuzzi (2006) Ciesla, F. J. & Cuzzi, J. N. 2006, Icarus, 181, 178
  * Dipierro & Laibe (2017) Dipierro, G. & Laibe, G. 2017, MNRAS, 469, 1932
  * Dra̧żkowska & Alibert (2017) Dra̧żkowska, J. & Alibert, Y. 2017, A&A, 608, A92
  * Dubrulle et al. (1995) Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
  * Elbakyan et al. (2020) Elbakyan, V. G., Johansen, A., Lambrechts, M., Akimkin, V., & Vorobyov, E. I. 2020, A&A, 637, A5
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  * Gárate et al. (2020) Gárate, M., Birnstiel, T., Drażkowska, J., & Stammler, S. M. 2020, A&A, 635, A149
  * Gonzalez et al. (2017) Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984
  * Gressel et al. (2015) Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84
  * Hyodo et al. (2021) Hyodo, R., Guillot, T., Ida, S., Okuzumi, S., & Youdin, A. N. 2021, arXiv e-prints, arXiv:2012.06700
  * Hyodo et al. (2019) Hyodo, R., Ida, S., & Charnoz, S. 2019, A&A, 629, A90
  * Ida & Guillot (2016) Ida, S. & Guillot, T. 2016, A&A, 596, L3
  * Ida et al. (2021) Ida, S., Guillot, T., Hyodo, R., Okuzumi, S., & Youdin, A. N. 2021, arXiv e-prints, arXiv:2011.13164
  * Kanagawa et al. (2018) Kanagawa, K. D., Muto, T., Okuzumi, S., et al. 2018, ApJ, 868, 48
  * Kretke & Lin (2007) Kretke, K. A. & Lin, D. N. C. 2007, ApJ, 664, L55
  * Lin (2020) Lin, M.-K. 2020, arXiv e-prints, arXiv:2011.12300
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
  * Mori et al. (2019) Mori, S., Bai, X.-N., & Okuzumi, S. 2019, ApJ, 872, 98
  * Mori et al. (2017) Mori, S., Muranushi, T., Okuzumi, S., & Inutsuka, S.-i. 2017, ApJ, 849, 86
  * Okuzumi & Hirose (2011) Okuzumi, S. & Hirose, S. 2011, ApJ, 742, 65
  * Okuzumi et al. (2016) Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82
  * Okuzumi et al. (2012) Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106
  * Pinilla et al. (2016) Pinilla, P., Flock, M., Ovelar, M. d. J., & Birnstiel, T. 2016, A&A, 596, A81
  * Ros & Johansen (2013) Ros, K. & Johansen, A. 2013, A&A, 552, A137
  * Schoonenberg & Ormel (2017) Schoonenberg, D. & Ormel, C. W. 2017, A&A, 602, A21
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
  * Simon et al. (2015) Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117
  * Stevenson & Lunine (1988) Stevenson, D. J. & Lunine, J. I. 1988, Icarus, 75, 146
  * Stoll & Kley (2016) Stoll, M. H. R. & Kley, W. 2016, A&A, 594, A57
  * Taki et al. (2016) Taki, T., Fujimoto, M., & Ida, S. 2016, A&A, 591, A86
  * Ueda et al. (2019) Ueda, T., Flock, M., & Okuzumi, S. 2019, ApJ, 871, 10
  * Weidenschilling (1977) Weidenschilling, S. J. 1977, MNRAS, 180, 57
  * Whipple (1972) Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
  * Yang et al. (2018) Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27
  * Youdin & Goodman (2005) Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459
  * Youdin & Lithwick (2007) Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588
  * Youdin & Shu (2004) Youdin, A. N. & Shu, F. H. 2004, in IAU Symposium, Vol. 202, Planetary Systems in the Universe, ed. A. Penny, 250
  * Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57

##  Appendix A The case of $K=0$

Here, we show the case of $K=0,$ where the diffusivity of pebbles is a fixed
constant: $D_{\rm p}=\alpha_{\rm mid}c_{\rm s}^{2}\Omega_{\rm
K}^{-1}/(1+\tau_{\rm s}^{2})$ (i.e., the back-reaction onto the diffusivity is
neglected). Figure 4 is the same as Fig. 2 but for the case of $K=0$.

As seen in the case of $K=1$ (Fig. 2), a runaway pile-up of pebbles within a
dead zone occurs as long as $F_{\rm p/g}$ meets the conditions discussed in
Section 3.2.2 (the analytical arguments were made for $K=0$). Compared to the
case of $K=1$, the radial width and the timescale of the pile-up are in
general narrower and longer in the case of $K=0$ because the diffusivity is a
constant for $K=0,$ while it becomes smaller for $K=1$ as pile-up proceeds.
These additional simulations with $K=0$ further show the robustness of the
physical mechanism presented in this letter and demonstrate that Drift-BKR,
which slows the drift velocity of pebbles, is a fundamental ingredient for the
ND mechanism.

Figure 4: Same as Fig. 2 but for the case of $K=0$.