© ESO, 2012

1. Introduction

Much of the information about the gas structure of protoplanetary discs is inferred from the emission by the dust component and an assumed dust-to-gas ratio. Interpretations of recent observations in the (sub)millimetre domain (Andrews & Williams 2005, 2007; Lommen et al. 2007) show that observed discs typically have masses between 10^-4 and 10^-1  M_⊙ and a spatial extent of a few hundred AU. Their radial surface density and temperature profiles are approximated by power laws (Σ ∝ r ^− p, ), whose respective exponents p and q have positive values typically of order unity.

Seminal studies describe the dust motion in protoplanetary discs, which depends strongly on the gas structure. Weidenschilling (1977a, hereafter W77) and Nakagawa et al. (1986, hereafter NSH86) demonstrated that dust grains from micron sizes to pre-planetesimals (a few metres in size) experience a radial motion through protoplanetary discs. This motion is called radial drift or migration. Due to its pressure gradient, the gas orbits the central star at a sub-Keplerian velocity. Grains therefore have a differential velocity with respect to the gas. The ensuing drag transfers linear and angular momentum from the dust to the gas. Thus, dust particles can not sustain the Keplerian motion they would have without the presence of gas and as a result migrate toward the central star.

This migration motion depends strongly on the grain size, which sets the magnitude of the drag, as well as the nature of the drag regime. Specifically, as shown by W77 and NSH86, grains of a critical size pass through the disc in a fraction of the disc lifetime. This catastrophic outcome is called the “radial-drift barrier” of planet formation. More precisely, we will adopt the subsequent definition for the “radial-drift barrier” in this study: “the ability of grains of be accreted onto the central star/depleted from the disc within its lifetime”. Historically, this process was first studied in a minimum mass solar nebula (MMSN, see Weidenschilling 1977b; Hayashi 1981; Desch 2007; Crida 2009), in which the critical size corresponds to metre-sized bodies and thus was called the “metre-size barrier”. However, planets are frequently observed (besides the 8 planets in our solar system, more than 700 extra-solar planets have been discovered to date¹): some solid material must therefore have overcome this barrier and stayed in the disc to form larger bodies. Moreover, if the small grains of every disc were submitted to the radial-drift barrier, we would barely detect them since their emission via optical/IR scattering and IR thermal radiation is due to small grains. As discs are frequently observed, the grains from which the emission is detected cannot be strongly depleted for a substantial fraction of discs.

From a theoretical point of view, such a discrepancy between the observations and the theoretical predictions imply that the seminal theory has to be extended (some physical element is lacking) or that it has not been fully exploited. This second option has been investigated by Youdin & Shu (2002, hereafter YS02). They highlight the fact that, contrary to the primary hypothesis of W77, observed dusty discs are drastically different from the MMSN prototype. As the radial surface density and temperature profiles fix both the radial pressure gradient and the magnitude of the gas drag, different values for the power-law exponents p and q affect the optimal grain size of migration and thus induce different radial motions for the dust through the disc. Specifically, YS02 showed that for steep surface density profiles and smooth temperature profiles, the grains radial velocity decreases when the grains reach the inner discs regions. Grains in such discs therefore experience a “pile-up”. However, while important, the work of YS02 does not provide a precise conclusion on the outcome of the grains nor any quantitative criterion for the “pile-up” process to be efficient enough to avoid the radial-drift barrier. Furthermore, YS02 restricts their study to the special case of a gas phase with a low density (e.g. the grain size smaller than the gas mean free path, called the Epstein regime). This hypothesis is not valid anymore when considering the radial drift of pre-planetesimals, whose grain sizes are larger than the gas mean free path and are submitted to the Stokes drag regime. Although the radial drift of pre-planetesimals has already been studied in different situations with numerical or semi-analytical methods – see e.g. Haghighipour & Boss (2003); Birnstiel et al. (2009); Youdin (2011) – its rigorous theory for the standard case of a simple disc has not yet been derived.

Within this context, we see that (i) the seminal theory describing the radial motion of dust grains has been developed within the limits of the Epstein regime but does not treat the Stokes regime; (ii) there exists no clear theoretical criterion to predict the impact of the “pile-up effect” on the outcome of the dust radial motion; (iii) there exists no criterion to predict whether a given disc will be submitted to the “radial-drift barrier” phenomenon. To answer these three points, we re-visit in this study the work of W77 and NSH86 and extend the developments of YS02 for both the Epstein and the Stokes regime. Performing rigorous perturbative expansions, we find two theoretical criteria (one for each regime) which predict when the “pile-up” effect is sufficient for the grains not to be accreted onto the central star. We then test when these theoretical criteria can be applied in real discs.

Additionally, our work is motivated by the recent observational results of Ricci et al. (2010a,b). From their observations they claim that “a mechanism halting or slowing down the inward radial drift of solid particles is required to explain the data”. In this work we aim to show that contrary to what is usually invoked, local pressure maxima due to turbulent vortices or spiral density waves may help but are not necessarily required to explain the observations. Ricci et al. (2010a) also mention that “the observed flux of the fainter discs are instead typically overpredicted even by more than an order of magnitude”. Here, we also aim to provide a quantitative criterion to determine which discs are faint and which one are not. Thus, revisiting the seminal theory of the radial drift is timely, all the more so than an important quantity of new data is soon to be provided by ALMA, the Atacama Large Millimeter/submillimeter Array.

In this paper, we first recall some general properties of grain motion in protoplanetary discs for both the Epstein and Stokes regime in Sect. 2. We then focus on the radial motion of non-growing grains in the Epstein regime. We expand the radial motion equations assuming a weak pressure gradient in Sect. 3 and detail the two different modes of migration which grains may experience in Sects. 3.1 and 3.2. This allows us to derive an analytic criterion which determines the asymptotic dust behaviour in the Epstein regime in Sect. 3.3. We transpose these derivations for the Stokes regime at low Reynolds number in Sect. 4 and obtain a similar criterion for this regime. We also discuss the grains outcome for large Reynolds numbers. In Sect. 5, we discuss the relevance of these criteria and study their implications for observed protoplanetary discs and planet formation in Sect. 6. Our conclusions are presented in Sect. 7.

2. Dynamics of dust grains

To reduce the parameter space for this study, we assume the following:

• 1.

  The disc is a thin, non-magnetic, non-self-graviting, inviscid perfect gas disc which is vertically isothermal. Its radial surface density and temperature are described by power-law profiles. Notations are described in Appendix A. The flow is laminar and in stationary equilibrium. Consequently, the gas velocity and density are described by well-known relations, which we present in Appendix B.

• 2.

  The grains are compact homogeneous spheres of fixed radius. The collisions between grains and the collective effects due to large dust concentrations are neglected. When the grains are small compared to the mean free path of the gas (λ_g > 4s/9, where s is the grain size), their interactions with the gas are treated by the Epstein drag force for diluted media (Epstein 1924; Baines et al. 1965; Stepinski & Valageas 1996). This drag is caused by the transfer of momentum by individual collisions with gas molecules at the grains surface. Assuming specular reflections on the grain and when the differential velocity with the gas is negligible compared to the gas sound speed, the now common expression of the drag force is $\{\begin{array}{c}\\ {F}\mathrm{D}& \mathrm{=}& \mathrm{-}\frac{{\mathit{m}}_{\mathrm{d}}}{{\mathit{t}}_{\mathrm{s}}}\hspace{0.17em}\mathrm{\Delta }{v}\\ {\mathit{t}}_{\mathrm{s}}& \mathrm{=}& \frac{{\mathit{\rho }}_{\mathrm{d}}\mathit{s}}{{\mathit{\rho }}_{\mathrm{g}}{\mathit{c}}_{\mathrm{s}}}\hspace{0.17em}\mathit{,}\end{array}$(1) where m_d is the dust grain’s mass, t_s the stopping time, ρ_g the gas density, c_s the local gas sound speed, ρ_d the intrinsic dust density, and Δv = v − v_g the differential velocity between dust and the mean gas motion. In classical T Tauri star (CTTS) protoplanetary discs, drag forces for particles smaller than  ~10 m are well described by the Epstein regime (Garaud et al. 2004, see also Sect. 6.1). Small grains which produce the emission of observed protoplanetary discs satisfy this criterion. The interactions between large dust particles (λ_g < 4s/9) and the gas are treated by the Stokes drag force (Whipple 1972; Stepinski & Valageas 1996). In this case, the gas mean free path is small and the dust particle is locally surrounded by a viscous fluid. Depending on the local Reynolds number of the flow around the grains ${\mathit{R}}_{\mathrm{g}} \mathrm{=} \frac{\mathrm{2}\mathit{s}\mathrm{\left|}\mathrm{\Delta }{v}\mathrm{\right|}}{\mathit{\nu }}$, where ν is the microscopic kinematic viscosity of the gas, the drag force takes the following expression: ${F}\mathrm{D}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{C}}_{\mathrm{D}}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{g}}\left|\mathrm{\Delta }{v}\right|\mathrm{\Delta }{v}\mathit{,}$(2)where the drag coefficient C_D is given by ${\mathit{C}}_{\mathrm{D}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{24}{\mathit{R}}_{\mathrm{g}}^{-1}& \mathrm{for}{\mathit{R}}_{\mathrm{g}}\mathit{<}\mathrm{1}\\ \mathrm{24}{\mathit{R}}_{\mathrm{g}}^{-0.6}& \mathrm{for}\mathrm{1}\mathit{<}{\mathit{R}}_{\mathrm{g}}\mathit{<}\mathrm{800}\\ \mathrm{0.44}& \mathrm{for}\mathrm{800}\mathit{<}{\mathit{R}}_{\mathrm{g}}\hspace{0.17em}\mathit{.}\end{array}$(3)If R_g < 1, the drag force remains linear in Δv.

In this work, the physical relations are written in cylindrical coordinates (r, θ, z). The related unit vector system is given by $\left({e}{\mathit{r}}_{}\mathit{,}{e}{\mathit{\theta }}_{}\mathit{,}{e}{\mathit{z}}_{}\right)$. As the system is invariant by rotation around the vertical axis e_z, the physical quantities depend only on r and z. The physical quantities of the gas, designated by subscript g, are first determined in a general way. Then, the limit z = 0 is taken to study the restricted radial motion.

Dust dynamics depends on both the magnitude of the drag (driven by the differential velocity) and on its relative contribution with respect to the gravity of the central star. Seminal studies of dust dynamics were conducted by Whipple (1972), W77, Weidenschilling (1980) and NSH86, and extended by others (YS02; Takeuchi & Lin 2002; Haghighipour & Boss 2003; Garaud et al. 2004; Youdin & Chiang 2004). Here we recall the major points of those studies. We consider two forces acting on the grain: the gravity of the central star and gas drag. (We assume that the momentum transferred by drag from a single grain on the gas phase is negligible.) Thus ${\mathit{m}}_{\mathrm{d}}\frac{\mathrm{d}{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}{F}\mathrm{D}\mathrm{+}{\mathit{m}}_{\mathrm{d}}{g}\mathit{,}$(4)where F_D is the drag force. As shown by Eqs. (1)–(2), a general expression of the ratio $\frac{{F}{\mathrm{D}}_{}}{{\mathit{m}}_{\mathrm{d}}}$ is of the form $\frac{{F}\mathrm{D}}{{\mathit{m}}_{\mathrm{d}}}\mathrm{=}\mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ \mathrm{𝒞}\end{array}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}}{{\mathit{s}}^{\mathit{y}}}\mathrm{|}{v}\mathrm{-}{v}\mathrm{g}{\mathrm{|}}^{\mathit{\lambda }}\mathrm{(}{v}\mathrm{-}{v}{\mathrm{g}}^{\mathrm{)}}\mathit{,}$(5)where the quantities $\begin{array}{}{\mathrm{˜}}_{}\\ \mathrm{𝒞}\end{array}$, y and λ are defined for both the Epstein and the Stokes regime in Appendix C.

3. Radial motion in the Epstein regime: perturbation analysis at small pressure gradients

Considering the Epstein (small grains) regime, Eq. (4) reduces to: ${\mathit{m}}_{\mathrm{d}}\frac{\mathrm{d}{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{{\mathit{m}}_{\mathrm{d}}}{{\mathit{t}}_{\mathrm{s}}}\mathrm{(}{v}\mathrm{-}{v}{\mathrm{g}}^{\mathrm{)}}\mathrm{+}{\mathit{m}}_{\mathrm{d}}{g}\mathit{.}$(6)Writing Eq. (6) in ${\mathrm{(}}^{}\mathit{r,\theta ,}{\mathit{z}}^{\mathrm{)}}$ coordinates leads to $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$(7)To highlight the important parameters involved in the grains dynamics, we introduce dimensionless quantities (see Appendix C). It is crucial to note that the ratio $\frac{{\mathit{t}}_{\mathrm{s}}}{{\mathit{t}}_{\mathrm{k}}}$ of the two timescales related to the physical processes acting on the grain is given by $\frac{{\mathit{t}}_{\mathrm{s}}}{{\mathit{t}}_{\mathrm{k}}}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{s}\mathrm{0}}}{{\mathit{t}}_{\mathrm{k}\mathrm{0}}}{\mathit{R}}^{\mathit{p}}{\mathrm{e}}^{\frac{{\mathit{Z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}^{\mathrm{3}\mathrm{-}\mathit{q}}}}\mathrm{=}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}}{\mathrm{e}}^{\frac{{\mathit{Z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}^{\mathrm{3}\mathrm{-}\mathit{q}}}}\mathit{.}$(8)With Eq. (1) and noting Ω_k the Keplerian angular velocity, this ratio can be written as $\frac{{\mathit{t}}_{\mathrm{s}}}{{\mathit{t}}_{\mathrm{k}}}\mathrm{=}\frac{\mathit{s}}{\left(\frac{{\mathit{\rho }}_{\mathrm{g}}{\mathit{c}}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{d}}{\mathrm{\Omega }}_{\mathrm{k}}}\right)}\mathrm{=}\frac{\mathit{s}}{{\mathit{s}}_{\mathrm{opt}}}\mathrm{=}\mathit{S,}$(9)where ${\mathit{s}}_{\mathrm{opt}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{g}}{\mathit{c}}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{d}}{\mathrm{\Omega }}_{\mathrm{k}}}$. This timescale ratio therefore corresponds to a dimensionless size $\mathit{S}\mathrm{=}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}}{\mathit{e}}^{\frac{{\mathit{Z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{R}}^{\mathrm{3}\mathrm{-}\mathit{q}}}}$ for the grain. If S ≪ 1 (resp. S ≫ 1), the effects of drag will occur much faster (resp. slower) than gravitational effects. If S ≃ 1, both gravity and drag will act on the same timescale. Interestingly, s_opt varies in the disc midplane as r ^− p, as does surface density.

Then using Eqs. (8), (B.2) and (C.1), we obtain for (e_r,e_θ,e_z): $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$

These equations depend on five control parameters: the initial dimensionless grain size, S₀, the radial surface density and temperature exponents, p and q, the square of the disc aspect ratio, ${\mathit{\phi }}_{\mathrm{0}}\mathrm{=}{{}^{\mathrm{(}}\mathit{H}_{\mathrm{0}}\mathit{/}{\mathit{r}}_{\mathrm{0}}{}^{\mathrm{)}}}^{\mathrm{2}}$, and the subkeplerian parameter, η₀, given by Eq. (C.2). The equations can be simplified in some cases, e.g. if the vertical motion is considered to occur faster than the radial motion, R ≃ 1 and $\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{Z}}{\mathrm{d}{\mathit{T}}^{\mathrm{2}}}$ simplifies to the damped harmonic oscillator equation. If we consider only the radial motion (for a 2D disc), we have Z = 0, and $\begin{array}{ccc}\{\begin{array}{c}\\ \\ \end{array}& & \end{array}$(11)Even for discs in two dimensions, Eq. (11) is not analytically tractable. However, as some of the parameters involved in the equation are small, approximations of the solution can be found by performing perturbative expansions. Some of the classical results detailed below have been studied in W77 and NSH86, but are here properly justified. The principle of those expansions is described on Fig. 1.

At t = 0, r = r₀ which implies that R⁽T = 0⁾ = 1. Because of gas drag, a grain dissipates both its energy and angular momentum and therefore, experiences a radial inward motion, i.e. R < 1. The first parameter with respect to which a perturbative expansion can be performed is η₀ (linked to the pressure gradient by Eq. (C.3)) as it takes values of approximately 10^-3–10^-2 in real protoplanetary discs (see NSH86), and thus η₀ ≪ 1. We consider that this inequality also implies that ${\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}\mathrm{\ll }\frac{\mathrm{1}}{\mathit{R}}\mathrm{·}$(12)This inequality is always verified when q ≤ 1 and thus applies to observed discs (see Sect. 6). For q > 1, there is a region where this inequality is not verified. However, in this case, the pressure gradient has the same order of magnitude as the gravity of the central star and the model of a power-law profile for the radial temperature is not accurate enough to model realistic discs. We thus consider that for real discs Eq. (12) is always justified. Then, following NSH86, we consider the system of equations given by Eq. (11). We set $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{+}{\mathit{\eta }}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{+}{\mathit{\eta }}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\end{array}$(13)and look at the orders $\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$, ,... of the expansion – see Appendix E. We find that: $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}& \mathrm{=}& {\mathit{\eta }}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\\ & \mathrm{=}& \end{array}$(14)The pressure gradient term has been retained to keep the generality, however since we assume that η₀ ≪ 1, we equivalently have $\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}\mathrm{-}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{=}\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{2}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}$(15)Thus, to order , $\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}$(16)or equivalently, using Eqs. (C.4) and (C.1), ${\mathit{v}}_{\mathit{r}}\mathrm{=}\frac{\mathit{r}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{k}}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\frac{\mathrm{(}{\mathit{t}}_{\mathrm{s}}\mathit{/}{\mathit{t}}_{\mathrm{k}}\mathrm{)}}{\mathrm{1}\mathrm{+}{\mathrm{(}{\mathit{t}}_{\mathrm{s}}\mathit{/}{\mathit{t}}_{\mathrm{k}}\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(17)S₀R^p is the dimensionless expression of the ratio t_s/t_k. Eq. (16) shows that ${\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\ll }\mathrm{1}$ or ${\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\gg }\mathrm{1}$, and thus R^pS₀ ≪ 1 or R^pS₀ ≫ 1, resulting in asymptotic behaviours for the radial grain motion. These asymptotic regimes were first described by W77 for the particular case p = 0. They correspond physically to two limiting cases: where the gas drag dominates, which we call the A-mode, and where gravity dominates, which we call the B-mode. In the next sections, we study and describe these two so-called “regimes of migration” or “modes of migration” before treating the global evolution of grains given by Eq. (16).

3.1. A-mode (radial differential migration)

The A-mode corresponds to the regime R^pS₀ ≪ 1 (or equivalently t_s/t_k ≪ 1). In the A-mode, Eq. (16) reduces to $\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{T}}\mathrm{=}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}$(18)or equivalently ${\mathit{v}}_{\mathit{r}}\mathrm{=}\frac{\mathit{r}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{k}}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\frac{{\mathit{t}}_{\mathrm{s}}}{{\mathit{t}}_{\mathrm{k}}}\mathit{,}$(19)where the $\mathrm{𝒪}{}^{\mathrm{(}}\mathit{R}^{\mathit{p}}{\mathit{S}}_{\mathrm{0}}{}^{\mathrm{)}}$ has been neglected. In this mode of migration, the stopping time is much smaller than the Keplerian time scale. Considering one grain’s orbit around the central star, its orbital velocity is forced by the gas drag to become sub-Keplerian in just a few stopping times, i.e. almost instantaneously. Thus, the centrifugal acceleration is not efficient enough to counterbalance the gravitational attraction of the central star, and the grain feels an inward radial differential acceleration. The gas drag counterbalances this radial motion and the grain reaches a local limit velocity in a few stopping times. We call the physical process of the A-mode of migration “Radial Differential Migration”.

The A-mode of migration originates first from a perturbative expansion for η₀ ≪ 1 (rigorously for η₀R ^− q + 1 ≪ 1) and second from a perturbative expansion for S₀ ≪ 1 (rigorously for S₀R^p ≪ 1). Formally, we have performed: $\begin{array}{}\lim \\ {\mathit{S}}_{\mathrm{0}}\mathrm{\ll }\mathrm{1}\end{array}\begin{array}{}\lim \\ {\mathit{\eta }}_{\mathrm{0}}\mathrm{\ll }\mathrm{1}\end{array}\left[\mathit{...}\right]$. Historically, the A-mode had been derived by W77 to explain the radial motion of small grains. In his study, he neglected the radial dependence of the stopping time and assumed S₀ ≪ 1 (this approximation also implies that ${\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\ll }\mathrm{1}$, as R < 1, see Appendix F).

It is straightforward to integrate the differential equation Eq. (18) by separating the R and T variables. Noting that R⁽T = 0⁾ = 1, we have:

• If $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\ne \mathrm{0}$: $\{\begin{array}{c}\\ \mathit{R}\mathrm{=}{\left[\mathrm{1}\mathrm{-}\left(\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right){\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}\mathit{T}\right]}^{\frac{\mathrm{1}}{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}\\ \mathit{T}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{(}\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}}\mathrm{·}\end{array}$(20)

• If $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0}$: $\{\begin{array}{c}\\ \mathit{R}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}\mathit{T}}\\ \mathit{T}\mathrm{=}\mathrm{-}\frac{\ln \mathrm{(}{\mathit{R}}^{\mathrm{)}}}{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}}\mathrm{·}\end{array}$(21)

The outcome of the dust radial motion comes from a competition between two effects. As the grain reaches smaller radii, (1) gas drag increases, which slows down the radial drift, and (2) the differential acceleration due to the pressure gradient increases which enhances the migration efficiency. Point (1) is related to s_opt, which scales as the surface density profile, while the acceleration due to the pressure gradient in (2) is related to the temperature profile (see Eq. (C.1)). Depending on which process is dominant, the grain’s dynamics can lead to two regimes:

• If $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$: then as $\mathit{T}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{~}\frac{\mathrm{-}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{(}\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}}\mathrm{=}\mathrm{𝒪}\left({\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathit{,}$(22)the time it takes the grain to reach smaller and smaller radii increases drastically, according to the diverging power-law. Importantly, this behaviour constitutes our definition of the grain “pile-up”. Mathematically speaking, accretion onto the central star occurs in an infinite time, i.e. $\begin{array}{}\lim \\ \mathit{T}\mathrm{\to }\mathrm{+}\mathrm{\infty }\end{array}\mathit{R}\mathrm{=}\mathrm{0.}$(23)

• If $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$: the grain is accreted onto the central star in a finite migration time given by ${\mathit{T}}_{\mathrm{m}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}\mathrm{(}\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}$(24)which increases as S₀ and η₀ decrease, so that $\begin{array}{}\lim \\ \mathit{T}\mathrm{\to }{\mathit{T}}_{\mathrm{m}}\end{array}\mathit{R}\mathrm{=}\mathrm{0.}$(25)

The presence or absence of a physical grain pile-up is therefore demonstrated considering the asymptotic behaviour of R(T) at large times. It is important to realise that the pile-up is a cumulative effect that can not arise from velocities only (which however provide qualitative information on the grain’s motion) but can only be found by integrating the equation of motion. This rigorously allows us to distinguish two different behaviours for the outcome of the grain’s radial motion, and thus two classes of discs with respect to the A-mode.

3.2. B-mode (drift forced by a resistive torque)

Returning to Eq. (16), the B-mode corresponds to the other asymptotic regime, where R^pS₀ ≫ 1 (or equivalently t_s/t_k ≫ 1). In this case, $\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{T}}\mathrm{=}\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}$(26)or equivalently ${\mathit{v}}_{\mathit{r}}\mathrm{=}\frac{\mathit{r}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{k}}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\frac{{\mathit{t}}_{\mathrm{k}}}{{\mathit{t}}_{\mathrm{s}}}\mathrm{·}$(27)In this mode of migration, the stopping time is much larger than the Keplerian time scale. Hence, the orbital velocity of a grain around the central star is almost the Keplerian velocity. However, because of the pressure gradient, the gas orbits around the central star at a sub-Keplerian velocity. Thus, the azimuthal differential velocity between the gas and the grain generates an azimuthal drag force whose torque extracts angular momentum from the grain. Given that the Keplerian angular momentum increases with radius ($\mathit{l}\mathrm{\propto }\sqrt{\mathit{r}}$), this torque results in the inward migration of the grain. We call the physical process of this B-mode of migration “drift forced by a resistive torque”.

As for the A-mode, the B-mode of migration can also be derived first from an expansion in ${{}^{\mathrm{(}}\mathit{S}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}}{\mathrm{)}}^{}}^{-2}$ and then from an expansion in η₀. Historically, W77 found an expression while only assuming that S₀ ≫ 1 since he considered a flat density profile. To find the expression derived by W77 for large grains, we must assume that S₀R^p ≫ 1. It is crucial to see that this expression does not imply that S₀ ≫ 1 when p > 0 and R → 0.

It is straightforward to integrate the differential equation Eq. (26) by separating the R and T variables. Noting that R⁽T = 0⁾ = 1, we have:

• If $\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\ne \mathrm{0}$: $\{\begin{array}{c}\\ \mathit{R}\mathrm{=}{\left[\mathrm{1}\mathrm{-}\left(\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right)\frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{S}}_{\mathrm{0}}}\mathit{T}\right]}^{\frac{\mathrm{1}}{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}\\ \mathit{T}\mathrm{=}\frac{{\mathit{S}}_{\mathrm{0}}}{{\mathit{\eta }}_{\mathrm{0}}}\frac{\mathrm{1}\mathrm{-}{\mathit{R}}^{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}}\mathit{.}\end{array}$(28)

• If $\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0}$: $\{\begin{array}{c}\\ \mathit{R}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{S}}_{\mathrm{0}}}\mathit{T}}\\ \mathit{T}\mathrm{=}\mathrm{-}\frac{{\mathit{S}}_{\mathrm{0}}}{{\mathit{\eta }}_{\mathrm{0}}}\ln \mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathit{.}\end{array}$(29)

As for the A-mode, two kinds of behaviour appear, depending on the p and q exponents:

• If $\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$: the grain migrates inwards, piles-up in the disc’s inner regions and falls onto the star in an infinite time: $\begin{array}{}\lim \\ \mathit{T}\mathrm{\to }\mathrm{+}\mathrm{\infty }\end{array}\mathit{R}\mathrm{=}\mathrm{0.}$(30)However, the negative exponents required to be in this regime do not correspond to physical discs. Therefore, the grain dynamics in the B-mode in real discs belong to the second case:

• If $\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$: the grain is accreted onto the central star in a finite time ${\mathit{T}}_{\mathrm{m}}\mathrm{=}\frac{{\mathit{S}}_{\mathrm{0}}}{{\mathit{\eta }}_{\mathrm{0}}\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}$(31)which increases as S₀ increases and η₀ decreases and so that $\begin{array}{}\lim \\ \mathit{T}\mathrm{\to }{\mathit{T}}_{\mathrm{m}}\end{array}\mathit{R}\mathrm{=}\mathrm{0.}$(32)

As for the A-mode, considering the limit of R(T) at large times also proves the existence of two classes of discs with respect to the B-mode of migration. The radial motion of a grain in the B-mode of migration is also driven by a competition between the increase of both the drag and the acceleration due to the pressure gradient. However, for real discs, p > 0 and q > 0, and therefore, $\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$. Grains migrating in B-mode in such discs fall in a finite time onto the central star.

3.3. Radial evolution and asymptotic behaviour of single grains

As we have seen, the grains behaviour is divided into two asymptotic regimes, called the A-mode and the B-mode, which come from two different physical origins. However, the two criteria determining if the grains are accreted onto the central object in a finite or infinite time differ for the A- and the B-mode. It is thus crucial to determine in which mode a dust grain ends its motion to predict if the grain is ultimately accreted or not. Returning again to Eq. (16), we have $\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{T}}\mathrm{=}\frac{\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{·}$(33)We can separate the R and T variables and integrate to obtain an expression for T(R): $\mathit{T}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}}\left[{\mathit{T}}_{\mathrm{A}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{+}{\mathit{T}}_{\mathrm{B}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\right]\mathit{,}$(34)where ${\mathit{T}}_{\mathrm{A}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{1}\mathrm{-}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}& \mathrm{if}\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\ne \mathrm{0}\\ \mathrm{-}\ln \mathrm{(}{\mathit{R}}^{\mathrm{)}}& \mathrm{if}\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0}\end{array}$(35)and ${\mathit{T}}_{\mathrm{B}}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{=}\{\begin{array}{c}\\ {\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\frac{\mathrm{1}\mathrm{-}{\mathit{R}}^{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}& \mathrm{if}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\ne \mathrm{0}\\ \mathrm{-}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\ln \mathrm{(}{\mathit{R}}^{\mathrm{)}}& \mathrm{if}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0.}\end{array}$(36)Equation (34) provides the asymptotic behaviour of the grains at large times. Interestingly, as p ≥ 0 for realistic discs, the contribution of the B-mode becomes negligible when $\mathit{R}\mathrm{\ll }{\mathit{S}}_{\mathrm{0}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathit{p}}$. Hence, grains initially migrating in the A-mode stay in the A-mode, but grains initially migrating in the B-mode end their radial motion in the A-mode. This behaviour is summarized on Fig. 2 and detailed in Appendix G. This result was not predicted by W77, as he neglected the radial dependence of the stopping time. Mathematically speaking, it comes from the fact that the perturbative expansion of W77 has been performed with respect to powers of S₀ and not powers of S₀R^p.

To illustrate the radial motion of dust grains in protoplanetary discs, we numerically integrate the equations of motion for different values of the parameters η₀, S₀, p, q. We set η₀ = 10^-2 to mimic a realistic disc and vary the order of magnitude of S₀ from 10^-4 to 10² for two sets of (p,q) values. First, we choose (p = 0, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$); according to the NSH86 expansion, the grains are accreted by the central star in a time T_m. Second, we set ($\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}$, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$); the grains fall onto the central star in an infinite time from the same approximation. This set of (p,q) values is taken to mimic discs profiles that are commonly used and for which $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$ can take a positive or a negative value. Consequently, we interpret the radial motion $\mathit{R}{}^{\mathrm{\right(}}\mathit{T}^{\mathrm{\left)}}$ of dust grains plotted in Fig. 3.

• The top panel of Fig. 3 shows the results for (p = 0, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$): Grains fall onto the central star, initially in the A-mode for the small grains and in the B-mode for the large ones. The radial-drift process is long for small and large grains but is optimal for grains with S = S_m = 1 for which the accretion time is T_m = 1.6/η₀, or T_m = 160 with η₀ = 10^-2.

  Fig. 3 Radial motion $\mathit{R}{}^{\mathrm{\right(}}\mathit{T}^{\mathrm{\left)}}$ of dust grains in the Epstein regime for η0 = 10-2. S0 varies from 10-4 to 102. Top: p = 0, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$, here $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$ and the grain is accreted onto the central star in a finite time. Bottom: $\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}$, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$, here $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{<}\mathrm{0}$ and the grain piles up and is consequently accreted onto the central star in an infinite time.

• The bottom panel of Fig. 3 shows the results for ($\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}$, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$): In this case, the radial density profile is steep enough to ensure that the grains are not accreted onto the central star. To reach a given radius (for example R_f = 0.1), the optimal size is ${\mathit{S}}_{\mathrm{m}\mathit{,}\mathrm{f}}\mathrm{\simeq }\mathrm{2.9}\mathrm{=}\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$ (see Eq. (G.1)). Hence, in this case grains efficiently reach the disc inner regions without ever being accreted onto the central star. The transition from the B-mode to the A-mode (for which R ∝ T^-4 in this case) for the large grains is visible in this plot.

4. Radial motion in the Stokes regime

Radial migration of large particles occurs in the Stokes drag regime, which depends on the dynamical viscosity μ of the gas. For hydrogen molecules: $\mathit{\mu }\mathrm{=}\frac{\mathrm{5}\mathit{m}\sqrt{\mathit{\pi }}}{\mathrm{64}{\mathit{\sigma }}_{\mathrm{s}}}\sqrt{\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{m}}}\mathit{,}$(37)where m = 2m_H = 3.347446922 × 10^-27 kg and σ_s = 2.367 × 10^-19 m² is the molecular cross section of the molecule (Chapman & Cowling 1970). The kinematic viscosity ν is then defined by μ = ρ_gν and the gas collisional mean free path is given by ${\mathit{\lambda }}_{\mathit{g}}\mathrm{=}\sqrt{\frac{\mathit{\pi }}{\mathrm{2}}}\frac{\mathit{\nu }}{{\mathit{c}}_{\mathrm{s}}}\mathrm{·}$(38)We now generalise the procedure outlined in Sects. 2 and 3 to the three Stokes regimes of Eq. (3). Using the dimensionless coordinates described above, we have $\begin{array}{ccc}\mathit{\mu }& \mathrm{=}& \\ {\mathit{\mu }}_{\mathrm{0}}& \mathrm{=}& \end{array}$where c_s0 is given by Eq. (B.23). Thus, the expression of the kinematic viscosity ν is $\begin{array}{ccc}\mathit{\nu }& \mathrm{=}& \\ {\mathit{\nu }}_{\mathrm{0}}& \mathrm{=}& \end{array}$First, if R_g < 1, the drag force is linear in v − v_g and thus has the same structure as for the Epstein regime. Comparing the expressions of for the Epstein and the linear Stokes regime (see Appendix C), all the results found for the radial motion in the Epstein regime can therefore be directly transposed by setting q′ = q and ${\mathit{p}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{q}\mathrm{-}\mathrm{3}}{\mathrm{2}}$. In this case, the grain radial motion does not depend on p anymore and the NSH86 expansion of the radial motion for small pressure gradients provides (see Eq. (33)) $\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{T}}\mathrm{=}\frac{\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}\mathrm{-}\mathrm{1}}}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathit{q}\mathrm{-}\mathrm{3}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{4}}}\mathrm{·}$(43)These crucial results follow:

• In the A-mode (${\mathit{S}}_{\mathrm{0}}\mathrm{\ll }{\mathit{R}}^{\frac{\mathrm{3}\mathrm{-}\mathit{q}}{\mathrm{4}}}$), grains experience a pile-up and migrate onto the central star in an infinite time if $\mathrm{-}{\mathit{p}}^{\mathrm{\prime }}\mathrm{+}{\mathit{q}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$, i.e. if q ≤  − 4 (which never occurs in real discs).

• In the B-mode (${\mathit{S}}_{\mathrm{0}}\mathrm{\gg }{\mathit{R}}^{\frac{\mathrm{3}\mathrm{-}\mathit{q}}{\mathrm{4}}}$), grains migrate onto the central star in an infinite time if ${\mathit{p}}^{\mathrm{\prime }}\mathrm{+}{\mathit{q}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$, i.e. if $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$.

Thus, similar to the Epstein regime, we derived one criterion for each mode and need to determine in which mode the grain ends its motion. For observed discs, q − 3 < 0 (see Sect. 6), and as the particle migrates inward, R becomes smaller than ${\mathit{S}}_{\mathrm{0}}^{\frac{\mathrm{4}}{\mathrm{3}\mathrm{-}\mathit{q}}}$ and grains end their radial motion in the B-mode (see Fig. 4). This result is fundamentally different to the one we obtained for the Epstein regime. Indeed, for grains migrating in the A-mode in the Stokes regime at low Reynolds numbers, the criterion obtained for a pile-up in the A-mode is never satisfied for real discs. However, after migrating inside a critical radius, grains switch to the B-mode, for which the pile-up can potentially occur, depending on the value of q. The corollary is that in discs having $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$, i.e. a shallow enough temperature profile, large grains in the Stokes regime at small Reynolds numbers remain in the disc. Such a criterion is applicable for real protoplanetary discs.

Second, if R_g > 800, the drag force is quadratic in v − v_g. Assuming that the radial motion is decoupled from the vertical motion, we perform the NSH expansion at small pressure gradient (cf. Eq. (13)). We find that whatever the integer j, (η₀R^−p)^j(v − v_g) → 0 at the limit η₀ → 0. This means that both v_r and ${\mathit{v}}_{\mathit{\theta }}\mathrm{-}\sqrt{\mathrm{1}\mathit{/}\mathit{R}}$ are flat functions as their Taylor series expansion equals zero at each order. Consequently, they can not be determined by perturbation analysis. This property comes from the quadratic dependency of the drag with respect to the differential velocity and thus is not related to the grain size. Consequently, in this drag regime, the drag force is extremely efficient and the corrections to the Keplerian motion are negligible at every order of the perturbative expansion. The particles are very well coupled to the gas and do not migrate significantly.

Third, for the intermediate case, we could not manage to perform the expansion at small pressure gradients. However, we expect an intermediate behaviour between the two Stokes regime at small and large Reynolds numbers. Consequently, if R_g > 1, the migration motion becomes less efficient as the drag force is no longer linear with respect to the differential velocity between the gas and the dust particles. Thus, the main constraint for the radial-drift barrier due to the Stokes drag comes from the low Reynolds number regime for which the migration motion is the most efficient.

Finally, confusion often arises when defining the “radial-drift barrier” as the difficulty a grain has of “overcoming s = s_opt” (i.e.) reaching the B-mode. Indeed, as we have shown, grains can survive their migration motion in the Epstein regime when they are in the A-mode whenever  − p + q + 1/2 < 0, and grains can start their migration motion in the B-mode but be accreted in a finite time if  − p + q + 1/2 > 0. This study also shows for the Stokes regime that a grain ends its migration motion in the B-mode. However, as demonstrated in this work, the ability of the grain to overcome the radial drift barrier is only linked to the value of q. If q > 2/3, the grain will be accreted onto the central star in a finite time, even if it has s > s_opt. Thus, we would argue that the definition of the radial-drift barrier has to remain the ability of the grain to be accreted onto the central star or depleted from the disc within its lifetime.

5. Limitations of the model

We have demonstrated that the time it takes for grains to reach smaller and smaller radii increases dramatically under certain conditions. Specifically, if $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{<}\mathrm{0}$ (resp. $\mathit{q}\mathit{<}\frac{\mathrm{2}}{\mathrm{3}}$), grains experience a pile-up in the Epstein (resp. Stokes) regime. However, the model developed for the radial evolution of dust grains in this paper remains simple in that we neglect several important physical processes: turbulence, grain growth, collective motion of dust grains, dust feedback on the gas surface density and temperature profiles. We now discuss how those processes can modify the criteria derived above.

• 1.

  The local pressure maxima created by turbulence (Cuzzi et al. 2001, 2008) and the collective effects due to the dust drag onto the gas phase (Youdin & Goodman 2005) are known to slow down the dust particles. However, the efficiency of these processes – such as the non linearity of the streaming instability in global disc models and the life time of the pressure maxima – in real discs remains difficult to quantify. Omitting these phenomena constitutes therefore an upper limit for the grain migration efficiency, which will be slowed by these additional processes.

• 2.

  In this study, we assume that changing the dust distribution does not change the thermal profile of the disc. We also neglect the viscous evolution of the disc, assuming that the viscous timescales are larger than the characteristic timescales of the initial dust evolution. It implies that we assume that p is constant during the whole grains evolution. We can expect that for long term evolution, the surface density profile will flatten, leading to a smaller value of p. This makes the Epstein criterion harder to be met, while the Stokes criterion is not affected.

• 3.

  We have shown that even if the velocity of the grain’s inward motion depends on their sizes, their outcome only depends on the surface density and temperature profiles. Thus, if we now consider growing (or fragmenting) grains, we expect that (i) the intensity of the inward motion depends on the growth efficiency (this point will be discussed in detail in a forthcoming paper), but that (ii) the grains outcome remains determined by our criteria, regardless of the growth regime.

Following this discussion, our simple model deals with processes that are optimized for the grains to be depleted on the central object. Consequently, our Epstein and Stokes criteria for the radial-drift barrier constitute the least favourable limit for grain survival. Thus, we are confident when claiming that the radial-drift barrier does not occur in some classes of discs. One may however expect that more discs are retaining their grains due to the complementary processes mentioned above.

It should be noted that the criterion  − p + q + 1/2 quantifies the outcome of the radial-drift motion of the grains, but not their kinematics (which depends on the grain size, the grain density, etc.). Thus, we provide predictions for which discs will retain the largest mass of solid particles, but do not predict for which discs the radial migration to the inner disc regions is the fastest. Full simulations like those developed in Brauer et al. (2008) are required to make predictions of the dust kinetics, even more so when the grain size evolution is driven by a complex model of growth and fragmentation. However, this study suggests that while complex simulations are useful to study the details of the dust dynamics, they are not required to determine the grains outcome.

As a conclusion of this section, we have mentioned that the physics treated by our model is not exhaustive. In real discs, the limits $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0}$ and $\mathit{q}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}$ may be softened by the effects of additional physical phenomena. However, these neglected processes (such as turbulence and grain growth) tend to decrease the efficiency of the dust radial motion. Our predictions of when the “radial-draft barrier” does not occur therefore remain valid. Our model represents a powerful indicator for predicting the dust behaviour in discs with given power law profiles: we expect that (i) discs satisfying $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$ retain their small grain population and that (ii) discs satisfying $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$ keep their large solids. On the contrary, discs for which $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$ (resp. $\mathit{q}\mathit{>}\frac{\mathrm{2}}{\mathrm{3}}$) likely lose their small (resp. large) particles.

6. Application to observed discs and planet formation

6.1. Validity of the criteria in real protoplanetary discs

We now study how the criteria we derived can be applied when considering the physical evolution of grains in observed discs, with finite inner radii and finite lifetimes. The analytic expressions of the previous sections have been derived using dimensionless quantities. We now provide the physical timescales of the radial dust motion estimating the parameters involved in real protoplanetary discs. We consider in this section a typical CTTS disc, of mass M_disc = 10^-2  M_⊙ around a 1 M_⊙ star, extending from r_in = 10^-2 AU to r_out = 10³ AU. The disc inner edge is chosen to correspond to the dust sublimation radius for a 1 M_⊙ star, whereas its outer boundary is representative of the largest observed discs. Its vertical extent is set by the choice of the temperature scale. We take T(1 AU) = 150 K, a typical value obtained by Andrews & Williams (2007) in their disc observations.

The transition from the Epstein to the Stokes regime occurs when ${\mathit{\lambda }}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{4}\mathit{s}}{\mathrm{9}}$, or $\mathit{s}\mathrm{=}\frac{\mathrm{9}}{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{45}{\mathit{\pi }}^{\frac{\mathrm{3}}{\mathrm{2}}}}{\mathrm{256}}\frac{\mathit{m}{\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}}{{\mathit{\sigma }}_{\mathrm{s}}{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{r}}^{\mathit{p}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}}\mathit{,}$(44)and is represented in the (r,s) plane in Fig. 5 for this typical disc for different values of the surface density and temperature power-law exponents p and q. The Stokes regime is seen to apply to large bodies in the disc inner regions and for large values of both p and q.

Disc lifetimes are generally thought to be a few Myr (Haisch et al. 2001; Carpenter et al. 2005), and thus we take t_disc ~ 10⁶ yr. For a grain starting at a distance r₀ from a 1 M_⊙ star, the dimensionless value T_disc is therefore ${\mathit{T}}_{\mathrm{disc}}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{disc}}}{{\mathit{t}}_{\mathrm{k}\mathrm{0}}}\mathrm{=}\frac{\sqrt{\mathrm{𝒢}{\mathit{M}}_{\mathrm{\odot }}}\hspace{0.17em}{\mathit{t}}_{\mathrm{disc}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{~}\frac{\mathrm{6}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}}{\mathrm{[}{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathrm{AU}\mathrm{\right)}{\mathrm{]}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}$(45)The dimensionless value of the dust disc inner radius (r_in ~ 0.01 AU) for a grain starting at r₀ is ${\mathit{R}}_{\mathrm{in}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{in}}}{{\mathit{r}}_{\mathrm{0}}}\mathrm{~}\frac{{\mathrm{10}}^{-2}}{{\mathit{r}}_{\mathrm{0}}\mathrm{\left(}\mathrm{AU}\mathrm{\right)}}\mathrm{·}$(46)The link between dimensionless and real grain sizes is made through the optimal size for radial migration. Considering first the Epstein regime, its midplane value has a radial dependence given by ${\mathit{s}}_{\mathrm{opt}}\mathrm{\left(}\mathit{r,}\mathrm{0}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\sqrt{\mathrm{2}\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{d}}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{r}}^{\mathrm{-}\mathit{p}}}{\sqrt{\mathrm{2}\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{d}}}\mathrm{·}$(47)The dimensionless size of a particle of size s starting its migration at position r₀ is therefore ${\mathit{S}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{s}}{{\mathit{s}}_{\mathrm{opt}\mathit{,}\mathrm{0}}}\mathrm{=}\frac{\sqrt{\mathrm{2}\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{d}}}{{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}}\hspace{0.17em}\mathit{s}\hspace{0.17em}{\mathit{r}}_{\mathrm{0}}^{\mathit{p}}\mathit{.}$(48)Values of S₀ are plotted in the (r₀,s) plane in Fig. 6 for a disc of total mass M_disc = 0.01  M_⊙ extending from r_in = 10^-2 AU to r_out = 10³ AU, with grains of intrinsic density ρ_d = 1000 kg m^-3 and η₀ = 10^-2, for both p = 0 and $\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}$.

The dimensionless time T_in for a grain to reach R_in is given by Eq. (34). In combination with Eq. (46), this gives an expression of T_in as a function of S₀ and r₀, whereas Eq. (45) gives an expression of T_disc as a function of r₀. Equating them yields a second order equation in S₀ as a function of r₀, which can be solved to determine under which conditions a grain reaches the disc inner edge at the end of its lifetime. Using Eq. (48) gives the corresponding relationship between the grain size and its initial position: $\begin{array}{ccc}\mathit{s}& \mathrm{=}& \frac{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{2}\sqrt{\mathrm{2}\mathit{\pi }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{d}}}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{r}}_{\mathrm{0}}^{\mathrm{-}\mathit{p}}}{\mathrm{1}\mathrm{-}{\mathrm{(}{\frac{{\mathit{r}}_{\mathrm{in}}}{{\mathit{r}}_{\mathrm{0}}}}^{\mathrm{)}}}^{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}[\frac{\sqrt{\mathrm{𝒢}\mathit{M}}\hspace{0.17em}{\mathit{t}}_{\mathrm{disc}}\hspace{0.17em}{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\\ & & \end{array}$(49)which is plotted as a thick line in Fig. 6. It separates the (r₀,s) plane into regions in which grains reach r_in and leave the disc before it dissipates (t_in < t_disc) or survive in the disc throughout its lifetime (t_in > t_disc). We call this curve the survival limit.

For $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$, illustrated by the case $\mathrm{(}\mathit{p}\mathrm{=}\mathrm{0}\mathit{,q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\mathrm{)}$, Fig. 6 shows as expected that most of the grains are lost during the disc lifetime. However, small and large grains initially in the outer disc survive, therefore even with this profile, the disc retains a fraction of its grain population before it dissipates. Moreover, one may expect growing grains that reach S₀ = 1 to be inevitably accreted onto the star (unless the growth process is fast enough for grains to outgrow the fast migrating sizes before they leave the disc, see Laibe et al. 2008). Such discs may not form planets, but their remaining dust content may still make them observable, although they would likely be faint.

For $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$, illustrated by the case $\mathrm{(}\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{,q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\mathrm{)}$, even though all grains fall on the star in an infinite time, some of them reach the disc’s inner edge before it dissipates. On the contrary, for r₀ > 350 AU, grains of all sizes remain in the disc. This is also the case for all (sub)micron-sized grains, whatever their initial location, as well as most of the grains up to 0.1 mm. These grains likely make the disc bright and easy to observe, since they are the grains contributing most to the disc emission at IR and submm wavelengths. A large reservoir of grains is available to participate in the planet formation process, however a firm conclusion on their survival as their size evolves would require incorporating a treatment of grain growth, as discussed in Sect. 5.

It should be noted that the disc used in these examples represents a lower limit, as it is low mass and very extended. A more massive disc with a smaller outer radius would have a larger surface density, and the corresponding survival limit in Fig. 6 would be shifted vertically towards larger sizes and more and more grains of larger sizes would survive.

The left panels of Fig. 7 show the influence of the surface density and temperature profiles on the location and shape of the survival limit curve in the (r₀,s) plane in the Epstein regime. Increasing p from 0 to 2 shifts the curve towards smaller radii and larger grain sizes, as well as slightly tilts it clockwise. The outer disc region in which grains of all sizes survive extends inwards, as well as the surviving population of small grains as the curve’s lower branch shifts upwards and becomes flatter. On the contrary, the steepening of the curve’s upper branch, confining the population of surviving large solids to the disc outer regions, is less dramatic. Increasing q from $\frac{\mathrm{1}}{\mathrm{4}}$ to $\frac{\mathrm{3}}{\mathrm{4}}$ also tilts the curve clockwards, but shifts it towards larger radii and smaller grain sizes. However, its effect is more limited than that of changing p. A disc with a steeper surface density profile and a shallower temperature profile is therefore more efficient at retaining a larger quantity of small grains and up to larger sizes. Indeed, large p and small q values are required to meet the $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{<}\mathrm{0}$ criterion introduced in Sect. 3.

An equation very similar to Eq. (49) can be obtained for the linear Stokes regime by replacing p and q by ${\mathit{p}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{q}\mathrm{-}\mathrm{3}}{\mathrm{2}}$ and q′ = q (since the equation of motion has the same structure for both drag regimes, see Sect. 4), and using the expression of s_opt,0 for that regime, given in Table C.1. Here ${\mathit{s}}_{\mathrm{opt}\mathit{,}\mathrm{0}}\mathrm{\propto }{\mathit{T}}^{\frac{\mathrm{1}}{\mathrm{4}}}$: the weak temperature dependence results in very little change for a large range of temperatures below or above our adopted value of T(1 AU) = 150 K. The right panels of Fig. 7 show the survival limit in the Stokes regime for the different values of q (note that it no longer depends on p). When q increases, the curve’s lower branch slides towards larger radii, making the survival of particles in the Stokes regime less and less favourable.

Given the form of Eq. (49) and its different expressions for each drag regime, it is not possible to compute analytically the survival limit for a grain transitioning from the Epstein to the Stokes regime as it migrates inwards. However, large values of p and q, for which the Stokes region is the largest, are not observed in protoplanetary discs (see Sect. 6.2), and in practical cases the Stokes regime only applies to a small area of the (r₀,s) plane. Small grains, which are detected in disc observations at IR and sub-millimetre (submm) wavelengths, are mostly subject to the Epstein drag. We therefore focus on that regime in the following.

Equation (49) can give quantitative information about the outcome of the grain population. Replacing t_disc by any time t gives the location in the (r₀,s) plane of grains reaching the disc inner edge (at r = r_in) in that time t, which is therefore the survival time of those grains. Its isocontours are shown in Fig. 8 for different values of p and for $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$. Only one value of q is shown as the q dependence is moderate, as can be seen from the left panels of Fig. 7. The fate of particular dust grains can easily be obtained from these figures. For example, in the context of disc observations, 1 mm grains initially at 100 AU fall on the 10⁵ yr contour for p = 0. Their survival time decreases to a few 10⁴ yr for p ~ 0.8, and increases again to values larger than 10⁶ yr as p increases. At an initial position of a few hundred AU, 1 mm grains survive longer than 10⁶ yr for any p, therefore long enough to contribute to the disc emission over its entire lifetime. As noted above, such grains have longer survival times for higher p and lower q values. As another example, in the context of planetesimal formation, the survival time of a 1 m particle initially at 1 AU is  ~10⁵ yr for p = 0, decreases to  ~10² yr for p ~ 1, and increases again to  ~10³ − 10⁴ yr for p = 2. The ability of such particles to remain for long enough in the disc to grow to larger sizes therefore strongly depends on the surface density profile. As a general rule, the survival of pre-planetesimals in the inner disc is favoured by small values of p.

Similarly, replacing now r_in with any radius r in Eq. (49) gives the locus in the (r₀,s) plane of grains reaching that radius r at any time t. Alternatively, one can plot isocontours of the initial position r₀ of grains in the (r,s) plane at various times, thus showing the radial evolution of grains with the same initial position but different sizes. This is shown in Fig. 9 for (p = 0, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$) and Fig. 10 for ($\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}$, $\mathit{q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}$). These plots make it easy to compare the radial evolution of any particle to any physical timescale of interest in the disc. In particular, they show that the disc still contains particles at all radii at the end of its evolution (t = 10⁶ yr). No grains are found to the right of the r₀ = 10³ AU contour, since this was the initial outer disc radius. In the disc with $\mathrm{(}\mathit{p}\mathrm{=}\mathrm{0}\mathit{,q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\mathrm{)}$, no grains between  ~0.06 and  ~0.2 mm remain, and grains of other sizes still present were initially in the outer disc. Given that the grains of sizes which contribute to IR and submm emission have come from a small fraction of the initial disc, this disc is likely faint. In the disc with $\mathrm{(}\mathit{p}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\mathit{,q}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\mathrm{)}$, only grains with s ~ 0.1 mm are absent from the very outer regions, and the observable grains come from a larger portion of the disc, likely making the disc brighter than in the previous case.

As a conclusion, the analytic criteria derived above apply even when taking into account the finite lifetime (or inner radius) of the disc. For most CTTS discs, the dust is in the Epstein drag regime (except for some extreme values for the grains sizes and discs profiles). Therefore, the grain’s radial outcome is given by the value of $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$. However, the transition between discs for which the radial-drift barrier occurs or not consists more of a continuum around the value $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0}$ than in the sharp transition predicted by the analytic model. Therefore, the radial motion of the grains has to be studied on a case-by-case basis for discs close to the transition $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0}$, using the figures shown above in this section.

6.2. Constraining physical systems

We now turn to observed discs and check if they meet our Epstein and Stokes criteria to determine whether the radial-drift barrier is constraining for planet formation. To estimate the values of the p and q exponents for real discs, we use the results of disc modeling obtained by Andrews & Williams (2005, 2007) from data on 63 discs in ρ Ophiuchi, Taurus and Aurigae. Using sub-millimetre fluxes measured at several wavelengths, they fit a range of disc parameters assuming a geometrically thin irradiated disc with opacities from Beckwith et al. (1990), a gas-to-dust ratio of 100, a disc radius of 100 AU and zero disc inclination. The temperature exponent q is well constrained by the observational data set: the histogram of most probable q values is shown in Fig. 11.

However, p is not well constrained and is usually assumed to be $\frac{\mathrm{3}}{\mathrm{2}}$. Very flat profiles with $\mathit{p}\mathit{<}\frac{\mathrm{1}}{\mathrm{2}}$ and very steep profiles with $\mathit{p}\mathit{>}\frac{\mathrm{3}}{\mathrm{2}}$ seem to be excluded (Dutrey et al. 1996; Wilner et al. 2000; Kitamura et al. 2002; Testi et al. 2003; Isella et al. 2009; Andrews & Williams 2007; Andrews et al. 2009). We represent the disc distribution modeled by Andrews & Williams from observations in the (p,q) diagram of Fig. 12: the histogram of Fig. 11 is represented by the gray-shaded area and spread over a range of p values, taking into account that extreme values of p are less probable. The dashed line ( − p + q + 1/2 = 0) represents the border between migration in an infinite time and accretion onto the central star for the A-mode of migration in the Epstein regime, while the thick dotted line (q = 3/2) represents that same border for the B-mode of migration in the Stokes regime at low Reynolds number. The two black circles indicate the discs used as examples in Sect. 3.3.

We have split the disc distribution in four regions in the (p,q) plane:

• 1.

  region 1: $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$ and $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$: both small and large grains experience the pile-up effect. Those discs are potentially observable and may favour planet formation.

• 2.

  region 2: $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$ and $\mathit{q}\mathit{>}\frac{\mathrm{2}}{\mathrm{3}}$: only small grains experience the pile-up effect: even though such discs retain their small grains, the population of pre-planetesimals in the disc inner regions may efficiently be accreted onto the central star (at least until they reach the high-R_g Stokes regime).

• 3.

  region 3: $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$ and $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$: if the pre-planetesimals can form before the entire distribution of small grains has been accreted onto the central object, they will remain in the disc and may constitute planet embryos.

• 4.

  region 4: $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$ and $\mathit{q}\mathit{>}\frac{\mathrm{2}}{\mathrm{3}}$: both small and large grains are accreted onto the central star.

The Epstein criterion indicates that for q values in the range constrained by observations, discs which keep their small grain population, and are therefore likely to be bright in the IR and submm, should have p values approximately in the $\mathrm{\left[}\mathrm{1}\mathrm{;}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{\right]}$ range. This is indeed what is found in most disc surveys (Ricci et al. 2010a,b). On the contrary, smaller p values should correspond to discs which lose most of their small grains, and are therefore more difficult to detect. This is what is found by Andrews et al. (2010), who pushed their previous observations of the Ophiuchus star forming region (Andrews et al. 2009) down to fainter discs, finding for this new sample a median p value of 0.9, lower than for brighter discs. The criterion we derive in this paper for small grains in the Epstein regime provides therefore the correct behaviour for explaining the range of p values of observed discs. However, this result has to be considered carefully for two reasons. Firstly, the p and the q exponents determined from the observations have to be considered with their respective errors. Given these uncertainties, one may not be able to distinguish between a strict negative or positive value for $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$. Second, the boundary between the different zones of the (p,q) plane consists more of a continuum rather than a strict limit due to the finite lifetime/inner radii of the discs. The outcome of the grains may thus not be predicted when the value of $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$ is close to zero.

Now turning to the Stokes criterion for large solids, Figs. 11 and 12 show that the vast majority of observed protoplanetary discs have shallow temperature profiles ($\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$) and are thus able to retain their population of pre-planetesimals. These discs are therefore relevant places to find evidence of planet formation, provided small grains can efficiently grow to form pre-planetesimals. For the remainder of the disc population, the outcome of pre-planetesimals will likely depend on their ability to reach the high Reynolds number Stokes regime. However, the case of a steep radial temperature profile can be encountered in at least one particular situation: circumplanetary discs which typically have temperature profiles with q = 1 (Ayliffe & Bate 2009). In this environment, we predict from our Stokes criterion that planetesimals will be accreted onto the planet. The timescale of the planet formation by the core-accretion process, which usually corresponds to the time required to release the gravitational energy of the accreted bodies (Pollack et al. 1996), is thus increased as the drag from the gas onto the planetesimals releases an additional thermal contribution.

7. Conclusion and perspectives

In this study, we have generalised the radial grain motion studies of W77, NSH86 and YS02 for both the Epstein and the Stokes regimes, taking into account the effects of both the surface density and temperature profiles in the disc. As observations do not provide direct information about the three dimensional structure of discs, radial profiles of surface density and temperature are often described by power laws: $\mathrm{\Sigma }{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{p}}$ and $\mathrm{𝒯}{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}\mathrm{=}{\mathrm{𝒯}}_{\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{q}}$, where both p and q take positive values. The radial dust behaviour in those discs is governed by the competition between gravity and gas drag. The final outcome of the radial motion is set by two counterbalancing effects. First, the temperature increases when the radius decreases. Consequently, the deviation from the Keplerian velocity increases, which accelerates the dust’s radial inward motion. At the same time, the surface density also increases, which increases the gas drag efficiency and slows down the dust motion. The competition between these two effects fixes the ultimate mode of migration of the grain (A-mode, where the drag dominates or B-mode, where the gravity dominates) and thus the final outcome for the dust motion. In this work, we have shown that it can be represented by an analytical criterion which depends on the drag regime. For the Epstein drag regime (in which the ultimate radial motion is in the A-mode), if $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$, the dust particle is accreted onto the central star in a finite time, and if $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$, the grain pile-up results in an infinite accretion time and small dust grains remain in the disc. We have shown that, as expected, these conclusions are somewhat mitigated when taking into account the finite disc lifetime and finite disc size. However, the outcomes still remain similar: the bulk of the small grain population is lost to the star in the first case, whereas in the second case the disc keeps most of its small grains. A similar criterion is found for the Stokes regime at low Reynolds number: if $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$, the accretion time is infinite and large pre-planetesimals remain in the disc and can constitute the primary material for planet formation. However, the Stokes radial motion differs from the Epstein regime as the ultimate radial motion occurs in the B-mode.

The observational consequence is that discs with a large population of small grains should be strong emitters in the infrared and sub-millimetre and should be easier to observe, and that those having lost most of their small dust should be fainter and harder to detect. This is indeed what is found: a large fraction of the observed discs have large p values whereas fainter discs tend to have lower p values (Andrews et al. 2010), in agreement with this Epstein criterion. In addition, most of the observed discs have $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$, allowing them to retain also their large pre-planetesimals. As noted by Ricci et al. (2010a,b), explaining the data requires a mechanism halting or slowing down the radial migration of dust grains. We show here that local pressure maxima need not be invoked, but rather that the combination of adequate surface density and temperature profiles is sufficient. The p and q exponents used to reach our conclusions are of course strongly dependent on the model used to fit the data. However, even varying the fitting models, a large majority of discs still satisfy both conditions $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$ and $\mathit{q}\mathrm{\le }\frac{\mathrm{2}}{\mathrm{3}}$. Consequently, the radial-drift barrier (or the so-called metre-size barrier when considering an MMSN disc) does not appear to constitute a problem for planet formation for the discs that we do observe.

Our conclusions presented in this study assumed that the grain size remains constant during its motion. However, observations tell us that grains do grow (Testi et al. 2003; Wilner et al. 2003; Apai et al. 2005; Lommen et al. 2007, 2009). Grain growth is studied in various theoretical studies (Schmitt et al. 1997; Stepinski & Valageas 1997; Suttner et al. 1999; Tanaka et al. 2005; Dullemond & Dominik 2005; Klahr & Bodenheimer 2006; Garaud 2007; Brauer et al. 2008; Laibe et al. 2008; Birnstiel et al. 2009). In a forthcoming paper, we will generalise the formalism developed here to explain the radial and vertical behaviour of growing dust grains.

Online material

Appendix B: Disc structure

B.1. Hydrostatic equilibrium

At stationary equilibrium ($\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{=}\mathrm{0}$), gas velocities v_gr, v_gθ, v_gz and the gas density ρ_g obey mass conservation and the Euler equation: $\{\begin{array}{c}\\ \frac{\mathrm{1}}{\mathit{r}}\frac{{\mathit{\rho }}_{\mathrm{g}}\mathit{\partial r}{\mathit{v}}_{\mathrm{g}\mathit{r}}}{\mathit{\partial r}}\mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{g}}\mathit{\partial }{\mathit{v}}_{\mathrm{g}\mathit{z}}}{\mathit{\partial z}}& \mathrm{=}& \mathrm{0}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{g}}{v}\mathrm{g}\mathit{.}\mathrm{\nabla }{v}\mathrm{g}& \mathrm{=}& \mathrm{-}\mathrm{\nabla }\mathit{P}\mathrm{+}{\mathit{\rho }}_{\mathrm{g}}{g}\mathit{.}\end{array}$(B.1)The solution of Eq. (B.1) requires: ${\mathit{v}}_{\mathrm{g}\mathit{r}}\mathrm{=}{\mathit{v}}_{\mathrm{g}\mathit{z}}\mathrm{=}\mathrm{0}\mathit{,}$(B.2)which ensures mass conservation. Projecting the Euler equation on e_z: $\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\mathit{\partial P}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathrm{𝒢}\mathit{Mz}}{{\left({\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}$(B.3)Assuming that: $\mathit{P}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}{\mathit{\rho }}_{\mathrm{g}}\mathit{,}$(B.4)and dividing both sides of Eq. (B.3) by ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, we have: $\frac{\mathit{\partial }\hspace{0.17em}\ln \mathrm{\left(}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\right)}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathrm{𝒢}\mathit{Mz}}{{\left({\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{·}$(B.5)Integrating Eq. (B.5) between 0 and z provides: ${\mathit{\rho }}_{\mathrm{g}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{=}\frac{\mathit{P̅}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}}{\mathrm{e}}^{\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\frac{\mathrm{𝒢}\mathit{M}{\mathit{z}}^{\mathrm{\prime }}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}}{{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{\prime }\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{r,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}}}\mathit{.}$(B.6)This expression can be simplified by the following approximations:

• In the special vertically isothermal case, wherethe sound speed depends only on the radial co-ordinate, Eq. (B.6) simplifies to:${\mathit{\rho }}_{\mathrm{g}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{=}\mathit{\rho ̅}\mathrm{g}\mathrm{(}{\mathit{r}}^{\mathrm{)}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{𝒢}\mathit{M}}{\mathit{c̅}\begin{array}{}\mathrm{2}\\ \mathrm{s}\end{array}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}\left[\frac{\mathrm{1}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\right]}\mathit{.}$(B.7)

• Further, assuming a thin disc (${{}^{\mathrm{\right(}}\frac{\mathit{z}}{\mathit{r}}^{\mathrm{\left)}}}^{\mathrm{2}}\mathrm{\ll }\mathrm{1}$), a Taylor series expansion of Eq. (B.7) leads to: ${\mathit{\rho }}_{\mathrm{g}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{=}\mathit{\rho ̅}\mathrm{g}\mathrm{(}{\mathit{r}}^{\mathrm{)}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}\mathit{H}{\mathrm{(}{\mathit{r}}^{\mathrm{)}}}^{\mathrm{2}}}}\mathit{,}$(B.8)with: $\mathit{H}\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}\frac{\mathit{r}\mathit{c̅}\mathrm{s}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}{{\mathit{v}}_{\mathrm{k}}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}\mathit{,}$(B.9)which is the classical scale height for vertically isothermal thin discs.

B.2. Azimuthal velocity

The radial component of the Euler equation is given by: $\mathrm{-}\frac{{\mathit{v}}_{\mathrm{g}\mathit{\theta }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\mathit{\partial P}}{\mathit{\partial r}}\mathrm{-}\frac{\mathrm{𝒢}\mathit{Mr}}{{\left({\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(B.10)where ρ_g is given by Eq. (B.6). Thus: ${\mathit{\rho }}_{\mathrm{g}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{=}\frac{\mathit{P̅}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}}{\mathrm{e}}^{\mathrm{-}{\mathit{I}}_{\mathrm{1}}\left(\mathit{r,z}\right)}\mathit{,}$(B.11)with: $\{\begin{array}{c}\\ {\mathit{I}}_{\mathrm{1}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}& \mathrm{=}& {\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathrm{𝒢}\mathit{M}{\mathit{c}}_{\mathrm{s}}^{-2}\left(\mathit{r,}{\mathit{z}}^{\mathrm{\prime }}\right){\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}\mathit{f}\left(\mathit{r,}{\mathit{z}}^{\mathrm{\prime }}\right)\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{,}\\ \mathit{f}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\mathit{.}\end{array}$(B.12)To simplify Eq. (B.10), we first use the following identity: $\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\hspace{0.17em}\mathit{\partial P}}{\mathit{\partial r}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathit{\partial }\hspace{0.17em}\ln \mathrm{\left(}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\right)}}{\mathit{\partial r}}\mathit{,}$(B.13)which becomes with Eq. (B.11): $\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\mathit{\partial P}}{\mathit{\partial r}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathrm{\left(}\mathit{P̅}\mathrm{\right)}}{\mathrm{d}\mathit{r}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathit{\partial }\hspace{0.17em}{\mathit{I}}_{\mathrm{1}}}{\mathit{\partial r}}\mathrm{·}$(B.14)Noting that f⁽r,z = 0⁾ = 0 and integrating I₁ by parts provides: ${\mathit{I}}_{\mathrm{1}}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{=}\mathrm{𝒢}\mathit{M}{\mathit{c}}_{\mathrm{s}}^{-2}\mathit{f}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathrm{𝒢}\mathit{Mf}\left(\mathit{r,}{\mathit{z}}^{\mathrm{\prime }}\right){\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}{\mathit{c}}_{\mathrm{s}}^{-2}\left(\mathit{r,}{\mathit{z}}^{\mathrm{\prime }}\right)\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{,}$(B.15)and: $\begin{array}{ccc}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}}\frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{r}}& \mathrm{=}& {\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathrm{(}\mathit{P̅}\mathrm{)}}{\mathrm{d}\mathit{r}}\mathrm{-}(\mathrm{-}\frac{\mathrm{𝒢}\mathit{M}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{𝒢}\mathit{Mr}}{{\left({\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\\ & & \end{array}$(B.16)Then, Eq. (B.10) becomes: $\begin{array}{ccc}\frac{{\mathit{v}}_{\mathrm{g}\mathit{\theta }}^{\mathrm{2}}}{\mathit{r}}& \mathrm{=}& {\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathrm{(}\mathit{P̅}\mathrm{)}}{\mathrm{d}\mathit{r}}\mathrm{+}\frac{\mathrm{𝒢}\mathit{M}}{{\mathit{r}}^{\mathrm{2}}}\\ & & \end{array}$(B.17)Noting that : $\frac{\mathit{\partial }}{\mathit{\partial r}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathit{f}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\frac{\mathit{\partial f}}{\mathit{\partial r}}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathit{f}{\mathit{\partial }}_{\mathit{r}}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{,}$(B.18)and integrating the last term of the right hand side of Eq. (B.18) by parts provides: ${\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathit{f}{\mathit{\partial }}_{\mathit{r}}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathrm{=}\mathit{f}{\mathit{\partial }}_{\mathit{r}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}\mathit{f}{\mathit{\partial }}_{\mathit{r}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{.}$(B.19)Therefore, Eq. (B.17) reduces to: $\frac{{\mathit{v}}_{\mathrm{g}\mathit{\theta }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\frac{\mathrm{𝒢}\mathit{M}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\mathrm{+}\mathrm{𝒢}\mathit{M}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathrm{(}{\mathit{\partial }}_{\mathit{r}}\mathit{f}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}{\mathit{\partial }}_{{\mathit{z}}^{\mathrm{\prime }}}\mathit{f}{\mathit{\partial }}_{\mathit{r}}{{\mathit{c}}_{\mathrm{s}}^{-2}}^{\mathrm{)}}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{,}$(B.20)which can be more elegantly written as: $\frac{{\mathit{v}}_{\mathrm{g}\mathit{\theta }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\frac{\mathrm{𝒢}\mathit{M}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\mathrm{+}\mathrm{𝒢}\mathit{M}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\mathrm{[}\mathrm{\nabla }\mathit{f}\mathrm{\times }\mathrm{\nabla }{{\mathit{c}}_{\mathrm{s}}^{-2}}^{\mathrm{]}}\mathit{.}{e}\mathit{\theta }\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{.}$(B.21)Thus, the expression of the azimuthal velocity of such a disc can be separated in three terms, called the Keplerian, the pressure gradient and the baroclinic terms respectively. This last term is neglected in most studies. For a three dimensional disc, this term rigorously cancels for c_s =  constant. In this case, the flow is inviscid and derives from a potential, and isobars and isodensity surfaces coincide: thus, there is no source of vorticity and the azimuthal velocity depends only on the radial coordinate. This terms also cancels out for flat discs in two dimensions. If the disc is vertically isothermal, Eq. (B.21) becomes: $\frac{{\mathit{v}}_{\mathrm{g}\mathit{\theta }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\frac{\mathrm{𝒢}\mathit{M}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\mathrm{-}\mathrm{𝒢}\mathit{M}\left[\frac{\mathrm{1}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\right]{\mathit{\partial }}_{\mathit{r}}\ln {\mathit{c}}_{\mathrm{s}}^{-2}\mathit{.}$(B.22)

B.3. Radial profiles of surface density and temperature

In this section, we consider that the disc surface density and the temperature (and thus the sound speed) depend only on the radial coordinate and are given by the following power-law profiles: $\{\begin{array}{c}\\ \mathrm{\Sigma }\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}}^{\mathrm{)}}}^{\mathrm{-}\mathit{p}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{p}}\mathit{,}\\ \mathit{T}\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathit{T}}_{\mathrm{0}}{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}}^{\mathrm{)}}}^{\mathrm{-}\mathit{q}}\mathrm{=}{\mathit{T}}_{\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{q}}\mathit{,}\\ {\mathit{c}}_{\mathrm{s}}\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{0}}{\mathrm{(}{\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}}^{\mathrm{)}}}^{\mathrm{-}\mathit{q}\mathit{/}\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{q}\mathit{/}\mathrm{2}}\mathit{.}\end{array}$(B.23)For vertically isothermal thin discs, the vertical density is therefore given by Eq. (B.8) with the scale height given by Eq. (B.9), which can be expressed as: $\mathit{H}\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathit{H}}_{\mathrm{0}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}}\mathrm{=}{\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}}\mathit{,}$(B.24)with: ${\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}\mathrm{0}}^{\mathrm{\prime }}}{\sqrt{\mathrm{𝒢}\mathit{M}}}\mathrm{·}$(B.25)The expression of ρ_g compatible with the vertical hydrostatic equilibrium and providing the power-law profile set by Eq. (B.23) is written as: ${\mathit{\rho }}_{\mathrm{g}}\mathrm{=}{\mathit{\rho }}_{\mathrm{g}\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{x}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}}\mathit{.}$(B.26)Indeed: $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\rho }}_{\mathrm{g}\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{x}}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}^{\mathrm{)}}}}\mathrm{d}\mathit{z}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{g}\mathrm{0}}^{\mathrm{\prime }}\sqrt{\mathrm{2}\mathit{\pi }}\mathit{H}\mathrm{(}{\mathit{r}}^{\mathrm{)}}{\mathit{r}}^{\mathrm{-}\mathit{x}}\mathit{,}\\ & \mathrm{=}& {\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}{\mathit{r}}^{\mathrm{-}\mathit{p}}\mathit{.}\\ & & \end{array}$(B.27)Hence, with ${\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\rho }}_{\mathrm{g}\mathrm{0}}^{\mathrm{\prime }}{\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}$ and $\mathit{x}\mathrm{=}\mathit{p}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}$, ${\mathit{\rho }}_{\mathrm{g}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}}{\mathit{r}}^{\mathrm{-}\left(\mathit{p}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}{\mathrm{e}}^{\mathrm{-}\left[\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }\mathrm{2}}{\mathit{r}}^{\mathrm{3}\mathrm{-}\mathit{q}}}\right]}\mathit{,}$(B.28)which gives the correct surface density profile when integrated with respect to z. With this expression of ρ_g, $\mathit{P̅}{}^{\mathrm{\right(}}\mathit{r}^{\mathrm{\left)}}$ is given by: $\mathit{P̅}\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}{\mathit{r}}^{\mathrm{)}}{\mathit{\rho }}_{\mathrm{g}}\mathrm{(}\mathit{r,z}\mathrm{=}{\mathrm{0}}^{\mathrm{)}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{0}}^{\mathrm{\prime }\mathrm{2}}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{\prime }}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{0}}^{\mathrm{\prime }}}{\mathit{r}}^{\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}\mathit{,}$(B.29)which ensures that: $\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{(}\mathit{p}\mathrm{+}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}{\frac{\mathrm{3}}{\mathrm{2}}}^{\mathrm{)}}}{\mathit{r}}\mathit{,}$(B.30)and, using Eq. (B.22), we find: ${\mathit{v}}_{\mathrm{g}\mathit{\theta }}\mathrm{=}\sqrt{\frac{\mathrm{𝒢}\mathit{M}}{\mathit{r}}\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right){\mathit{c}}_{\mathrm{s}}^{\mathrm{\prime }\mathrm{2}}{\mathit{r}}^{\mathrm{-}\mathit{q}}\mathrm{-}\mathrm{𝒢}\mathit{Mq}\left(\frac{\mathrm{1}}{\mathit{r}}\mathrm{-}\frac{\mathrm{1}}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\right)}\mathit{.}$(B.31)

Appendix C: Dimensionless quantities and equations of motion

To highlight the important physical parameters involved, we set ${\mathit{v}}_{\mathrm{k}\mathrm{0}}\mathrm{=}\sqrt{\frac{\mathrm{𝒢}\mathit{M}}{{\mathit{r}}_{\mathrm{0}}}}$ and introduce dimensionless quantities given by the following expressions: $\{\begin{array}{c}\\ \mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}& \mathrm{=}& \mathit{R,}\\ \mathrm{𝒯}\mathit{/}{\mathrm{𝒯}}_{\mathrm{0}}& \mathrm{=}& {\mathit{R}}^{\mathrm{-}\mathit{q}}\mathit{,}\\ \mathrm{\Sigma }\mathit{/}{\mathrm{\Sigma }}_{\mathrm{0}}& \mathrm{=}& {\mathit{R}}^{\mathrm{-}\mathit{p}}\mathit{,}\\ {\mathit{v}}_{\mathrm{k}}\mathit{/}{\mathit{v}}_{\mathrm{k}\mathrm{0}}& \mathrm{=}& {\mathit{R}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\mathit{,}\\ \mathit{H}\mathit{/}{\mathit{H}}_{\mathrm{0}}& \mathrm{=}& {\mathit{R}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\frac{\mathit{q}}{\mathrm{2}}}\mathit{,}\end{array}$with: ${\mathit{\eta }}_{\mathrm{0}}\mathrm{=}\left(\mathit{p}\mathrm{+}\frac{\mathit{q}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right){\left(\frac{{\mathit{c}}_{\mathrm{s}\mathrm{0}}}{{\mathit{v}}_{\mathrm{k}\mathrm{0}}}\right)}^{\mathrm{2}}\mathit{.}$(C.2)The dimensionless parameter η₀ gives the order of magnitude of the relative discrepancy between the Keplerian motion and the gas azimuthal velocity. We note that: $\frac{\mathit{r}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{k}}}\frac{\mathrm{d}\ln \mathit{P̅}}{\mathrm{d}\mathit{r}}/{\mathit{v}}_{\mathrm{k}\mathrm{0}}\mathrm{=}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(C.3)Then, we set ${\mathit{t}}_{\mathrm{k}\mathrm{0}}\mathrm{=}\sqrt{\frac{{\mathit{r}}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{𝒢}\mathit{M}}}$ and define $\{\begin{array}{c}\\ \frac{\mathit{t}}{{\mathit{t}}_{\mathrm{k}\mathrm{0}}}& \mathrm{=}& \mathit{T}\\ \frac{{\mathit{v}}_{\mathit{r}}}{{\mathit{v}}_{\mathrm{k}\mathrm{0}}}& \mathrm{=}& \frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{T}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\\ \frac{{\mathit{v}}_{\mathit{\theta }}}{{\mathit{v}}_{\mathrm{k}\mathrm{0}}}& \mathrm{=}& \mathit{R}\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}\mathit{T}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\end{array}$Writing the coefficient $\begin{array}{}{}_{\mathrm{˜}}\\ \mathrm{𝒞}\end{array}{}^{\mathrm{(}}\mathit{r,}{\mathit{z}}^{\mathrm{)}}$ of the drag force of Eq. (5) as: $\begin{array}{}\mathrm{˜}\\ \mathrm{𝒞}\end{array}\mathrm{(}\mathit{r,}{\mathit{z}}^{\mathrm{)}}\mathrm{=}{\mathrm{𝒞}}_{\mathrm{0}}\mathrm{𝒞}\mathrm{(}\mathit{R,}{\mathit{Z}}^{\mathrm{)}}\mathit{,}$(C.5)and using dimensionless coordinates, we have: $\frac{{F}\mathrm{D}\mathit{/}{\mathit{m}}_{\mathrm{d}}}{{\mathit{v}}_{\mathrm{k}\mathrm{0}}\mathit{/}{\mathit{t}}_{\mathrm{k}\mathrm{0}}}\mathrm{=}\mathrm{-}\frac{\mathrm{𝒞}\mathrm{(}\mathit{R,}{\mathit{Z}}^{\mathrm{)}}}{{\left[\frac{\mathit{s}}{{\left({\mathit{v}}_{\mathrm{k}\mathrm{0}}^{\mathit{\lambda }}{\mathit{t}}_{\mathrm{k}\mathrm{0}}{\mathrm{𝒞}}_{\mathrm{0}}\right)}^{\frac{\mathrm{1}}{\mathit{y}}}}\right]}^{\mathit{y}}}\mathrm{|}\begin{array}{}\mathrm{˜}\\ {v}\end{array}\mathrm{-}{\begin{array}{}\mathrm{˜}\\ {v}\end{array}}_{\mathrm{g}}{\mathrm{|}}^{\mathit{\lambda }}\mathrm{(}\begin{array}{}\mathrm{˜}\\ {v}\end{array}\mathrm{-}{\begin{array}{}\mathrm{˜}\\ {v}\end{array}}_{\mathrm{g}}\mathrm{)}\mathit{.}$(C.6)We also introduce: ${\mathit{s}}_{\mathrm{opt}\mathit{,}\mathrm{0}}\mathrm{=}{\mathrm{\left(}{\mathit{v}}_{\mathrm{k}\mathrm{0}}^{\mathit{\lambda }}{\mathit{t}}_{\mathrm{k}\mathrm{0}}{\mathrm{𝒞}}_{\mathrm{0}}\mathrm{\right)}}^{\frac{\mathrm{1}}{\mathit{y}}}\mathit{,}$(C.7)and ${\mathit{S}}_{\mathrm{0}}\mathrm{=}\frac{\mathit{s}}{{\mathit{s}}_{\mathrm{opt}\mathit{,}\mathrm{0}}}\mathit{,}$(C.8)so that Eq. (C.6) becomes: $\frac{{F}\mathrm{D}\mathit{/}{\mathit{m}}_{\mathrm{d}}}{{\mathit{v}}_{\mathrm{k}\mathrm{0}}\mathit{/}{\mathit{t}}_{\mathrm{k}\mathrm{0}}}\mathrm{=}\mathrm{-}\frac{\mathrm{𝒞}\mathrm{(}\mathit{R,}{\mathit{Z}}^{\mathrm{)}}}{{\mathit{S}}_{\mathrm{0}}^{\mathit{y}}}\mathrm{|}\begin{array}{}\mathrm{˜}\\ {v}\end{array}\mathrm{-}{\begin{array}{}\mathrm{˜}\\ {v}\end{array}}_{\mathrm{g}}{\mathrm{|}}^{\mathit{\lambda }}\mathrm{(}\begin{array}{}\mathrm{˜}\\ {v}\end{array}\mathrm{-}{\begin{array}{}\mathrm{˜}\\ {v}\end{array}}_{\mathrm{g}}\mathrm{)}\mathit{.}$(C.9)Physically, s_opt,0 corresponds to the grain size at which the drag stopping time equals the Keplerian time at r₀. In Table C.1, we give the expressions of y, λ, s_opt,0, for the Epstein and the three Stokes drag regimes. The dimensionless equations of motion for a dust grain are then: $\begin{array}{ccc}\{\begin{array}{c}\\ & & \\ & & \\ & & \end{array}& & \end{array}$

Appendix D: Lemma for the different expansions

Lemma: Let x be either r or θ and i the order of the perturbative expansion. If:

• $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{=}\mathrm{0}$, and

• $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{xi}}\end{array}$ can be written as a function of R ($\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{xi}}\end{array}\mathrm{=}\mathit{f}{}^{\mathrm{\right(}}\mathit{R}^{\mathrm{\left)}}$) with $\mathit{f}\mathrm{=}\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$ of the expansion in η₀,

then, $\frac{\mathrm{d}\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{xi}}\end{array}}{\mathrm{d}\mathit{T}}$ is of order .

Proof: $\frac{\mathrm{d}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{xi}}\end{array}}{\mathrm{d}\mathit{T}}\mathrm{=}\frac{\mathrm{d}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{xi}}\end{array}}{\mathrm{d}\mathit{R}}\frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{T}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}{\mathit{f}}^{\mathrm{\prime }}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{=}{\mathit{\eta }}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}{\mathit{f}}^{\mathrm{\prime }}\mathrm{(}{\mathit{R}}^{\mathrm{)}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{=}\mathrm{𝒪}\left({\mathit{\eta }}_{\mathrm{0}}\right)\mathit{.}$(D.1)

Appendix E: Epstein regime: perturbation analysis

• Order $\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$: At this order of expansion,η₀R ^− q is negligible compared to $\frac{\mathrm{1}}{\mathit{R}}$. Thus, substituting Eq. (13) into Eq. (11) provides $\{\begin{array}{c}\\ \frac{\mathrm{d}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}}{\mathrm{d}\mathit{T}}\mathrm{=}\frac{\stackrel{}{{{\mathit{v}}_{\mathit{\theta }\mathrm{0}}^{\mathrm{2}}}_{\mathrm{˜}}}}{\mathit{R}}\mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}\mathrm{-}\frac{\mathrm{1}}{{\mathit{R}}^{\mathrm{2}}}\\ \frac{\mathrm{d}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}}{\mathrm{d}\mathit{T}}\mathrm{=}\mathrm{-}\frac{\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}}{\mathit{R}}\mathrm{-}\frac{\left(\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{-}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\right)}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}\mathit{.}\end{array}$(E.1)At this stage, we do not know the order of $\frac{\mathrm{d}\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}}{\mathrm{d}\mathit{T}}$. We show in the lemma of Appendix D that $\frac{\mathrm{d}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}}{\mathrm{d}\mathit{T}}\mathrm{=}\mathrm{𝒪}\left({\mathit{\eta }}_{\mathrm{0}}\right)$. Applying this lemma, we see that at the order $\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$, taking $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{=}\mathrm{0}$ and $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}$ (which ensures that $\frac{\mathrm{d}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}}{\mathrm{d}\mathit{T}}\mathrm{=}\mathrm{𝒪}\left({\mathit{\eta }}_{\mathrm{0}}\right)$) is a relevant solution for the equations of motion (which corresponds to circular Keplerian motion). Thus, $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{=}\mathrm{0}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathit{.}\end{array}$(E.2)

• Order : Applying the lemma in this order of expansion and noting that $\frac{\mathrm{d}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}}{\mathrm{d}\mathit{T}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{R}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{2}}{\mathit{R}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}$(E.3)Eq. (11) becomes $\{\begin{array}{c}\\ \mathrm{0}\mathrm{=}\mathrm{-}\frac{{\mathit{R}}^{\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}}{{\mathit{S}}_{\mathrm{0}}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{+}\mathrm{2}{\mathit{R}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\\ \mathrm{0}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{R}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{-}\frac{\left(\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{0}}}\left(\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}\mathrm{-}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\right)\right)}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\left(\mathit{p}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}\mathit{.}\end{array}$(E.4)Solving the linear system Eq. (E.4) for ${\mathrm{(}}^{}\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathit{,}{\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}}^{\mathrm{)}}$ provides $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{0}}}\frac{\mathrm{2}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}}}\right)}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{0}}}\frac{{\mathit{R}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}}}\right)}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{.}\end{array}$(E.5)

In addition to the expression of $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{r}}\end{array}$ given in Sect. 3, we also note that $\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}& \mathrm{=}& \sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{+}{\mathit{\eta }}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\\ & \mathrm{=}& \sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{-}\frac{{\mathit{R}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}}}\right)}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}\end{array}$(E.6)which provides with Eq. (15) $\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{2}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{.}$(E.7)

Appendix F: Link with W77’s original derivation

Following W77’s historic reasoning for small grains (see Sect. 3.1), we perform a perturbative expansion of the radial equation of motion with the S₀ variable. We verify that taking the limit at small η₀ provides the expression found for the A-mode in the NSH86 expansion. (Formally, we will show that $\begin{array}{}\lim \\ {\mathit{S}}_{\mathrm{0}}\mathrm{\ll }\mathrm{1}\end{array}\begin{array}{}\lim \\ {\mathit{\eta }}_{\mathrm{0}}\mathrm{\ll }\mathrm{1}\end{array}\left[\mathit{...}\right]\mathrm{=}\begin{array}{}\lim \\ {\mathit{\eta }}_{\mathrm{0}}\mathrm{\ll }\mathrm{1}\end{array}\begin{array}{}\lim \\ {\mathit{S}}_{\mathrm{0}}\mathrm{\ll }\mathrm{1}\end{array}\left[\mathit{...}\right]$ ). Hence, we set: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{+}{\mathit{S}}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{+}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{2}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{3}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{+}{\mathit{S}}_{\mathrm{0}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{+}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{2}}\end{array}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{3}}}^{\mathrm{)}}\mathit{,}\end{array}$(F.1) where we have used for convenience the same formalism as for the expansion in η₀ – see Eq. (13) (noting of course that $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{v}\end{array}$ represents different functions). An important point is that the lemma of Appendix D holds when substituting S₀ to η₀. Therefore, substituting Eq. (F.1) into Eq. (11) provides the equations of motion for different orders of $\mathrm{𝒪}{}^{\mathrm{(}}\mathit{S}_{\mathrm{0}}{}^{\mathrm{)}}$:

• Order $\mathrm{𝒪}{}^{\mathrm{\right(}}\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}}^{\mathrm{\left)}}$: Eq. (11) provides$\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$ expressions for the velocities: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{=}& \mathrm{0}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{=}& \sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}\mathit{.}\end{array}$(F.2) In this order of expansion, the azimuthal velocity corresponds to the sub-Keplerian velocity of the gas. There is no radial motion.

• Order $\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{=}& \mathrm{-}\mathrm{2}{\mathit{R}}^{\mathit{p}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}}}\right)\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{=}& \mathrm{0.}\end{array}$(F.3)

• Order $\mathrm{𝒪}{}^{\mathrm{(}}\mathit{S}_{\mathrm{0}}{}^{\mathrm{)}}$: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{2}}\end{array}\mathrm{=}& \mathrm{0}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{2}}\end{array}\mathrm{=}& {\mathit{R}}^{\mathit{p}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}}[{\mathit{\eta }}_{\mathrm{0}}\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\\ & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{-}\frac{\mathrm{1}}{{\mathit{R}}^{\mathrm{2}}}\mathrm{+}{\mathit{\eta }}_{\mathrm{0}}\mathit{q}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{-}\mathrm{1}}}{\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}]\mathit{.}\end{array}$(F.4)

Finally, we obtain expressions for $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{r}}\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}$: $\begin{array}{ccc}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$We now compare the NSH86 expansion at R^pS₀ ≪ 1 (A-mode) and the W77 expansion at η₀ ≪ 1.

• NSH86: From Eqs. (16) and Eq. (E.7):$\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}& \mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\\ \mathrm{=}& \mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}& \sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{2}}\frac{{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{1}\mathrm{+}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{2}\mathit{p}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\end{array}$

• W77 small grains: From Eq. (F.5): $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}& \mathrm{-}\mathrm{2}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}}}\right)\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\\ \mathrm{=}& \mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}& \sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}\mathrm{+}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathit{p}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}}[{\mathit{\eta }}_{\mathrm{0}}\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\end{array}$

Clearly, Eqs. (F.6) and (F.7) are identical, demonstrating that the theories of W77 and NSH86 are consistent. We also note that if the simplification of Eq. (15) is not performed, the two W77 expansions directly appear as the expansion of NSH86 in $\mathrm{𝒪}{}^{\mathrm{(}}\mathit{S}_{\mathrm{0}}{}^{\mathrm{)}}$ or $\mathrm{𝒪}{}^{\mathrm{\right(}}{\mathit{S}}_{\mathrm{0}}^{-1}^{\mathrm{\left)}}$.

Now, in the case of large grains, we perform a perturbative expansion of the radial equation of motion with respect to $\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}}$ while assuming that S₀R^p ≫ 1, and verify that taking the limit at small η₀ provides the expression found for the B-mode in the NSH86 expansion.

Taking the same precautions as for the previous expansions, we write: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}& \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{+}\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}& \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{+}\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{.}\end{array}$(F.8) Following the same method as for the small grain sizes expansion, we obtain:

• Order $\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{0}}\end{array}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{0}}\end{array}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathit{.}\end{array}$(F.9) It this order of expansion, the azimuthal velocity of the grain is the standard Keplerian velocity.

• Order $\mathrm{𝒪}{}^{\mathrm{\right(}}\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}}^{\mathrm{\left)}}$: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}\mathrm{1}}\end{array}\mathrm{=}\mathrm{-}\mathrm{2}\left(\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{-}\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}\right){\mathit{R}}^{\mathrm{-}\mathit{p}}\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }\mathrm{1}}\end{array}\mathrm{=}\mathrm{0.}\end{array}$(F.10)

The expansion at higher order is more complicated and will not be used for further developments. At the order $\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)$, we have for $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{r}}\end{array}$ and $\begin{array}{}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}$: $\{\begin{array}{c}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{2}}{{\mathit{S}}_{\mathrm{0}}}\left(\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{-}\sqrt{\frac{\mathrm{1}}{\mathit{R}}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}}}\right){\mathit{R}}^{\mathrm{-}\mathit{p}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{,}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{\theta }}\end{array}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{R}}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{.}\end{array}$(F.11) We now compare the expressions provided by the NSH86 expansion at ${\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\gg }\mathrm{1}$ (B-mode) and the W77 expansion at η₀ ≪ 1 for the radial velocity:

• NSH86: $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}& \mathrm{=}& \mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}{\mathit{R}}^{\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{1}\mathrm{+}{\mathit{R}}^{\mathrm{2}\mathit{p}}{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}\\ & \mathrm{=}& \mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{.}\\ & & \end{array}$(F.12)

• W77 large grains: $\begin{array}{ccc}\begin{array}{}\mathrm{˜}\\ {\mathit{v}}_{\mathit{r}}\end{array}& \mathrm{=}& \mathrm{-}\frac{\mathrm{2}}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{-}\mathit{q}\mathrm{+}\mathrm{1}}}\right)\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{,}\\ & \mathrm{=}& \mathrm{-}\frac{{\mathit{\eta }}_{\mathrm{0}}}{{\mathit{S}}_{\mathrm{0}}}{\mathit{R}}^{\mathrm{-}\mathit{p}\mathrm{-}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\mathrm{𝒪}\left(\frac{\mathrm{1}}{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathrm{+}\mathrm{𝒪}\mathrm{(}{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}\\ & & \end{array}$(F.13)

Once again, we show that the W77 and NSH86 theories are consistent.

Appendix G: Asymptotic radial behaviour of a single grain

Noting ${\mathit{R}}_{\mathrm{1}}{}^{\mathrm{\right(}}\mathit{T}^{\mathrm{\left)}}$ the position of a grain integrated directly from the equation of motion (Eq. (11)) and ${\mathit{R}}_{\mathrm{2}}{}^{\mathrm{\right(}}\mathit{T}^{\mathrm{\left)}}$ the position integrated from the NSH86 approximation (Eq. (34)), we highlight (Fig. G.1) that the discrepancy between the motion from the exact equations and its NSH86 approximation is negligible (the relative error is lower than 10^-3 for all the considered sizes). It is therefore justified to use the analytical results derived in Sect. 3 to interpret the grain behaviour.

Thus, from Eq. (34), we see that the time for a grain starting at R = 1 to reach some final radius R_f is minimized for an optimal grain size S_m,f given by ${\mathit{S}}_{\mathrm{m}\mathit{,}\mathrm{f}}\mathrm{=}\sqrt{\frac{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\times }{\mathit{k}}_{\mathrm{f}}\mathrm{(}{\mathit{R}}_{\mathrm{f}}\mathrm{)}}\hspace{0.17em}\mathit{,}$(G.1)with ${\mathit{k}}_{\mathrm{f}}\mathrm{\left(}{\mathit{R}}_{\mathrm{f}}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \frac{\left(\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{f}}^{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\right)}{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{R}}_{\mathrm{f}}^{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}^{\mathrm{)}}}& \mathrm{if} \mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\ne \mathrm{0}\\ \mathrm{-}\ln \mathrm{(}{\mathit{R}}_{\mathrm{f}}\mathrm{)}& \mathrm{if} \mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{=}\mathrm{0.}\end{array}$(G.2)As shown in Eq. (34), the outcome of the grain radial motion depends on the value of $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$:

• If $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\le }\mathrm{0}$: $\begin{array}{}\lim \\ \mathit{T}\mathrm{\to }\mathrm{+}\mathrm{\infty }\end{array}\mathit{R}\mathrm{=}\mathrm{0.}$(G.3)For such disc profiles, all grains pile-up and fall onto the central star in an infinite time. Indeed, the surface density profile given by $\mathit{p}\mathrm{\ge }\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$ is steep enough to counterbalance the increase of the acceleration due to the pressure gradient. Therefore, grains fall onto the central star in an infinite time, whatever their initial size. Such an evolution happens because the grains always end migrating in the A-mode when they reach the disc’s inner regions. A crucial consequence is that grains are not depleted on the central star and therefore stay in the disc where they can potentially form planet embryos.

• If $\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{>}\mathrm{0}$: $\begin{array}{}\lim \\ \mathit{T}\mathrm{\to }{\mathit{T}}_{\mathrm{m}}\end{array}\mathit{R}\mathrm{=}\mathrm{0}\mathit{,}$(G.4)where ${\mathit{T}}_{\mathrm{m}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{0}}{\mathit{S}}_{\mathrm{0}}}\left(\frac{\mathrm{1}}{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{+}\frac{{\mathit{S}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}\right)\mathit{.}$(G.5)In this case, grains fall onto the central star in a finite time. The surface density profile given by $\mathit{p}\mathit{<}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}$ is now too flat to counterbalance the increasing acceleration due to pressure gradient. We note that:

  • For small sizes(S₀ ≪ 1), .

  • For large sizes (S₀ ≫ 1), ${\mathit{T}}_{\mathrm{m}}\mathrm{=}\mathrm{𝒪}\left(\frac{{\mathit{S}}_{\mathrm{0}}}{{\mathit{\eta }}_{\mathrm{0}}}\right)$.

  • T_m reaches a minimal value for a size S_m given by ${\mathit{S}}_{\mathrm{m}}\mathrm{=}\sqrt{\frac{\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}}}\mathit{.}$(G.6)Therefore ${\mathit{T}}_{\mathrm{m}}\mathrm{\left(}{\mathit{S}}_{\mathrm{m}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{{\mathit{\eta }}_{\mathrm{0}}\sqrt{\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}\mathrm{(}\mathrm{-}\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}}}\mathit{.}$(G.7)S_m is of order unity and corresponds to an optimal size of migration. Values of S_m and ${\mathit{\eta }}_{\mathrm{0}}{\mathit{T}}_{\mathrm{m}}{}^{\mathrm{(}}\mathit{S}_{\mathrm{m}}{}^{\mathrm{)}}$ in the (p,q) plane are shown in Fig. G.2. When S ≃ S_m, both the A- and B-modes contribute in an optimal way to the grains radial motion.

  In this case, grains can be efficiently accreted by the central star if $\mathit{S}\mathrm{\simeq }{\mathit{S}}_{\mathrm{m}}\mathrm{=}\mathrm{𝒪}{}^{\mathrm{\right(}}\mathrm{1}^{\mathrm{\left)}}$. Thus, they can not contribute to the formation of pre-planetesimals. This process is called the radial-drift barrier for planet formation.

---

Acknowledgments

This research was partially supported by the Programme National de Physique Stellaire and the Programme National de Planétologie of CNRS/INSU, France, and the Agence Nationale de la Recherche (ANR) of France through contract ANR-07-BLAN-0221. The authors want to thank C. Terquem, L. Fouchet, S. Arena and E. Crespe for useful comments and discussions. We also thank the referee for greatly improving the quality of this work by suggesting we include the important Stokes regime and present real numbers (time scales and grain sizes) for our criteria.

References

Appendix A: Notations

The notations and conventions used throughout this paper are summarized in Table A.1.

All Tables

All Figures