Dynamical demixing of a binary mixture under sedimentation
A. S. Nunes R. C. V. Coelho V. C. Braz M. M. Telo da Gama N. A. M. Araujo

	
	 
	
	
	Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa,
		P-1749-016 Lisboa, Portugal
    Departamento de Física, Faculdade
		de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal
	Electric Ant Lab, Science Park 106, 1098 XG, Amsterdam, The Netherlands
	
    
	Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa,
		P-1749-016 Lisboa, Portugal
    Departamento de Física, Faculdade
		de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal	
	
	
	Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa,
		P-1749-016 Lisboa, Portugal
    Departamento de Física, Faculdade
		de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal
		
	
	Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa,
		P-1749-016 Lisboa, Portugal
    Departamento de Física, Faculdade
		de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal
	
	
	nmaraujo@fc.ul.pt
	Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa,
		P-1749-016 Lisboa, Portugal
    Departamento de Física, Faculdade
		de Ciências, Universidade de Lisboa, P-1749-016 Lisboa, Portugal
	
	
		
		We investigate the sedimentation dynamics of a binary mixture, the species of which differ by their Stokes coefficients but are identical otherwise. We analyze the sedimentation dynamics and the morphology of the final deposits using Brownian dynamics simulations for mixtures with a range of sedimentation velocities of both species. In addition, we use the lattice Boltzmann method to study hydrodynamic effects. We found a threshold in the difference of the sedimentation velocities above which the species in the final deposit are segregated. The degree of segregation increases with the difference in the Stokes coefficients or the sedimentation velocities above the threshold. We propose a simple analytical model that captures the main features of the simulated deposits. 
		
	
	
	Dynamical demixing of a binary mixture under sedimentation
    Nuno A. M. Araújo
    26.7.2018
==========================================================



§ INTRODUCTION

	
	
	The process of sedimentation, where particles in suspension settle in the presence of a gravitational field, is ubiquitous over a wide range of length scales <cit.>. For example, sedimentation plays a relevant role in natural water transport, affecting the chemical composition of the seabed <cit.> and the water quality in reservoirs <cit.>. At the other end of the scale, sedimentation by ultracentrifugation is used as an analytical tool in medical, biological and pharmaceutical applications, where the constituents of a suspension are separated by molecular weight <cit.>. At the fundamental level, sedimentation experiments were developed and used extensively in statistical physics and colloidal science to evaluate the equation of state of hard spheres <cit.> and to study the phase diagram of colloidal particles <cit.>. 
	
	Studies of the sedimentation of mixtures of particles that differ in their buoyant mass revealed a rich phase stacking diagram under thermodynamic equilibrium conditions <cit.>. The structure of the final deposit depends not only on the difference in buoyant masses but also on the particle-particle interactions <cit.>. The roughness of the particle surface is known to affect the hydrodynamics of the surrounding fluid, e.g., alters the lubrication film thickness <cit.>. These conditions could be realized in an experiment with particles composed by a rigid core covered by an elastic surface layer, which would affect the hydrodynamics while the pairwise interactions remain dominated by the rigid cores. In order to shed light on the role of the hydrodynamic radius on the sedimentation dynamics, we consider a binary mixture of particles that differ only through their Stokes coefficient when using molecular dynamics. Following a method developed previously <cit.>, we consider that the particles differ by their Stokes coefficients only, and are identical otherwise. Thus, the thermodynamic phase is perfectly mixed and demixing, if it occurs, is dynamically driven. Hydrodynamic interactions are know to be relevant in different limits, triggering, for example, a number of different instabilities <cit.>. To account for hydrodynamic effects, we consider a complementary set of simulations by modelling the fluid-particle interaction with the lattice Boltzmann method (LBM). 
	
	When thermal fluctuations are negligible, colloidal particles in solution are expected to sediment with a sedimentation velocity that depends only on the strength of the gravitational field and their Stokes coefficient. Thus, distinct Stokes coefficients imply different sedimentation velocities. In what follows, we show that the morphology of the final deposit depends crucially on the ratio of the sedimentation velocities. Above a certain threshold, which will be quantified below, the particles are segregated in the final deposit, as they arrive at the substrate at different rates and do not have time to relax to the thermodynamic equilibrium mixed state. We investigate this segregation and discuss its dependence on the model parameters.
	
	The paper is organized in the following way. In Section <ref>, we describe the model and the details of the simulations. Results from the particle-based simulations (Brownian dynamics), LBM and an analytical model are discussed in Section <ref>. Finally, we draw some conclusions in Section <ref>. 
	
	

§ MODEL AND SIMULATIONS 

	
	
	

 §.§ Molecular dynamics

	We consider a binary mixture of identical spherical particles where the two species are characterized by distinct Stokes coefficients. The particles are in a uniform gravitational field along the vertical direction (y-direction) and inside a rectangular two-dimensional box of width L_x and height L_y. The boundary conditions are periodic in the x-direction and are rigid walls in the y-direction. The trajectory of each particle i is obtained by solving the Langevin equation, in the overdamped regime,
	
	
    γ_i dr⃗_i/dt = -∇ _i[∑_j^N_t V_ij(r) ] +mg⃗+ ξ⃗_i,   j≠ i

	
	where, r⃗_i is the position of particle  i, V_ij is the pairwise potential, N_t the total number of particles, m is a parameter that takes into account the effects of mass and buoyancy of the particles, g⃗=-ge⃗_y the gravitational field, ξ⃗_i a stochastic force, and γ_i is the Stokes coefficient. The two species differ through the values of γ_i: γ_f for fast particles and γ_s for slow ones, such that γ_f<γ_s. The diffusion coefficient of each species is also different, as given by the Stokes-Einstein relation D_i=k_BT/γ_i. As a result, the two species have different sedimentation velocities, v⃗_i=m/γ_ig⃗ since the the particles have the same mass. The fluid is in thermodynamic equilibrium at a thermostat temperature T and hydrodynamic effects are neglected. Thus the time series of the stochastic force is drawn from a Gaussian distribution with zero mean and uncorrelated second moments in time and space, given by ⟨ξ_i^k(t)ξ_i^l(t')⟩ = 2k_BTγ_iδ_klδ(t-t'), where k and l refer to the coordinates of the vector ξ⃗_i. 
	
	In order to focus on purely dynamical effects, we consider that the particle-particle interactions are identical for all the particles. We describe this interaction through a (repulsive) Lennard-Jones potential, truncated at a cut-off distance r_cut=2^1/6σ,
	
	
    V_ij(r)=ϵ[ ( σ/r) ^12-( σ/r) ^6],

	
	where ϵ sets the energy scale and σ the size of the particles. Thus, the potential depends only on the distance r=|r⃗_i-r⃗_j| between the particles i and j.
	
	Hereafter, σ sets the unit of length. The energy is expressed in units of k_BT, time is defined in units of the Brownian time τ=σ^2γ (k_BT)^-1 and the strength of the external field, g, is given in units of k_BT/(m σ). Equation <ref> is integrated using a second-order stochastic Runge-Kutta numerical scheme, proposed by Brańka and Heyes  <cit.>, with a time-step of Δ t=10^-4τ. Initially, the particles are distributed uniformly at random (without overlapping) in the simulation box. Unless stated otherwise, we set ϵ=1 and g=12. The box size is L_y=200 and L_x=37.5 and the binary mixture consists of N=3000 particles, with N/2 particles of each species. The initial number density is ρ_0 = 0.4.  
	
	
	

 §.§ Lattice Boltzmann method

	
	
To account for hydrodynamic effects, we employ a Lattice Boltzmann method <cit.> to describe the motion of the suspended medium and couple it to a discrete element simulation of the particles. For the fluid, the position x⃗ is discretized on a regular lattice with lattice spacing Δ x and the velocity of the particles in the suspending fluid is discretized according to the D3Q19 lattice (with 19 vectors in three dimensions): the weights for the different velocity vectors are w_α=1/3 for |𝐜|^2=0, w_α=1/18 for |𝐜|^2=1Δ x and w_α=1/36 |𝐜|^2=2Δ x, where α=1, …, 19 labels the velocity vector. The speed of sound in this lattice is c_s=1/√(3). The LBM is based on solving the dynamics of the distribution function f_α, which is governed by the Boltzmann equation. Its discretized version with the single relaxation time collision operator is:

    f_α (x⃗+c⃗_α, t+Δ t)  -  f_α(x⃗,  t) =  - f_α - f^eq_α/τ_fΔ t ,

where Δ t is the time step and τ_f is the relaxation time which sets the kinematic viscosity of the fluid ν=c_s^2(τ_f-1/2). The equilibrium distribution function is given by:

    f^eq_α = ρ w_α[ 1 + c⃗_α·u⃗/c_s^2 - u^2/2c_s^2 + (c⃗_α·u⃗)^2/2c_s^4],

where ρ is the density and u⃗ is the macroscopic velocity of the fluid, calculated as follows:

    ρ = ∑_α f_α ,   u⃗ = 1/ρ∑ _αc⃗_α f_α  .

The simulation results using LBM are expressed in lattice units, in which ρ=1, Δ t =1 and Δ x=1.


At the solid surfaces moving with velocity u⃗_s, the no-slip condition is applied though the bounce-back boundary condition:

    f_α̅(x⃗_s,t+Δ t) = f_α^∗(x⃗_s,t) - 2w_αρ_s c⃗_α·u⃗_s/c_s^2,
	
where α̅= -α represents the opposite velocity vector, f_α^∗ is the distribution propagated towards the solid surface and ρ_s is the average density of the fluid node neighbors. The drag force on the solid is calculated using the momentum exchange method <cit.>:

    F⃗_⃗d⃗ = Δ x ^3/Δ t∑_x⃗_s, α[   2f_α^∗(x⃗_s,t) - 2w_αρ_s c⃗_α·u⃗_s/c_s^2] c⃗_α,

which uses the reflected distributions in the bounce-back f_α^∗. When two particles approach each other and there are few or no fluid nodes between their surfaces, Eq. <ref> becomes inaccurate since the space between the particles is interpreted as vacuum. To correct for this, we consider a hydrodynamic lubrication force, which is repulsive and points in the direction connecting the center of the particles of strength <cit.>:

    F_l = 3πρν R^2 v_r/2h,

where v_r is the relative velocity between the particles, h is the minimum distance between their surfaces and R is the radius of the particles. This force diverges as h→ 0; thus we introduce a cutoff at h=0.01R. All the forces acting on the particles are combined to describe its trajectory using the Verlet method. As the particle moves in the fluid, liquid nodes are destroyed and created. For the destroyed fluid nodes, the fluid information is erased and the fluid momentum is transferred to the solid. For the new fluid nodes as in Ref. <cit.>, the density is set to the average density of the neighbouring fluid nodes in the first belt ρ̅ (considering all nodes reached by the lattice vectors of D3Q19) and the velocity field to the solid velocity u⃗_s. The distribution function of these new nodes is the corresponding equilibrium distribution, f_α = f^eq_α(ρ̅, u⃗_s).

The collisions between particles i and j are mediated by an elastic force <cit.>: 

    F⃗^el_ij =  
    
        0,   ε_ij=0
    
        k_lε_ij,   dε_ij/dt≥ 0  
    
        k_uε_ij,   dε_ij/dt < 0  ,

where ε_ij is the overlapping vector, defined as:

    ε_ij =  
    
        0,   if  2R - |r⃗_ij|≤ 0
    
        2R - |r⃗_ij|   otherwise.

Here, r⃗_ij is the vector connecting the centers of the two particles
from j to i and k_l and k_u are the elastic coefficients when the
colliding particles are approaching or moving away from each other, respectively. The interaction with the walls also follows Eq. <ref>, but considering the distance of the particle to the wall.

In the LBM simulations, we use the following parameters, all in lattice units. The dimensions of the simulation box are L_x × L_y = 113 × 413 and the relaxation time is τ_f=0.7. We simulate 100 circular particles with radius R=5, the same density as the fluid, k_u=0.1 and k_s=0.2, up to t=4× 10^6, which is sufficient for all the particles to sediment and stop moving.  In order to consider particles with different terminal velocities, we apply different external forces to each species (fast and slow). To the slow particles we apply an acceleration (or force density) a_s = 3× 10^-6 and, to the fast ones, a_f = a_s/β, where β≤ 1 is the ratio of the forces. As will be discussed in Sec. <ref> we find that β corresponds to the ratio between the velocities of the slow v_s and fast particles v_f for single particles: β = v_s/v_f . The particles of each species are chosen uniformly at random and we consider 10 samples for each value of β.

	
	
	

§ RESULTS 

	
	
	
		
	In the overdamped regime and neglecting thermal fluctuations, one expects that single particles move with a constant sedimentation velocity given by v_i=mg/γ_i. The rate of particle accumulation at the bottom depends on the flux density of each species at the growing front of the deposit, 𝒥_i=ρ_iv_i, where ρ_i is the particle number density of species i. Differences in the flux density of each species result in demixing along the vertical direction as particles accumulate at different rates on the bottom and thermal fluctuations are not strong enough to promote mixing. In what follows, we set the particle densities to be the same (equimolar mixture) and vary their velocities only. The demixing that occurs during sedimentation is, therefore, purely dynamical in nature. The relevant control parameter is the ratio between the sedimentation velocities of the two species, v=v_s/v_f=γ_f/γ_s. In order to investigate the dependence on this parameter, in the results that follow, we fix γ_f and vary v by changing γ_s, i.e., by changing v_s, with v_f constant. The degree of demixing depends on v as seen from the final deposit density profiles ρ_f(y) and ρ_s(y) in Fig. <ref>. When  v=1 the density profiles are identical as the particles are indistinguishable. When v ≠ 1, the final deposit can be divided into two regions: one, at the bottom, where the density of the fast particles is higher than the density of the slow ones and another, at the top of the deposit, composed essentially by slow particles. The difference in the particle densities in the first region and the thickness of the second region increase as v decreases (see Figs. <ref> (b) and (c)). 
		
	
	
	
	
	In order to characterize the segregation along the vertical direction, we define a parameter 
	
    Φ  = 1/L_y'∫_0^L_y'|ρ_f-ρ_s|/ρ_f+ρ_s dy,
 
	
	where L_y' corresponds to the height at which the last moving (slow) particle is located, given by L_y'= L_y-v_st. 
		
	
	
	
	

 §.§ Numerical results 

	
	To evaluate numerically the integral in Eq. (3), we divided the simulation domain into horizontal slices of height Δ y=1.5 and width L_x. The integral is then converted into a sum,
	
	
    Φ_n(t)= 1/N_b∑_i^N_b| N_f - N_s|/N_f + N_s,

	
	where N_b is the number of slices and, N_f and N_s are the number of fast and slow particles in each slice, respectively. This parameter is one if the species are completely segregated and zero if they are perfectly mixed. 
	
	The time evolution of Φ_n is shown in Fig. <ref> (bottom plot). Initially, Φ_n grows until it reaches a maximum, Φ_max, at a time defined as t^*. For t>t^*, this parameter decreases until it saturates asymptotically. Note that Φ_n is not zero at t=0. Since the particles are initially distributed uniformly at random, the average absolute difference | N_f-N_s| for a given slice can be estimated from a binomial distribution. This difference follows a half-normal distribution with mean μ=√(2N_t/π), where N_t=N_f+N_s is the total number of particles in the slice. Therefore, the value of Φ_n at the starting configuration is Φ_n(0)=Φ_0=√(2/ρ_0L_x Δ y π), where ρ_0 is the initial total density of particles and it vanishes only in the thermodynamic limit. The initial increase in Φ_n corresponds to the dynamical demixing regime where there is a rapid accumulation of the fast particles in the deposit with a fraction of the slow particles dragged along while the remaining slow particles lag behind. The peak is reached when all the fast particles deposit at t=t^* (see Fig. <ref>(c)). For t>t^*, Φ_n decreases until the remaining slow particles deposit, on a top layer consisting (almost) exclusively of slow particles, if the difference in the velocities is sufficiently large (see Fig. <ref>(d)). We define this instant as the saturation time, t_final. Note that, the level of segregation of the particles remains almost the same from t^* to t_final, as seen from the snapshots (c) and (d). This parameter depends on the level of segregation and on L'_y(t) that decreases with time. The region where only slow particles are present occupies a larger area at t^* than at t_final and the same number of particles contribute to the integral. 
	
	To characterize the structure of the deposit we measured the six-fold bond order parameter, ⟨ϕ_6 ⟩, defined  as
	
	
    ⟨ϕ_6 ⟩ = 1/N∑_l^N1/6| ∑_j^N_be^i6θ_lj| ,

	
	where N is the total number of particles, N_b is the number of neighbors (within a cutoff of 2.5) and θ_lj is the angle between the line that connects the particles j and l with the x-axis. This is a continuous order parameter that is one when the particles are arranged in a perfect hexagonal structure and it decreases when the order decreases. At the low temperatures considered, the particles in the deposit form a nearly perfect hexagonal structure with ⟨ϕ_6 ⟩ = 0.9 and the level of segregation remains the same until the end of the simulation (Φ_n is constant). This does not correspond to the configuration expected at thermodynamic equilibrium where the particles form a completely mixed phase, as they are indistinguishable <cit.>. Obviously, the deposits observed at the end of the simulations are transient but their relaxation towards equilibrium occurs on much longer timescales than the observation time.
	
	The solid lines with open symbols in Fig. <ref> show the time dependence of Φ obtained numerically for different values of the velocity ratio v. As expected, Φ_max increases as the ratio v decreases from one, revealing that, as the difference of the sedimentation velocities increases, higher levels of segregation are attained in the deposit. 
	
	Reducing v, decreases the sedimentation velocity of the slow particles (keeping the velocity of the fast particles constant). However, the effective velocity of the fast particles also decreases due to the interactions with the slow ones (which act as obstacles) and, as a consequence, the time when the peak in Φ is reached, t^*, also increases. The saturation time, t_final, is also affected by v since the difference between t^* and t_final is the time taken by the remaining slow particles to deposit and this depends only on their sedimentation velocity. 
	
	The parameter Φ_n measures the segregation in the entire system but we can also measure the segregation in the deposit using the following definition <cit.>
    
	
    Φ_d= 1/N_d∑_i^N_d( N_f - N_s)^2/(N_f + N_s)^2,

	
	where N_d is the total number of particles with ϕ_6>0.9. N_f and N_s are the number of fast and slow particles within the cutoff distance 2.5 of particle i. The segregation in the deposit increases monotonically with time (see Fig. <ref>). At t^*, the slope of the curve changes signalling the phase where only slow particles arrive.
	
	
	
	
	

 §.§ Analytical model 

	
	
	
	We consider now a simple analytical model. We assume that the particles move with a constant sedimentation velocity v_i = m g/γ_i that depends only on the particle species, and we neglect particle-particle interactions and thermal fluctuations. We define the thickness of the packed deposit as l^*(t). The number of particles of a species in the region y<l^*(t) at a given time is the number of particles initially at y<l^*(t) plus those of that species that entered into that region. The latter can be estimated considering that the fraction of particles of species i that entered into the region is  v_i/v_f+v_s. Despite the fact that particle interactions are soft, we assume an upper bound for the density, given by the packing fraction of disks with diameter σ=1, i.e.,  ρ_max= 2√(3)/3. We can then estimate the density of particles by
	
	
    ρ_i(y<l^*,t)= ρ_0/2 + (ρ_max-ρ_0)v_i/v_f+v_s,

	
	where (ρ_max-ρ_0) is an estimate of the increase in the local density.
	
	The integral in Eq. <ref> for t<t^* can be replaced by the sum of three terms, corresponding to the contribution of three different regions, as shown in Fig. <ref>. Region I corresponds to the deposit, where the density of each species is given by Eq. <ref>. We also account for the non-zero segregation factor, Φ_0, which arises from the initial uniform distribution of the particles and the discrete nature of the numerical integration. Region II, where the two types of particles are perfectly mixed, is delimited by the surface of the deposit, y ≈ l^*(t) and  y ≈ L_y-v_ft, the height of the last fast particle. Here, we assume that the number of particles of either species that leave region II is approximately the same as the number that enters it and, therefore, the integrand is the constant Φ_0. Region III contains particles of only one type and is delimited by y ≈ L_y-v_st ≡ L'_y, the position of the last slow particle. In this region, the integrand is one as only slow particles are present. The number of fast, N_f, and slow particles, N_s,  in regions I and II can then be regarded as the number of successes and failures in a binomial distribution of N_t trials. The probability that a particle in each of these regions is a fast particle is
	
    p_f,α= ρ_0/2+(ρ_α-ρ_0)1/1+v/ρ_α,

	
	and the probability that it is a slow particle is
	
    p_s,α= ρ_0/2+(ρ_α-ρ_0)v/1+v/ρ_α,

	where ρ_α is the the total density of each region (ρ_max in region I and ρ_0 in region II). The term |  N_f - N_s| corresponds to the absolute difference between successes and failures of N_t trials that follow a folded normal distribution whose expected value is given by: 
	  
    μ'_α = σ_α√(2/π)exp(-μ_α^2/2 σ_α^2) + μ_α(μ_α/√(2 σ_α)).

	Here,  μ_α = N_t(p_f,α-p_s,α) and σ_α^2 =4 N_tp_f,α p_s,α + μ_α correspond to the expected value and to the variance of the Gaussian distribution (not “folded”) of the difference between successes and failures for the probabilities p_f,α and p_s,α. 
	The integral for t ≤ t^* is then approximated by
	
	
    Φ(t ≤ t^*) = 1/L_y'   [ ∫_0^l^*(t)μ'_I/ρ_maxΔ_y L_x dy +.
       . + ∫_l^*^L_y-v_ftμ'_II/ρ_0Δ_y L_x dy  +  ∫_L_y-v_ft^L_y' dy ].

	
	For t^*<t<t_final the thickness of region II is zero and then
	
	
    Φ(t > t^*) = 1/L_y'   [  ∫_0^L^*μ'_I/ρ_maxΔ_y L_x dy  +∫_L^*^L_y' dy ],

	
	where L^*=l^*(t^*) is the length of the structure at t=t^*.
	
	We estimate l^*(t) in the following way. The number of deposited particles is given by the flux of particles through the line defined by y=l^*(t) plus the number of particles that is initially present in the region below this height. We consider that particles travel with constant velocity, v_i, for y>l^*(t) and that the flux through that line for each species can be written as j_i(t)=ρ_0/2v_iL_xt. The number of particles in the deposit, N_s, is
	
	
    N_s=ρ_max L_x l^*(t) = L_x l^*(t)  ρ_0+ (ρ_0/2 v_f+ρ_0/2 v_s ) t L_x.

	
	We can now rearrange the terms and the expression for l^* is
	
	
    l^*(t) = ρ_0/2(1+ v)/ρ_max-ρ_0v_ft.

	For t>t^*,  all the fast particles are deposited and
	
	
    L^* = ρ_0L_y (1+v)/2(ρ_max-ρ_0+ ρ_0/2(1+v)).

	
	Since L^*=l^*(t^*),
	
	
    t^* = (ρ_max-ρ_0)L_y/v_f[1/2ρ_0(1+v)+(ρ_max-ρ_0)   ].

	
	The height of the final deposit is l^*(t_final)=N/ρ_maxL_x, where t_final is the time when the last slow particle deposits, which is estimated as t_final=(L_y-N/ρ_maxL_x)/v_s. At later times the parameter Φ remains constant as no more particles are added to the deposit.
	
	Thus, the parameter Φ is
	
    Φ(t)=
    	1/L_y-v_st[ ρ_0(1+v)/2 (ρ_max-ρ_0)(μ'_Iρ_0-μ'_IIρ_max/ρ_maxρ_0 Δ y L_x-
     -μ'_II)/ρ_0 Δ y L_x+1 - v) v_f t +
    	  L_y μ'_II/ρ_0 Δ t L_x],  t ≤ t^* 
    
    
    	
    	1+L^*/L_y-v_st[μ'_I/ρ_maxΔ y L_x -1],  t^* ≤ t ≤ t_final
    
    
    	1+L^*/L_y-v_st_final[μ'_I/ρ_maxΔ y L_x -1],    t ≥ t_final .

	
	
	The dashed lines in Fig. <ref> show the time dependence of Φ obtained from Eq. <ref> for different values of the ratio v. The comparison with the numerical results reveals that the analytical calculation reproduces the main features of the simulations. The highest levels of segregation are achieved for the lowest v, when the difference between the velocities is largest. In this limit, the slow particles travel much slower than the fast ones and are, on average, the last to deposit forming a thick layer on top of the first deposit consisting of slow particles only. For v=1, the particles are indistinguishable and the maximum value of Φ is the (finite-size) initial value. 
	
	
	
	
	
	So far, we discussed segregation along the y-direction. However, the snapshot in Fig. <ref>(d) suggests that at the bottom of the deposit there are linear-like clusters along the y-direction, which may promote segregation along the x-direction. To further investigate this question, we identified all the clusters of particles in the final deposit below the top layer of slow particles and calculated their inertia tensor. We defined the parameter r_I = 1/N_c∑ r_c, where r_c is the ratio of the smallest and largest (non-zero) eigenvalues of the inertia tensor of cluster c and the sum is over all clusters of the same species with N_c the number of such clusters. This parameter is one when all the clusters are symmetric and falls below one as the clusters become asymmetric. We measure r_I for clusters larger than two particles excluding the largest cluster, where finite size effects may be significant. A deposit with randomly distributed particles is obtained for v=1, when the particle are indistinguishable. In this limit, r_I is low, which suggests the prevalence of symmetric clusters. In the limit of low v, the number of fast particles is much larger than the number of slow ones in the region under consideration and they form a single large cluster. Accordingly, as shown in Fig. <ref> (top panel), r_I increases for fast particles as v decreases. The opposite occurs for the slow particles where r_I decreases, showing that, as the ratio of velocities decreases, the clusters of slow particles become asymmetric. In the bottom panel of Fig. <ref>, we plot the average angle, ⟨θ⟩, with the y-direction of the eigenvector corresponding to the largest eigenvalue of the inertia tensor for each cluster of slow particles. As v decreases, ⟨θ⟩ also decreases revealing that the clusters in the deposit tend to extend along the y-direction as the difference of the particles velocities increases. Accordingly, the distribution of the angles of the clusters with the y-direction (inset of Fig. 6) the number of clusters with small angles increases for low v. This is likely a consequence of the laning phenomenon observed in binary mixtures of species moving at different velocities in region II <cit.>.  
	
		
	
	Although we have neglected the effect of particle-particle collisions during sedimentation in the analytical model, the simulation and analytic results are in good agreement for intermediate and low densities. In the limit of high densities, however, the simulation and analytic results disagree, as particle-particle interactions becomes relevant. In Fig. <ref> we plot the time-dependence of the order parameter for v=0.5 for the system depicted in Fig. <ref> (in blue) with density ρ_0 = 0.4 and for a system with a density 60% higher, for two system sizes (red and green). In spite of the fact that at higher densities the time dependence of the order parameter obtained from Eq. <ref> is not in quantitative agreement with the numerical simulations, the maximum of the order parameter and its final value are similar.
	

	

 §.§ Hydrodynamic effects

	
	
We consider now the role of hydrodynamics in the sedimentation of particles with different velocities. As described in Sec. <ref>, we simulate the fluid flow and fluid-particle interaction using the LBM and the particle sedimentation is driven by an external force a⃗, with different values for the two species (fast and slow). 

    
	
	
    
	
Figure <ref> shows the initial configuration of one of the samples and its final state for three different force ratios β=a_s/a_f. From visual inspection, we notice that for β=0.1, where the fast particles have the highest velocity, particles in the final deposit are more segregated than for the other two values of β. This behaviour is in line with what we found without hydrodynamics. When the time evolution of the particles is analyzed, as shown in Fig. <ref>, some major differences are revealed, which result from hydrodynamic effects. For example, we observed Rayleigh-Taylor instabilities <cit.>  as shown in Fig. <ref>. Initially, in order for the particles to go down, the fluid below them has to go up generating upcurrents. Despite the sedimentation force acting on the particles, these upcurrents can move the particles up, by contrast to the results of the MD simulations and the analytical model. Thus, quantitative differences due to hydrodynamics are to be expected.

	
	
		
	

In Fig. <ref>, we show the profile of the number of particles as a function of the height for fast and slow particles for different values of β. 
This corresponds to the final states, shown in Fig. <ref>, but now averaged over 10 independent samples. The behaviour of the profiles is qualitatively the same as those shown in Fig. <ref>. Segregation is observed for the system with β=0.1, and no segregation occurs for β=1 (equal sedimentation velocities for both species). In Fig. <ref>, we compare the time evolution of the order parameter Φ (Eq.  (<ref>)) from the LBM simulations with the analytical model. To obtain the sedimentation velocity using the analytical model, we simulated the sedimentation of a single particle for different forces and measured the terminal velocity. In this case, the velocity increases linearly with the applied force, as expected for the values of the Reynolds numbers considered. Thus, the ratio of the forces is equal to the ratio of the velocities: β = v. The analytical curves using this velocity are represented by the solid lines in Fig. <ref>. We find that the sedimentation process for this velocity is much faster without hydrodynamics. This difference is due the upwards fluid displacement when the particles are sedimenting, which reduces their velocity (and even inverts the direction of motion). Thus, the average velocity of all particles is much lower due to the retarding effects of hydrodynamics. Numerically, we observe that the average velocity of particles of one specie (fast or slow) is nearly one fifth the velocity of a single particle of the same specie sedimenting in the fluid. 

We also show in Fig. <ref> the analytical results using this average velocity (dashed lines). Changing the velocity changes 
the time scale of the sedimentation, which becomes closer to the time scale of the LBM simulations. The final value of the order parameter is different from that predicted analytically as the velocity is highly non-uniform, which is not considered in the analytical model. However, the existence of a maximum of the order parameter is captured by the analytical model.

	
	

§ CONCLUSIONS 

	
	We investigated the dynamics of sedimentation of a simple binary mixture, and observed purely dynamical demixing. 
	
	We started by considering that the two species differ in their Stokes coefficients only, i.e., they have different sedimentation velocities. Since the species travel at different velocities, they demix dynamically as they move towards the bottom of the container in a gravitational field. We measured the degree of demixing without hydrodynamics using Brownian dynamics simulations for different ratios of the velocities and proposed a simple analytical description in the low density limit. We found that the analytical model captures the dynamics and the degree of segregation in the system even though particle-particle interactions are not taken into account. We also considered the sedimenting of particles in a hydrodynamic setting using the lattice Boltzmann method and found similar qualitative behavior, in particular, the existence of purely dynamical demixing. Quantitative differences, however, were observed resulting from the currents set up during the sedimentation which strongly affect the particles velocities, making them non-uniform.
	
	We focused on equimolar binary mixtures but the same demixing mechanism will occur for mixtures with any other composition. In fact, the composition of the deposit is determined by the ratio of the fluxes of sedimenting of the two species and, as a result, the initial composition of the mixture will affect the composition of the final deposit. The mechanism described here could also be used to obtain fully mixed deposits, since by tuning v an equimolar deposit may be formed from mixtures poor in one of the two components.
	
	As a final note, we focused on colloidal suspensions but our conclusions can be extended to systems with larger particles, as we considered the limit of high Péclet number, where the dominant mechanism of mass transport is advection and thermal fluctuations are negligible. Finally, since in this limit the relevant mechanisms depend only on the ratio of sedimentation velocities, we expect the same behavior for particles with the same shape but different buoyant masses. The nature of the field is also irrelevant, and thus similar results are to be expected for other (constant) external fields (e.g., electromagnetic field).
	
	

§ ACKNOWLEDGEMENTS

	
	We acknowledge financial support from the Portuguese Foundation for Science and Technology
	(FCT) under the contracts no. EXCL/FIS-NAN/0083/2012, SFRH/BD/119240/2016, UIDB/00618/2020 and UIDP/00618/2020.

	
	
	
