Sedimentation of particles in Stokes flow

In this paper, we consider $N$ identical spherical particles sedimenting in a uniform gravitational field. Particle rotation is included in the model while inertia is neglected. Using the method of reflections, we extend the investigation of [R. M. H\"ofer, Sedimentation of inertialess particles in Stokes flows, arXiv:1610.03748, (2016)] by discussing the optimal particle distance which is conserved in finite time. We also prove that the particles interact with a singular interaction force given by the Oseen tensor and justify the mean field approximation of Vlasov-Stokes equations in the spirit of [M. Hauray and P. E. Jabin, Particle approximation of Vlasov equations with singular forces : propagation of chaos, Ann. Sci. Ec. Norm. Super. (4), (2015)] and [M. Hauray, Wasserstein distances for vortices approximation of Euler-type equations, Math. Models Methods Appl. Sci. 19, (2009), pp. [1357,1384]].

1. Introduction. In this paper, we consider a system of N spherical particles (B i ) 1≤i≤N with identical radii R immersed in a viscous fluid satisfying the following Stokes equation: completed with the no-slip boundary conditions where (V i , Ω i ) ∈ R 3 × R 3 , 1 ≤ i ≤ N represent the linear and angular velocities, We describe the intertialess motion of the rigid spheres (B i ) 1≤i≤N by adding to the instantaneous Stokes equation the classical Newton dynamics for the particles where m denotes the mass of the identical particles adjusted for buoyancy, g the gravitational acceleration, F i (resp. T i ) the drag force (resp. the torque) applied 996 AMINA MECHERBET by the the fluid on the i th particle B i defined by with n the unit outer normal to ∂B i and σ(u N , p N ) = 2D(u N ) − p N I, the stress tensor where 2D(u N ) = ∇u N + ∇u N . Note that the constant velocities (V i , Ω i ) of each particle are unknown and are determined by the prescribed force and torque F i = mg and T i = 0. In [16], the author shows that the linear mapping on R 6N is bijective for all N ∈ N * . This ensures existence and uniqueness of (u N , p N ) and the velocities.
Remark 1 (About the modeling and nondimensionalization). Equations (1)- (3) describe suspensions sedimenting in a uniform gravitational field. Equations (1), (2) are derived starting from the Navier-Stokes equations and neglecting the inertial terms by means of the Reynolds and Stokes number, see [6, Chapter 1 Section 1], [1], [16] and all the references therein. Analogously, the ODE system (3) is obtained by neglecting particle inertia. We refer also to [4] where a formal derivation taking into account the slow motion of the system is performed. When considering one spherical particle sedimenting in a Stokes flow, the linear relation between the drag force F and the velocity V is given by the Stokes law F = −6πRV , see Section 2.1 for more details. Stokes law leads to the well-known formula for the fall speed of a sedimenting single particle under gravitational force denoted by κg := m 6πR g .
It is important to point out that in our model, a scaling with respect to the velocity fall κg has been performed. This means that the drag forces (F i ) 1≤i≤N and the gravitational force mg are terms of order R. Consequently, in this paper, κg is a constant of order one. For more details on the derivation of the model, we refer to [11,Section 1.1] where a nondimensionalization including physical units is provided. Moreover, as in [4], the particle radius R is assumed to be proportional to 1 N so that the collective force applied by the particles on the fluid is of order one. This will be made precise in the presentation of the main assumptions.
Given initial particle positions x i (0) := x 0 i , 1 ≤ i ≤ N , we are interested in the asymptotics of the solution when the number of particles N tends to infinity and the radius R tends to zero. The main motivation is to justifiy the representation of the motion of a dispersed phase inside a fluid using Vlasov-Stokes equations in spray theory [7], [2].
The analysis of the dynamics is done in [13] in the dilute case i.e. when the minimal distance between particles is at least of order N −1/3 . The authors prove that the particles do not get closer in finite time. Moreover, in the case where the minimal distance between particles is much larger than N −1/3 the result in [13] shows that particles do not interact and sink like single particles. We refer finally to [11] where the author considers a particle system with minimal distance of order N −1/3 and proves that, under a relevant time scale, the spatial density of the cloud converges in a certain averaged sense to the solution of a coupled transport-Stokes equation (15).
Since the desired threshold for the minimal distance is of order N −2/3 , which allows to tackle randomly distributed particles, we are interested in extending the results for lower orders of the minimal distance. Therefore, in this paper, we continue the investigation of [11] by looking for a more general set of particle configurations that is conserved in time and prove the convergence to the kinetic equation (15). Also, we include particle rotation in the modeling.
1.1. Main assumptions and results. In this Section, we describe the configuration of particles that we consider and present the main results : Theorem 1.2 and Theorem 1. 3.
We recall that the particles B i are spherical with identical radii R with r 0 a positive constant satisfying a smallness assumption (see Theorem 1.2). Due to the quasi-static modeling, the velocities (V i (t), Ω i (t)) 1≤i≤N at time t ≥ 0 depend only on the prescribed force (F i ) 1≤i≤N , torque (T i ) 1≤i≤N and the particles position (x i (t)) 1≤i≤N at the same time t. Consequently, we drop the dependence with respect to time in the definition of the set of particle configurations. Keeping in mind that the idea is to start from a configuration of particles that lies in the set and show that it remains in it for a finite time interval. Definition 1.1 (Definition of the set of particle configuration). Let (X N ) N ∈N * be a configuration of particles, where X N := (x 1 , · · · , x N ). We define the minimal distance d N min by d N min := min i =j 1≤i,j≤N We introduce the particle concentration M N defined for each positive sequence (λ N ) N ∈N * by Given two positive constantsM , E and a sequence (λ N ) N ∈N * , we define X (M , E) as the set of configurations for which (d N min ) N ∈N * and (M N ) N ∈N * satisfy the following assumptions: λ N must satisfy the following compatibility conditions: 998 AMINA MECHERBET Remark 2. Note that, according to the definition of M N , assumption (5) ensures which yields thanks to assumption (6) Since R ∼ 1 N , this leads also lim which ensures that the particles do not overlap.
Furthermore, for the proof of the second Theorem 1.3, the following assumption must be satisfied initially: Finally, we define ρ N the spatial density of the cloud by In the rest of this paper, if needed, we make clear the dependence with respect to time by writing for all N ∈ N * , X N (t) = (x 1 (t), · · · , x N (t)) for the particles configuration, d N min (t) for the minimal distance and M N (t) the particle concentration at time t ≥ 0.
The main results of this paper are the two following theorems. The first one ensures that the particle configurations considered herein are preserved in a short time interval depending only on the data r 0 ,M , E, κ|g|. Theorem 1.2. Let (X N (0)) N ∈N * be the initial position of the particles. Assume that there existsM , E and a sequence (λ N ) N ∈N * such that (X N (0)) N ∈N * lies in the set X (M , E) i.e. assumptions (5), (6), (7) hold true initially. IfM 1/3 r 0 is small enough, there exists N * ∈ N * depending on (r 0 ,M , E) and T > 0 depending on (r 0 , E,M , κ|g|) such that for all t ∈ [0, T ] and N ≥ N * The second part of the result is the justification of the convergence of ρ N when N tends to infinity. Theorem 1.3. Consider the maximal time T > 0 introduced in Theorem 1.2 and the additional assumption (11). Let ρ 0 be a positive regular density such that R 3 ρ 0 = 1. We denote by (ρ, u) the unique solution to the coupled equation (15).
There exists some positive constants C 1 , C 2 depending on (r 0 ,M , E, ρ 0 L ∞ , κ|g|) and N * ∈ N * depending on (r 0 ,M , E, ρ 0 L ∞ , κ|g|, T ) such that for all N ≥ N * and t ∈ [0, T ] This shows that if the initial particle distribution ρ N 0 converges to ρ 0 then the particle distribution ρ N converges toward the unique solution ρ of equation (15) for all time 0 ≤ t ≤ T . Moreover, Theorem 1.3 provides a quantitative convergence rate in terms of the initial Wasserstein distance W 1 (ρ 0 , ρ N 0 ). Remark 3. The regularity assumption on the initial density ρ 0 is the one introduced by Höfer in [11] which is ρ 0 , ∇ρ 0 ∈ X β , for some β > 2. See Section 5.1 for the definition of X β . In particular, the assumption is satisfied if ρ 0 is compactly supported and C 1 .
The idea of proof of Theorem 1.3 is to formulate the problem considered as a mean-field problem. The mean-field theory consists in approaching equations of motion of large particles systems (X 1 , · · · , X N ) when the number of particles N tends to infinity. In mean-field theory, the ODE governing the particle motion is known and is given by where the kernel F is the interaction force of the particles. The limit model describing the time evolution for the spatial density ρ(t, x) is given by In our case, the first difficulty is to extract a system similar to (12) for the particle motion and to identify the interaction force F . A key step is then a sharp expansion of the velocities for large N . We obtain for each 1 ≤ i ≤ N where Φ is the Green's function for the Stokes equations, also called the Oseen tensor (see formula (26) for a definition). κg is the fall speed of a sedimenting single particle under gravitational force and is of order one in our model, see Remark 1. This shows that the particle system satisfies approximately equation (12) with an interaction force given by the Oseen tensor. Since the convolution term Kρ appearing in (13) corresponds to the solution of a Stokes equation in our case, the limiting model describing (1), (2), (3) is a coupled transport-Stokes equation The proof of Theorem 1.3 is based on the two papers [9], [8] where, in the first one, the authors justify the mean field approximation and prove the propagation of chaos for a system of particles interacting with a singular interaction force and where the ODE governing the particle motion is second order. In [8] the author considers a different mean-field equation where the particle dynamics is a first order ODE. The results obtained hold true for a family of singular kernels and applies to the case of vortex system converging towards equations similar to the 2D Euler equation in vorticity formulation. The associated kernel in this case is the Biot-Savard kernel.
In order to extract the first order terms for the velocities (V i , Ω i ) we apply the method of reflections. This method is introduced by Smoluchowski [17] in 1911. The main idea is to express the solution u N of N separated particles as superposition of fields produced by the isolated N particle solutions. We refer to [14,Chapter 8] and [6,Section 4] for an introduction to the method. A convergence proof based on orthogonal projection operators is introduced by Luke [16] in 1989. We refer also to the method of reflections developped in [12] which is used by Höfer in [11]. In this paper, we design a modified method of reflections that takes into account the particle rotation and relies on explicit solutions of Stokes flow generated by a translating, rotating and straining sphere. To obtain the convergence of the method of reflections we need to identify particle configuration that can be propagated in time. The particle configuration considered herein is the one introduced in [10] to study the homogenization of the Stokes problem in perforated domains. The novelty is that the author considers the minimal distance d N min together with the particle concentration M N as parameters to describe the cloud. The result in [10] extends in particular the validity of the homogenization problem for randomly distributed particles i.e. particle configurations having a minimal distance of order at least N −2/3 . Note that the notion of particle concentration appears also in [9] to describe the cloud.
1.2. Discussion about the particle configuration set. As stated above, the assumptions introduced in Definition 1.1 are based on [10]. Assumptions (5) and (7) means that there exists a uniformly bounded discrete spatial density that approximates ρ N . Indeed, if we defineρ N bỹ one can show that W 1 (ρ N , ρ N ) ≤ λ N . Assumption (5) ensures that there exists a sequence λ N for which the infinite norm ofρ N is bounded byM , see formula (95). This suggests that ρ ∞ andM are equivalent.
We recover the result of [13] in the case where λ N = N −1/3 and the minimal distance d N min is much larger than N −1/3 , the explicit formula for the velocities (14) implies in this case which is in accordance with the "non-interacting scenario" explained in [13]. In our case, the smallness assumption on r 0M 1/3 means that we consider a density of particles such that ρ ∞ is small but of order one. Indeed, the second term in the velocity formula (14) can be seen as a perturbation of order one of the velocity fall κg in the case whereM (or the particle density ρ ∞ ) is small. This can be also seen in the coupled equation (15) where the velocity term u is proportional to ρ ∞ .
The second assumption (6) ensures the conservation of the minimal distance, see Proposition 7. In particular, for λ N = N −1/3 , Theorem 1.2 extends the previous known results to configurations having minimal distance at least of order N −1/2 , see assumption (6). This lower bound for the minimal distance appears naturally in our analysis and is closely related to the properties of the Green's function for the Stokes equations. We emphasize that this critical minimal distance appears also in the mean-field analysis due to [8]. Precisely, computations in the proof of [8,Theorem 2.1] show the convergence for a short time interval under the assumption that Definition 5.4 for the definition of the infinite Wasserstein distance W ∞ . Standard measure-theory arguments show that the infinite Wasserstein distance ensures assumption (5). In other words, one can take λ N to be the infinite Wasserstein distance, which yields finally the same assumption (6).
The first assumption in formula (7) means that we are interested in cases where there is more than one particle per cube of length λ N . As pointed out by Hillairet in [10], one can choose a larger sequence (λ N ) N ∈N * such that the compatibility assumption holds true. Note also that, in the case where λ N is the infinite Wasserstein distance, this compatibility assumption is satisfied by definition. Finally, assumption (11) is needed for the control of the Wasserstein distance.
1.3. Outline of the paper and main notations. The remaining Sections of this paper are organized as follows.
In Section 2 we recall the classical results for the existence and uniqueness of the Stokes solution u N . We recall also the definition of the drag force F i , torque T i and stresslet S i and present in Section 2.1 the particular solutions to a Stokes flow generated by a translating, a rotating or a straining sphere. Finally, the end of Section 2 is devoted to the approximation of the stresslets S i . In Section 3 we present and prove the convergence of the method of reflections in order to compute the first order terms for the velocities (V i , Ω i ) 1≤i≤N . Section 4 is devoted to the proof of Theorem 1.2. In Section 5 we recall some definitions associated to the Wasserstein distance. We present then the strong existence, uniqueness and stability theory for equation (13). In the second part of Section 5 we show that the discrete density ρ N satisfies weakly a transport equation. Section 6 is devoted to the proof of the second Theorem 1.3. Finally, some technical Lemmas are presented in the appendix.
Notation. In this paper, n always refer to the unit outer normal to a surface.
The following shortcut will be often used where we drop the dependence with respect to N in order to simplify the notation. Given an exterior domain Ω with smooth boundaries, we set and the following norm for all u ∈ C ∞ (Ω) we define then the homogeneous Sobolev space D(Ω) as the closure of C ∞ (Ω) for the norm · 1,2 (see [5,Theorem II.7.2]). We also use the notation D σ (Ω) for the subset of divergence-free D(Ω) fields Which is also the closure of the subset of divergence-free C ∞ (Ω) fields for the · 1,2 norm. Analogously, if Ω = R 3 we use the notatioṅ We denote by × the cross product on R 3 and by ⊗ the tensor product on R 3 which associates to each couple (u, v) ∈ R 3 × R 3 the 3 × 3 matrix defined as For all 3 × 3 matrices A, B, we use the classical notation In R 3 , | · | stands for the Euclidean norm while | · | ∞ represents the l ∞ norm. We use the notation B ∞ (x, r) for the ball with center x and radius r for the l ∞ norm. Finally, in the whole paper we use the symbol to express an inequality with a multiplicative constant independent of N and depending only on r 0 ,M , E and eventually on κ|g| which is uniformly bounded according to Remark 1.

2.
Reminder on the Stokes problem. In this Section we recall some results concerning the Stokes equations. We remind that for all N ∈ N we denote by (u N , p N ) the solution to (1) - (2). Keeping in mind that the linear mapping, that associates to the linear and angular velocities the forces and torques, is bijective (see [16]) together with the classical theory for the Stokes equations yields: The velocity field u N can be extended to V i + Ω i × (x − x i ) on each particle B i . This extension denoted also u N is inḢ 1 σ (R 3 ). We recall the definition of the force F i ∈ R 3 , torque T i ∈ R 3 and stresslet S i ∈ M 3 (R) applied by the particle B i on the fluid (see [6,Section 1.3]) The matrix M i represents the first momentum which is decomposed into a symmetric and skew-symmetric part M i = T i + S i , the symmetric part S i is called stresslet, see [6,Section 2.2.3]. Since the skewsymmetric part of a 3 × 3 matrix M has only three independent components, it can be associated to a unique vector T such that In this paper, we allow the confusion between the skew-symmetric matrix ssym(M ) and the vector T . Hence, we define the torque T i ∈ R 3 as being the skew-symmetric part of the first momentum M i which satisfies 2.1. Particular Stokes solutions. The linearity of the Stokes problem allows us to develop powerful tools that will be used in the method of reflections. In particular, we investigate in what follows the analytical solution to a Stokes flow generated by a translating, a rotating or a straining sphere. The motivation in considering these cases is that the fluid motion near a point x 0 may be approximated by hence, if we replace the boundary condition on each particle by its Taylor series of order one, we can use these special solutions to approximate the flow u.
completed by the boundary condition On the other hand, the force F , torque T and stresslet S exerted by a translating sphere B as defined in (18) read We recall now an important formula that links the solution to the Green's function of the Stokes problem. For all x ∈ R 3 \ B(a, R) we have where Φ is the Green's function for Stokes flow also called Oseen-tensor Φ(x) = 1 8π The 3 × 3 matrix ∆Φ represents the Laplacian of Φ and is given by The first term in the right-hand side of (25) is the point force solution also called stokeslet, see [6,Section 3.1]. In this paper we use the term stokeslet to define the whole solution U a,r [V ] which can bee seen as an extension of the point force solution.
Remark 4. Formula (25) is closely related to the Faxén law which represents the relations between the force F and the velocity V . We refer to [6, Section 2.3] and [14, Section 3.5] for more details on the topic.
Remark also that in (25) the point force solution retains the most slowly decaying portion, which is of order R |x| . This property is useful in order to extract the first order terms for the velocities (V i ) 1≤i≤N , see Lemma 3.3.
Moreover, we recall a Lipschitz-like inequality satisfied by the Oseen tensor , ∀ x , y = 0.
Finally, in this paper, the velocity field U a,R [V ] is extended by V on B(a, R).
completed with the boundary conditions a,R [ω] represents the flow generated by a sphere rotating with angular velocity ω. In particular we have P On the other hand, the hydrodynamic torque resulting from the fluid traction on the surface defined in (19) is given by Finally, there exists C > 0 such that for all x ∈ R 3 \ B(a, R) completed with the boundary conditions The velocity field A (2) a,R [E] is the flow generated by a sphere submitted to the strain E(x − a). In this case, the drag force and torque vanishes On the other hand, the symmetric part of the first momentum S as defined in (19) is given by completed by the boundary conditions We set then D = E + ω with E = sym(D) and ω = ssym(D). As stated in the definition (19), ω represents also a 3D vector. Hence, the boundary condition (37) reads We have, thanks to the linearity of the Stokes equation, that Since the two solutions have the same decay-rate, this yields for all x ∈ R 3 \ B(a, R) Analogously, for the second derivative we have 2.2. Approximation result. In this part we consider the unique solution (v, p) of the following Stokes problem: completed with the boundary conditions with V ∈ R 3 and D a trace-free 3 × 3 matrix. We set We aim to show that the velocity field v 1 is a good approximation of the unique solution v.
Lemma 2.1. For N sufficiently large, we have the following error bound: as v and v 1 satisfy the same boundary condition on ∂B 1 this yields . In order to bound the first term we construct an extensionṽ of the boundary conditions of v and apply the variational principle. We definẽ where χ is a truncation function such that χ = 1 on B(0, 1) and χ = 0 out of B(0, 2). Thanks to formula (10), for N sufficiently large we have R < d N min /4 and thus suppṽ ⊂ B(x 1 , d N min /2).f is defined as follows . We refer to [5, Theorem III.3.1] for a complete definition of the Bogovskii operator. In particular, from [10, where With this constructionṽ is a divergencefree field satisfying the same boundary conditions as v. Moreover, applying formula (43), there exists (another) constant C > 0 such that Thanks to (22) and (38) we have: Reproducing an analogous computation for the first term we obtain finally This yields the expected result.
Proof. We fix i = 1. Let E be a trace-free symmetric 3 × 3 matrix. We define v as the unique solution to the following Stokes equation completed with the boundary conditions We also denote by (v 1 , p 1 ) the special solution (A We have thanks to the symmetry of E Using an integration by parts we have for the second term in the right hand side since v 1 corresponds to a flow submitted only to a strain, see (33). For the first term in the right hand side of (47), using Lemma 2.1 we have It remains to estimate ∇u N L 2 (R 3 \ i Bi) . One can reproduce the same arguments as for the proof of Lemma 2.1 or follow the same proof as [10, Lemma 10] to get Gathering all the estimates we obtain this being true for all symmetric trace-free 3 × 3 matrix E, we obtain the desired result.
3. Analysis of the stationary Stokes equation. This Section is devoted to the analysis of a method of reflections and computation of the unknown velocities (V i , Ω i ) 1≤i≤N . We remind that, for fixed time, u N is the unique solution to the stationary Stokes problem completed with the no-slip boundary conditions In this Section, we show that at each fixed time t ≥ 0, the convergence of the method of reflections toward the unique solution u N holds true in the case where (X N (t)) N ∈N * ∈ X (M , E) and under the assumption that r 0M 1/3 is small enough.
3.1. The method of reflections. In this part, we present and prove the convergence of a modified method of reflections for the velocity field u N for arbitrary N ∈ N * , we remind that u N is the unique solution to the stationary Stokes problem (1), (2), with unique velocities (V i , Ω i ) satisfying (48). The main idea is to express u N as the superposition of N fields produced by the isolated N particle. Thanks to the superposition principle, we know that the velocity field which represents the error committed on the boundary conditions when approaching u N by the sum of the particular Stokes solutions. In this paper, for all u * ∈ 1010 AMINA MECHERBET C ∞ ( i B i ) we use the notation U [u * ] to define the unique solution of the Stokes completed by the boundary conditions hence, we can write Note that the boundary condition u (1) * is not constant on each particle B i , thus, the idea is to approach u on each particle B i and write U [u (1) * ] as follows: We set then U [u (2) * ] the new error term satisfying We iterate then the process by setting for all and for p ≥ 1, for the error term we set and define for all With this construction the following equality holds true for all k ≥ 1 Remark 5. This method of reflection is obtained by expanding the error term u * up to the first-order which leads us to formula (56). If one consider an expansion of u * up to the zerothorder then one obtain only a stokeslet development: The main difference between these two expansions is that the first one allows us to tackle the particle rotation. It also helps us to obtain a converging method of reflections for a more general assumption on the minimal distance.
We emphasize that we only need to show that the series for all 1 ≤ i ≤ N converge to obtain the convergence of the expansion (56), see Lemma 3.1 and Proposition 3. Precisely, the only assumptions needed to obtain the convergence of the series are the smallness ofM 1/3 r 0 , assumption (5) and the fact that The second step is to show that the expansion is a good approximation of the unique solution u N . This is ensured by Proposition 5. Precisely, in addition of the previous assumptions, we need the following uniform bound One can show that this assumption is less restrictive than (6) and allows us to consider smaller minimal distance. To reach lower bound for the minimal distance, one may develop u * at higher orders.

Preliminary estimates.
Recall that the dependence in time is implicit in this Section. All the following estimates hold true under the assumption that there exists a sequence (λ N ) N ∈N * and two positive constantsM , E such that (X N ) N ∈N * ∈ X (M , E), see Definition 1.1 andM 1/3 r 0 is small enough.
Proof. Using formulas (53) and (55) we get and This yields, for all 1 ≤ i ≤ N , using the decay-rate of the special solutions (38), (22) and Lemma A.1 with k = 1 and k = 2 . For the second term on the right hand side we have which vanishes when N tends to infinity according to (6) and (7). Moreover, if r 0M 1/3 is small enough, this ensures the existence of a positive constant K < 1/2 such that for N large enough and depending on r 0 ,M and E.
We have the following estimate.
Proof. Considering only the stokeslet expansion we have The first term in the right hand side of (59) can be computed using the fact that

AMINA MECHERBET
For the second term in the right hand side of (59) we write for all i = j According to the decay properties of the stokeslet (22) we have where δ ij is the Kronecker symbol. On the other hand, we recall that the triangle inequality Using (60), formula (61) twice and Lemma A.1 we obtain for all i = j N l=1 where we kept only the largest terms using the fact that d N min vanishes according to (7) for N large enough. Hence, the second term in the right hand side of (59) yields using Lemma A.1 where we used the fact that R d N min 1 thanks to (10) and 1 according to (6) and (7). The term involving rotating and straining solutions A xi,R [D i ] is handled analogously.
Since the series are convergent, we denote their limit by Thanks to the linearity of the Stokes equation and Proposition 3, the expansion term uniformly in N to the expansion where we replace the series by their limit. This shows that the error term U [u (k) * ] has a limit when k → ∞. In order to quantify the error term, we begin by the following estimate Proposition 4. For all k ≥ 1 we set Under the same assumptions as Lemma 3.1, there exists N (r 0 ,M , E) ∈ N * such that for all N ≥ N (r 0 ,M , E) and Let x ∈ B i , using formula (10) we recall that for i = j Applying this, formula (10), the decay properties of the second gradient of single particle solutions (23), (39) and the iteration formula (55) together with Lemma A.1 for k = 3 yields hence, we iterate the formula and use the fact that ∇ 2 u (0) * = 0 according to formula (54) to get which yields the expected result by applying Lemma 3.1.

Estimate of ∇u
again, according to (10), note that for N large enough, 1 + R d N min ≤ 2. We conclude using assumption (6) to bound the right hand side by η (k) up to constants depending onM , E, r 0 .

Approximation result.
We can now state the main result of this Section. Proposition 5. Assume that there existsM , E and a sequence (λ N ) N ∈N * such that (X N ) N ∈N * ∈ X (M , E). Assume moreover thatM 1/3 r 0 is small enough.
There exists a positive constant C = C(r 0 ,M , E) and N (r 0 ,M , E) ∈ N * satisfying for all N ≥ N (r 0 ,M , E) Proof. The aim is to estimate ∇U [u , with χ a truncation function such that χ = 1 on [0, 1] and χ = 0 outside [0, 2]. We have then In what follows we introduce the notation A(x, r, R) := B(x, R) \ B(x, r) for r < R.
For the second term we set: i ), where B is the Bogovskii operator, see [10, Appendix A Lemma 15 and 16] for more details.
The construction satisfies: and thanks to the variational formulation we have where we used the fact that the v i have disjoint support.
Thanks to the properties of the Bogovskii operator B xi , R , 2R we get Thanks to Proposition 4 we have with K < 1 according to Lemma 3.1, we get Since K < 1 for N large enough the second term on the right hand side, which is uniformly bounded with respect to N , vanishes when k → ∞. This yields The second term on the right hand side can be bounded using assumptions (6), (8) and (10) which is the desired result.
Remark 6. According to Proposition 4 we have for all as for the proof of Proposition 5 the second term vanishes when k → ∞ and we obtain lim 3.1.3. Some associated estimates. We recall that we aim to compute the velocities (V i , Ω i ) associated to the unique solution u N of the Stokes equation: completed with the no-slip boundary conditions The method of reflections obtained in this Section helps us to describe the velocity field u N in terms of explicit flows In order to extract a formula for the unknown velocities (V i , Ω i ), 1 ≤ i ≤ N we need to compute first the velocities V where κg is defined thanks to formula (4).
Proof. For the sake of clarity we fix i = 1 and the same result holds for all 1 ≤ i ≤ N . Let V ∈ R 3 , D a trace-free 3 × 3 matrix. 1 and v = 0 on the other ∂B j , j = 1. We choose v the unique solution to the Stokes equation:

The main idea is to apply an integration by parts with a suitable test function
completed by the boundary conditions We extend u N and v by their boundary values on all B i , 1 ≤ i ≤ N . We set E = sym(D), Ω = ssym(D). An integration by parts yields see (18) and (19) for the definition of the force F 1 , torque T 1 and stresslet S 1 . On the other hand, we apply the method of reflections to get For the first term we integrate by parts to get for all 1 Recall that v vanishes on ∂B i , i = 1 and hence, the sums above are reduced to the first term. Applying (34) (30) and (24) there holds for all where δ 1j is the Kronecker symbol. For the second term on the right hand side of formula (67), we consider v 1 : To bound the last term we apply Lemma 2.1 and Proposition 5 We focus now on the first term on the right hand side of formula (68), we have using the decay properties (22), (38) we have According to Remark 6 , we have for all Finally we get Identifying formula (66) and (67) and gathering all the inequalities above we have for all V , Ω ∈ R 3 and symmetric trace- with F 1 + mg = 0, T 1 = 0. Note that the value of the stresslet S i , see (19) for the definition, is unknown. However, we only need to approximate its value using Proposition 2. We conclude by identifying the terms involving V ∈ R 3 to obtain for the skew-symmetric part we get , and for the symmetric part using Proposition 2 Hence, according to Lemma 3.1 we have K 1−K < 1. Moreover, assumption (9) ensures that R |d N min | 3/2 which vanishes when N goes to infinity.

3.2.
Extraction of the first order terms for the velocities (V i , Ω i ). In order to control the motion of the particles, we want to provide a good approximation of the unknown velocities (V i , Ω i ). Thanks to the method of reflections, the velocity field u N can be approached by a superposition of analytical solutions to a Stokes flow generated by a translating, a rotating and a straining sphere (See Proposition 5) with the associated velocities (V ∞ i , ∇ ∞ i ). This allows us to compute the first order terms for (V i , Ω i ) applying Lemma 3.2 and Corollary 1. Keeping in mind that all the computations are done for a fixed time t ≥ 0, the main result of this Section is the following Proposition.
Proposition 6. Assume that, for a fixed time, we have the existence of a sequence (λ N ) N ∈N * and two positive constantsM , E such that (X N ) N ∈N * ∈ X (M , E). Assume moreover thatM 1/3 r 0 is small enough. Then, there exists N (r 0 ,M , E) ∈ N * such that for all N ≥ N (r 0 ,M , E), for all 1 ≤ i ≤ N we have We begin by the following lemma: Proof of Lemma 3.3. Thanks to formula (25) we have for i = j We have thanks to assumptions (6), (8) and (10) Analogously, we obtain the second bound by applying A.1 with k = 2 this time.
We can now prove the main result.
Proof of Proposition 6. Let fix 1 ≤ i ≤ N . According to Lemma 3.2 and Corollary 1 we have we apply Lemma 3.3, Lemma 3.2 and Corollary 1 together with (9) and (10) to get: min . Now, we rewrite the sum as follows: and we bound the error term using the decay rate (22), Lemma 3.2 and Lemma A.
We conclude by replacing the stokeslets by the Oseen tensor thanks to Lemma 3.3. Finally we have for all 1 ≤ i ≤ N For the angular velocities we obtain thanks to Lemma 3.2 and formula (58) for ∇ As before, using Lemma 3.2 we bound the first term by where we used the fact that R|λ N | 3 |d N min | 3 is uniformly bounded according to (6) and (10). 4. Control of the particle distance and concentration. In this Section, we make precise the particle behavior in time. Precisely we want to prove that if initially there exists two positive constantsM , E and a sequence (λ N ) N ∈N * such that (X N (0)) N ∈N * ∈ X (M , E) (see Definition 1.1), then the same holds true for a finite time. Recall that the initial distribution of particles satisfies: • The minimal distance is at least of order |λ N | 3/2 . • The maximal number of particles concentrated in a cube of width λ N satisfies assumption (5).
We aim to show that there exists a small interval of time [0, T ] independent of N on which the particle distance and concentration stay at the same order. The idea is to use a Gronwall argument and the computation of the velocities (V i ) 1≤i≤N at each fixed time t ≥ 0.

4.1.
Proof of Theorem 1.2. We assume that initially there exists two positive constantsM , E and a sequence (λ N ) N ∈N * such that (X N (0)) N ∈N * ∈ X (M , E). Let T > 0 be such that This maximal time T > 0 exists and we aim to prove that it is independent of N . As long as t < T we have a control on the particle concentration.
Proof. We recall the definition of M N (t) We introduce the following quantity: One can show that the two definitions of concentration L N (t) and M N (t) are equivalent in the sense that see Lemma A.2 for the proof. We also need to introduce the following notation for all β > 0: The idea is to show that the concentration L N is controlled in time and hence, the same applies to M N according to Lemma A.2. Recall that we have for all Which means that We obtain Hence taking the maximum over 1 ≤ i ≤ N we obtain L N This shows that as long as t < T we have (X N (t)) N ∈N * ∈ X (8 4M , 4E). This implies the following control.
Proof. For the sake of clarity we fix i = 1 and j = 2. The computations below are independent of this choice. Thanks to Proposition 6 we obtain : Hence, according to assumption (6), formula (27) and using Lemma A.1 for k = 2 we obtain We set then C > 0 the universal constant implicit in the above estimate.
We have the following control.
Proof. Thanks to (70) and Lemma 4.1 we have for all t < T that Hence, all computations from Proposition 7 hold true up to time T . In other words, there exists a positive constant C = C(r 0 ,M , E, κ|g|) such that for all indices This entails which is the desired result.

Conclusion. Thanks to Lemma 4.2 and Lemma 4.1 we have for all 1
this shows that T is independent of N and is at least of order log (2) C where C depends on (r 0 ,M , E, κ|g|).

5.
Reminder on Wasserstein distance and analysis of the limiting equation. In this part we recall some important results of existence, uniqueness, regularity and stability concerning the mean-field equation (13). We recall also the definition of the Monge-Kantorovich-Wasserstein distance of order one and infinite. We refer to [18,Part I,chapter 6] for definition and properties of the order one distance W 1 . To define the infinite Wasserstein distance we start with some associated notions. We refer to [3] for more details.
Recall that for all probability measure λ ∈ P(R 3 × R 3 ) we have We recall also the definition of the support for a (non-negative) measure. With this definition for the support one can show that there holds λ − esssup |x − y| := sup{|x − y| : (x, y) ∈ supp λ}).
We can now define the infinite Wasserstein distance W ∞ : Definition 5.4 (Infinite Wasserstein distance). The infinite Wasserstein distance between two probability measures µ and ν is defined as follows: A transference plan π * ∈ Π(µ, ν) satisfying is called an optimal transference plan.
We recall also the definition of a transport map.
Definition 5.5 (Transport map). Given two probability measures µ and ν, a transport map T is a measurable mapping T : supp µ → R 3 such that We emphasize that T (R 3 ) ⊂ supp ν µ -almost everywhere. Indeed Remark 7. Note that, for all transport map T from µ to ν one may associate a transference plan (Id, T )#µ ∈ Π(µ , ν) i.e. the pushforward of µ by the map x → (x, T (x)) and we have (Id, T )#µ − esssup |x − y| , Note that this yields It is then natural to investigate in which conditions one has the existence of a transport map T associated to an optimal transference plan. As in [9] we refer to [3] for the following existence result.
Theorem 5.6 (Champion, De Pascale, Juutinen). Assume that µ is absolutely continuous with respect to the Lebesgue measure. Then there exists optimal transference plans, and at least one of them is given by a transport map T . If moreover ν is a finite sum of Dirac masses, this optimal transport map is unique.

5.1.
Existence, uniqueness and stability for the mean-field equation. Consider the following problem ∂ρ ∂t + div((κg + Kρ)ρ) = 0 , ρ(0, ·) = ρ 0 , recall the definition of the kernel K . We refer to the existence and uniqueness result due to Höfer [11,Theorem 9.2] in the case where the initial data ρ 0 and its gradient ∇ρ 0 are in the space X β for some β > 2 where Theorem 5.7 (Höfer). Assume that ρ 0 , ∇ρ 0 ∈ X β for β > 2. There exists a unique solution ρ ∈ W 1,∞ ((0, T ), X β ) to equation (73) for all T > 0 and a unique well defined flow X satisfying such that Remark 8. The flow X is measure-preserving i.e. for a test function φ ∈ C b (R 3 ) we have φ(y)ρ(s, y)dy = φ(X(s, t, y))ρ(t, y)dy, for all s , t ∈ [0, T ]. This allows us to separate the dependence of time s in the integral with respect to the measure ρ(t, ·).
Remark 9. Note that for all η ∈ L ∞ (R 3 )∩L 1 (R 3 ), the velocity field Kη is Lipschitz Moreover, if one assume that ρ 0 is only Lipschitz and compactly supported, then one can show the existence and uniqueness of the solution ρ to equation (73) in the space L ∞ ((0, T ); L ∞ (R 3 ) ∩ L 1 (R 3 )). The method of proof is related to the stability result due to G. Loeper in [15] which gives a stability estimate in terms of Wasserstein distance for the Vlasov-Poisson equation. This result is adapted by M. Hauray in [8, Theorem 3.1] for a more general class of kernels K satisfying a (C α ) see [8]. This condition being satisfied by the Oseen tensor Φ we have the following result.
Theorem 5.8 (Hauray-Loeper). Given T > 0, consider two solutions ρ 1 , ρ 2 ∈ L ∞ ((0, T ), L ∞ (R 3 ) ∩ L 1 (R 3 )) of equation (73) associated to two initial data ρ 1 0 , ρ 2 0 ∈ L ∞ (R 3 ) ∩ L 1 (R 3 ). There holds We refer to [8,Theorem 3.1] for a complete proof which introduces the main ideas used also in [9] for the mean field approximation result. 5.2. ρ N as a weak solution to a transport equation. According to Theorem 1.2, there exists a time T > 0 independent of N for which the particles do not overlap. This shows that the empirical measure is well defined on [0, T ]. Recall that we are interested in the limiting behaviour of ρ N ∈ P([0, T ] × R 3 ) when N → ∞. According to Proposition 6, particles (x i ) 1≤i≤N satisfy the following system: In order to prove Theorem 1.3 we want to compare the particle system to the continuous density ρ which is solution to equation (73). Hence, we need to express ρ N as a weak solution to a transport equation. The remainder of this Section is devoted to establish such a formulation.
Analogously to the continuous case, we are interested in giving a sense to the quantity which is not well defined because Φ is singular. On the other hand, as the only values of Φ that matters are the terms Φ(x i − x j ), i = j we define the following regularization and ψ is a truncation function such that ψ = 0 on B(0, 1/4) and ψ = 1 outside B(0, 1/2). We can now define the operator K N Since Theorem 1.2 ensures that the particles satisfy Now, it remains to check that ρ N is a weak solution of a transport equation. We recall that ρ N is a weak solution of a transport equation

AMINA MECHERBET
Note that this integral yields ) .
In particular if we choose V such that V (t, x i (t)) = V i (t) one has On the other hand, we recall that from Proposition 6 we can write for all 1 ≤ i ≤ N we can define V as Construction of E N . We fix χ a truncation function such that χ = 1 on B(0, 1) and χ = 0 on c B(0, 2). For all i we set By construction, E i is a divergence-free compactly supported vector field satisfying Furthermore, E i is supported in B(x i (t), 2R). Thanks to Theorem 1.2, this entails that supp(E i ) ∩ supp(E j ) = ∅ for i = j. We set then By construction, this velocity field is divergence-free and is regular The only statement that needs further explanation is (77). For all and for all x ∈ B(x i , 2R) \ B(x i , R), direct computations yields Therefore We can now state the following proposition.
Proposition 8. For arbitrary N we have that κg . Moreover, the velocity field satisfies for some positive constant C independent of N .
Proof. As the kernel is regularized, the two first properties are satisfied by construction. For all (t, x) ∈ [0, T ] × R 3 we have Reproducing the arguments of Lemma A.1 for k = 1 together with assumptions (6), (7) and Theorem 1.2 yields Furthermore, the velocity field E N is uniformly bounded according to (77).
This allows us to state the following result.
on [0, T ] × R 3 . Moreover, the characteristic flow defined for all s, t ∈ [0, T ] by is of class C 1 for all N ≥ 1 and the following classical formula holds true: thus, ρ N is a weak solution for (79). According to Proposition 8, the ode governing the characteristic flow satisfies the assumptions of the Cauchy-Lipschitz theorem. Therefore, the ode admits a unique maximal solution X N ∈ C 1 ([0, T ] × [0, T ] × R 3 ) thanks to formula (78). Equality (81) holds true thanks to the classical theory for transport equations.
6. Proof of Theorem 1.3. At this point, we proved that the particles interact two by two with an interaction force given by the Oseen-tensor with an additional error term.
We want to estimate the Wasserstein distance W 1 (ρ N (t, ·), ρ(t, ·)) for all time 0 ≤ t ≤ T . To this end, we follow the ideas of [8] and [9] and show that the additional error term E N can be controled. As in [9], we introduce an intermediate densitȳ ρ N . 6.1.
Step 1. Estimate of the distance between ρ andρ N . We defineρ N 0 as the regularized density of ρ N 0 :ρ N 0 := ρ N 0 * χ λ N where χ λ N (x) := 1 |λ N | 3 χ x λ N a mollifier compactly supported. Note that the support of χ is not important, we consider for instance χ such that supp χ = B(0, 1). We emphasize that the regularized density is uniformly bounded according to assumption (5). Moreover, we have Now, we defineρ N as the unique solution to equation (73) associated to the initial dataρ N 0 . The stability Theorem 5.8 allows us to compare ρ andρ N : where C = C( χ ∞ ,M , ρ 0 ∞ ). We split the distance W 1 (ρ 0 ,ρ N 0 ) as follows , and use the fact that to get 6.2.
Step 2. Estimate of the distance betweenρ N to ρ N . It remains to estimate W 1 (ρ N (t, ·),ρ N (t, ·)). We have the following result.
Lemma 6.1. For N large enough, there exists a positive constant C such that for all t ∈ [0, T ] W 1 (ρ N (t, ·),ρ N (t, ·)) λ N + td N min e Ct . Theorem 1.3 is a consequence of estimate (86) and Lemma 6.1. The rest of this Section is devoted to proving the above lemma.
Proof of Lemma 6.1. According to Theorems 5.7 and 5.9 we have the explicit formulas for all s , t ∈ [0, T ]ρ N (t, ·) = X(t, s, ·)#ρ N s , ρ N (t, ·) = X N (t, s, ·)#ρ N s . At t = 0 we have the existence of an optimal transport map T 0 fromρ N 0 to ρ N 0 thanks to Theorem 5.6 We construct then a transport map T t fromρ N to ρ N at all time t ∈ [0, T ] by following T 0 along the two flows X and X N T t = X N (t, 0, ·) • T 0 • X(0, t, ·).
One can remark that for all 0 ≤ s ≤ t: As in [9] we set then , and thanks to (85) we have We reproduce the same steps as in [9] and introduce the following notation for a generic "particle" of the continuous system with position x t at time t such that we fix in what follows 0 ≤ t 2 ≤ t 1 and recall the following formula We aim now to estimate |T t1 (x t1 ) − x t1 | for all test particle x t1 where we used the fact that ρ N s = T s #ρ N s to get We set then t 1 = t and t 2 = t 1 − τ = t − τ , τ > 0. We obtain for almost every x t here we used Remark 8 with y s = X(s, t, y t ). In addition we defined This being true for almost every x t we obtain Hence, it remains to control the last quantity. We split the integral on R 3 into two terms: the first one denoted J 1 is the integral over the subset I and the second one denoted J 2 the integral over R 3 \ I where where L will be defined later.
Using Remarks 8 and 9, formula (75) and the uniform bounds (83), (84), the Lipschitz constant of Kρ N is uniformly bounded. This allows us to define the constant L as We can make precise now the constant L := Lip (Kρ N ) which is uniformly bounded with respect to N and t ∈ [0, T ].
We have for all 0 ≤ t − τ ≤ s ≤ t and τ small enough Analogously, for almost all x s and y s where we used the fact that f N (t) ≥ f N (s). According to the definition of I = {y t : |x t − y t | ≥ 4f N (t)e τ L }, this yields for τ small enough Moreover, recall that T s (x s ) and T s (y s ) are in the support of ρ N (s, ·) i.e. there exists i , j such that T s (x s ) = x i (s) and T s (y s ) = x j (s). In addition, estimate (90) and the definition of I ensures that i = j. We have then Finally, using estimates (89), (90), formula (91) and the Lipschitz-like estimate (27) for Φ we obtain where we used Remark 8, formula (75) and the uniform bounds (83), (84).
We focus now on For the remaining term we apply Theorem 1.2 and get for all t − τ ≤ s ≤ t

Conclusion.
Gathering these bounds, there exists a constant K > 0 independent of N such that for τ small enough and 0 < t ≤ T .
We can now apply a discrete Gronwall argument: Note that at time t = 0, assumption (11) and formula (87) ensures the existence of a positive constant C 1 > 1 such that hence, we define T * ≤ T as the maximal time for which which vanishes according to assumption (7) and (11). This shows that we can take N large enough and depending on T , K and C 1 such that T * → T and formula (93) holds true up to time T . Hence, for N large enough we have for all t ∈ [0, T ] f N (t) ≤ f N (0)e 2C1t + te 2C1t E N ∞ . Using (87) and the fact that W 1 (ρ N ,ρ N ) ≤ W ∞ (ρ N ,ρ N ) ≤ f N , this implies Lemma 6.1.
Appendix A. Technical lemmas. We state here an important lemma which is the extension of [13, Lemma 2.1] to the new assumptions on the dilution regime introduced in [10]. We introduceρ N an approximation of ρ N defined as ρ N is L ∞ and using (5), one can check that Moreover,ρ N is L 1 and we have ρ N L 1 = 1 by construction. Lemma A.1. For all k ∈ [0, 2], under assumptions (5), (7), if N is large enough, there exists a positive constant C > 0 such that for all fixed 1 ≤ i ≤ N : Moreover, if k = 3 we have Proof. We fix i = 1 and the same holds true for all 1 ≤ i ≤ N . We use the following shortcut I 1 := {j ∈ {1, · · · , N } such that |x 1 − x j | ∞ ≤ λ N }.
The sum can be written as follows: .
For the second term in the right hand side, note that, for all y ∈ B ∞ (x j , λ N /3), j ∈ I 1 we have Hence, we have for all constant L > 2/3λ N One can show that the optimal constant L > 2/3λ N is L = 1 M 1/3 . Since lim N →∞ λ N = 0, this choice of L is possible for N large enough such that λ N < 3 2M 1/3 . Hence, we If k = 3, we integrate the term r −1 keeping the same value for L as before The following results are used for the control of the particle concentration M N : We recall the definition of L N introduced in (71): The following lemma shows that the two definitions are equivalent.
Indeed, for all x ∈ R 3 there existsx k , k = 1, · · · , 8 such that this yields i ∈ {1, · · · , N } such that x i ∈ B ∞ (x, λ N ) Taking the supremum in the right hand side and then in the left one we obtain Moreover, we remark that the supremum in the right hand side over all x ∈ R 3 can be reduced to the supremum over which means that for all x ∈ i B ∞ (x i , λ N 2 ) there exists 1 ≤ i 0 ≤ N such that 1 ≤ j ≤ N, such that x j ∈ B ∞ (x, λ N /2) Taking the maximum over all i 0 in the right hand side, and then the supremum over all Gathering inequality (97) and (98) concludes the proof.
More generally we define for all β > 0: Proof. For sake of clarity we set β = 1 and the proof remains the same for all β > 0. The idea is to show an equivalent formula for M N and use Lemma A.2. Analogously to the proof of Lemma A.2, for all x ∈ R 3 there existsx k , k = 1, · · · , λ 3 such that This yields, with the definition of M N λ : M N α ≤ α 3 M N (t). Finally, we apply Lemma A.2 to get , which completes the proof.