Filippo Pantellini

Simone Landi

Species segregation in one-dimensional granular system simulations

Keywords: 45, 70, Mg Granular flow: mixing, segregation and stratification -45, 50, Tn Collisions -02, 70, Ns Molecular dynamics and particle methods

We present one dimensional molecular dynamics simulations of a two species, initially uniform, freely evolving granular system. Colliding particles swap their relative position with a 50% probability allowing for the initial spatial ordering of the particles to evolve in time and frictional forces to operate. Unlike one dimensional systems of identical particles, two species one dimensional systems of quasi-elastic particles are ergodic and the particles' velocity distributions tend to evolve towards Maxwell-Boltzmann distributions. Under such conditions, standard fluid equations with merely an additional sink term in the energy equation, reflecting the non elasticity of the interparticle collisions, provide an excellent means to investigate the system's evolution. According to the predictions of fluid theory we find that the clustering instability is dominated by a non propagating mode at a wavelength of the order 10πL/N ε, where N is the total number of particles, L the spatial extent of the system and ε the inelasticity coefficient. The typical fluid velocities at the time of inelastic collapse are seen to be supersonic, unless N ε 10π. Species segregation, driven by the frictional force occurs as a result of the strong temperature gradients within clusters which pushes the light particles towards the clusters' edges and the heavy particles towards the center. Segregation within clusters is complete at the time of inelastic collapse.

Introduction

Granular materials are ubiquitous in the macroscopic world [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF]. The consequence of the inelasticity of the collisions between grains is that the standard fluid equations used to describe a gas of atoms or molecules need to be modified. Even for the simplest systems, the mathematical expressions for the transport coefficients turn out to be much more involved than the expressions for a system of elastic particles. In addition, granular materials are often polydisperse. For example sand grains are generally characterized by a broad range of particle sizes and shapes. Thus, theoretical models and numerical simulations based on granular systems of identical particles may be too limited to describe a real system, even on a purely qualitative level. On the other hand, systems of identical particles have the non negligible advantage of being mathematically simple while retaining many important features of real dissipative systems.

One of the most spectacular consequences of the inelastic nature of collisions in a granular system is the so called inelastic collapse. In the absence of energy injection statistical fluctuations in an initially uniform system can lead some regions to cool faster than others producing linearly growing density fluctuations which can be described in the frame of continuum model theories [START_REF] Mcnamara | Kinetics of a onedimensional granular medium in the quasielastic limit[END_REF][START_REF] Goldhirsch | Clustering instability in dissipative gases[END_REF][START_REF] Efrati | Hydrodynamic Singularities and Clustering in a Freely Cooling Inelastic Gas[END_REF]. In the non linear regime the high density regions collapse producing clumps of particles where the collision frequency grows to infinity [START_REF] Bernu | One-dimensional bounce of inelastically colliding marbles on a wall[END_REF][START_REF] Mcnamara | Inelastic collapse and clumping in a one-dimensional granular medium[END_REF]. The instability is a long wavelength instability which is hardly avoidable in freely evolving systems provided the number of particles in the system exceeds some minimum value which very much depends on the system's characteristics [START_REF] Mcnamara | Inelastic collapse and clumping in a one-dimensional granular medium[END_REF][START_REF] Goldhirsch | Clustering instability in dissipative gases[END_REF] It is often acceptable to describe granular systems by means of fluid equations derived from the inelastic Boltzmann equation using the Chapman-Enskog procedure [START_REF] Jenkins | Kinetic theory for binary mixtures of smooth, nearly elastic spheres[END_REF][START_REF] Brey | Hydrodynamics for granular flow at low density[END_REF][START_REF] Garzó | Hydrodynamics for a granular binary mixture at low density[END_REF]. Hydrodynamic models have been applied to the case of inelastic particles in a constant gravitational field [START_REF] Eggers | Sand as Maxwell's Demon[END_REF][START_REF] Brey | Hydrodynamic Maxwell demon in granular systems[END_REF][START_REF] Ramírez | Temperature inversion in granular fluids under gravity[END_REF]. These systems often behave non intuitively. Hence, it is shown in [START_REF] Ramírez | Temperature inversion in granular fluids under gravity[END_REF] that a collection of inelastic particles confined to a vertical box with permanent injection of energy from the base (e.g. by shaking the base of the box vertically) the granular temperature profile is a non monotonic function of height. The temperature decreases near the bottom, goes through a minimum, and finally rises indefinitely with height above the minimum.

It has been recently shown in [START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF] that during the homogeneous regime of the temporal evolution of a one dimensional system of identical particles the spatial ordering of the particles appears to play a crucial role in shaping the particles' velocity distribution. Thus, works [START_REF] Ben-Naim | Multiscaling in inelastic collisions[END_REF], based on the so-called pseudo-Maxwell model, a homogeneous model where no spatial ordering is assumed by construction, have been found to be unable to describe the evolution of a system of inelastically interacting particles, even during the homogeneous regime. The late evolution of systems of identically spatially ordered particles has been shown to be strongly affected by the development of spatial correlations between the particles' velocities eventually leading to an inelastic collapse [e.g. [START_REF] Mcnamara | Kinetics of a onedimensional granular medium in the quasielastic limit[END_REF][START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF].

One dimensional systems of identical particles behave in a very peculiar way. Thus, in the elastic limit, the overall velocity distribution of the particles' is not modified by collisions. Colliding particles do merely exchange their velocities and the system is not ergodic. Adjunction of just one anomalous particle forces the velocity distribution to relax towards a Maxwell-Boltzmann distribution.

In this paper we use molecular dynamics (MD) simulations to investigate the evolution of a two species onedimensional, periodic, system of N point particles. N/2 particles have mass m 1 = 1 (species 1) and N/2 particles have mass m 2 = 4 (species 2). The system is a prototype for more complex systems with broader mass distributions. One important ingredient of the model is that colliding particles swap their relative position with a 50% probability allowing for the initial spatial ordering of the particles to evolve in time. This introduces the possibility of species to stream with respect to eachother giving frictional forces the opportunity to operate.

The restitution coefficient r, appearing in the collision rules [START_REF] Ramírez | Temperature inversion in granular fluids under gravity[END_REF] and [START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF], is close to elastic such that the thermalization time is always much shorter than the cooling time. As a consequence, the velocity distribution function of both species are always close to Maxwellian, at least as long as the species do not segregate. Such dissipative systems are always linearly unstable with respect to the so-called clustering instability. In one dimension the clustering instability is a fluid instability which is found to triggers the formation of non propagating spatial density inhomogeneities at scales of the order 10πL/N ε, where ε ≡ 1 -r 2 and L is the size of the system. The clusters of particles which form due to the clustering instability are characterized by a temperature profile which increases exponentially as a function of the particle index l away from the center of the cluster (cf e.g. Figure 6). These strong temperature gradients within clusters drive the light particles away from their centers towards the edges of the clusters while the opposite happens to the heavy particles which tend to become concentrated in a very small region near the center of the cluster. (cf Figures 4 and5). The species segregation is driven by the frictional force [START_REF] Jenkins | Kinetic theory for binary mixtures of smooth, nearly elastic spheres[END_REF] which points in opposite directions for the light and heavy particles. At the time of inelastic collapse, fluid motions may or may not be supersonic depending on whether or not the number of particles N exceeds a critical number of order 10π/ε. In order to avoid confusion, we do loosely define the time of inelastic collapse in our simulations as the time for which the molecular dynamic simulation is essentially "frozen". By "frozen" we mean that the total energy and, consequently the fluid velocity profiles of the system do not change significantly by doubling the number of collisions.

The paper is organized as follows. In Section 2 we introduce the set of fluid equations whereon we base our analysis of the simulation results. The numerical model is presented in Section 3. Simulations with N = 19600 particles and 3 different values of the restitution coefficient r are presented in Section 4, In Appendix A we develop the linear theory of the clustering instability for a twospecies, one dimensional system. The collision frequency for such a system in the homogeneous limit and close to the thermodynamic equilibrium, is given in Appendix B.

Fluid equations for a mixture of two species

In the general case [see e.g. 15], the expressions for the hydrodynamic transport coefficients describing a granular gas with two species of inelastically interacting particles are so involved that they are of little practical use to describe real systems. It is not even clear under which conditions fluid equations are relevant for real granular systems. However, over the years it has been shown that fluid equations are extremely useful to help understand results from granular system simulations. Over the few last decades, explicit expressions for the transport coefficients in 2 and 3 spatial dimensions have been published by various authors. To our knowledge, expressions for the 1d case have not been published until now, certainly because of the non-ergodicity of most 1d systems. Yet, in the case of a mixture of unequal particles, collisions among particles of different mass do inevitably drive the velocity distribution functions for all species towards Maxwell-Boltzmann. The tendency to thermalization in such one dimensional mixtures indicates that "standard" fluid equations may apply to such systems as well.

The fluid equations of this section all stem from a Chapman-Enskog type treatment of the Boltzmann equation by various authors [START_REF] Jenkins | Kinetic theory for binary mixtures of smooth, nearly elastic spheres[END_REF][START_REF] Brey | Hydrodynamics for granular flow at low density[END_REF][START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF][START_REF] Garzó | Hydrodynamics for a granular binary mixture at low density[END_REF] which implies that their validity is at best limited by the assumptions on which the Chapman-Enskog procedure is based on. One of the most stringent conditions, in a system characterized by a mean free path l mfp and a spatially varying temperature T (x), is the requirement of the Knudsen number K ≡ l mfp ∂ log T /∂x being much smaller than unity. The K ≪ 1 condition means that the particles' distribution functions do not vary substantially over distances of the order of the mean free path l mfp and that the velocity distributions are close to Maxwellian. A related condition is that fluid time scales must be long compared to the collisional time scale, i.e. the time between successive collisions of a thermal particle.

We write the continuity equation for two species α = 1, 2 in terms of the species mean velocities u α and mass densities ̺ α = m α n α :

∂̺ α ∂t + ∂ ∂x (̺ α u α ) = 0. ( 1
)
Defining the total density ̺ ≡ ̺ 1 + ̺ 2 and the center of mass flow velocity u ≡ (u 1 ̺ 1 +u 2 ̺ 2 )/̺ one can add the two continuity equations ( 1) and write a one fluid continuity equation:

∂̺ ∂t + ∂ ∂x (̺u) = 0. (2)
According to [START_REF] Jenkins | Segregation in Binary Mixtures under Gravity[END_REF] we write the momentum equation for both species in the form

̺ α ∂u α ∂t + u α u α ∂x = - ∂[p α (1 + K α )] ∂x + φ α (3)
where φ α represents the frictional force on the species α due to the collisions with the other species in the system and where K α is the contribution to the pressure tensor from collisions with all species (shear and bulk viscosity). Momentum conservation implies φ 1 = -φ 2 . The sign difference between φ 1 and φ 2 indicates that if the frictional force is dominant, it can act as a species filter. The partial pressures p α in equation ( 3) are defined with respect to the center of mass flow velocity u, i.e.

p α ≡ m α (v -u) 2 f α (v)dv = n α T α . (4)
For the remaining of this section, and throughout the paper, we assume that collisions are efficient enough to ensure energy equipartition among species, i.e. T = T 1 = T 2 . This restriction implies the characteristic cooling time for each species being long with respect to the thermalization time. The energy equation for the mixture in the onedimensional case is then given by [START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF] ∂p ∂t

+ u ∂p ∂x + 3p ∂u ∂x + 2 ∂q ∂x - 8 3 η ∂u ∂x 2 = -ζp (5)
where p ≡ p 1 +p 2 is the total pressure, q the total heat flux, η the shear viscosity and ζ the cooling rate due to inelastic collisions. The bulk viscosity is zero by construction [see 9]. Similarly, the momentum equation for the mixture can be written in the conventional form [START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF] ∂u ∂t

+ u ∂u ∂x = - 1 ̺ ∂ ∂x p - 4 3 η ∂u ∂x . (6)
We note that this equation is equivalent to equations (3) up to terms which are quadratic in the species relative drift velocity δu 1 ≡ (u 1 -u). In the quasi-elastic limit, for small relative drift velocities, and assuming energy equipartition among species the following relations do conveniently close the above set of fluid equations. According to [START_REF] Jenkins | Segregation in Binary Mixtures under Gravity[END_REF] the frictional force may be cast into the form

φ 1 = a(n 1 , n 2 ) 2T n 1 n 2 n ∂ ∂x ln n 1 n 2 T m12 (7)
with

m 12 ≡ (m 2 -m 1 )/m, m ≡ m 1 + m 2 and n ≡ n 1 + n 2 .
In equation [START_REF] Jenkins | Kinetic theory for binary mixtures of smooth, nearly elastic spheres[END_REF] a is a dimensionless function of n 1 and n 2 which takes into account the species volume fraction and the radial distribution function of the contacting pairs. Explicit forms have been given for the 2 and 3 dimensional case [START_REF] Jenkins | Segregation in Binary Mixtures under Gravity[END_REF]. In this paper we do not care about the exact form of equation 7. We are just interested in the qualitative dependence of φ 1 on the temperature T and on the relative species concentration n 1 /n 2 . Following [START_REF] Garzó | Hydrodynamics for a granular binary mixture at low density[END_REF][START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF] we write the constitutive relation for the heat flux as

q = - nT ν 0 √ m 1 m 2 ∂T ∂x b 1 n 1 n m m 1 + b 2 n 2 n m m 2 (8)
where we have defined the frequency ν 0 ≡ V 12 nσ and the characteristic velocity V 2 12 = 2T m/(m 1 m 2 ), which, for species with Maxwellian distributions, is just the thermal velocity for the relative velocities between the two species (cf Appendix B). Again, both, b 1 and b 2 , are rather involved functions of n 1 /n, m 1 , m 2 and even V 12 , which we do not need to care about in this paper. In the above equation σ has the dimension of a surface and may be interpreted as a collisional cross section. Formally, we express σ in terms of the total number of particles N , the box length L and the mean number density n:

σ = N nL . (9)
With this definition, in a uniform system one has ν 0 = V 12 N/L. Starting from equation ( 27) of Jenkins and Mancini [START_REF] Jenkins | Kinetic theory for binary mixtures of smooth, nearly elastic spheres[END_REF] we find that ζ must depend on the species densities via an expression of the kind

ζ = εν 0 c 11 n 2 1 n 2 m 2 m + c 22 n 2 2 n 2 m 1 m + c 12 n 1 n 2 n 2 (10)
where the numerical coefficients c 11 , c 12 , and c 22 are unknown for the 1d case and where ε is the fractional energy loss per collision (see Section 3). Equation ( 10) is correct up to terms of order two in the gradients of hydrodynamic quantities [START_REF] Brey | Hydrodynamics for granular flow at low density[END_REF][START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF]. Following Garzo et al [START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF] we write the shear viscosity coefficient η as

η = nT ν 0 d 1 n 1 n m m 2 + d 2 n 2 n m m 1 (11)
where, once more, d 1 , and d 2 are order unity coefficients for which we do not care about the exact form. As a final remark, we note that in the set of fluid equations ( 2),( 5) and ( 6) the inelasticity coefficient ε does only appear in the energy equation ( 5) through the cooling rate ζ ∝ εν 0 .

In particular, neither the friction φ 1 , nor the heat flux q depend on ε, which is only acceptable in the quasi-elastic limit ε ≪ 1 [START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF].

Simulation model

We use a one-dimensional model which is close to the traditional point-like hard particles on a ring model [START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF] with the notable difference that N/2 particles have mass m 1 and N/2 particles have mass m 2 . In the strict one dimensional case the post collision velocities v ′ i and v ′ j for two colliding particles i and j = i + 1 of mass m i and m j and center of mass velocity

v 0ij = (m i v i + m j v j )/(m i + m j ) are given by v ′ i = ± m j v ij m i + m j r + v 0ij (12)
v ′ j = ∓ m i v ij m i + m j r + v 0ij ( 13
)
where r is the restitution coefficient (r = 1 for elastic collisions), v k the velocity of particle k = i, j before the collision and v ij ≡ v j -v i < 0 (we implicitly suppose a one dimensional axis x with coordinates increasing from left to right x i < x j ). With this definition of the restitution coefficient r the kinetic energy of two colliding particles is reduced by a factor ε ≡ 1 -r 2 in the center of mass frame. The + sign in [START_REF] Ramírez | Temperature inversion in granular fluids under gravity[END_REF] and the -sign in ( 13) correspond to the normal case where particle i and j do not exchange their relative positions during collision. Thesign in [START_REF] Ramírez | Temperature inversion in granular fluids under gravity[END_REF] and the + sign in (13) correspond to the case where particles i and j swap their relative positions during collision. In the latter case, taking r = 1, one has v ′ i = v i and v ′ j = v j , i.e. particles just ignore each other. In principle one should reject this possibility as it is a non physical event for the case of a system of beads on a ring. We shall take advantage of this possibility as it allows for the different particles to become distributed spatially independently of the initial ordering with the possibility of species segregation to operate as in 2 and 3 dimensional systems.

Simulation results

In all simulations N particles are initially distributed uniformly over the spatial domain [0, L[. We consider a simple two species system of N/2 particles with mass m 1 = 1 and N/2 particles with mass m 2 = 4. Particles of both species are initially disposed in alternating order. The initial velocity of a particle with mass m α is randomly selected in the interval [-0.5, 0.5]/m 1/2 α following a uniform probability distribution as in the [START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF] where the reference case of N identical particles is discussed. The mean energy per particle is therefore the same for all species. Bulk velocities of all species are initially set to zero. The parameters for the 3 runs which will be discussed in the paper are listed in Table 1. The somewhat arbitrary mass ratio m 2 /m 1 = 4, which will be used throughout the paper, is sufficiently small to ensure efficient interspecies energy thermalization, so that we can assume T 1 = T 2 , and sufficiently large for the frictional force φ 1 (see equation ( 7)) to play a non negligible role.

Homogeneous phase

The evolution of a one-dimensional system with not all particles having the same mass is very different from that of a one dimensional system of identical particles [e.g 6, 2, 13]. The difference is substantial, even in the elastic limit r = 1. Indeed, while in the case of elastically interacting identical particles the velocity distribution does not change in time, the asymptotic velocity distribution is a Maxwell-Boltzmann velocity distribution when different species coexist in the system. The global thermalization rate strongly depends on the particles' mass distribution, each species α thermalizing at different rates depending on the species number density N α and on the mass m α of the species particles. For example in the case of a system with N ≫ 1 and just one anomalous particle, the evolution towards a Maxwell-Boltzmann distribution is much slower than in the case of a two species system with N/2 particles for each species, since only collisions among unequal particles can modify the velocity distribution functions in one dimensional systems.

In Figure 1 a snapshot of the velocity distribution function of both species are shown for Run 1 after a total number of 8 × 10 7 collisions . The distribution for both, light and heavy particles, which should be compared to the double peaked distribution of Figure 1 in [START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF] for the one species case, do closely follow a Maxwell-Boltzmann distribution disregarding of the fact that the initial distribution was strongly non thermal. This tendency of the velocity distributions to evolve towards Maxwell-Boltzmann 00000000 00000000 00000000 11111111 11111111 11111111 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 allows the use of standard fluid theories [START_REF] Jenkins | Kinetic theory for binary mixtures of smooth, nearly elastic spheres[END_REF][START_REF] Brey | Hydrodynamics for granular flow at low density[END_REF][START_REF] Garzó | Mass and heat fluxes for a binary granular mixture at low density[END_REF][START_REF] Garzó | Hydrodynamics for a granular binary mixture at low density[END_REF] to analyze the system's behavior.

Figure 2 shows the evolution of the energy in Run 1. The energy decreases exponentially as a function of the number of collisions I during the first 5 10 8 collisions. We shall define this phase as the homogeneous phase as it is perfectly well described by the homogeneous version of the fluid energy equation ( 5

Inhomogeneous phase

After the initial homogeneous phase the system enter a new, inhomogeneous, regime where cluster formation becomes the dominant effect. The transition between the homogeneous and the inhomogeneous regime is clearly visible as a departure from the exponential energy decay in Figure 2. In this section we discuss the cluster structures observed in the simulations in the light of the fluid equations exposed in Section 2.

Species segregation

As in the one species case, the inhomogeneous regime in the multi species case is characterized by the presence of clusters of particles. Given the smallness of the inelasticity coefficient ε and the initial homogeneity of the system, one expects the pressure to be spatially uniform for both species, except near the center of the collapsing cluster where the cooling rate ζ goes to infinity. Indeed, the temperature and density profiles, plotted for the light particles at the end of Run 1 are seen to vary in antiphase (cf Figure 3) ensuring an approximate pressure balance. However, light and heavy particles are not evenly distributed within the cluster. The relative concentration of light particles as a function of the particle index l is shown in Figure 4 which clearly shows that the heavy particles are concentrated in the central part of the cluster, leaving an excess of light particles in the wings of the 6000 particles cluster.

In order to appreciate the spatial scales associated with the concentrations of light and heavy particles in the cluster we compute the integrated heavy particles excess func- tion Y defined as

Y (x) = N l=1 (δ m l m2 -δ m l m1 ) x 0 δ(ξ -x l )dξ (15)
Contribution to the integral is -1 for light particles (m l = 1) and +1 for heavy particle (m l = 4). The function Y (x) is plotted in Figure 5 and illustrates how disparate the spatial repartition of light and heavy particles are in the system. The core of the cluster contains some 3×10 3 heavy particles concentrated in a region of size 10 -8 . From Figure 4 we know that the 3 × 10 3 light particles of the cluster are largely excluded from this small central region dominated by the heavy particles but Figure 5 shows that the light particles dominate over a much larger region of size 10 -5 , the characteristic size of the cluster at the end of the simulation. The species segregation is the visible effect of the frictional force φ 1 (cf equation ( 7)) of the two fluid momentum equations [START_REF] Mcnamara | Inelastic collapse and clumping in a one-dimensional granular medium[END_REF]. In case of equal densities n 1 = n 2 and a spatially varying temperature profile T (x) the frictional force φ 1 ∝ (m 2 -m 1 )∂T /∂x points up the temperature gradient for the light species m 1 and down the temperature gradient for the heavy species m 2 . Segregation appears as a ineluctable consequence of the temperature gradients within clusters.

Temperature profiles in clusters

Figure 3 shows that the temperature in a cluster strongly decreases towards its center. Figure 6 shows the temperature profiles in the inhomogeneous regime for Run 1 and 2. The figure illustrates the fact that the temperature drops exponentially as a function of the particles index l towards the center of the cluster. In these particular cases the temperature decreases by several order of magnitudes from edge to center. As expected, the thermalization time being Fig. 6. Fluid temperature (circles) and light particles' temperature (plus) as function of the particle index for Run 1 (top panel) and Run 2 (bottom panel) during the inhomogeneous regime. In both simulations the temperature is found to decrease exponentially towards the center of the cluster according to a characteristic scale of order ε -1/2 . much shorter than the cooling time, the fluid and light particle temperature profiles closely follow each other. The interesting point is that the granular temperature decreases exponentially towards the center of the cluster with a scale which seems to be of order ε -1/2 . Again, the ε -1/2 scaling of the temperature profile in Fig. 6 can be deduced easily from the fluid equations. In order to do so, we evaluate the temperature profiles within clusters in the static limit. Using the energy equation ( 5) we can write

2 ∂q ∂x = -ζp (16)
Let us suppose that in a given region of the system species 1 is the dominating one, i.e. n ≈ n 1 ≫ n 2 . In this case the heat flux reads (cf equation ( 8))

q = - b 1 σ T 2m 1 1/2 ∂T ∂x ( 17
)
whereas the right hand side of Eq. ( 16) becomes

-ζp = -εc 11 σ 2T m 1 1/2 T n 2 (18)
and the static energy equation ( 16) writes

∂ ∂x T 1/2 ∂T ∂x = ε c 11 b 1 σ 2 T -1/2 p 2 0 ( 19
)
where we have used the static limit p = p 0 = const of the momentum equation [START_REF] Mcnamara | Inelastic collapse and clumping in a one-dimensional granular medium[END_REF]. In order to compute the temperature profile as a function of the particle index l we have to introduce the transformation from the spatial coordinate x to the particle index l.

l(x) = A x 0 n(ξ)dξ ( 20
)
where A is a surface such that An(ξ)dξ is the number of particles in the interval [ξ, ξ + dξ]. Using the transformation (20) we write the energy equation ( 19) in terms of the particle index l

∂ ∂l T -1/2 ∂T ∂l = ε 4 c 11 b 1 T 1/2 (21)
where we have used the fact that the particles' cross section σ is equal to A/2 given the 50% swap probability for colliding particles. The solution of ( 21) is

T ∝ exp ±l ε c 11 2b 1 1/2 (22)
If c 11 and b 1 are numerical constants of order unity, it follows, as already suggested by Figure 6, that the typical scale of variation for the temperature profile is of order ε -1/2 .

Cluster formation and linear fluid theory

As in one species simulations [START_REF] Mcnamara | Inelastic collapse and clumping in a one-dimensional granular medium[END_REF], Figure 7 shows that at the time of inelastic collapse, the number of clusters is an increasing function of the N ε. Even the Mach number of the fluid motions is seen to depend on N ε, the simulation with the smallest value of N ε being characterized by subsonic motions and the two other simulations being characterized by supersonic motions with the evident formation of shocks. As pointed out by Ben-Naim et al [START_REF] Ben-Naim | Shocklike Dynamics of Inelastic Gases[END_REF] the piecewise linear velocity profiles observed in Figure 7 are a direct consequence of the momentum equation [START_REF] Mcnamara | Inelastic collapse and clumping in a one-dimensional granular medium[END_REF]. Indeed, as already discussed in the context of Figure 3, the pressure in the system is essentially constant everywhere except within clusters (shocks in the supersonic regime). In this case the momentum equation reduces to ̺(∂u/∂t + u∂u/∂x) = 0.75η∂ 2 u/∂x 2 with solutions of the type u(x, t) = u 0 + (x -x 0 )/t, where u 0 and x 0 are integration constants. The difference between the velocity profiles at the time of inelastic collapse for the three runs shown in Figure 7 can be explained on the basis of a linear analysis of the fluid equations of Section 2 developed in Appendix A. The linearized system of equations has three eigenmodes, two of which correspond to the right and left propagating sound wave and the non propagating (zero real frequency) entropy mode. The main difference with respect to classical gas dynamics is that all modes are unstable at long wave length. The other important prediction of linear theory is that in general, but more specifically for wave vectors k such that k ≡ kL/(N ε) 1, the non propagating entropy mode grows substantially faster than the sound mode. Indeed, the growth rate of the entropy mode always peaks at k = 0 whereas the sound mode's growth rate vanishes as k → 0. A careful analysis of the growth rates does actually shows that the entropy mode is always the fastest growing mode (see Appendix A).

The dominant mode

What determines the number of clusters in a one dimensional system at the time of inelastic collapse? Knowledge of the initial conditions for the velocity fluctuations of all modes makes it possible to estimate different characteristic time scales in terms of the dimensionless (pseudo) time variable τ (defined in A.5) using the approximations of linear theory.

Initially, in all systems, particles are distributed uniformly in the spatial interval [0, 1[ with N/2 particles of mass m 1 = 1 and N/2 particles of mass m 2 = 4. The velocity of a particle of mass m is sorted randomly in the velocity interval [-0.5, 0.5]/m 1/2 so that the average kinetic energy per particle is the same for all species corresponding (formally) to an initial temperature T 0 = 1/12 and thermal velocity v 2 0 = 2/(6m) = 1/15. Thermalization of such a system takes only a few collisions per particle so that all species can be assumed to be initially distributed according to a Maxwellian distribution at temperature T 0 = 1/12, ε being much smaller than unity. The initial amplitudes of the velocity fluctuations can be evaluated using the Wiener-Khintchine theorem which states that the power spectrum of the fluid velocity field for a set particles with uncorrelated velocities is flat (white noise). The immediate consequence of the Wiener-Khintchine theorem is that energy is evenly distributed among all Fourier modes, i.e. all modes have equal amplitude |δu 0 |. Requiring the total kinetic energy of all particles in the system being N T /2 one must have

δu 2 0 = T m 4 N = 2 v 2 0 N (23)
In all simulations, at t = 0, we have T = 1/12, N = 19600, m 1 = 1 and m 2 = 4. The initial amplitude of the velocity fluctuations are therefore of order δu 0 = 2.6 × 10 -3 with v 0 = 0.258, which corresponds to

δU 0 = (2/N ) 1/2 = 1.01× 10 -2 .
From linear theory we know that at small values of the wave vector k 1/2 the growth rate is approximately given by γN/ε = 1 -k2 /2 indicating that the mode with the longest possible wavelength is always the fastest growing mode. It is then expected that clustering always occurs at the largest possible scale allowed by the system with only one cluster growing in the system. This is effectively the case in Run 1 only but not in Run 2 (with 2 clusters) and particularly in Run 3 where several clusters occur. This can be understood by noting that linear theory predicts that the density fluctuation of the entropy mode is proportional to the wave vector k:

δn k n = -i kδU k (24)
where δU k ≡ δu k /v 0 is essentially the fluctuation's Mach number. Thus, for two modes k 1 and k 2 with equal initial velocity fluctuation amplitude δU k1 = δU k2 , the one with the largest wave vector is the one with the largest associated density fluctuation and therefore a more likely candidate for triggering the inelastic collapse, provided the growth rate of the two modes are not too dissimilar. For example, in Run 3, there are some 10 visible density clusters corresponding to an average distance between clusters of the order λ = L/10. Therefore, the mode which triggers the inelastic collapse has k = 0.16 and a linear growth rate which is only roughly 1% less than the growth rate of the smallest possible wave vector kmin ≡ 2π/N ε, but, according to equation 24, with a 10 times larger density fluctuation. Thus, as long as the velocity fluctuation amplitude of the mode k = 0.16 is not much smaller than the velocity fluctuation amplitude of the fastest growing mode, i.e. potentially during the first 100γ -1 ≈ 100N/ε = 9.8 × 10 7 collisions, the strongest density fluctuations in the system are due to the mode with k = 0.16. Indeed, in Run 3 particle clustering and supersonic motions are already present after I ≈ 4 × 10 6 collisions, which, using the definition for τ given by (B.7) in Appendix B corresponds to τ = 2.06I ≈ 8 × 10 6 . This time is too short by more than one order of magnitude for the fastest growing mode to become dominant in terms of density fluctuations. We conclude that, for quasi elastic restitution coefficients ε ≪ 1, the mode which dominates the inelastic collapse has a wave vector k ∼ 1/5 corresponding to a wavelength λ ∼ 10πL/N ε. If such a long wavelength does not fit into the system, the dominant mode is simply the one with the longest possible wavelength, i.e. λ = L.

Supersonic fluid motions

Let us now address the question of the magnitude of the characteristic fluid velocities at the time of the inelastic collapse. Velocities are subsonic in Run 1 and supersonic in Runs 2 and 3.

We introduce the sonic time scale τ s , defined as the time for the amplitude of the velocity fluctuation δU k to reach unity assuming linear growth, i.e.

τ s ( k) = - ln(δU 0 ) γ( k) . ( 25
)
For a system of uniformly distributed particles with uncorrelated velocities, the initial fluctuation amplitude is δU 0 ≈ (2/N ) 1/2 . Similarly, using equation ( 24) we define τ n as the time for the density fluctuation δn(k)/n to reach an amplitude of order unity, i.e.

τ n = - ln( k δU 0 ) γ( k) . (26)
Both, τ s and τ n are listed in Table 2 for the mode with the longest wavelength kmin .

In Run 1 the motions remain subsonic at the time of the inelastic collapse despite the simulation lasting significantly longer than τ s , meaning that supersonic motions had enough time to develop. We conclude that an inelastic collapse is not necessarily associated with supersonic motions. The explanation at hand can be found in Table 2 which shows that in Run 1 the characteristic times τ s and τ n are very much the same for the dominant mode k = k min ≡ 2π/N ε while the sonic time scale τ s is significantly shorter than τ n for the dominant mode in Runs 2 and 3. From Table 2 it follows that supersonic motions develop before the occurrence of an inelastic collapse provided kmin 1/5, i.e. 10π N ε.

Conclusions

We have performed simulations of a one dimensional two species periodic system of N point particles undergoing inelastic collisions. In order to keep the possibility for the system to arrange particles independently of the initial configuration, particles are allowed to change their relative positions during collisions with a 50% probability. Initial velocities for particles of species α are selected randomly according to a constant distribution in the range [-0.5, 0.5]/m 1/2 α . Unlike one dimensional systems of identical particles, where the homogeneous regime is characterized by double peaked velocity distributions [START_REF] Baldassarri | Influence of correlations on the velocity statistics of scalar granular gases[END_REF][START_REF] Benedetto | A Kinetic equation for granular media[END_REF], one dimensional multi-species systems tend to evolve towards a Maxwell-Boltzmann distribution. This tendency justifies the use of "standard" fluid equations to describe such systems. As for one dimensional systems of identical particles, the inhomogeneous regime is characterized by the formation of clusters which eventually make an inelastic collapse. It has been shown in [START_REF] Mcnamara | Kinetics of a onedimensional granular medium in the quasielastic limit[END_REF] that clustering in one species systems is driven by a beam instability. We show that in a two species system clustering is triggered by a fluid instability of the non propagating entropy mode, for which density and temperature fluctuations vary in antiphase ensuring spatial pressure balance. The sound mode is also found to be unstable but its growth rate is always substantially smaller than the growth rate of the entropy mode. Because of the non propagating nature of the entropy mode, the one dimensional system does evolve naturally towards a series of clusters with symmetric temperature and density profiles and antisymmetric velocity profiles.

The number of particles in the clusters is determined by the number of particles in one wavelength of the dominant mode of the clustering instability, i.e. the mode which first reaches a relative density fluctuation of order unity. The typical wavelength of this mode has been found to be of the order λ ≈ 10πL/N ε corresponding to 10π/ε particles per cluster. Clustering is generally characterized by supersonic fluid motions unless the total number of particles N 10π/ε.

A peculiar aspect of the late evolution of a multispecies systems is the appearance of species segregation within the collapsing clusters, heavy particles becoming concentrated in an extremely small region around the center of the cluster and the light particles filling the space between the central region and the edges of the cluster. Segregation is driven by the temperature gradient inside the cluster [START_REF] Hsiau | Granular thermal diffusion in flows of binary-sized mixtures[END_REF] via the frictional force φ in the one-species fluid momentum equation (3).

We thank both referees for the many pertinent and helpful comments which strongly helped us to improve the original manuscript.

A Linear theory for a one-dimensional system

In this appendix we develop the linear theory for the clustering instability in a one-dimensional and spatially uniform system of inelastically colliding particles using the set of fluid equations of Section 2. Assuming particles of mean mass m, we write the set of equations pertinent to the one-dimensional one-fluid case as

∂n ∂t + u ∂n ∂x + n ∂u ∂x = 0 (A.1) m n ∂u ∂t + nu ∂u ∂x + ∂ ∂x nT - 4 3 η ∂u ∂x = 0 (A.2) ∂ ∂t (nT ) + u ∂ ∂x (nT ) + 3nT ∂u ∂x +2 ∂q ∂x - 8 3 η ∂u ∂x 2 + c ε N L v 0 T n = 0 (A.3)
where c is a constant of order unity, which will be computed explicitly for a two species system in Appendix B. The other quantities in equation (A.3) are defined as

v 0 ≡ 2T m 1/2 η ≡ d nT ν 0 = d nT v 0 L N q ≡ -b nT ν 0 m ∂T ∂x = - b 2 nv 0 L N ∂T ∂x
with b and d being two more constants of order unity. We note that for a two species system the normalization velocity v 0 is related to the thermal velocity V 12 of the species relative velocity of Section 2 via

V 2 12 = 2 m m 1 m 2 v 2 0 . (A.4)
The linear analysis of the system is easier if performed in terms of the dimensionless variables τ and U :

dτ = N 2 v 0 L dt (A.5) U = u v 0 . (A.6)
where dτ /dt is a measure of the number of collisions per time unit in the system during the homogeneous phase and U is essentially the flow's Mach number. We assume a static equilibrium u(x, τ ) = 0 with a uniform density n = n 0 , and temperature T (x, τ ) = T 0 (τ ) decreasing according to equation (A.3), i.e. The fact that the entropy mode is a non propagating mode implies that density and temperature fluctuations vary spatially in antiphase, ensuring pressure balance. Indeed, for k → 0, one has -iω = 1 -k 2 /2 (see Figure 8) from where, using the second line of the matrix in equation (A.11), one immediately deduces the linear pressure equilibrium condition δn/n 0 = -δθ which is valid up to order one in k. To same order in k, the first equation in (A.11) implies the δn/n 0 = -ikδU . This latter relation shows that in the entropy mode density and velocity fluctuations are out of phase by π/2 corresponding to the fluid flowing from the low density to the high density parts of the wave. These linear motions, eventually enforced by non linear effects (i.e. velocity profile steepening) at late times, ineluctably drive the system toward the final inelastic collapse.

B Collision frequency for Maxwellian distributions

The number of collisions per time unit of a particle of species α with the particles of species β for relative velocities in the range [u, u + du] is given by dν αβ = |u|f αβ (u)du, u < 0.

(B.1)

In equation (B.1) f αβ (u) denotes the distribution of the relative velocities. If one assumes Maxwellian distributions for both species, the distribution of the relative velocities is also Maxwellian. In one dimensions, with N β particles uniformly distributed in a system of length L it must be that

f αβ (u) = 1 V αβ π 1/2 N β L e -u 2 /V 2 αβ (B.2)
where V αβ is the thermal velocity defined as where m ≡ (m α + m β )/2 and v 2 0 ≡ 2T /m. Using the expression for the collision frequency given by (B.6) it becomes possible to specify the constant c which establishes a relation between the time variable t and the collision index I (cf equations (A.3) and (A.7)). Of course, given the above assumptions, the relation between t and I only holds as long as the inhomogeneities in the system are weak. By comparing (A.7) and ( 14) one obtains the following relation between the pseudo time variable τ and the collision index I: We note that c does not depend on the temperature, as long as energy equipartition is a valid approximation.

V 2 αβ ≡ V 2 α + V 2 β = 2 T α m α + 2 T β m β . (B.

Fig. 1 .

 1 Fig. 1. Homogeneous regime for the case of r 2 = 0.9995 and N = 19600 particles half of which have mass m1 = 1 and the other half m2 = 4 (Run 1). Shown are the velocity distribution histograms for both species after 8 × 10 7 collisions. The solid line represents the Maxwell-Boltzmann distribution corresponding to the measured particles' mean energy.

Fig. 2 .

 2 Fig. 2. Energy versus collision index I for Run 1. The exponential part of the profile closely follows the E = E0 exp(-I/2 × 10 7 ) profile of the fluid theory.

  ) (see equation (A.7) together with the relation (B.7)) which predicts a per collision energy loss r 2 = 0.9995 and N = 19600 one has δE/E ≈ 5.1 × 10 -8 which corresponds to the value 2 × 10 7 observed in Figure 2.

Fig. 3 .

 3 Fig.3. Run 1: Temperature and density profiles for the light particles (i.e. m = 1) normalized to their average value at the end of the simulation. Note the high concentration and low temperature within the cluster near x = 0.416. Also note the presence of a second cluster forming near x = 0.8.

Fig. 4 .

 4 Fig. 4. Light particles histogram (m=1) at the end of Run 1. Bin values represent the fractional number of light particles n1/n. Each bin contains 400 particles. The main cluster includes all particles with index l in the range 5000 to 11000.

Fig. 5 .

 5 Fig. 5. Function Y (x) (see equation (15)) for Run 1. The slope of the function Y is a measure of the heavy (positive slope) or light (negative slope) particle excess.

x 3 :Fig. 7 .

 37 Fig. 7. Inhomogeneous regime: scatter plot for the m = 4 particles for all runs at the time of inelastic collapse.

Fig. 8 .

 8 Fig.8. Growth rate for the both the sound mode and the entropy mode for a particular value of the thermal diffusion coefficient B and the shear viscosity coefficient D. Plotted is also the real frequency of the sound mode ω(k) and the asymptotic values of the growth rates for both modes in the limit k 1, and D = B = 0.

3 ) 4 ) 2 N 5 )

 3425 Integration of equation (B.1), assuming the Maxwellian distribution (B.2), gives the collision frequency of a particle α with the particles of species β:The total collision frequency ν in the system is obtained by adding the contributions from all kinds of collisions, i.e.ν = N α ν αβ + N β ν βα + N α ν αα + N β ν ββ =Assuming energy equipartition T = T α = T β and an equal number of particles for both species N α = N β = N/2, one ends up with ν(m α , m β , T, N ) =

2δI = cδτ. (B. 7 )

 7 Given the relation between τ and the time variable t (see equation (A.5)) and using the fact that the collision frequency is just ν = δI/δt it follows that

Table 1 .

 1 Basic parameters for the 3 simulations discussed in the paper. The total number of particles in each simulation is N = 19600 with half of the particles with mass m1 = 1 and the other half with mass m2 = 4.

	Run ε	r 2	N ε
	1	5 × 10 -4 0.9995 9.8
	2	2 × 10 -3 0.998	39.2
	3	2 × 10 -2 0.98	392

Table 2 .

 2 Characteristic time scales of the entropy mode for the smallest possible wave vector k = kmin ≡ 2π/N ε. The column labelled Supersonic indicates whether or not the fluid motions are supersonic at the time of inelastic collapse. γ refers is the linear growth rate of the mode, τs and τn are defined in equations (25) and (26), respectively. The parameters of the runs are given in Table1.

	Run Supersonic kmin	γ -1	τs	τn
	1	no	0.64	4.9 × 10 7 2.3 × 10 8 2.5 × 10 8
	2	yes	0.160 1.0 × 10 7 4.6 × 10 7 6.4 × 10 7
	3	yes	0.016 1.0 × 10 6 4.5 × 10 6 8.6 × 10 6

During the linear phase of the instability c is a constant which only depends on the species relative densities and on the species particle mass (see Appendix B). For example, if N/2 particles have mass m 1 = 1 and N/2 particles m 2 = 4, one has c = 0.97. Thus, c = 0.97 must be used for all runs discussed in the paper except Run 1, where all particles are identical and for which the linear theory presented here is not applicable.

Let us now investigate the response of the equilibrium characterized by n = n 0 and T = T 0 (τ ) to small perturbations δn, δθ and δU by setting

The coefficients of the resulting linear system are independent of τ and x so that we can solve it using standard Fourier techniques. We therefore assume perturbations of the type δn, δθ, δU ∝ exp[i(kx-ωτ )] which, when plugged into the linearized system, lead to

where M is the 3 × 3 matrix given by

with k ≡ kL/(N ε), ω ≡ ωN/ε, B ≡ εb and D ≡ εd. Non trivial solutions of the system correspond to a vanishing determinant of the matrix M. Splitting the complex frequency ω = ω r + iγ into real and imaginary parts ω r and γ, and assuming real wave vectors k leads to two equations for the real and the imaginary part of the condition det M = 0 (in the remaining of this appendix we drop the tildes to ease readability):

and either

for the modified sound wave, or

for the modified, non propagating, entropy wave. The dispersion relation for sound waves, ω 2 r = 3k 2 /2, is immediately obtained by setting γ = B = D = 0 in equation (A.14) (we shall remember that frequencies are measured in terms of the pseudo time τ and not the real time t). Switching back to real time using equation (A.5), the dispersion relation reduces to the more familiar relation ω 2 r = 3T k 2 /m for a one dimensional system. The growth rate γ(k) for both, sound and entropy mode, can be obtained by substituting ω 2 r from equation (A.14) or (A.15), respectively, into (A.13). After some juggling one ends up with

for the sound mode, and

for the entropy mode. For small values of B and D, i.e. at sufficiently large spatial scales such that both thermal conduction time scale 1/Bk 2 and the viscous time scale 1/Dk 2 both exceed the linear time scale 1/ω, one may use the asymptotic expressions for the imaginary part of both the sound wave and the entropy wave, viz.

and

We note that f s,e (0, k, 0, 0) < 0 for all values of k. Since f s,e (γ → ∞) = ∞, a positive real solutions γ > 0 must exist for both equations (A.18) and (A.19). Thus, both the sound and the entropy mode are always unstable for k → 0. There is indeed, only one real solution of (A.18) which goes as k 2 for |k| → 0 and tends towards the asymptotic value 1/6 for |k| → ∞. Similarly, for |k| → ∞, the entropy mode's growth rate tends towards 2/3. Of course, for sufficiently large values of k both modes must be damped by diffusive effects. However, unlike the entropy mode, which is always unstable for |k| → 0, the sound mode is completely stabilized for thermal diffusivity values B > 1/2.


Molecular and cellular basis of depigmentation in vitiligo patients

Véronique Delmas 1-3 and Lionel Larue 1-3 (1) Institut Curie, PSL Research University, INSERM U1021, Normal and Pathological Development of Melanocytes, Orsay, France,

(2) Univ Paris-Sud, Univ Paris-Saclay, CNRS UMR 3347, Orsay, France, (3) Equipe Labellisée Ligue Contre le Cancer, Orsay, France

Abstract

Vitiligo is a chronic skin disease characterized by the appearance of zones of depigmentation. It is mostly described as an autoimmune disease in which the immune system destroys the melanocytes. Consistent with this origin, genetic studies have implicated genes encoding proteins mediating the immune response targeting melanocytes in the etiology of this disease, together with proteins specific to these cells. However, the destruction of melanocytes by the immune system is neither global nor complete, because the patients do not display total depigmentation. The etiopathology of vitiligo is clearly complex and cannot be simply reduced to an autoimmune reaction directed against pigmented cells. Intrinsic changes have been observed in the melanocytes, keratinocytes and dermal cells of vitiligo patients.

Identification of the molecular and cellular changes occurring in normally pigmented skin in vitiligo patients, and an understanding of these changes, are essential to improve the definition of trigger events for this disease, with a view to developing treatments with long-term efficacy. This review focuses on the early events identified to date in the non-lesional regions of the skin in vitiligo patients, and discusses the process of repigmentation from melanocyte stem cells.

Changes to the epidermis before depigmentation

1.a. Changes in epidermal structure

The structure of the epidermis, including its thickness and the tegumental folds, differs between anatomic sites. At a given site, the structure of the epidermis and the number of melanocytes are similar in all individuals, regardless of their phototype and ethnic origin. The epidermis in lesional regions of the skin is abnormally thick in patients with non-segmental vitiligo. Interestingly, this thickening of the epidermis is also observed in non-lesional regions of the skin located some distance away from the depigmented zones. Histological analyses of biopsy specimens from undamaged areas normally pigmented skin in patients without vitiligo, patients with nonsegmental vitiligo (NSV, characterized by bilateral depigmented areas) and patients with segmental vitiligo (SV, characterized by highly localized, unilateral depigmentation) have revealed that the epidermis has different characteristics specific to these groups. Patients with NSV display a thickening of the spinous layer (acanthosis) of the epidermis, whereas SV patients have a larger number of tegumental folds than patients without vitiligo ((1) and Figure 1). However, the molecular basis of these histological differences remains unknown. The thickening of the epidermis is observed in normally pigmented skin of NSV patients compared to the skin of control patients. The thickening is only observed at the level of the spinous layer (stratum spinosum). It is likely to involve changes in the balance between keratinocyte proliferation and differentiation. One explanation could be an abnormal production of cytokine/growth factors in the entire skin of vitiligo patients and/or modifications of keratinocyte genes involved in proliferation and differentiation.

A transcriptomic analysis comparing normal appearing non-lesional skin of NSV patients with healthy patient skin reveals an up-regulation of Epiregulin (EREG) encoding a secreted peptide hormone, member of the epidermal growth factor (EGF) family of proteins (2). Our hypothesis is that Epiregulin promote cell proliferation of the keratinocytes of the spinous layer. However, the overall anatomy of the epidermis is preserved.

The epidermal rete ridges in the skin of SV patients probably reflect a specific modification (or modifications) at the basal level, between the dermis and epidermis.

Changes in dermal fibroblasts have been reported in patients with NSV, but few studies have considered the nature of the dermal compartment involved in patients with SV [START_REF] Bastonini | Skin Pigmentation and Pigmentary Disorders: Focus on Epidermal/Dermal Cross-Talk[END_REF][START_REF] Kovacs | Vitiligo Skin: Exploring the Dermal Compartment[END_REF]. A recent study showed that altered cross-talk between epidermal and dermal component can be detected in normally pigmented areas of vitiligo skins (4).

These changes include the presence of myofibroblastic cells and increase secretion of various factors such as IL-6 and HGF ligands, and the ECM-protein fibronectin.

Unfortunately this study did not address the differences between dermal components of SV, NSV and control skins. This hypothesis would be important to test in the future. Since the histological modifications in NSV and SV are observed in biopsies taken from clinically normal pigmented looking skin, it most likely not a compensatory response to pigment loss in non-lesional skin. These observations suggest that the epidermal and dermal microenvironment of the melanocytes is different between NSV and SV patients, even before depigmentation. These differences may have consequences for the initiation and progression of the disease, which is likely to be extensive and unpredictable in NSV patients, but rapid and very limited in SV patients.

1.b. Adhesion defect

Within the epidermis, melanocytes and keratinocytes interact with each other directly, through cell-to-cell adhesion molecules. These interactions are favored by the secretion into the microenvironment of the various factors involved in proliferation and differentiation, in a free form or in vesicles [START_REF] Hirobe | How are proliferation and differentiation of melanocytes regulated? Pigment cell & melanoma[END_REF]. The integrity of this association of melanocytes and keratinocytes, known as the epidermal melanization unit, is essential for pigmentation. In vitiligo patients, the epidermal melanization unit is no longer present in the depigmented zones, but is also altered in the normally pigmented zones. Indeed, melanocytes are absent from the depigmented zones and there are about 25% fewer melanocytes in the basal area of the normally pigmented zones of affected patients than in that of unaffected individuals (6). E-cadherin, a cellular adhesion molecule with calcium-dependent activity, is considered to be an epithelial marker. It is responsible for adhesion between melanocytes and keratinocytes and, thus, for the integrity of the epidermal melanization structure. Ecadherin plays an essential role in embryonic development and in the progression of many carcinomas [START_REF] Larue | E-cadherin null mutant embryos fail to form a trophectoderm epithelium[END_REF][START_REF] Larue | A role for cadherins in tissue formation[END_REF][START_REF] Takeichi | Cadherin cell adhesion receptors as a morphogenetic regulator[END_REF]. There is less E-cadherin on the surface of melanocytes than on epithelial cells, including keratinocytes in particular. Melanocyte-keratinocyte adhesion is, therefore, much weaker than keratinocyte-keratinocyte adhesion. We have shown that E-cadherin is absent, or present in only very small amounts, at the membrane of melanocytes from the pigmented skin of patients with vitiligo, regardless of its type (NSV or SV; Figure 2). This absence or low level of E-cadherin on the surface of melanocytes is linked to the smaller number of melanocytes in the basal layer of the epidermis, and the larger number of these cells in the upper layers (spiny and granular). The change in the amount of E-cadherin at the melanocyte membrane is similar in the two types of vitiligo (NSV and SV). E-cadherin deficiency is not, therefore, sufficient in itself to account for the difference in the progression of depigmentation in these two types of vitiligo. The decrease in the amount of Ecadherin, or its total loss at the surface of melanocytes is a change that does not, in itself modify pigmentation. Nevertheless, the decrease/loss of E-cadherin weakens melanocyte adhesion to keratinocytes, and to the basal layer.

Mechanical and oxidative stresses aggravate depigmentation in vitiligo [START_REF] Picardo | Vitiligo[END_REF]. We have shown, both in vitro and in vivo, that if mechanical and oxidative stresses are associated with a decrease/loss of E-cadherin in melanocytes, these cells detach from the basal layer and move to the suprabasal level of the epidermis. In the long term, melanocytes are completely eliminated from the treated zone. We can hypothesize that the detached melanocytes are periodically replaced by melanocytes produced from melanocyte stem cells (McSCs). This replacement is accelerated by chronic or continuous stresses, but remains effective until the McSC reservoir is exhausted. Furthermore, half the suprabasal melanocytes of the pigmented epidermis of vitiligo patients are apoptotic or necrotic and, thus, potentially antigenic.

An autoimmune reaction directed against abnormal and apoptotic/necrotic suprabasal melanocytes may then be induced, to complete the destruction and elimination of these cells from the epidermis of the affected area.

1.c. Increase in oxidative stress

Oxidative stress is increasingly being identified as a major determinant in the etiology of vitiligo [START_REF] Dell'anna | Membrane lipid defects are responsible for the generation of reactive oxygen species in peripheral blood mononuclear cells from vitiligo patients[END_REF][START_REF] Glassman | Vitiligo, reactive oxygen species and T-cells[END_REF]. An increase in reactive oxygen species (ROS) has been observed in the lesional and non-lesional epidermis of vitiligo patients. These ROS may be produced by exposure of the skin to UV irradiation, for example. An accumulation of hydrogen peroxide, together with changes in the antioxidant machinery, such as a decrease in catalase levels or an increase in superoxide dismutase levels, has been demonstrated in vitiligo [START_REF] Schallreuter | In vivo and in vitro evidence for hydrogen peroxide (H2O2) accumulation in the epidermis of patients with vitiligo and its successful removal by a UVB-activated pseudocatalase[END_REF][START_REF] Maresca | Increased sensitivity to peroxidative agents as a possible pathogenic factor of melanocyte damage in vitiligo[END_REF][START_REF] Dammak | Antioxidant enzymes and lipid peroxidation at the tissue level in patients with stable and active vitiligo[END_REF]. Oxidative stress can be detected as the presence of lipid peroxides in the membrane of epidermal cells from normally pigmented skin in patients with vitiligo (Figure 3). In the limited number of samples from NSV patients, SV patients and controls analyzed for lipid peroxidation (≈20), this peroxidation (visualized by labeling with 4-HNE antibodies) affected all cells (keratinocytes and melanocytes) in the epidermis of NSV patients, whereas less intense, discontinuous peroxidation was observed in the epidermis of SV patients, and almost no peroxidation was observed in the skin of controls. The structure of the epidermis and the oxidation state of the epidermis from normally pigmented skin in NSV and SV patients merit more detailed analysis, at a large scale and at different time points in the progression of the disease, to provide solid, factual information about these aspects of the disease.

1.d. Adhesion and oxidative stress: two closely linked processes

The increase in oxidative stress and the adhesion defect in the epidermis are not independent events. Instead, they are closely linked through multiple mechanisms.

Indeed, an increase in oxidative stress effectively decreases the amount of Ecadherin at the cell membrane, due to activation of the Src kinase, which destabilizes the E-cadherin/β-catenin complex at the membrane by phosphorylation [START_REF] Guntaka | Epidermal growth factor protects the apical junctional complexes from hydrogen peroxide in bile duct epithelium[END_REF][START_REF] Chan | Major role of epidermal growth factor receptor and Src kinases in promoting oxidative stress-dependent loss of adhesion and apoptosis in epithelial cells[END_REF]. Furthermore, lipid peroxidation modifies the organization of lipids within cell membranes. This disturbance of the plasma membrane lipids alters the distribution of membrane proteins, such as E-cadherin. Finally, a decrease in E-cadherin induces the activation of nuclear factor E2-related factor 2 (Nrf2), which protects against oxidative stress, and the expression of which is induced in the skin of vitiligo patients [START_REF] Kim | E-cadherin inhibits nuclear accumulation of Nrf2: implications for chemoresistance of cancer cells[END_REF][START_REF] Natarajan | Transcriptional upregulation of Nrf2dependent phase II detoxification genes in the involved epidermis of vitiligo vulgaris[END_REF]. Oxidative stress probably decreases the amount of E-cadherin at the membrane in all epidermal cells, but this decrease is particularly marked for melanocytes, which already have smaller amounts of E-cadherin at their surface than keratinocytes. As oxidative stress directly affects the amount of E-cadherin at the melanocyte membrane and consequently the number of melanocytes at the base of the epidermis, the peroxidation state of membrane lipids may be considered a predictive marker for the progression of depigmentation in vitiligo patients. Finally, the entry of suprabasal melanocytes into apoptosis or necrosis may be induced jointly by the loss of melanocyte adhesion to the basal membrane rich in growth and survival factors and the presence of oxidative stress in the epidermis, affecting the integrity of the endoplasmic reticulum and favoring the accumulation of incorrectly folded proteins and cell apoptosis.

Melanocyte stem cells and repigmentation

2.a. Melanocyte stem cells

Understanding the molecular and cellular mechanisms underlying the establishment of melanocytes in the epidermis and their renewal is a fundamental issue for the development of effective treatment strategies for depigmentation diseases. It is very difficult to perform fundamental research studies on human McSCs. The mouse model, in which genetic modification is easy, is the most widely used system for investigating the mechanisms of mammalian stem cell activation and the de novo production of melanocytes in the epidermis. The principal reservoir of McSCs is the bulge of the hair follicle [START_REF] Nishimura | Melanocyte stem cells: a melanocyte reservoir in hair follicles for hair and skin pigmentation[END_REF]. Other reservoirs of melanocyte stem cells have been identified in the interfollicular epidermis and dermis of glabrous skin [START_REF] Seleit | Immunohistochemical study of melanocytemelanocyte stem cell lineage in vitiligo; a clue to interfollicular melanocyte stem cell reservoir[END_REF][START_REF] Li | Human dermal stem cells differentiate into functional epidermal melanocytes[END_REF] [START_REF] Goldstein | Narrow Band Ultraviolet B Treatment for Human Vitiligo Is Associated with Proliferation, Migration, and Differentiation of Melanocyte Precursors[END_REF], and in the sweat glands [START_REF] Okamoto | A melanocyte--melanoma precursor niche in sweat glands of volar skin[END_REF]. In the follicular bulge region, the MsSCs interact closely with hair follicle stem cells (HFSCs), which provide them with their functional niche. This interaction facilitates the growth of the pigmented hair through crosstalk between various signaling pathways, including the TGFβ and Wnt pathways [START_REF] Rabbani | Coordinated activation of Wnt in epithelial and melanocyte stem cells initiates pigmented hair regeneration[END_REF][START_REF] Tanimura | Hair follicle stem cells provide a functional niche for melanocyte stem cells[END_REF].

The hair cycle has three phases: anagen (active hair growth), catagen (regression of the transitory zone of the hair) and telogen (rest phase before the next hair cycle).

McSCs are activated during anagen, by the Shh or Wnt signaling pathways, or through β-catenin activation [START_REF] Lowry | Defining the impact of beta-catenin/Tcf transactivation on epithelial stem cells[END_REF]. The entry of McSCs into anagen may also be triggered by environmental factors. For example, in response to UVB irradiation, McSCs leave quiescence, proliferate and concomitantly differentiate into melanocytes, migrating from the hair follicles to the interfollicular epidermis [START_REF] Chou | Direct migration of follicular melanocyte stem cells to the epidermis after wounding or UVB irradiation is dependent on Mc1r signaling[END_REF].

Consequently, narrow-band UVB irradiation is often used to stimulate melanocyte production in repigmentation therapies for vitiligo patients. β-catenin is both involved in vitiligo depigmentation and required for repigmentation. Indeed, like E-cadherin, βcatenin is absent from the membrane of the melanocytes in the pigmented skin of vitiligo patients, and a decrease in the transcription induced by the WNT/β-catenin pathway is observed in the depigmented epidermis of vitiligo patients (Figure 2 and [START_REF] Wagner | Altered E-Cadherin Levels and Distribution in Melanocytes Precede Clinical Manifestations of Vitiligo[END_REF][START_REF] Regazzetti | Transcriptional Analysis of Vitiligo Skin Reveals the Alteration of WNT Pathway: A Promising Target for Repigmenting Vitiligo Patients[END_REF]. UVB irradiation enhances β-catenin activity, by promoting the translocation of this molecule to the nucleus, leading to the production of differentiated melanocytes in the epidermis [START_REF] Bertrand | UVB represses melanocyte cell migration and acts through beta-catenin[END_REF]. However, β-catenin activity should be induced with caution, because this protein is oncogenic, and melanocyte production may lead to transformation into melanoma [START_REF] Delmas | Beta-catenin induces immortalization of melanocytes by suppressing p16INK4a expression and cooperates with N-Ras in melanoma development[END_REF][START_REF] Gallagher | Beta-catenin inhibits melanocyte migration but induces melanoma metastasis[END_REF].

2.b. Exhaustion of the McSC reservoir in the depigmented zones

Whatever the mechanism by which melanocytes are eliminated from the epidermis in vitiligo patients, they are not replaced in the depigmented epidermis. Two nonmutually exclusive hypotheses can be put forward: melanocytes may be produced

from

  McSCs and then systematically destroyed by the cells of the immune system, or the self-renewal capacity of the McSCs is exhausted. Obviously, these two mechanisms may operate together during the various phases of vitiligo, with one predominating over the other at a particular point in disease progression. In elderly humans and mice, pigmented melanocytes are present in the bulge region and are associated with exhaustion of the McSC reservoir (33). Moreover, it has been shown that the amount of McSC is not affected in the bulge and the precursors found in the infundibulum was decreased (23). In other words, the presence of McSC in the bulge or in the infundibulum can be used as an indicator of the state of exhaustion of McSC reservoir and of treatment efficacy. In conclusion, we propose one of the potential sequential events for vitiligo pathogenesis: 1) E-cadherin is involved in triggering vitiligo, 2) melanocyte detachment and exhaustion of the McSC reservoir are facilitated by oxidative/mechanical stress, and 3) activation of the immune response induced by the presence of apoptotic/necrotic melanocytes in the suprabasal epidermis during disease progression. program launched by the French Government and implemented by ANR Labex CelTisPhyBio (ANR-11-LBX-0038 and ANR-10-IDEX-0001-02 PSL).

Figure 1 :

 1 Figure 1: Epidermal structures are different in (A) control (CTRL), (B) non-segmental

Figure 2 :

 2 Figure 2: E-cadherin (A,B,E,F,I,J) and β-catenin (C,D,G,H,K,L) membranous

Figure 3 :

 3 Figure 3: Oxidative stress is more abundant in NSV (B) than in SV (C), it is quasi-

Figure 4 :

 4 Figure 4: Schematic representation of the normal, NSV and SV normally pigmented

  Figure 1

Acknowledgments

We thank the members of the laboratory who participated to the work. This work was supported by ARC, INCa, ITMO Cancer, and part of the "Investissements d'Avenir"