# Linear growth of streaming instability in pressure bumps

Jérémy Auffinger1, Guillaume Laibe1,
1Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis-Laval, France E-mail:
<EMAIL_ADDRESS>
###### Abstract

Streaming instability is a powerful mechanism which concentrates dust grains
in protoplanetary discs, eventually up to the stage where they collapse
gravitationally and form planetesimals. Previous studies inferred that it
should be ineffective in viscous discs, too efficient in inviscid discs, and
may not operate in local pressure maxima where solids accumulate. From a
linear analysis of stability, we show that streaming instability behaves
differently inside local pressure maxima. Under the action of the strong
differential advection imposed by the bump, a novel unstable mode develops and
grows even when gas viscosity is large. Hence, pressure bumps are found to be
the only places where streaming instability occurs in viscous discs. This
offers a promising way to conciliate models of planet formation with recent
observations of young discs.

###### keywords:

protoplanetary discs — planets and satellites: formation

††pagerange: Linear growth of streaming instability in pressure
bumps–B††pubyear: 2017

##  1 Introduction

The main challenge of planet formation consists in figuring out how solids
originating from the interstellar medium concentrate and grow over orders of
magnitude up to form planetary cores (Chiang & Youdin, 2010). Up to decimetric
sizes, surface forces are strong enough for dust grains to grow by hit-and-
stick collisions (Blum & Wurm, 2008). This is not the case anymore for larger
pebbles, and solid aggregates are expected instead to undergo bouncing or
fragmentation (e.g. Güttler et al. 2010; Zsom et al. 2010). On the other hand,
rocky structures should typically reach hundreds of metres in size to be glued
by their own gravity. Thus, a third mechanism must bridge the gap and collect
pebbles up the the stage where their local weight becomes sufficient. Goodman
& Pindor (2000) suggested that such a concentration may originate from an
hydrodynamical instability. Youdin & Goodman (2005); Youdin & Johansen (2007)
demonstrated that the flow resulting from the radial drift of dust particles
in weakly viscous discs is actually linearly unstable. In the non-linear
regime, this so-called streaming instability develops dust over-concentrations
(Johansen & Youdin, 2007), which may ultimately form planetesimals by
gravitational instabilities in discs of sufficient metallicities (Johansen et
al., 2007; Johansen et al., 2009; Bai & Stone, 2010b, c; Carrera et al.,
2015). Hence, streaming instability may be responsible for the Initial Mass
Function of planetesimals in discs (Simon et al., 2016; Schäfer et al., 2017).
The robustness of the streaming instability has been tested against several
numerical schemes (Balsara et al., 2009; Miniati, 2010; Tilley et al., 2010;
Bai & Stone, 2010a; Johansen et al., 2012, 2014), towards the aim of
simulating its effect in a global disc (Lyra & Kuchner, 2013; Kowalik et al.,
2013; Yang & Johansen, 2014). Other physical processes such as vortices
(Raettig et al., 2015), photo-evaporation (Carrera et al., 2017), presence of
small grains (Laibe & Price, 2014), grain growth (Dra̧żkowska & Dullemond,
2014) or snow lines (Schoonenberg & Ormel, 2017) may reinforce the ability of
streaming instability to concentrate dust.

Discs may therefore contain only a moderate amount of dust grains, since up to
$50\%$ of their retained solid material may be converted into planetesimals
(e.g. Johansen et al. 2015; Dra̧żkowska et al. 2016). However, the emission
from a continuous dense phase of millimetre grains is commonly detected in
young discs (e.g. ALMA Partnership et al. 2015; Andrews et al. 2016), except
at some specific locations. Dark rings are often associated to ongoing planet
formation (e.g. Zhang et al. 2015; Gonzalez et al. 2015; Okuzumi et al. 2016)
or even to planets (e.g. Dipierro et al. 2015; Dong et al. 2015; Picogna &
Kley 2015; Rosotti et al. 2016). To explain the persistence of the dust
population almost everywhere, one may invoke the turbulent viscosity of the
gas, which damps efficiently the small-scale perturbations at which streaming
instability develops, but this would prevent planetesimal formation (see
however Johansen et al. 2007; Dittrich et al. 2013).

Such an interpretation is based on properties of the streaming instability
derived for discs with monotonically decreasing pressure profiles. On the
other hand, local pressure maxima may be created in the disc at some locations
by internal processes (e.g. Ruge et al. 2016; Estrada et al. 2016; Béthune et
al. 2016; Gonzalez et al. 2017). These pressure bumps are privileged locations
for planetesimal formation since they concentrate dust (e.g. Nakagawa et al.
1986; Haghighipour 2005). Finding a way for streaming instability to develop
specifically in these pressure bumps would conciliate current scenarii of
planetesimal formation and recent observations of discs. Simulations have been
recently performed to show that streaming instability may deform the bump
(Taki et al., 2016). So far, the resilience of the streaming instability
against viscosity in pressure bumps has not been investigated.

In this study, we show that the development of the streaming instability in
local pressure maxima is more complex than for discs with monotonic pressure
profiles. We address the problem analytically by performing a linear
perturbation analysis in a shearing box centred around a pressure maximum. In
Sect. 2, we derive new solutions for the steady state, since pressure
curvature provides additional advection compared to the usual case. The
analysis of the unstable modes of the system as a function of the steepness of
the bump is performed using a WKB approximation. The apparition of a second
unstable mode for the streaming instability and its resilience against viscous
damping are analysed in Sect. 3. Properties of this mode are brought back into
the context of planet formation in Sect. 4.

##  2 Equations of motion

###  2.1 Evolution in a global disc

####  2.1.1 Hypothesis

We consider non self-gravitating, non magnetic, vertically isothermal discs
made of perfect gas. The background surface densities and temperatures of the
gas are modelled by power-laws of decreasing exponents, $\Sigma\propto r^{-p}$
and $T\propto r^{-q}$. We use $p=1$ and $q=0.4$ to be consistent with models
that include detailed radiative transfer (e.g. Pinte & Laibe 2014). Under
these assumptions, the gas density $\rho_{\mathrm{g}}$ and the pressure $P$ in
the midplane of the disc scale as $\rho_{\rm g}\propto r^{-\xi}$ and $P\propto
r^{-\xi-q}$, where $\xi\equiv p-q/2+3/2$. The effective viscosity of the disc
is parametrised using an alpha prescription (Shakura & Sunyaev, 1973).
Numerical simulations of visco-turbulent discs exhibit values of $\alpha\sim
10^{-3}-10^{-2}$ (e.g. Meheut et al. 2015). Dust grains are assumed to be
compact, spherical, uncharged, of constant size and density. In observed discs
(e.g. Williams & Best 2014), local gas surface densities are low enough for
dust grains to be in the dilute Epstein drag regime (Epstein, 1924; Baines et
al., 1965). The drag stopping time of the particles is denoted
$t_{\mathrm{stop}}$. For millimetre-in-size grains, $t_{\mathrm{stop}}$ is of
order of the orbital period at $\sim 50\,$AU (Laibe et al., 2012). The ratio
between the drag and the orbital times, often called the Stokes number of the
flow, is denoted $\tau_{\rm s}$, consistently with the notations of Youdin &
Goodman (2005). The dust phase is modelled by a continuous viscousless and
pressureless fluid (Saffman, 1962; Garaud et al., 2004). The local dust-to-gas
ratio $\epsilon\equiv\rho_{\rm p}/\rho_{\rm g}$ is larger than the typical
$1\%$ of the interstellar medium since dust concentrates vertically and
radially in the disc. We follow Youdin & Goodman (2005) and neglect the
vertical stratification of the disc. Our study is therefore relevant for
grains with typical sizes $\gtrsim 10\mu$m that have settled close enough to
the midplane (Dubrulle et al., 1995; Fromang & Nelson, 2009).

####  2.1.2 Pressure maximum

Figure 1: Zoomed-in density profiles around pressure bumps. Solid and dashed-
dotted lines correspond to relative width of the bump
$\sigma/r_{\mathrm{u}}=1$ and $\sigma/r_{\mathrm{u}}=10$, while grey and black
lines correspond to relative amplitudes $\tilde{A}=1$ and $\tilde{A}=10$
respectively. The gas density profile of exponents $p=1$, $q=0.4$ in absence
of perturbation is given as reference (dashed line/light grey).

We model a local pressure maximum by superimposing a Gaussian perturbation to
the usual gas density profile, i.e.

$\rho_{\mathrm{g}}(r)=\rho_{0}\left[\left(\dfrac{r}{r_{\mathrm{u}}}\right)^{-\xi}+A\mathrm{e}^{-\dfrac{(r-r_{0})^{2}}{2\sigma^{2}}}\right]\,.$
(1)

Hence, the amplitude and the width of the density maximum are parametrised by
$A$ and $\sigma$ respectively. The radial coordinate is scaled with a radius
$r_{\rm u}$ and the Gaussian bump is centred around a position $r_{0}$. Fig. 1
shows different shapes of density profiles obtained when varying $A$ and
$\sigma$. Note that $A$ should be large enough for the Gaussian perturbation
to dominate locally over the decreasing background and thus, for the pressure
maximum to exist. The relative amplitude $\tilde{A}$ of the maximum
respectively to the background is

$\tilde{A}\equiv A\left(r_{0}/r_{\rm u}\right)^{\xi}.$ (2)

$\tilde{A}$ varies from $\sim 0.1$ for a perturbation due to a Neptune-like
mass planet (e.g. Dipierro & Laibe 2017) up to $\sim 10$ for self-induced dust
traps (Gonzalez et al., 2017). The width of the bump $\sigma$ is of the order
of $\sim H$, the pressure scale height. In absence of any pressure maximum, or
when the pressure perturbation is negligible, the orbital correction with
respect to a pure Keplerian rotation is of order
$\mathcal{O}\left(H^{2}/r^{2}\right)$. In a pressure bump, the orbital
correction is of order $\mathcal{O}\left(\tilde{A}\,H^{2}/\sigma^{2}\right)$
and should remain small enough for the disc to be supported by rotation.

###  2.2 Shearing box approximation

For simplicity, the evolution of gas and dust is studied in a local frame
corotating with the disc at a location $\hat{r}_{0}$ and a frequency
$\Omega_{0}$. In this shearing box, the coordinates are expanded to the linear
order (Goldreich & Lynden-Bell, 1965), and $x\equiv r-\hat{r}_{0}$,
$y\equiv\hat{r}_{0}\left(\theta-\Omega_{0}t\right)$ and $z$ denote the radial,
azimuthal and vertical directions respectively.

Figure 2: Comparison between the gas density and its expansion used in the
shearing box approximation (dashed and solid lines respectively). Errors are
of order $1\%$ in $\mathrm{L}^{1}$ norm, peaking at $5\%$ at the edges of the
box of width $\pm\sigma$ delimited by a grey zone. The amplitude and width of
the bump are $\tilde{A}=1$ and $\sigma/r_{\rm u}=1$.

In the shearing box approximation, the large-scale contribution of the
background pressure gradient comes under the form of a constant force (Youdin
& Goodman, 2005). In a box centred around a pressure bump, the large-scale
background pressure gradient is a linear function of the distance from the
pressure maximum. For convenience, we centre the box around
$\hat{r}_{0}=r_{0}$, the maximum of the Gaussian perturbation. Expanding the
pressure gradient term to the second order in $x$ provides

$-\dfrac{\bm{\nabla}P}{\rho_{\mathrm{g}}}\simeq
2r_{0}\Omega_{0}^{2}\left(\eta+\dfrac{\Gamma}{2r_{0}}x\right)\mathrm{\textbf{u}}_{x}\,,$
(3)

with

$\displaystyle\eta\equiv\dfrac{1}{2}\left.\frac{\mathrm{d}\ln\rho_{\mathrm{g}}}{\mathrm{d}\ln
r}\right|_{r_{0}}\left(\dfrac{H_{0}}{r_{0}}\right)^{2}\,,$ (4)
$\displaystyle\Gamma\equiv
r_{0}^{2}\left.\frac{\mathrm{d}^{2}\ln\rho_{\mathrm{g}}}{\mathrm{d}r^{2}}\right|_{r_{0}}\left(\dfrac{H_{0}}{r_{0}}\right)^{2}\,,$
(5)

where $H_{0}$ denotes the scale height of the disc at $r_{0}$. Following
Youdin & Goodman (2005) notations, we denote
$\eta_{0}=\left.\eta\right|_{A=0}$ the pressure gradient term in absence of
bump. Young protoplanetary discs are denser and warmer in the inner regions,
which implies $\eta_{0}>0$. In typical discs, $\eta_{0}\simeq 10^{-2}$. Note
that discs with power-law profiles have in general non-zero values for
$\Gamma$. This contribution from the curvature of the density profile is
neglected in Youdin & Goodman (2005), which is an excellent approximation. A
pressure bump is defined by $\Gamma>0$. The pressure maximum position $x_{\rm
max}\equiv-2r_{0}\eta/\Gamma$ is slightly shifted with respect to the centre
of the box as a result of the small contribution of the decreasing unperturbed
pressure profile. Using Eq. 1 in Eqs. 4 – 5 gives

$\displaystyle\eta\equiv$
$\displaystyle\dfrac{\xi}{2(1+\tilde{A})}\left(\dfrac{H_{0}}{r_{0}}\right)^{2},$
(6) $\displaystyle\Gamma\equiv$
$\displaystyle\left\\{-\dfrac{\xi(\xi+1)}{1+\tilde{A}}+\dfrac{\xi^{2}}{(1+\tilde{A})^{2}}\right\\}\left(\dfrac{H_{0}}{r_{0}}\right)^{2}+\dfrac{\tilde{A}}{1+\tilde{A}}\left(\dfrac{H_{0}}{\sigma}\right)^{2}\,.$
(7)

The first term of the right-hand side of Eq. 7 corresponds to the contribution
of the background profile, which increases with $\xi$ (steeper and more curved
density profiles). The second term corresponds to the contribution added by
the Gaussian perturbation and scales like
$\tilde{A}\left(H_{0}/\sigma\right)^{2}$. The relative contribution between
the maximum and the background is therefore of order
$\tilde{A}\left(r_{0}/\sigma\right)^{2}$.

Figure 3: Similar to Fig. 2: comparison between the normalised pressure force
(dashed line) and its linearisation expansion (solid line) around $r_{0}$, the
maximum of the Gaussian perturbation. Errors are of order $10\%$ in
$\mathrm{L}^{1}$ norm, peaking at $20\%$ at the edges of the box.

Fig. 2 shows that the size of the box is constrained by the validity of the
linear approximation. For Eq. 3 to remain valid, the size of the box should
not exceed $\sim\min\left(\sigma,H\right)$. Note that on the other hand, very
narrow bumps may be Rayleigh unstable (Yang & Menou, 2010). In practice, we
use $\sigma=1\,r_{\rm u}$. The discrepancy between the real pressure force and
its linear approximation is of order $\sim 20\%$ at most at the edges of the
box (Fig. 3). Viscous forces are similarly decomposed in large and small
scales components. In a typical $\alpha-$disc, the relative contribution
between the viscosity and the radial pressure gradient is of order $\alpha\ll
1$ at large scale. This contribution is therefore neglected. On the other
hand, the small scale viscous forces, which damp local gas fluctuations, are
treated as usual. To avoid unnecessary complications, the local gas sound
speed $c_{\rm s}$, the viscosity $\nu=\alpha c_{\rm s}H$, the stopping time of
dust grains $t_{\mathrm{stop}}$ and the background dust-to-gas ratio
$\epsilon$ are assumed to be constant over the size of the box. This implies
that there is more dust in the centre of the bump than at the edges. The
equations of motion for the gas and the dust are therefore

$\displaystyle\dfrac{\partial\rho_{\mathrm{g}}}{\partial
t}+\bm{\nabla}\cdot(\rho_{\mathrm{g}}\mathrm{\textbf{V}}_{\mathrm{g}})=0\,,$
(8) $\displaystyle\dfrac{\partial\rho_{\mathrm{p}}}{\partial
t}+\bm{\nabla}\cdot(\rho_{\mathrm{p}}\mathrm{\textbf{V}}_{\mathrm{p}})=0\,,$
(9) $\displaystyle\left(\dfrac{\partial}{\partial
t}+\mathrm{\textbf{V}}_{\mathrm{g}}\cdot\bm{\nabla}\right)\mathrm{\textbf{V}}_{\mathrm{g}}=$
$\displaystyle-
r_{0}x\left.\dfrac{\mathrm{d}\Omega_{\mathrm{K}}^{2}}{\mathrm{d}r}\right|_{r_{0}}\mathrm{\textbf{u}}_{x}-2\Omega_{0}\mathrm{\textbf{u}}_{z}\times\mathrm{\textbf{V}}_{\mathrm{g}}$
$\displaystyle+2r_{0}\Omega_{0}^{2}\left(\eta+\dfrac{\Gamma}{2r_{0}}x\right)\mathrm{\textbf{u}}_{x}+\nu\Delta\mathrm{\textbf{V}}_{\mathrm{g}}\,$
$\displaystyle+\dfrac{\rho_{\mathrm{p}}}{\rho_{\mathrm{g}}}\dfrac{\mathrm{\textbf{V}}_{\mathrm{p}}-\mathrm{\textbf{V}}_{\mathrm{g}}}{t_{\mathrm{stop}}}\,,$
(10) $\displaystyle\left(\dfrac{\partial}{\partial
t}+\mathrm{\textbf{V}}_{\mathrm{p}}\cdot\bm{\nabla}\right)\mathrm{\textbf{V}}_{\mathrm{p}}=$
$\displaystyle-
r_{0}x\left.\dfrac{\mathrm{d}\Omega_{\mathrm{K}}^{2}}{\mathrm{d}r}\right|_{r_{0}}\mathrm{\textbf{u}}_{x}-2\Omega_{0}\mathrm{\textbf{u}}_{z}\times\mathrm{\textbf{V}}_{\mathrm{p}}$
$\displaystyle-\dfrac{\mathrm{\textbf{V}}_{\mathrm{p}}-\mathrm{\textbf{V}}_{\mathrm{g}}}{t_{\mathrm{stop}}}\,.$
(11)

In particular, gas and dust have different advection velocities
$\mathrm{\textbf{V}}_{\mathrm{g,p}}$ and the drag from the dust onto the gas
is not neglected (Youdin & Goodman, 2005). For simplicity, physical quantities
are used in a dimensionless form, i.e.
$\tilde{\omega}\equiv\omega/\Omega_{0}$, $\tilde{x}\equiv x/\eta_{0}r_{0}$,
$\tilde{k}\equiv k\eta_{0}r_{0}$, $\tilde{v}\equiv v/\eta_{0}r_{0}\Omega_{0}$.
As a remark, Eqs. 8 – 11 differ from the system studied by Taki et al. (2016),
where the large scale pressure gradient is similar to the constant background
introduced in Youdin & Goodman (2005), and on top of which a small local
Gaussian perturbation of width half of the box is superimposed. This situation
corresponds to a pressure perturbation developing at small scales that is weak
enough for not affecting background radial velocities of the gas and the dust.

###  2.3 Steady state

Figure 4: Radial and azimuthal (top/bottom) velocities for the gas and the
dust (left/right) in the box. Bumps with relative amplitudes
$\tilde{A}=0.1-1-10$ are considered (dashed, dot-dashed and dotted black
lines), other parameters being those of the linB problem. Velocities in
absence of pressure bump are indicated for reference (grey dotted lines).

The steady state velocities for gas and dust in the pressure maximum are
determined by seeking for solutions of the form
$\left(\overline{U}_{\mathrm{g}},\overline{V}_{\mathrm{g}}\right)=\left(a_{\mathrm{g}}x+b_{\mathrm{g}},\alpha_{\mathrm{g}}x+\beta_{\mathrm{g}}\right)$
and
$\left(\overline{U}_{\mathrm{p}},\overline{V}_{\mathrm{p}}\right)=\left(a_{\mathrm{p}}x+b_{\mathrm{p}},\alpha_{\mathrm{p}}x+\beta_{\mathrm{p}}\right)$.
The linear dependancy of the velocities with respect to $x$ is enforced to
ensure consistency between the presence of the pressure maximum and the
shearing box formalism. Using this Ansatz in Eqs. 8 – 11 provides the non-
linear system of equations

$\displaystyle
a_{\mathrm{p}}^{2}-3\Omega_{0}^{2}-2\Omega_{0}\alpha_{\mathrm{p}}+\dfrac{1}{t_{\mathrm{stop}}}(a_{\mathrm{p}}-a_{\mathrm{g}})=0\,,$
(12) $\displaystyle
a_{\mathrm{p}}b_{\mathrm{p}}-2\Omega_{0}\beta_{\mathrm{p}}+\dfrac{1}{t_{\mathrm{stop}}}(b_{\mathrm{p}}-b_{\mathrm{g}})=0\,,$
(13) $\displaystyle
a_{\mathrm{p}}\alpha_{\mathrm{p}}+2\Omega_{0}a_{\mathrm{p}}+\dfrac{1}{t_{\mathrm{stop}}}(\alpha_{\mathrm{p}}-\alpha_{\mathrm{g}})=0\,,$
(14) $\displaystyle
b_{\mathrm{p}}\alpha_{\mathrm{p}}+2\Omega_{0}b_{\mathrm{p}}+\dfrac{1}{t_{\mathrm{stop}}}(\beta_{\mathrm{p}}-\beta_{\mathrm{g}})=0\,,$
(15) $\displaystyle
a_{\mathrm{g}}^{2}-3\Omega_{0}^{2}-2\Omega_{0}\alpha_{\mathrm{g}}-\dfrac{\epsilon}{t_{\mathrm{stop}}}(a_{\mathrm{p}}-a_{\mathrm{g}})-\Gamma\Omega_{0}^{2}=0\,,$
(16) $\displaystyle
a_{\mathrm{g}}b_{\mathrm{g}}-2\Omega_{0}\beta_{\mathrm{g}}-\dfrac{\epsilon}{t_{\mathrm{stop}}}(b_{\mathrm{p}}-b_{\mathrm{g}})-2r_{0}\Omega_{0}^{2}\eta=0\,,$
(17) $\displaystyle
a_{\mathrm{g}}\alpha_{\mathrm{g}}+2\Omega_{0}a_{\mathrm{g}}-\dfrac{\epsilon}{t_{\mathrm{stop}}}(\alpha_{\mathrm{p}}-\alpha_{\mathrm{g}})=0\,,$
(18) $\displaystyle
b_{\mathrm{g}}\alpha_{\mathrm{g}}+2\Omega_{0}b_{\mathrm{g}}-\dfrac{\epsilon}{t_{\mathrm{stop}}}(\beta_{\mathrm{p}}-\beta_{\mathrm{g}})=0\,.$
(19)

These non-linear terms originate from non-trivial advection terms specific to
the pressure bump and require care for numerical root finding (see Appendix B
for technical details). Neglecting this additional advection provides a crude
approximation of the solution of Eqs. 12 – 19 by taking the solution of
Nakagawa et al. (1986) given in Appendix B and replacing $\eta$ by
$\eta+\frac{\Gamma}{2r_{0}}x$. We find errors of order $\sim 10\%$ up to $\sim
100\%$ between the two approaches, a discrepancy becoming important for the
radial velocity of the gas. Fig. 4 illustrates the dust and the gas motion
inside the bump for various relative amplitudes of the pressure maximum. As a
reminder, pure Keplerian shear is $v_{y}=-\frac{3}{2}\Omega_{0}x$. In absence
of a bump, dust (resp. gas) drifts inwards (resp. outwards) by conservation of
angular momentum. Both the gas and the dust are sub-Keplerian. Inside a bump,
dust drifts towards the pressure maximum, while gas drifts outwards. This
requires for gas and dust to orbit at super-Keplerian frequency in the inner
edge of the bump. Both the gas and the dust radial velocities are rigorously
zero at $x_{\rm max}<0$ the location of the pressure maximum. Although the
different velocities appear to cross each other at the same location in Fig.
4, this is actually not the case. Mathematically, the intersection between the
line corresponding to $\tilde{A}=0$ and the other lines depends slightly on
the different physical parameters. In particular, for increasing values of
$\Gamma$, the intersecting point becomes closer to the centre of the box.

###  2.4 Perturbation

The linear stability of the system Eqs. 8 – 11 is investigated by looking for
perturbations of the form

$\displaystyle\rho_{\mathrm{g}}=\rho_{\mathrm{g}}^{0}(1+\delta_{\mathrm{g}})\,,$
(20)
$\displaystyle\rho_{\mathrm{p}}=\rho_{\mathrm{p}}^{0}(1+\delta_{\mathrm{p}})\,,$
(21)
$\displaystyle\mathrm{\textbf{V}}_{\mathrm{g}}=\overline{\mathrm{\textbf{V}}}_{\mathrm{g}}+\mathrm{\textbf{v}}_{\mathrm{g}}=\overline{\mathrm{\textbf{V}}}_{\mathrm{g}}+u_{\mathrm{g}}\mathrm{\textbf{u}}_{x}+v_{\mathrm{g}}\mathrm{\textbf{u}}_{y}+w_{\mathrm{g}}\mathrm{\textbf{u}}_{z}\,,$
(22)
$\displaystyle\mathrm{\textbf{V}}_{\mathrm{p}}=\overline{\mathrm{\textbf{V}}}_{\mathrm{p}}+\mathrm{\textbf{v}}_{\mathrm{p}}=\overline{\mathrm{\textbf{V}}}_{\mathrm{p}}+u_{\mathrm{p}}\mathrm{\textbf{u}}_{x}+v_{\mathrm{p}}\mathrm{\textbf{u}}_{y}+w_{\mathrm{p}}\mathrm{\textbf{u}}_{z}\,,$
(23)

where

$\delta(x,z,t)=\Delta(x)e^{i(k_{x}x+k_{z}z-\omega t)}\,.$ (24)

In absence of a pressure bump, the perturbation must develop in both the $x$
and the $z$ direction to become unstable (Youdin & Goodman, 2005; Jacquet et
al., 2011). This property originates from local conservation of the gas mass.
With a pressure bump, advection terms enforce the amplitude of the
perturbation $\Delta$ to depend on $x$. For simplicity, we focus on cases
where this amplitude varies slowly compared to the phase ($Hk_{x}\gg 1$) and
use a WKB approximation to compute spatial derivatives, i.e.

$\frac{\partial\delta(x,z,t)}{\partial x}\simeq
ik_{x}\Delta(x)e^{i(k_{x}x+k_{z}z-\omega t)}.$ (25)

We obtain the following set of 8 linear equations for the perturbation

$\displaystyle-i\omega\delta_{\mathrm{g}}$
$\displaystyle+i\delta_{\mathrm{g}}k_{x}\overline{U}_{\mathrm{g}}+\delta_{\mathrm{g}}a_{\mathrm{g}}+ik_{x}u_{\mathrm{g}}+ik_{z}w_{\mathrm{g}}=0\,,$
(26) $\displaystyle-i\omega\delta_{\mathrm{p}}$
$\displaystyle+i\delta_{\mathrm{p}}k_{x}\overline{U}_{\mathrm{p}}+\delta_{\mathrm{p}}a_{\mathrm{p}}+ik_{x}u_{\mathrm{p}}+ik_{z}w_{\mathrm{p}}=0\,,$
(27) $\displaystyle-i\omega\mathrm{\textbf{v}}_{\mathrm{g}}$
$\displaystyle+ik_{x}\overline{U}_{\mathrm{g}}\mathrm{\textbf{v}}_{\mathrm{g}}+(u_{\mathrm{g}}a_{\mathrm{g}}-2v_{\mathrm{g}}\Omega_{0})\mathrm{\textbf{u}}_{x}+u_{\mathrm{g}}(\alpha_{\mathrm{g}}+2\Omega_{0})\mathrm{\textbf{u}}_{y}$
$\displaystyle-\dfrac{\epsilon}{t_{\mathrm{stop}}}(\mathrm{\textbf{v}}_{\mathrm{p}}-\mathrm{\textbf{v}}_{\mathrm{g}}+(\delta_{\mathrm{p}}-\delta_{\mathrm{g}})(\overline{\mathrm{\textbf{V}}}_{\mathrm{p}}-\overline{\mathrm{\textbf{V}}}_{\mathrm{g}}))$
$\displaystyle+ic_{s}^{2}\delta_{\mathrm{g}}\mathrm{\textbf{k}}+\nu(k_{x}^{2}+k_{z}^{2})\mathrm{\textbf{v}}_{\mathrm{g}}=0\,,$
(28) $\displaystyle-i\omega\mathrm{\textbf{v}}_{\mathrm{p}}$
$\displaystyle+ik_{x}\overline{U}_{\mathrm{p}}\mathrm{\textbf{v}}_{\mathrm{p}}+(u_{\mathrm{p}}a_{\mathrm{p}}-2v_{\mathrm{p}}\Omega_{0})\mathrm{\textbf{u}}_{x}+u_{\mathrm{p}}(\alpha_{\mathrm{p}}+2\Omega_{0})\mathrm{\textbf{u}}_{y}$
$\displaystyle+\dfrac{1}{t_{\mathrm{stop}}}(\mathrm{\textbf{v}}_{\mathrm{p}}-\mathrm{\textbf{v}}_{\mathrm{g}})=0\,.$
(29)

where the four background velocities
$\overline{U}_{\mathrm{g}},\overline{U}_{\mathrm{p}},\overline{V}_{\mathrm{g}},\overline{V}_{\mathrm{p}}$
are the linear functions determined in Sect. 2.3 and not the solutions derived
by Nakagawa et al. (1986). In practice, the eigenmodes $\omega_{1,8}$ of the
system Eqs. 26 – 29 are determined by finding zeros of the determinant of the
perturbation matrix numerically, using a sufficient precision. Note that the
eigenmodes obtained by this procedure depend on $x$, i.e. $\omega_{1,8}(x)$,
which may sound inconsistent with Eq. 24. However, at small times,

$k_{x}^{-1}\frac{\mathrm{d}\omega}{\mathrm{d}x}t\simeq\frac{\Omega_{0}t}{k_{x}H}.$
(30)

In this case, the solution is consistent with the initial Anzatz and the WKB
approximation over a number $n=k_{x}H\gg 1$ of orbital periods. To compare our
results with the case of an inviscid disc with no pressure bump, we use the
test cases linA, linB, linC, linD studied in Youdin & Johansen (2007); Bai &
Stone (2010a) (see parameters in Appendix A). Our procedure provides the
expected coefficients with similar precision.

##  3 Results

###  3.1 Unstable modes

Figure 5: Growth rates of the two epicyclic modes and the secular mode of the
streaming instability for increasing values of the relative amplitude
$\tilde{A}$ of the bump. In absence of any pressure maximum, the only unstable
mode is the usual secular mode (dot-dashed/black). In pressure maxima, a novel
unstable mode develops (dashed/black). A bifurcation between the secular mode
and the other epicycle (dotted/black) occurs at $\tilde{A}\simeq 0.2$.
Parameters correspond to the linA problem.

The growth rates $s$ obtained from the procedure described in Sect. 2.4 are
shown in Fig. 5. Three modes are considered: two modified epicyclic modes and
the secular mode of the streaming instability identified in Youdin & Goodman
(2005), the other modes playing no particular role in this problem. For tiny
perturbations of the pressure profile, the only unstable mode is the secular
mode, similarly to what happens in a disc with no bump. Novel features appear
when increasing progressively the relative amplitude of the bump. Fig. 5 shows
that one of the two epicyclic modes becomes unstable. The related growth rate
may become larger than the growth rate of the secular mode as the amplitude of
the bump increases. The possible instability of the epicyclic modes was
mentioned in Youdin & Goodman (2005). When increasing further the amplitude of
the bump, a transition between two distinct regimes occurs for relative
amplitudes of order $\tilde{A}\gtrsim 0.2$, a conservative value for typical
bumps. Fig. 5 shows that a bifurcation between the secular and the stable
epicyclic modes that gives birth to two new modes. We find that this
bifurcation is universal and does in particular not depend on the dust-to-gas
ratio. The exact number of unstable modes depends on the distance to the
pressure maximum. Near $x_{\rm max}$, two modes are unstable, whereas at the
edges of the bump, only one mode is unstable. This mode corresponds to the
novel unstable mode originating from the bifurcation and is not the secular
mode of the streaming instability which grows in absence of maximum.
Physically, this novel instability originates from the strong differential
advection between the gas and the dust at the edges of the bump powered by the
background velocities in the bump.

Fig. 6 shows how the growth rates depend on the radial location in the bump
for increasing amplitudes, showing the consequences of the bifurcation
identified above. Rigorously, the exact location of the pressure maximum
depends on $\tilde{A}$ and is slightly offset from the centre of the shearing
box due to the local curvature of the pressure profile in absence of maximum
(see Sect. 2.3). For shallow bumps ($\tilde{A}\lesssim 0.2$), streaming
instability develops everywhere except close to the pressure maximum. Its
efficiency is maximum at the edges of the bump, where the local pressure
gradient is the greatest. This behaviour is consistent with the linear
analysis of Youdin & Goodman (2005), valid for discs with monotonic pressure
profiles. For larger amplitudes ($\tilde{A}\gtrsim 0.2$), the streaming
instability grows more efficiently in the centre of the bump than at the
edges. As shown in Fig. 6, its efficiency is slightly reduced at $x=x_{\rm
max}$. At the edges of the bump, the growth rate decreases up to eventually
reach zero. No instability develops in this particular case. Hence, for a
shallow bump in an inviscid disc, gas and dust are linearly unstable
everywhere except at the exact location of the pressure maximum. On the
opposite, when the disc contains a bump that is large enough (i.e. for
pressure perturbation to be of a least a few ten percents), dust concentration
may occur preferentially inside the maximum, and not at the edges, as one
would have expected with the classical linear stability analysis.

Figure 6: Growth rate of the streaming instability in the box for bumps of
relative amplitudes $\tilde{A}=0.1-0.25-1-10$ (dashed/grey, dot-dashed/grey,
dotted/black and dashed/black respectively). For shallow bumps, streaming
instability develops everywhere except near the maximum where the pressure
gradient is zero. For large bumps, the instability develops more efficiently
close to the maximum. Parameters are those of the linD problem. The growth
rate in absence of maximum is given by the dotted/grey line. Figure 7: Growth
rate of the instability for two same ratios $\lambda/l_{\rm stop}=1$. Growth
rates are almost similar, discrepancies of a few percent originating from
slightly different advection velocities. We use $\tilde{A}=10$, $\epsilon=2$,
$\tau_{\mathrm{s}}=0.01$, $\tilde{k}_{x,z}=100$ (black dashed line) and
$\tau_{\mathrm{s}}=0.001$, $\tilde{k}_{x,z}=1000$ (grey solid line).

The growth rates depend on the ratio between the wavelength $\lambda$ of the
perturbation and $l_{\rm
stop}~{}\equiv~{}\eta_{0}r_{0}\Omega_{0}t_{\mathrm{stop}}$, the length over
which the gas decouples from the dust, sometimes referred as the stopping
length. The ratio $\lambda/l_{\rm stop}$ measures the number of perturbations
over which the stopping length spreads. Fig. 7 shows that for similar values
of $\lambda/l_{\rm stop}$ or equivalently $\tilde{k}_{x,z}\,\tau_{\rm s}$, the
growth rates obtained are almost identical. Corrections of order a few
percents are due to slightly different values for the advection velocities.
Fig. 8 shows growth rates obtained for different wavelengths $\lambda$. For
$l_{\rm stop}\gg\lambda$, dust and gas experience the details of the pressure
profile before being coupled together by the drag, and the growth rate profile
is narrow. Instead, for $l_{\rm stop}<\lambda$, dust and gas are quickly
coupled by the drag, their differential velocity is proportional to the local
pressure gradient, and this information is carried away by the perturbation.
In this case the growth rate profile is wider. The ratio
$\lambda/l_{\mathrm{stop}}$ sets the width of the region where no instability
develops for shallow bumps, and the width of the central region where the
instability is weakened for large bumps.

Figure 8: Growth rates for a bump of relative amplitude $\tilde{A}=10$ and
different wavenumbers $\tilde{k}_{x,z}=50-100-300-1500$ (dotted/black, dot-
dashed/grey, dashed/grey, dotted/grey). The region where the instability
develops is more extended for larger wavelengths. Parameters are those of the
linC problem.

###  3.2 Viscosity

Figure 9: Growth rate of the instability in absence of bump for increasing
viscosities. The unstable mode is damped for $\alpha\gtrsim 10^{-5}$.
Parameters correspond to the linA problem.

The ability for viscosity to damp the instability is now investigated. Fig. 9
shows that in absence of a pressure bump, viscosity prevents the development
of the streaming instability for values of $\alpha$ as low as $\sim 10^{-5}$.
This value depends on the wavelength of the perturbation, smaller fluctuations
being damped more efficiently by viscosity. Streaming instability does
therefore not grow in typical visco-turbulent discs where $\alpha\sim
10^{-3}-10^{-2}$. However, when a significant bump is present in the disc, a
different behaviour is observed. Fig. 10 shows that for increasing values of
$\alpha$, one of the two unstable modes is suppressed, whereas the other one
is only weakened, but not damped. The growth rate is reduced by one order of
magnitude compared to the inviscid case (see Fig. 11). Still, this mode
develops in a time relevant for planetesimal formation (see Sect.4). Hence,
streaming instability is found to always develop in a pressure bump, even in
highly viscous discs. The fact that viscosity does not entirely damp the
instability in the bump may appear counterintuitive. This apparent conundrum
can be explained by noting that viscosity damps only the perturbations of the
gas velocity. Fig. 12 shows that these perturbations are actually suppressed
more and more efficiently when viscosity increases. However, in pressure
bumps, the required gradients of back-reaction are provided by the background
velocities resulting from the local pressure profile. Hence, the instability
can grow even if perturbations in the gas velocity are damped by viscosity.
This effect is not observed in a disc with monotonic pressure profiles, since
local gradients of back-reaction originate only from perturbations of the gas
velocity that are killed by viscosity. Consistently, growth rates in the
viscous regime do almost not depend on the wavelength, as shown in Fig. 13.
Fig. 14 and Fig. 15 show that similarly to the usual case, streaming
instability is most efficient for $\tau_{\mathrm{s}}\sim 1$ and $\epsilon\sim
1$, such parameters being typical protoplanetary discs.

Figure 10: Growth rates of the unstable modes for increasing viscosities. An
unstable mode exists even in highly viscous discs. Its growth rate is weakened
by one order of magnitude compared to the inviscid case. Parameters correspond
to the linA problem. Figure 11: Growth rates for viscous parameters
$\alpha=0-10^{-6}-10^{-5}-10^{-4}-10^{-3}$ (dashed/light, dotted, dashed,
dash-dotted grey and dashed/black lines respectively). Parameters are those of
the linA problem, the relative amplitude of the bump is $\tilde{A}=1$. Figure
12: Amplitude of the perturbation of the radial velocity of the gas
$|u_{\mathrm{g}}|$ for increasing viscosities. The fluctuations are damped
efficiently at large viscosities. The instability does not require
fluctuations in the gas velocity to grow, since differential back-reaction is
provided by the bump itself. Parameters are those of the linA problem. Figure
13: Growth rate of the instability for different wavenumbers $\tilde{k}_{x,z}$
in a viscous disc with $\alpha=10^{-2}$. No dependancy is found, consistently
with a mechanism powered by the background profile of the bump. Parameters are
those of the linA problem. Figure 14: Growth rate of the instability as a
function of the Stokes number $\tau_{\mathrm{s}}$ for $\alpha=10^{-2}$.
Maximal efficiency is reached for $\tau_{\mathrm{s}}\sim 1$ as expected.
Parameters are those of the linA problem. Figure 15: Growth rate of the
instability as a function of the dust-to-gas ratio $\epsilon$ for
$\alpha=10^{-2}$. Maximal efficiency is reached for $\epsilon\sim 1$ as
expected. Parameters are those of the linA problem.

Fig. 16 summarises where streaming instability develops in discs, whether a
pressure bump is present or not and whether the disc is inviscid or viscous. A
dusty disc with monotonic pressure profile is linearly unstable everywhere if
it is inviscid, and stable everywhere if it is viscous (Youdin & Goodman,
2005). The latter is not true anymore in presence of a pressure maximum:
streaming instability always grows in the bump, even if the disc is viscous.

Figure 16: Sketch illustrating where the streaming instability grows in
different discs. An inviscid disc is linearly unstable everywhere, with or
without a pressure bump. In visco-turbulent discs, the streaming instability
develops only inside pressure bumps.

##  4 Discussion

Figure 17: Timescales of streaming instability (dotted/black) and dust drift
(dashed/black) at the outer edge of the pressure bump. Top: inviscid disc,
instability grows faster than dust drifts. Bottom: viscous disc
($\alpha=10^{-3}$), the two timescales are comparable for significant pressure
bumps. Parameters correspond to the linA problem.

We now verify that the assumptions made in Sect. 3 are consistent with the
global evolution of the disc. The dust-to-gas ratio, the stopping time and the
background velocities have been assumed to be constant during the time $t_{\rm
growth}$ it takes for the instability to grow. Fig. 17 compares $t_{\rm
growth}$ to the drift time of dust into the bump $t_{\rm drift}$ and shows
that in the inviscid case, the disc remains in a steady state during the
growth of the perturbation since $t_{\mathrm{growth}}\ll t_{\mathrm{drift}}$.
In viscous discs, $t_{\rm growth}$ increases by one order of magnitude and the
two timescales become comparable. Long-lived global numerical simulations are
therefore required to determine how the instability develops in this case.
Moreover, the local depletion of gas caused by back-reaction empties the bump.
Its lifetime is compared to the growing time by estimating the ratio $t_{\rm
drift}/t_{\rm growth}\epsilon$. As long as $\epsilon\sim 1$, this condition is
similar to the precedent one. Hence, our assumptions are essentially valid in
real discs and streaming instability is expected to develop linearly as
described in Sect. 3.

Gonzalez et al. (2017) found that formation of pressure bumps can be triggered
by the sole action of dust drag onto the gas, with a spatial resolution
insufficient to capture streaming instability. Their Fig. 9 shows that the
relative amplitude of the bump is a least of order a few ten per cent, even
for viscous discs with $\alpha=10^{-2}$. Streaming instability in self-induced
dust traps offers therefore a promising way to build planetesimals at specific
locations. Millimetre grains are converted efficiently into planetesimals –
and potentially later on into planets – in the trap and only there. Such a
scenario is consistent with the dark rings probed recently by ALMA. Similarly,
streaming instability may also grow specifically at the edges of the gap
created by a massive planet. Since the tidal torque from the planet acts both
on dust and gas, it does not induce additional differential velocities and we
expect the results derived above to hold. Johansen et al. (2009) found that
for large enough metallicities ($Z\gtrsim 0.03$), dust clumps formed by
streaming instability collapse gravitationally. We expect a similar trend in a
pressure bump, although numerical simulations at high resolution are required
to investigate the non-linear stages of the instability.

##  5 Conclusion

In this study, we extend the linear theory of streaming instability in dusty
protoplanetary discs to include the eventual presence of local pressure
maxima, where solids drift in and pile-up. We find that for pressure bumps
with relative amplitude larger than $\sim 20\%$, a bifurcation occurs, giving
birth to a novel unstable mode. The instability is powered by the strong
differential advection locally imposed by the bump and as such, is resilient
against viscous damping from the gas. The growth rate of the instability is
typically reduced by one order of magnitude compared to the inviscid case.
Hence, in viscous discs, streaming instability is found to grow in and only in
local pressure maxima. In particular, streaming instability in self-induced
dust traps provides a scenario for the early stages of planet formation
consistent with the recent millimetre observations of dark rings in young
discs. Numerical simulations are required to investigate the further non-
linear development of the instability.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership et al. 2015, ApJ, 808, L3
  * Andrews et al. (2016) Andrews S. M., et al., 2016, ApJ, 820, L40
  * Bai & Stone (2010a) Bai X.-N., Stone J. M., 2010a, ApJS, 190, 297
  * Bai & Stone (2010b) Bai X.-N., Stone J. M., 2010b, ApJ, 722, 1437
  * Bai & Stone (2010c) Bai X.-N., Stone J. M., 2010c, ApJ, 722, L220
  * Baines et al. (1965) Baines M. J., Williams I. P., Asebiomo A. S., 1965, MNRAS, 130, 63
  * Balsara et al. (2009) Balsara D. S., Tilley D. A., Rettig T., Brittain S. D., 2009, MNRAS, 397, 24
  * Béthune et al. (2016) Béthune W., Lesur G., Ferreira J., 2016, A&A, 589, A87
  * Blum & Wurm (2008) Blum J., Wurm G., 2008, ARA&A, 46, 21
  * Carrera et al. (2015) Carrera D., Johansen A., Davies M. B., 2015, A&A, 579, A43
  * Carrera et al. (2017) Carrera D., Gorti U., Johansen A., Davies M. B., 2017, ApJ, 839, 16
  * Chiang & Youdin (2010) Chiang E., Youdin A. N., 2010, Annual Review of Earth and Planetary Sciences, 38, 493
  * Dipierro & Laibe (2017) Dipierro G., Laibe G., 2017, preprint (arXiv:1704.06664)
  * Dipierro et al. (2015) Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS, 453, L73
  * Dittrich et al. (2013) Dittrich K., Klahr H., Johansen A., 2013, ApJ, 763, 117
  * Dong et al. (2015) Dong R., Zhu Z., Whitney B., 2015, ApJ, 809, 93
  * Dra̧żkowska & Dullemond (2014) Dra̧żkowska J., Dullemond C. P., 2014, A&A, 572, A78
  * Dra̧żkowska et al. (2016) Dra̧żkowska J., Alibert Y., Moore B., 2016, A&A, 594, A105
  * Dubrulle et al. (1995) Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237
  * Epstein (1924) Epstein P. S., 1924, Physical Review, 23, 710
  * Estrada et al. (2016) Estrada P. R., Cuzzi J. N., Morgan D. A., 2016, ApJ, 818, 200
  * Fromang & Nelson (2009) Fromang S., Nelson R. P., 2009, A&A, 496, 597
  * Garaud et al. (2004) Garaud P., Barrière-Fouchet L., Lin D. N. C., 2004, ApJ, 603, 292
  * Goldreich & Lynden-Bell (1965) Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125
  * Gonzalez et al. (2015) Gonzalez J.-F., Laibe G., Maddison S. T., Pinte C., Ménard F., 2015, MNRAS, 454, L36
  * Gonzalez et al. (2017) Gonzalez J.-F., Laibe G., Maddison S. T., 2017, MNRAS, 467, 1984
  * Goodman & Pindor (2000) Goodman J., Pindor B., 2000, Icarus, 148, 537
  * Güttler et al. (2010) Güttler C., Blum J., Zsom A., Ormel C. W., Dullemond C. P., 2010, A&A, 513, A56
  * Haghighipour (2005) Haghighipour N., 2005, MNRAS, 362, 1015
  * Jacquet et al. (2011) Jacquet E., Balbus S., Latter H., 2011, MNRAS, 415, 3591
  * Johansen & Youdin (2007) Johansen A., Youdin A., 2007, ApJ, 662, 627
  * Johansen et al. (2007) Johansen A., Oishi J. S., Low M., Klahr H., Henning T., Youdin A., 2007, Nature, 448, 1022
  * Johansen et al. (2009) Johansen A., Youdin A., Mac Low M.-M., 2009, ApJ, 704, L75
  * Johansen et al. (2012) Johansen A., Youdin A. N., Lithwick Y., 2012, A&A, 537, A125
  * Johansen et al. (2014) Johansen A., Blum J., Tanaka H., Ormel C., Bizzarro M., Rickman H., 2014, Protostars and Planets VI, pp 547–570
  * Johansen et al. (2015) Johansen A., Mac Low M.-M., Lacerda P., Bizzarro M., 2015, Science Advances, 1, 1500109
  * Kowalik et al. (2013) Kowalik K., Hanasz M., Wóltański D., Gawryszczak A., 2013, MNRAS, 434, 1460
  * Laibe & Price (2014) Laibe G., Price D. J., 2014, MNRAS, 444, 1940
  * Laibe et al. (2012) Laibe G., Gonzalez J.-F., Maddison S. T., 2012, A&A, 537, A61
  * Lyra & Kuchner (2013) Lyra W., Kuchner M., 2013, Nature, 499, 184
  * Meheut et al. (2015) Meheut H., Fromang S., Lesur G., Joos M., Longaretti P.-Y., 2015, A&A, 579, A117
  * Miniati (2010) Miniati F., 2010, Journal of Computational Physics, 229, 3916
  * Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus3, 67, 375
  * Okuzumi et al. (2016) Okuzumi S., Momose M., Sirono S.-i., Kobayashi H., Tanaka H., 2016, ApJ, 821, 82
  * Picogna & Kley (2015) Picogna G., Kley W., 2015, A&A, 584, A110
  * Pinte & Laibe (2014) Pinte C., Laibe G., 2014, A&A, 565, A129
  * Raettig et al. (2015) Raettig N., Klahr H., Lyra W., 2015, ApJ, 804, 35
  * Rosotti et al. (2016) Rosotti G. P., Juhasz A., Booth R. A., Clarke C. J., 2016, MNRAS, 459, 2790
  * Ruge et al. (2016) Ruge J. P., Flock M., Wolf S., Dzyurkevich N., Fromang S., Henning T., Klahr H., Meheut H., 2016, A&A, 590, A17
  * Saffman (1962) Saffman P. G., 1962, Journal of Fluid Mechanics, 13, 120
  * Schäfer et al. (2017) Schäfer U., Yang C.-C., Johansen A., 2017, A&A, 597, A69
  * 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. (2016) Simon J. B., Armitage P. J., Li R., Youdin A. N., 2016, ApJ, 822, 55
  * Taki et al. (2016) Taki T., Fujimoto M., Ida S., 2016, A&A, 591, A86
  * Tilley et al. (2010) Tilley D. A., Balsara D. S., Brittain S. D., Rettig T., 2010, MNRAS, 403, 211
  * Williams & Best (2014) Williams J. P., Best W. M. J., 2014, ApJ, 788, 59
  * Yang & Johansen (2014) Yang C.-C., Johansen A., 2014, ApJ, 792, 86
  * Yang & Menou (2010) Yang C.-C., Menou K., 2010, MNRAS, 402, 2436
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
  * Youdin & Johansen (2007) Youdin A., Johansen A., 2007, ApJ, 662, 613
  * Zhang et al. (2015) Zhang K., Blake G. A., Bergin E. A., 2015, ApJ, 806, L7
  * Zsom et al. (2010) Zsom A., Ormel C. W., Güttler C., Blum J., Dullemond C. P., 2010, A&A, 513, A57

##  Appendix A Parameters of the linear tests

The sets of parameters for the linA, linB, linC and linD problems are provided
in Tab. 1.

Table 1: Parameters of the linA, linB, linC and linD problems, where $\eta_{0}r_{0}\Omega_{0}/c_{\mathrm{s}}=0.05$. Problem | $\tilde{k}_{x}$ | $\tilde{k}_{z}$ | $\tau_{\mathrm{s}}$ | $\epsilon$ | $\tilde{s}$
---|---|---|---|---|---
linA | 30 | 30 | 0.1 | 3 | 0.4190204
linB | 6 | 6 | 0.1 | 0.2 | 0.0154764
linC | $1\,500$ | $1\,500$ | 0.01 | 2 | 0.5980690
linD | $2\,000$ | $2\,000$ | 0.001 | 2 | 0.3154373

##  Appendix B Solving for background velocities

Solving Eqs. 12 – 19 by elimination leads to a 256th-order polynomial equation
in one of the variables (e.g. $a_{\mathrm{g}}$), making the numerical
resolution difficult. Simplifications can be made up to almost obtain an
analytical solution. Combining Eq. 12 $\pm$ Eq. 16, Eq. 13 $\pm$ Eq. 17, Eq.
14 $\pm$ Eq. 18, Eq. 15 $\pm$ Eq. 19, we obtain

$\displaystyle
a^{+}a^{-}-2\Omega_{0}\alpha^{-}+\gamma^{+}a^{-}+\Gamma\Omega_{0}^{2}=0\,,$
(31)
$\displaystyle\dfrac{1}{2}a^{+}b^{-}+\dfrac{1}{2}a^{-}b^{+}-2\Omega_{0}\beta^{-}+\gamma^{+}b^{-}+2r_{0}\eta\Omega_{0}^{2}=0\,,$
(32)
$\displaystyle\dfrac{1}{2}a^{+}\alpha^{-}+\dfrac{1}{2}a^{-}\alpha^{+}+2\Omega_{0}a^{-}+\gamma^{+}\alpha^{-}=0\,,$
(33)
$\displaystyle\dfrac{1}{2}b^{+}\alpha^{-}+\dfrac{1}{2}b^{-}\alpha^{+}+2\Omega_{0}b^{-}+\gamma^{+}\beta^{-}=0\,,$
(34)
$\displaystyle\dfrac{1}{2}(a^{+})^{2}+\dfrac{1}{2}(a^{-})^{2}-6\Omega_{0}^{2}-2\Omega_{0}\alpha^{+}+\gamma^{-}a^{-}-\Gamma\Omega_{0}^{2}=0\,,$
(35)
$\displaystyle\dfrac{1}{2}a^{+}b^{+}+\dfrac{1}{2}a^{-}b^{-}-2\Omega_{0}\beta^{+}+\gamma^{-}b^{-}-2r_{0}\eta\Omega_{0}^{2}=0\,,$
(36)
$\displaystyle\dfrac{1}{2}a^{+}\alpha^{+}+\dfrac{1}{2}a^{-}\alpha^{-}+2\Omega_{0}a^{+}+\gamma^{-}\alpha^{-}=0\,,$
(37)
$\displaystyle\dfrac{1}{2}b^{+}\alpha^{+}+\dfrac{1}{2}b^{-}\alpha^{-}+2\Omega_{0}b^{+}+\gamma^{-}\beta^{-}=0\,,$
(38)

where

$\begin{matrix}a^{\pm}\equiv a_{\mathrm{p}}\pm
a_{\mathrm{g}}\,,&&b^{\pm}\equiv b_{\mathrm{p}}\pm b_{\mathrm{g}}\,,\\\
\alpha^{\pm}\equiv\alpha_{\mathrm{p}}\pm\alpha_{\mathrm{g}}\,,&&\beta^{\pm}\equiv\beta_{\mathrm{p}}\pm\beta_{\mathrm{g}}\,,\\\
\end{matrix}$ (39)

and $\gamma^{\pm}\equiv\dfrac{1\pm\epsilon}{t_{\mathrm{stop}}}$. Eqs. 31, 33,
35 and 37 contain only $a^{\pm}$ and $\alpha^{\pm}$ and can be combined to
obtain $a^{+}$ as a root of a much simpler 5th-order polynomial. No analytic
root exists for such a polynomial. However, its single real root can easily be
found numerically. From there, the determination of all the other quantities
is algebraic, except choosing one root between the two of a second order
polynomial. Yet, one of the root is unphysical and can be easily discarded
when looking at the numerical values (e.g. $V_{\mathrm{p}}\simeq
10^{3}r_{0}\Omega_{0}$). In the case $\eta=0$ (pure maximum), we obtain
analytically $b_{\rm g,p}=\beta_{\rm g,p}=0$, i.e. the stationary velocities
are centred around $x=0$, as expected. The non-linearity of the system comes
from non-linear advection terms, even in the steady state. Although
alternative functional forms may solve the equations of motion, we restrain
ourselves to a linear form to remain consistent with the shearing box
approximation. In absence of any pressure maximum ($\Gamma=0$), we obtain the
expressions given in Youdin & Johansen (2007)

$\displaystyle\overline{U}_{\mathrm{g}}=\dfrac{2\epsilon\tau_{\mathrm{s}}\eta
r_{0}\Omega_{0}}{\tau_{\mathrm{s}}^{2}+(1+\epsilon)^{2}}\,,$ (40)
$\displaystyle\overline{V}_{\mathrm{g}}=-\dfrac{3}{2}\Omega_{0}x-\eta
r_{0}\Omega_{0}+\dfrac{\epsilon(1+\epsilon)\eta
r_{0}\Omega_{0}}{\tau_{\mathrm{s}}^{2}+(1+\epsilon^{2})}\,,$ (41)
$\displaystyle\overline{U}_{\mathrm{p}}=\dfrac{-2\eta\tau_{\mathrm{s}}r_{0}\Omega_{0}}{\tau_{\mathrm{s}}^{2}+(1+\epsilon)^{2}}\,,$
(42)
$\displaystyle\overline{V}_{\mathrm{p}}=-\dfrac{3}{2}\Omega_{0}x-\dfrac{(1+\epsilon)\eta
r_{0}\Omega_{0}}{\tau_{\mathrm{s}}^{2}+(1+\epsilon)^{2}}\,.$ (43)