Guillaume Laibe
email: guillaume.laibe@monash.edu

Jean-François Gonzalez
email: francois.gonzalez]@ens-lyon.fr

Sarah T Maddison
email: smaddison@swin.edu.au

J.-F Gonzalez

Revisiting the "radial-drift barrier" of planet formation and its relevance in observed protoplanetary discs

Keywords: planetary systems, protoplanetary discs -methods, analytical G. Laibe et al

Context. To form metre-sized pre-planetesimals in protoplanetary discs, growing grains have to decouple from the gas before they are accreted onto the central star during their phase of fast radial migration and thus overcome the so-called "radial-drift barrier" (often inaccurately referred to as the "metre-size barrier"). Aims. To predict the outcome of the radial motion of dust grains in protoplanetary discs whose surface density and temperature follow power-law profiles, with exponent p and q respectively. We investigate both the Epstein and the Stokes drag regimes which govern the motion of the dust. Methods. We analytically integrate the equations of motion obtained from perturbation analysis. We compare these results with those from direct numerical integration of the equations of motion. Then, using data from observed discs, we predict the fate of dust grains in real discs. Results. When a dust grain reaches the inner regions of the disc, the acceleration due to the increase of the pressure gradient is counterbalanced by the increase of the gas drag. We find that most grains in the Epstein (resp. the Stokes) regime survive their radial migration if -p + q + 1 2 ≤ 0 (resp. if q ≤ 2 3 ). The majority of observed discs satisfies both -p + q + 1 2 ≤ 0 and q ≤ 2 3 : a large fraction of both their small and large grains remain in the disc, for them the radial drift barrier does not exist.

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 , T ∝ r -q ), 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;[START_REF] Crida | [END_REF], 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 1 ): 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 radialdrift 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 radialdrift 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) here 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 nongrowing grains in the Epstein regime. We expand the radial mo-tion 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.

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

             F D = - m d t s ∆v t s = ρ d s ρ g c s , (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 [START_REF] Whipple | From Plasma to Planet[END_REF]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 R g = 2s|∆v| ν , where ν is the microscopic kinematic viscosity of the gas, the drag force takes the following expression:

F D = - 1 2 C D πs 2 ρ g |∆v| ∆v, (2)
where the drag coefficient C D is given by

C D =                  24R -1 g for R g < 1 24R -0.6 g for 1 < R g < 800 0.44 for 800 < R g .
(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 (e r , e θ , e z ). 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 [START_REF] Whipple | From Plasma to Planet[END_REF], 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

m d dv dt = -F D + m d g, (4)
where F D is the drag force. As shown by Eqs. ( 1)-( 2), a general expression of the ratio F D m d is of the form

F D m d = - C (r, z) s y |v -v g | λ v -v g , (5)
where the quantities C, y and λ are defined for both the Epstein and the Stokes regime in Appendix C.

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

Considering the Epstein (small grains) regime, Eq. ( 4) reduces to:

m d dv dt = - m d t s v -v g + m d g. (6)
Writing Eq. ( 6) in (r, θ, z) coordinates leads to

                               dv r dt - v 2 θ r + v r -v gr t s + GMr z 2 + r 2 3/2 = 0 dv θ dt + v r v θ r + v θ -v gθ t s = 0 dv z dt + v z -v gz t s + GMz z 2 + r 2 3/2 = 0. ( 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 t s t k of the two timescales related to the physical processes acting on the grain is given by

t s t k = t s0 t k0 R p e Z 2 2R 3-q = S 0 R p e Z 2 2R 3-q . ( 8
)
With Eq. ( 1) and noting Ω k the Keplerian angular velocity, this ratio can be written as

t s t k = s ρ g c s ρ d Ω k = s s opt = S , (9)
where s opt = ρ g c s ρ d Ω k . This timescale ratio therefore corresponds to a dimensionless size S = S 0 R p e Z 2 2R 3-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 ):

                                           dṽ r dT - ṽ2 θ R + ṽr S 0 R -(p+ 3 2 ) e -Z 2 2R 3-q + R R 2 + φ 0 Z 2 3/2 = 0 dṽ θ dT + ṽθ ṽr R +          ṽθ -1 R -η 0 R -q -q         1 R - 1 √ R 2 +φ 0 Z 2                  S 0 R -(p+ 3 2 ) e -Z 2 2R 3-q = 0 d 2 Z dT 2 + 1 S 0 dZ dT R -(p+ 3 2 ) e -Z 2 2R 3-q + Z R 2 + φ 0 Z 2 3/2 = 0. ( 10
)
These equations depend on five control parameters: the initial dimensionless grain size, S 0 , the radial surface density and temperature exponents, p and q, the square of the disc aspect ratio, φ 0 = (H 0 /r 0 ) 2 , and the subkeplerian parameter, η 0 , 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 d 2 Z dT 2 simplifies to the damped harmonic oscillator equation. If we consider only the radial motion (for a 2D disc), we have Z = 0, and

                   dṽ r dT = ṽ2 θ R - ṽr S 0 R -(p+ 3 2 ) - 1 R 2 dṽ θ dT = - ṽθ ṽr R - ṽθ -1 R -η 0 R -q S 0 R -(p+ 3 2 ) . (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 0 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 η 0 (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 η 0 ≪ 1. We consider that this inequality also implies that

η 0 R -q ≪ 1 R . ( 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 Fig. 1. Principle of the various perturbative expansions for the grain radial motion. Expanding first with respect to the small pressure gradient (η 0 R -q+1 ≪ 1) leads to NSH86 equations. Expanding first with respect to the grain sizes (S 0 R p ≪ 1 or S 0 R p ≫ 1) leads to W77 expressions for the particular case p = 0. Combining both leads to the A-and Bmode, respectively for small and large grains.

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

       ṽr = ṽr0 + η 0 ṽr1 + O η 2 0 ṽθ = ṽθ0 + η 0 ṽθ1 + O η 2 0 (13)
and look at the orders O (1), O (η 0 ),... of the expansion -see Appendix E. We find that:

ṽr = η 0 ṽr1 + O η 2 0 = - 2S 0 R p-1 2 1 -1 -η 0 R -q+1 1 + R 2p S 2 0 + O η 2 0 .
(14) The pressure gradient term has been retained to keep the generality, however since we assume that η 0 ≪ 1, we equivalently have

1 R -η 0 R -q - 1 R = - η 0 2 R -q+ 1 2 + O η 2 0 . (15)
Thus, to order O (η 0 ),

ṽr = - η 0 S 0 R p-q+ 1 2 1 + R 2p S 2 0 + O η 2 0 , (16)
or equivalently, using Eqs. (C.4) and (C.1),

v r = rc 2 s v k d ln P dr (t s /t k ) 1 + (t s /t k ) 2 . ( 17
)
S 0 R p is the dimensionless expression of the ratio t s /t k . Eq. ( 16) shows that R 2p S 2 0 ≪ 1 or R 2p S 2 0 ≫ 1, and thus R p S 0 ≪ 1 or R p S 0 ≫ 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).

A-mode (Radial differential migration)

The A-mode corresponds to the regime R p S 0 ≪ 1 (or equivalently t s /t k ≪ 1). In the A-mode, Eq. ( 16) reduces to

ṽr = dr dt = dR dT = -η 0 S 0 R p-q+ 1 2 , ( 18
)
or equivalently

v r = rc 2 s v k d ln P dr t s t k , ( 19
)
where the O (R p S 0 ) 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 η 0 ≪ 1 (rigorously for η 0 R -q+1 ≪ 1) and second from a perturbative expansion for S 0 ≪ 1 (rigorously for S 0 R p ≪ 1). Formally, we have performed: lim

S 0 ≪1 lim η 0 ≪1 [...].
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 0 ≪ 1 (this approximation also implies that

R 2p S 2 0 ≪ 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 -p + q + 1 2 0:                    R = 1 --p + q + 1 2 η 0 S 0 T 1 -p+q+ 1 2 T = 1 -R -p+q+ 1 2 -p + q + 1 2 η 0 S 0 . ( 20
) -If -p + q + 1 2 = 0:          R = e -η 0 S 0 T T = - ln (R) η 0 S 0 . ( 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 -p + q + 1 2 ≤ 0: then as T (R) ∼ -R -p+q+ 1 2 -p + q + 1 2 η 0 S 0 = O R -p+q+ 1 2 , ( 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.

lim

T →+∞ R = 0. ( 23
)
-If -p + q + 1 2 > 0: the grain is accreted onto the central star in a finite migration time given by

T m = 1 η 0 S 0 -p + q + 1 2 , ( 24
)
which increases as S 0 and η 0 decrease, so that lim

T →T m R = 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.

B-mode (Drift forced by a resistive torque)

Returning to Eq. ( 16), the B-mode corresponds to the other asymptotic regime, where R p S 0 ≫ 1 (or equivalently t s /t k ≫ 1).

In this case,

ṽr = dR dT = - η 0 S 0 R -p-q+ 1 2 , ( 26
)
or equivalently

v r = rc 2 s v k d ln P dr t k t s . (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 (l ∝ √ 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 (S 0 R p ) -2 and then from an expansion in η 0 . Historically, W77 found an expression while only assuming that S 0 ≫ 1 since he considered a flat density profile. To find the expression derived by W77 for large grains, we must assume that S 0 R p ≫ 1. It is crucial to see that this expression does not imply that S 0 ≫ 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 p + q + 1 2 0:                    R = 1 -p + q + 1 2 η 0 S 0 T 1 p+q+ 1 2 T = S 0 η 0 1 -R p+q+ 1 2 p + q + 1 2 . ( 28
) -If p + q + 1 2 = 0:          R = e - η 0 S 0 T T = - S 0 η 0 ln (R) . ( 29
)
As for the A-mode, two kinds of behaviours appear, depending on the p and q exponents:

-If p + q + 1 2 ≤ 0:
The grain migrates inwards, piles-up in the disc's inner regions and falls onto the star in an infinite time:

lim T →+∞ R = 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 p + q + 1 2 > 0:
The grain is accreted onto the central star in a finite time

T m = S 0 η 0 p + q + 1 2 , (31)
which increases as S 0 increases and η 0 decreases and so that lim

T →T m R = 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, p + q + 1 2 > 0. Grains migrating in B-mode in such discs fall in a finite time onto the central star.

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

dR dT = -η 0 S 0 R p-q+ 1 2 1 + R 2p S 2 0 . ( 33
)
We can separate the R and T variables and integrate to obtain an expression for T (R): where

T = 1 η 0 S 0 [T A (R) + T B (R)] , (34)
T A (R) =            1 -R -p+q+ 1 2 -p + q + 1 2 if -p + q + 1 2 0 -ln (R) if -p + q + 1 2 = 0. (35)
and

T B (R) =            S 2 0 1 -R p+q+ 1 2 p + q + 1 2 if p + q + 1 2 0 -S 2 0 ln (R) if p + q + 1 2 = 0. ( 36
)
Eq. ( 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 R ≪ S -1/p 0 . 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 0 and not powers of S 0 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 η 0 , S 0 , p, q. We set η 0 = 10 -2 to mimic a realistic disc and vary the order of magnitude of S 0 from 10 -4 to 10 2 for two sets of (p,q) values. First, we choose (p = 0, q = 3 4 ); according to the NSH86 expansion, the grains are accreted by the central star in a time T m . Second, we set (p = 3 2 , q = 3 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 -p + q + 1 2 can take a positive or a negative value. Consequently, we interpret the radial motion R (T ) of dust grains plotted in Fig. 3.

-The top panel of Fig. 3 shows the results for (p = 0, q = 3 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 . S 0 varies from 10 -4 to 10 2 . Top: p = 0, q = 3 4 , here -p + q + 1 2 > 0 and the grain is accreted onto the central star in a finite time. Bottom: p = 3 2 , q = 3 4 , here -p + q + 1 2 < 0 and the grain piles up and is consequently accreted onto the central star in an infinite time.

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/η 0 , or T m = 160 with η 0 = 10 -2 .

-The bottom panel of Fig. 3 shows the results for (p = 3 2 , q = 3 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 S m,f ≃ 2.9 = O (1) (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.

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:

µ = 5m √ π 64σ s k B T m , (37)
where m = 2m H = 3.347446922 × 10 -27 kg and σ s = 2.367 × 10 -19 m 2 is the molecular cross section of the molecule [START_REF] Chapman | The mathematical theory of non-uniform gases[END_REF]. The kinematic viscosity ν is then defined by µ = ρ g ν and the gas collisional mean free path is given by

λ g = π 2 ν c s . ( 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 where c s0 is given by Eq. (B.23). Thus, the expression of the kinematic viscosity ν is

µ = µ 0 R -q 2 , ( 39
)
µ 0 = 5m √ π 64σ s c s0 , (40)
ν = ν 0 R 3 2 +p-q e Z 2 2R 3-q , (41)
ν 0 = µ 0 ρ g0 . (42)
First, if R g < 1, the drag force is linear in vv g and thus has the same structure as for the Epstein regime. Comparing the expressions of C (R, 0) 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 p ′ = q-3 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))

dR dT = -η 0 S 2 0 R -q 2 -1 1 + R q-3 S 4 0 . (43)
These crucial results follow:

-In the A-mode (S 0 ≪ R 3-q 4 ), grains experience a pile-up and migrate onto the central star in an infinite time if -p ′ + q ′ + 1 2 ≤ 0, i.e. if q ≤ -4 (which never occurs in real discs). -In the B-mode (S 0 ≫ R 3-q 4 ), grains migrate onto the central star in an infinite time if p ′ + q ′ + 1 2 ≤ 0, i.e. if q ≤ 2 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 S 4 3-q 0 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 q ≤ 2 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 vv 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, (η 0 R -p ) j vv g → 0 at the limit η 0 → 0. This means that both v r and v θ -√ 1/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.

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 -p+q+ 1 2 < 0 (resp. q < 2 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(Cuzzi et al. , 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 radialdrift 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 -p + q + 1 2 = 0 and q = 2 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 -p + q + 1 2 ≤ 0 retain their small grain population and that (ii) discs satisfying q ≤ 2 3 keep their large solids. On the contrary, discs for which -p + q + 1 2 > 0 (resp. q > 2 3 ) likely lose their small (resp. large) particles.

Application to observed discs and planet formation

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 3 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 λ g = 4s 9 , or

s = 9 4 λ g = 45π 3 2 256 mH ′ 0 σ s Σ ′ 0 r p-q 2 + 3 2 , ( 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 6 yr. For a grain starting at a distance r 0 from a 1 M ⊙ star, the dimensionless value T disc is therefore

T disc = t disc t k0 = √ GM ⊙ t disc r 3/2 0 ∼ 6 × 10 6 [r 0 (AU)] 3/2 . ( 45
)
The dimensionless value of the dust disc inner radius (r in ∼ 0.01 AU) for a grain starting at r 0 is

R in = r in r 0 ∼ 10 -2 r 0 (AU) . ( 46
)
The link between dimensionless and real grain sizes is made through the optimal size for radial migration. Considering first Fig. 6. Values of S 0 as a function of grain size s and initial position r 0 for a disc of mass M disc = 0.01 M ⊙ extending from r in = 10 -2 AU to r out = 10 3 AU, with grains of intrinsic density ρ d = 1000 kg m -3 and η 0 = 10 -2 , for p = 0 and q = 3 4 (top) and p = 3 2 and q = 3 4 (bottom). The thick line shows the limit between grains that are accreted onto the star (t in < t disc ) and those that survive in the disc (t in > t disc ), i.e. the survival limit, in the Epstein regime.

the Epstein regime, its midplane value has a radial dependence given by

s opt (r, 0) = Σ(r) √ 2π ρ d = Σ ′ 0 r -p √ 2π ρ d . (47)
The dimensionless size of a particle of size s starting its migration at position r 0 is therefore

S 0 = s s opt,0 = √ 2π ρ d Σ ′ 0 s r p 0 . (48)
Values of S 0 are plotted in the (r 0 , 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 3 AU, with grains of intrinsic density ρ d = 1000 kg m -3 and η 0 = 10 -2 , for both p = 0 and p = 3 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 0 and r 0 , whereas Eq. ( 45) gives an expression of T disc as a function of r 0 . Equating them yields a second order equation in S 0 as a function of r 0 , 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:

s = p + q + 1 2 2 √ 2π ρ d Σ ′ 0 r -p 0 1 -r in r 0 p+q+ 1 2        √ GM t disc η 0 r 3/2 0 ± GM t 2 disc η 2 0 r 3 0 -4 1 -r in r 0 -p+q+ 1 2 1 -r in r 0 p+q+ 1 2 -p + q + 1 2 p + q + 1 2                , ( 49
)
which is plotted as a thick line in Fig. 6. It separates the (r 0 , 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 -p + q + 1 2 > 0, illustrated by the case (p = 0, q = 3 4 ), 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 0 = 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 -p + q + 1 2 ≤ 0, illustrated by the case (p = 3 2 , q = 3 4 ), 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 0 > 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 0 , 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 1 4 to 3 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 -p + q + 1 2 < 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 p ′ = q-3 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 s opt,0 ∝ T 1 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 0 , 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 0 , 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 q = 3 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 5 yr contour for p = 0. Their survival time decreases to a few 10 4 yr for p ∼ 0.8, and increases again to values larger than 10 6 yr as p increases. At an initial position of a few hundred AU, 1 mm grains survive longer than 10 6 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 5 yr for p = 0, decreases to ∼ 10 2 yr for p ∼ 1, and increases again to ∼ 10 3 -10 4 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 0 , s) plane of grains reaching that radius r at any time t. Alternatively, one can plot isocontours of the initial position r 0 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, q = 3 4 ) and Fig. 10 for (p = 3 2 , q = 3 4 ). These plots make it easy to compare the radial evolution of any particle to any physical Fig. 9. Time evolution (in nine snapshots from t = 10 -2 to 10 6 yr) of isocontours of the initial position of grains in the (r, s) plane for the disc with p = 0 and q = 3 4 . The label for each contour can be deduced from its abscissa in the upper left panel at t = 10 -2 yr. 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 6 yr). No grains are found to the right of the r 0 = 10 3 AU contour, since this was the initial outer disc radius. In the disc with (p = 0, q = 3 4 ), 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 (p = 3 2 , q = 3 4 ), 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 -p + q + 1 2 . However, the transition between discs for which the radial-drift barrier occurs or not consists more of a continuum around the value -p + q + 1 2 = 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 -p + q + 1 2 = 0, using the figures shown above in this section.

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 irradi-Fig. 11. Histogramm of the q parameters obtained from Andrews & Williams (2005, 2007) data of 63 observed discs. The distribution is roughly comprised between 0.4 and 0.8 and, centred around 0.55. Approximately 90 % of the discs satisfy q < 2/3. ated disc with opacities from Beckwith et al. (1990), a gas-todust 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 3 2 . Very flat profiles with p < 1 2 and very steep profiles with p > 3 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: -p + q + 1 2 ≤ 0 and q ≤ 2 3 : both small and large grains experience the pile-up effect. Those discs are potentially observable and may favour planet formation. 2. region 2: -p + q + 1 2 ≤ 0 and q > 2 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: -p + q + 1 2 > 0 and q ≤ 2 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: -p + q + 1 2 > 0 and q > 2 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 [1; 3 2 ] 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 -p + q + 1 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 -p + q + 1 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 (q ≤ 2 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 preplanetesimals. 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 in-creased as the drag from the gas onto the planetesimals releases an additional thermal contribution.

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: Σ (r) = Σ ′ 0 r -p and T (r) = T ′ 0 r -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 (Amode, 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 -p + q + 1 2 > 0, the dust particle is accreted onto the central star in a finite time, and if -p + q + 1 2 ≤ 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 q ≤ 2 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 q ≤ 2 3 , allowing them to retain also their large preplanetesimals. 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 -p + q + 1 2 ≤ 0 and q ≤ 2 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, obser-vations tell us that grains do grow (Testi et al. 2003;[START_REF] Wilner | [END_REF]Apai et al. 2005;Lommen et al. 2007Lommen et al. , 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.

G. Laibe et al.: Revisiting the "radial-drift barrier" of planet formation and its relevance in observed protoplanetary discs and dividing both sides of Eq. (B.3) by c 2 s , we have:

∂ ln c 2 s ρ g ∂z = - GMz z 2 + r 2 3/2 c 2 s . (B.5)
Integrating Eq. (B.5) between 0 and z provides:

ρ g (r, z) = P (r) c 2 s (r, z) e - z 0 GMz ′ dz ′ r 2 + z ′2 3/2 c 2 s (r, z ′ ) . (B.6)
This expression can be simplified by the following approximations:

-In the special vertically isothermal case, where the sound speed depends only on the radial coordinate, Eq. (B.6) simplifies to:

ρ g (r, z) = ρg (r) e - GM c2 s (r) 1 r - 1 √ r 2 + z 2 . (B.7)
-Further, assuming a thin disc ( z r 2 ≪ 1), a Taylor series expansion of Eq. (B.7) leads to:

ρ g (r, z) = ρg (r) e -z 2 2H(r) 2 , (B.8) with: H (r) = rc s (r) v k (r) , (B.9)
which is the classical scale height for vertically isothermal thin discs.

B.2. Azimuthal velocity

The radial component of the Euler equations is given by:

- v 2 gθ r = - 1 ρ g ∂P ∂r - GMr z 2 + r 2 3/2 , (B.10)
where ρ g is given by Eq. (B.6). Thus:

ρ g (r, z) = P (r) c 2 s (r, z) e -I 1 (r,z) , (B.11) with:                I 1 (r, z) = z 0 GMc -2 s r, z ′ ∂ z ′ f r, z ′ dz ′ , f (r, z) = 1 r - 1 √ r 2 + z 2 . (B.12)
To simplify Eq. (B.10), we first use the following identity: (B.14)

1 ρ g ∂P ∂r = c 2 s ∂ ln c 2 s ρ g ∂r , (
Noting that f (r, z = 0) = 0 and integrating I 1 by parts provides: (B.15) and:

I 1 (r, z) = GMc -2 s f (r, z) - z 0 GM f r, z ′ ∂ z ′ c -2 s r, z ′ dz ′ ,
1 ρ g dP dr = c 2 s d ln P dr -       - GM r 2 + GMr z 2 + r 2 3/2 (B.16) +c 2 s GM f ∂ r c -2 s -c 2 s ∂ ∂r z 0 GM f ∂ z ′ c -2 s dz ′ .
Then, Eq. (B.10) becomes:

v 2 gθ r = c 2 s d ln P dr + GM r 2 +c 2 s ∂ ∂r z 0 GM f ∂ z ′ c -2 s dz ′ -c 2 s GM f ∂ r c -2 s . (B.17) Noting that : ∂ ∂r z 0 f ∂ z ′ c -2 s dz ′ = z 0 ∂ f ∂r ∂ z ′ c -2 s dz ′ + z 0 f ∂ r ∂ z ′ c -2 s dz ′ , (B.18
) and integrating the last term of the right hand side of Eq. (B.18) by parts provides:

z 0 f ∂ r ∂ z ′ c -2 s dz ′ = f ∂ r c -2 s - z 0 ∂ z ′ f ∂ r c -2 s dz ′ . (B.19)
Therefore, Eq. (B.17) reduces to:

v 2 gθ r = GM r 2 + c 2 s d ln P dr + GMc 2 s z 0 ∂ r f ∂ z ′ c -2 s -∂ z ′ f ∂ r c -2 s dz ′ ,
(B.20) which can be more elegantly written as:

v 2 gθ r = GM r 2 + c 2 s d ln P dr + GMc 2 s z 0 ∇ f × ∇c -2 s .e θ dz ′ . (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:

v 2 gθ r = GM r 2 + c 2 s d ln P dr -GM 1 r - 1 √ r 2 + z 2 ∂ r lnc -2 s . (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:

             Σ (r) = Σ 0 r r 0 -p = Σ ′ 0 r -p , T (r) = T 0 r r 0 -q = T ′ 0 r -q , c s (r) = c s0 r r 0 -q/2 = c ′ s0 r -q/2 . (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:

H (r) = H 0 r r 0 3 2 -q 2 = H ′ 0 r 3 2 -q 2 , (B.24)
with:

H ′ 0 = c ′ s0 √ GM . (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: Hence, with Σ ′ 0 = √ 2πρ ′ g0 H ′ 0 and x = p -q 2 + 3 2 ,

ρ g = ρ ′ g0 r -x e
ρ g = Σ ′ 0 √ 2πH ′ 0 r -(p-q 2 + 3
2 ) e

v gθ = GM r -p + q 2 + 3 2 c ′2 s r -q -GMq 1 r - 1 √ r 2 + z 2 .
(B.31) 2R 3-q .

(C.1) with:

η 0 = p + q 2 + 3 2 c s0 v k0 2 . (C.2)
The dimensionless parameter η 0 gives the order of magnitude of the relative discrepancy between the Keplerian motion and the gas azimuthal velocity. We note that:

                                               t t k0 = T v r v k0 = dR dT = ṽr v θ v k0 = R dθ dT = ṽθ z H 0 = Z v z H 0 /t k0 = dZ dT .
(C.4)

Writing the coefficient C (r, z) of the drag force of Eq. ( 5) as:

C (r, z) = C 0 C (R, Z) , (C.5) and using dimensionless coordinates, we have:

F D /m d v k0 /t k0 = - C (R, Z) s (v λ k0 t k0 C 0 )
1 y y |ṽṽg | λ ṽṽg .

(C.6)

We also introduce: Physically, s opt,0 corresponds to the grain size at which the drag stopping time equals to the Keplerian time at r 0 . In Table C.1, we give the expressions of y, λ, s opt,0 , C (R, Z) for the Epstein and the three Stokes drag regimes. The dimensionless equations of motion for a dust grain are then: This is illustrated plotting the radial motion of dust grains for S 0 = 10 -2 , η 0 = 10 -2 and for p = 0, q = 3 4 (solid) and p = 3 2 , q = 3 4 (dashed). Top: R 1 (T ), middle: R 2 (T ), bottom: relative difference (R 2 (T ) -R 1 (T )) /R 1 (T ).

                                          
(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

S m,f = p + q + 1 2 -p + q + 1 2 × k f (R f ) , (G.1) with k f (R f ) =                1-R -p+q+ 1 2 f 1-R p+q+ 1 2 f if p + q + 1 2 0 -ln (R f ) if p + q + 1 2 = 0. (G.2)
As shown in Eq. ( 34), the outcome of the grain radial motion depends on the value of -p + q + 1 2 : -If -p + q + 1 2 ≤ 0: lim

T →+∞ R = 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 p ≥ q + 1 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 -p + q + 1 2 > 0: lim

T →T m R = 0, (G.4)
where

T m = 1 η 0 S 0       1 -p + q + 1 2 + S 2 0 p + q + 1 2       . (G.5)
In this case, grains fall onto the central star in a finite time.

The surface density profile given by p < q + 1 2 is now too flat to counterbalance the increasing acceleration due to pressure gradient. We note that: -For small sizes (S 0 ≪ 1), T m = O (S 0 η 0 ). -For large sizes (S 0 ≫ 1), T m = O S 0 η 0 . -T m reaches a minimal value for a size S m given by S m = p + q + 1 2 -p + q + 1 2 .

(G.6) Therefore T m (S m ) = 2

η 0 p + q + 1 2 -p + q + 1 2 . (G.7)
S m is of order unity and corresponds to an optimal size of migration. Values of S m and η 0 T m (S m ) in the (p, q) plane are shown in Fig. G.2. When S ≃ S m , both the Aand 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 S ≃ S m = O (1). Thus, they can not contribute to the formation of pre-planetesimals. This process is called the radial-drift barrier for planet formation.

Fig. 2 .

 2 Fig. 2. Radial evolution of grains in the (R,S 0 ) plane showing that a grain in the Epstein drag regime ends its radial motion in the A-mode. The solid curves represent R -p for various values of p, they separate the A-mode (below) from the B-mode (above) regions. The horizontal dashed lines show trajectories of grains as they migrate inwards from R = 1. The shaded area is a forbidden zone.

Fig. 3 .

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

Fig. 4 .

 4 Fig. 4. Radial evolution of grains in the (R,S 0 ) plane showing that a grain in the linear Stokes drag regime ends its radial motion in the Bmode. The solid curves represent R 3-q 4for various values of q, they separate the A-mode (below) from the B-mode (above) regions. The horizontal dashed lines show trajectories of grains as they migrate inwards from R = 1. Shaded area: forbidden.

Fig. 5 .

 5 Fig.5. Transition between the Epstein and the Stokes regimes in a protoplanetary disc of M disc = 10 -2 M ⊙ extending from r in = 10 -2 AU to r out = 10 3 AU for several values of p and q. Grains with sizes below (resp. above) than the curve experience the Epstein (resp. Stokes) drag regime.

Fig. 7 .

 7 Fig. 7.Survival limits of grains for different values of p and q. Left: Epstein regime, right: Stokes regime. Grains to the left of the curves (t in < t disc ) are accreted onto the star whereas those to the right (t in > t disc ) survive in the disc.

Fig. 8 .

 8 Fig. 8. Isocontours of the survival time (i.e. the time needed to reach the disc inner edge at r in = 0.01 AU) of grains of size s and initial position r 0 for q = 3 4 and different values of p.

Fig. 10 .

 10 Fig.10. Same as Fig.9for p = 3 2 and q = 3 4 .

Fig. 12 .

 12 Fig. 12. Location of the different outcomes of radial migration in the (p,q) plane. Dashed (resp. dotted) line: limit between accretion without or with grains pile-up resulting in a finite or infinite time in the A-mode of the Epstein regime (resp. B-mode of the Stokes regime at small Reynolds numbers). Shaded area: location of observed discs. Black dots: discs used as examples in Sect. 3.3.

  correct surface density profile when integrated with respect to z. With this expression of ρ g , P (r) is given by:P (r) = c 2 s (r) ρ g (r, z = 0) = c ′2s0 Eq. (B.22), we find:

  | λ ṽṽg . (C.9)

  , Z) |ṽṽg | λ + Z R 2 + φ 0 Z 2 3/2 = 0.(C.10)

Fig. G. 2 .

 2 Fig. G.2. Values of S m (left) and η 0 T m (S m ) (right) in the (p,q) plane.

Table C .

 C 1. Expressions of the coefficients y, λ, s opt,0 and C (R, Z) for different drag regimes.

	Drag regime Epstein	y 1	λ 0			s opt,0 Σ 0 √ 2πρ d					C (R, Z) R -(p+ 3 2 ) e -Z 2 2R 3-q
	Stokes (R g < 1)	2	0			9t k0 µ 0 2ρ d						R -q 2
	Stokes (1 < R g < 800) 1.6 0.4	      	v 0.6 k0	18r 0 ν 0.6 0 Σ 0 √ 2π2 1.6 ρ d H 0	      	1 1.6	R	-	2p 5 +	q 10 + 3 5 e -2 5	Z 2 2R 3-q
	Stokes (800 < R g )	1	1			1.32r 0 Σ 0 8ρ d √ 2πH 0			R -(p-q 2 + 3 2 ) e -Z 2 2R 3-q

http://exoplanet.eu

q 2 , z/H 0 = Z,

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

Symbol Meaning M

Mass of the central star g

Gravity field of the central star r 0

Initial distance to the central star ρ g Gas density ρg (r) ρ g (r, z = 0) c s Gas sound speed cs (r) c s (r, z = 0) c s0

Gas sound speed at r 0 T Dimensionless time T

Gas temperature (T 0 : value at r 0 ) Σ 0

Gas surface density at r 0 p Radial surface density exponent q Radial temperature exponent P Gas pressure v k Keplerian velocity at r v k0

Keplerian velocity at r 0 H 0 Gas scale height at r 0 φ 0 Square of the aspect ratio H 0 /r 0 at r 0 η 0 Sub-Keplerian parameter at r 0 s

Grain

Appendix A: Notations

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

Appendix B: disc structure

B.1. Hydrostatic equilibrium

At stationary equilibrium ( ∂ ∂t = 0), gas velocities v gr , v gθ , v gz and the gas density ρ g obey mass conservation and the Euler equation:

The solution of Eq. (B.1) requires:

Appendix C: Dimensionless quantities and equations of motion

To highlight the important physical parameters involved, we set v k0 = GM r 0 and introduce dimensionless quantities given by the following expressions:

Appendix D: Lemma for the different expansions

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

-ṽr0 = 0, and -ṽxi can be written as a function of

Appendix E: Epstein regime: perturbation analysis -Order O (1): At this order of expansion, η 0 R -q is negligible compared to 1 R . Thus, substituting Eq. ( 13) into Eq. ( 11) provides

At this stage, we do not know the order of dṽ θ0 dT . We show in the lemma of Appendix D that dṽ θ0 dT = O (η 0 ). Applying this lemma, we see that at the order O (1), taking ṽr0 = 0 and ṽθ0 = 1 R (which ensures that dṽ θ0 dT = O (η 0 )) is a relevant solution for the equations of motion (which corresponds to circular Keplerian motion). Thus,

-Order O (η 0 ): Applying the lemma in this order of expansion and noting that

Eq. ( 11) becomes

. (E.4) Solving the linear system Eq. (E.4) for (ṽ r1 , ṽθ1 ) provides

In addition to the expression of ṽr given in Sect. 3, we also note that

which provides with Eq. ( 15)

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 0 variable. We verify that taking the limit at small η 0 provides the expression found for the Amode in the NSH86 expansion. (Formally, we will show that lim

). Hence, we set:

where we have used for convenience the same formalism as for the expansion in η 0 -see Eq. ( 13) (noting of course that ṽ represents different functions). An important point is that the lemma of Appendix D holds when substituting S 0 to η 0 . Therefore, substituting Eq. (F.1) into Eq. ( 11) provides the equations of motion for different orders of O (S 0 ):

-Order O 1 S 0 : Eq. ( 11) provides O (1) expressions for the velocities:

In this order of expansion, the azimuthal velocity corresponds to the sub-Keplerian velocity of the gas. There is no radial motion. -Order O (1):

Finally, we obtain expressions for ṽr and ṽθ :

(F.5)

We now compare the NSH86 expansion at R p S 0 ≪ 1 (Amode) and the W77 expansion at η 0 ≪ 1.

-NSH86: From Eqs. ( 16) and Eq. (E.7):

(F.6) -W77 small grains: From Eq. (F.5):

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 O (S 0 ) or O S -1 0 . Now, in the case of large grains, we perform a perturbative expansion of the radial equation of motion with respect to 1 S 0 while assuming that S 0 R p ≫ 1, and verify that taking the limit at small η 0 provides the expression found for the B-mode in the NSH86 expansion.

Taking the same precautions as for the previous expansions, we write:

Following the same method as for the small grain sizes expansion, we obtain:

It this order of expansion, the azimuthal velocity of the grain is the standard Keplerian velocity. -Order O 1 S 0 :

(F.10)

The expansion at higher order is more complicated and will not be used for further developments. At the order O 1 S 2 0 , we have for ṽr and ṽθ :

Appendix G: Asymptotic radial behaviour of a single grain

Noting R 1 (T ) the position of a grain integrated directly from the equation of motion (Eq. ( 11)) and R 2 (T ) the position integrated from the NSH86 approximation (Eq. ( 34)), we highlight


C Poussot-Vassal

O Sename

♯ , P. Gáspár ♮ , Szabó ♮ L Dugard

J Bokor

Péter Gáspár
email: gaspar@sztaki.hu

The design of a chassis system based on multi-objective qLPV control

In this paper we compare LTI and qLPV H ∞ / H 2 controllers. The Pareto limit is used to show the compromise that has to be done when a mixed synthesis is achieved. Simulations on a nonlinear half vehicle model, with multiple objectives, are performed to show the efficiency of the method.

INTRODUCTION

The main role of suspensions is to improve comfort by isolating the vehicle chassis to an uneven ground and to provide a good road holding to ensure passenger safety. Suspension control of quarter vehicle have been widely explored the past few years to improve vertical movements either by applying LQ [START_REF] Hrovat | Survey of advanced suspension developments and related optimal control application[END_REF], Skyhook [START_REF] Sammier | Skyhook and H ∞ control of active vehicle suspensome practical aspects[END_REF], H ∞ control [START_REF] Gaspar | Active suspension design using LPV control[END_REF][START_REF] Zin | An LPV/H ∞ active suspension control for global chassis technology: Design and performance analysis[END_REF], LPV [START_REF] Fialho | Road adaptive active suspension design using linear parameter varying gain scheduling[END_REF] or mixed synthesis [START_REF] Abdellahi | On the design of active suspension system by H ∞ and mixed H 2 /H ∞ : An LMI approach[END_REF][START_REF] Takahashi | A multiobjective approach for H 2 and H ∞ active suspension control[END_REF]. Roll dynamic is catched by the half vehicle model and is directly linked to suspension behavior. Separated synthesis on each suspension can't guarantee global performances. The aim of the mixed H ∞ / H 2 control synthesis is to treat the standard H 2 and H ∞ optimal control problems as separate problems but in a unified state-space framework. This method yields a compensator that combines the H 2 quadratic performance criterion for disturbance rejection with the H ∞ performance criterion for maximum robustness against destabilizing uncertainties. The controller which minimizes the H 2 performance index is selected from the suitable H ∞ controllers, thus the desired criteria are met by creating a balance between H 2 and H ∞ norms [START_REF] Doyle | Mixed H 2 and H ∞ performance objectives: Optimal control[END_REF].

The mixed qLPV H ∞ / H 2 method is proposed here for the design of active suspension system, in which different optimization criteria are applied to guarantee the performance specifications and the nonlinearity of the suspension system. The nonlinearity in the suspension system is caused by the changes in the spring and damping coefficients. It is assumed that the nonlinear dynamics of road vehicles is approximated by LPV (qLPV) models, in which nonlinear terms are hidden with newly defined scheduling variables and they are available from calculated signals.

The active suspension based on the LPV model takes the nonlinear dynamics of the system into consideration. Performance limitations according to the importance given between the H ∞ and the H 2 criteria is shown with the Pareto limit.

The paper is organized as follows: in Section 2 we introduce a linear and a nonlinear model of the half vehicle. LTI and qLPV polytopic H ∞ / H 2 control, based on LMIs, are presented in Section 3. The Pareto limit, applied to the half vehicle provides smart indications in the way to choose H ∞ / H 2 attenuation parameters according to the desired performances. In Section 4, validations are done on the nonlinear model presented.

Vehicle model of suspension systems

Roll dynamic is the main movement that enters when a driver turns. The half vehicle model involved here is a chassis model that catches vertical and roll dynamics [START_REF] Sammier | Skyhook and H ∞ control of active vehicle suspensome practical aspects[END_REF] (Figure 1). It models the left/right vehicle load transfers that appears during a steering situation. The model is composed of two suspensions, each of them modeled by a spring (F k {f l,f r} ), a damper (u {f l,f r} = F c {f l,f r} , in the passive case) or an actuator (u {f l,f r} = u {f l,f r} H∞/H 2 + c {f l,f r} , in the active case) linked to a common suspended mass (m s ) and to a specific unsprung mass m us f l and m us f r . Tires k t f l and k t f r are linked to the ground and to the unsprung mass m us f l and m us f r respectively. The movements taken into account are the vertical displacement of the suspended mass (z s ), the unsprung masses (z us {f l,f r} ), the suspension deflections (z def {f l,f r} ) and the roll angle (θ) of the center of gravity of the suspended mass. The model is obtained by simply adding two suspensions and tires equations with the dynamical equation of the chassis as follows. First we derive the suspension (F sz i ) and tire (F tz i ) forces,

F k f l F kt f l z s z r f l m us f l m us f r z us f l z r f r z us f r m s , I x k t f r k f r θ 7 7 - t f t f - u f r u f l
         F tz f l = k t f l (z us f l -z r f l ) F tz f r = k t f r (z us f r -z r f r ) F sz f l = F k f l (z s f l -z us f l ) -k b θ 2l f l + u f l F sz f r = F k f r (z s f r -z us f r ) + k b θ
2l f r + u f r then the dynamic of the chassis and unsprung masses (bounce and roll) are given by

       m s zs = -(F sz f l + F tz f r + F dz ) m us f l zus f l = F sz f l -F tz f l m us f r zus f r = F sz f r -F tz f r I x θ = F sz f l t f -F tz f r t f + M dx (1)
where F k i and F c i , {i = f l, f r}, represent the force delivered by the spring and by the damper (either linear or nonlinear), k t i is the stiffness of the tire and k b models the influence of an antiroll bar. I x is the chassis inertia on the roll axis, t f is the distances of the unsprung masses to the center of gravity of the suspended mass. Finally, θ, z s , z us f l and z us f r represent the roll angle and the chassis, unsprung mass left and right bounce. Then z r f l and z r f r represent the road disturbances on the wheels. F dz , M dx represent the load and inertia disturbances. Note that when the passive system is considered,

u f r = F c f r and u f l = F c f l .
Then, the state space vector of the linear model is defined by

x = z us f l żus f l z us f r żus f r z s żs θ θ ,
the input are given by w = z r f l z r f r F dz M dx u f l u f r and the measured signal used for control y = z def f l z def f r .

3 MIXED H ∞ / H 2 LMI BASED SYNTHESIS

A LTI multi-objective controller

The multi-objective synthesis consists of giving different kind of constraints on the output of a system. With this formulation (for the case of H ∞ / H 2 ), let describe the system as follows

    ẋ z ∞ z 2 y     =     A B ∞ B 2 B C ∞ D ∞1 D ∞2 E ∞ C 2 D 21 D 22 E 2 C F ∞ F 2 0         x w ∞ w 2 u     (2)
the controller,

ẋc u = A c B c C c D c x y = S (3)
and the closed loop,

  ẋ z ∞ z 2   =   A B ∞ B 2 C ∞ D ∞1 D ∞2 C 2 D 21 0     x w ∞ w 2   (4)
The

H ∞ / H 2 synthesis consists of, imposing T ∞ = ||z ∞ /w ∞ || ∞ < γ ∞ and T 2 = ||z 2 /w 2 || 2 < γ 2 .
Hence the LMI based problem formulation is the following: minimize γ 2 and γ ∞ subject to K and Z. [2, 9]

                         A T K + KA KB ∞ C T ∞ B T ∞ K -γ 2 ∞ I D T ∞1 C ∞ D ∞1 -I   < 0 A T K + KA KB 2 B T 2 K -I < 0 , K C T 2 C 2 Z > 0 , T race(Z) < γ 2 , D 22 = 0 (5)
Then solving this problem gives the LTI controller that achieves the desired performances. Note that to relax BMIs (5) into LMIs we use the transformation given in [START_REF] Scherer | Multiobjective output-feedback control via LMI optimization[END_REF].

A qLPV multi-objective controller

Linear parameter varying theory is useful to tackle measurable and bounded nonlinearities. We talk about qLPV when the varying parameters only enter in the dynamic matrix A of the system. In the suspension system, the measure of the deflection (used as a controller input) can also be used to reconstruct the stiffness coefficient [START_REF] Zin | An LPV/H ∞ active suspension control for global chassis technology: Design and performance analysis[END_REF]. To build a qLPV controller, we use the polytopic approach which consists of building a controller to the k-corners of the polytope (formed by all the possible combinations of the upper and lower bounds of each varying parameters) and to schedule theses k-controllers by the measure of the varying variables. qLPV system is described as follows, with p a varying parameter,

   ẋ z ∞ z 2 y     =     A(p) B ∞ B 2 B C ∞ D ∞1 D ∞2 E ∞ C 2 D 21 D 22 E 2 C F ∞ F 2 0         x w ∞ w 2 u     (6)
the parameter dependent controller,

ẋc u = A c (p) B c (p) C c (p) D c (p) x c y = S(p) (7)
and the parameter dependent closed loop,

  ẋ z ∞ z 2   =   A(p) B ∞ B 2 C ∞ D ∞1 D ∞2 C 2 D 21 0     x w ∞ w 2   (8)
Then the corresponding mixed problem is similar to the LTI one: minimize γ 2 and γ ∞ subject to K and Z.

                         A(p) T K + KA(p) KB ∞ C T ∞ B T ∞ K -γ 2 ∞ I D ∞ 1 T C ∞ D ∞ 1 -I   < 0 A(p) T K + KA(p) KB 2 B T 2 K -I < 0 , K C T 2 C 2 Z > 0 , T race(Z) < γ 2 , D 22 = 0 (9)

Design characteristics and performances on the half vehicle model

In the case of a half vehicle, the measure is the suspension deflection and the selected varying variables are the stiffness of the suspension spring, i.e. k f l and k f r . The associated polytope is then formed by k = 4 corners (10) and k-controllers.

Θ =      k f l k f r k f l k f r k f l k f r k f l k f r      , k f l ∈ [k f l , k f l ] k f r ∈ [k f r , k f r ] (10)
According to the dissipative theory, each constraint can be expressed as a supply function, then translated into an LMI (5, 9) [START_REF] Scherer | Multiobjective output-feedback control via LMI optimization[END_REF][START_REF] Scherer | LMI in control[END_REF].

• H ∞ performance is used to enforce robustness to model uncertainties and to express frequency-domain performance specifications

• The H 2 performance can be used to minimize energy of the signal (note the equivalence of these norms in the frequency domain, but not in the time one)

Hence, coupled together, these specifications should improve the single H ∞ constraint. On a half vehicle model the performances we want to reach are multiple. As exposed in [START_REF] Poussot-Vassal | Optimal skyhook control for semi-active suspensions[END_REF][START_REF] Sammier | Skyhook and H ∞ control of active vehicle suspensome practical aspects[END_REF], some frequency specifications have to be specified concerning the suspension deflection, suspended mass and the unsprung masses (to reduce gain around sensitive low frequencies). Also, weight on the control signal prevent actuator saturation. To these frequency specifications, expressed by the H ∞ theory, the addition of H 2 constrain is used to minimize energy of time signals.

Then, the resulting generalized plant is (Figure 2), In the mixed synthesis the considered H ∞ and H 2 controlled output, noted z ∞ and z 2 respectively, are the following,

z ∞ = z s z θ z u f l z u f l z 2 = z s z θ z u f l z u f l
Note that when we will compare the mixed synthesis to the H ∞ one, the controlled outputs of the

H ∞ controller are, z ∞ = z s z θ z u f r z u f l . W zr {f l,f r} W u {f l,f r} - - - + ? W n {f l,f r} S(ρ) z s y u z r z u w 1 w 2 n u żdef - -z θ W θ θ - -z s W zs

Pareto limit

It is impossible to minimize both γ ∞ and γ 2 . In the literature, the mixed problem is generally solved by minimizing a convex combination of H ∞ and H 2 that represents a compromise between the two performances. Such a minimization can take the following form,

min{α 1 T ∞ + α 2 T 2 }, where {α 1 , α 2 } ∈ [0, 1] × [0, 1], α 1 + α 2 = 1
Hence a natural problem raises, how to choose in a smart way α 1 and α 2 . The concept of non-inferiority (also called Pareto optimality) is used here to characterize the objectives. A non-inferior solution is one in which an improvement in one objective requires a degradation of an other. In our case the objectives are H ∞ and H 2 . To plot the Pareto optimum, applied to our problem, we iteratively fix the γ ∞ and minimize the γ 2 . The corresponding results are given in Figure 3.

k {f l,f r} ∈ k nom × [1, 1.2], [1, 1.5], [1, 2]
The achievable combinations {γ ∞ ,γ 2 } is the set of couples located over the Pareto limit. The Pareto limit is also useful to measure the conservatism of a method and to exhibit how much one can decrease the performances with a qLPV approach compared to the LTI one. Such a Figure can also motivate researches on polytope reduction. In effect, the more you increase the size of your polytope (bounds of the parameters), the more far you go from the LTI Pareto optimum (Figure 3), and loose performance.

SIMULATION RESULTS

To validate the control design, first, simulations are done in order to show the advantages of mixed synthesis compared to single H ∞ objective, then, we study the influence of the choice of the couple {γ ∞ , γ 2 } on the reached performances. Finally, we compare the LTI mixed approach with the qLPV one. Note that on theses simulations, when a control law is considered, the damper is removed so that the considered suspension simulated is a real semi-active one i.e. a spring (nonlinear) plus an active actuator. In such a way we explicitly model the fact that the damper is replaced by the actuator. Such a control also justify the choice of {k f l , k f r } as varying parameters in the qLPV synthesis.

The LTI case

we show the advantages of the mixed H ∞ / H 2 compared to H ∞ synthesis (for the same γ ∞ ). In this simulation we generate a step disturbance on the first then on the second wheel, then a roll moment disturbance and we compare controllers performances according to the passive suspension.

By using the mixed synthesis instead of single H ∞ , we reduce the roll angle due to the roll energy minimization (Figure 4). Then, we compare the performances of the mixed synthesis for different couples {γ ∞ , γ 2 } (Figure 5). If one decrease the γ 2 attenuation value, then it increases the γ one (see Pareto limit Figure 3). The H 2 criteria's aim is to minimize the energy and variations of a signal (here, z s but also the control input are limited). Hence we observe that the z s variations are smoother by using a smaller attenuation gain on the H 2 criteria (less oscillations, i.e. ameliorate vertical comfort).

The qLPV case

LTI and the qLPV controllers are here (Figure 6 and7) investigated for parameters k {f l,f r} varying between k nom × [1, 1.95] (i.e. [k f l,f r , k f l,f r ]) and for a fixed γ ∞ = 0.25. Here, we assume bigger road step disturbance to reach the nonlinear area of the suspension deflection.

The qLPV synthesis improves the performances achieved by the LTI one. Then, such a control tackle the nonlinearities, hence enforce robustness. The α shows the scheduling done according to the parameters variations. Note that a qLPV approach, even if enforces robustness, LTI control signal qLPV control signal 7: Control signal u f l (up) and u f r (down) for the LTI (dashed) and qLPV (solid) controller exhibits more complexity than the LTI one because it increases the number of controllers to be synthesized (4 in our case) and requires to schedule them in real-time. Then, the control signal looks sensitive to the parameters variations (Figure 7). Nevertheless, we use in both synthesis (LTI and qLPV) the same number of measures.

CONCLUSION AND FUTURE WORKS

In this paper we investigate a multi-objective mixed qLPV H ∞ / H 2 control applied to a half vehicle model. A special interest is made on the advantages of such a synthesis and on the compromises that have to be done in multi-objective applications. By using the Pareto limit (non-inferior solution) we expose a smart way to select the objectives and show the influence on significative driving situations.

Figure 1 :

 1 Figure 1: Half vehicle model

Figure 2 :

 2 Figure 2: Generalized plant

Figure 3 :

 3 Figure 3: LTI (solid) vs qLPV (dashed) Pareto limit fork {f l,f r} ∈ k nom × [1, 1.2], [1, 1.5],[START_REF] Abdellahi | On the design of active suspension system by H ∞ and mixed H 2 /H ∞ : An LMI approach[END_REF][START_REF] Chilali | H ∞ design with pole placement constraints: an LMI approach[END_REF]

2 Figure 4 :

 24 Figure 4: Comparison between H ∞ and Mixed (solid) design with Passive (solid slim)

Figure 5 :

 5 Figure 5: Comparison of mixed synthesis performances according to different {γ ∞ , γ 2 }

Figure 6 :

 6 Figure 6: Comparison of LTI (dashed) and qLPV (solid) mixed synthesis performances

Acknowledgement: This work was supported by the Hungarian National Science Foundation (OTKA) under the grant T -048482 which are gratefully acknowledged.