A model for suspension of clusters of particle pairs

In this paper, we consider $N$ clusters of pairs of particles sedimenting in a viscous fluid. The particles are assumed to be rigid spheres and inertia of both particles and fluid are neglected. The distance between each two particles forming the cluster is comparable to their radii $\frac{1}{N}$ while the minimal distance between the pairs is of order $N^{-1/2}$. We show that, at the mesoscopic level, the dynamics are modelled using a transport-Stokes equation describing the time evolution of the position $x$ and orientation $\xi$ of the clusters. Under the additional assumption that the minimal distance is of order $N^{-1/3}$, we investigate the case where the orientation of the cluster is initially correlated to its position. In this case, a local existence and uniqueness result for the limit model is provided.


Introduction
We consider the problem of N rigid particles sedimenting in a viscous fluid under gravitational force. The inertia of both fluid and particles is neglected. At the microscopic level, the fluid velocity and the pressure satisfy a Stokes equation on a perforated domain. The mathematical derivation of models for suspensions in Stokes flow interested a lot of researches. One of the most investigated question is the effective computations of quantities such as the viscosity or the average sedimentation velocity, see for instance [2,7,8,10,13,12,20,30,33,35,37] and all the references therein. Regarding the analysis of the associated homogenization problem, it has been proved that the interaction between particles leads to the appearance of a Brinkman force in the fluid equation. This Brinkman force depends on the dilution of the cloud but also the geometry of the particles, see [1,3,5,8,18,19,34]. In the dynamic case, the justification of a mesoscopic model using a coupled transport-Stokes equation has been proved in [24] where authors show that the interaction between particles is negligible in the dilute case i.e. when the minimal distance between particles is larger than 1 N 1/3 . In [22,32] the justification has been extended to regimes that are not so dilute but where the minimal distance between particles is still large compared to the particles radii. The coupled equations derived are: (1)    ∂ t ρ + div((κg + u)ρ) = 0 −∆u + ∇p = 6πr 0 κgρ , div(u) = 0.
Here u is the fluid velocity, p its associated pressure, ρ is the density of the cloud. r 0 = RN, where R is the particles radii, g the gravity vector. The velocity κg = 2 9 R 2 (ρ p − ρ f )g represents the fall speed of a sedimenting single particle under gravitational force. The derivation of this model is a consequence of the method of reflections which consists in approaching the flow around several particles as the superposition of the flows associated to one particle at time, see [36], [27,Chapter 8], [31], [11,Section 4], [28], [23] for more details.
In this paper, we are interested in the case where the cloud is made up of clusters i.e. the case where the minimal distance between the particles is proportional to the particles radii R. The main motivation is to show the influence of the clusters configuration on the mean velocity fall. A first investigation in this direction is to consider clusters of pairs of particles where the minimal distance between the particles forming the pair is comparable to their radii. The cluster configuration is determined by the center x and the orientation ξ of the pair. Starting from the microscopic model and assuming the propagation in time of the dilution regime, the first result of this paper is the derivation of a fluid-kinetic model describing the sedimentation of the suspension at a mesoscopic scaling. The fluid-kinetic model obtained couples a Stokes equation for the fluid velocity and pressure (u, p) with a transport equation for the function µ(t, x, ξ) representing the density of clusters centered in x and having orientation ξ at time t, see Theorem 0.1. The mean velocity fall of clusters is formulated through the Stokes resistance matrices while the variation of its orientation involves the gradient of the fluid velocity. In particular, the presence of the gradient of the fluid velocity suggests a similarity with the model of suspension of rod-like particles where the density function depends on the center of the rods x and there orientations n ∈ S 2 , see [6,16,17] and the references therein.
Note that one can reproduce the same arguments for a cloud of clusters constituted of k ≥ 2 identical particles by considering the center of mass of the cluster x = 1 k k p=1 x p and , q = 1, · · · , k − 1. Extending the assumptions for the general case k ≥ 2, the limit density µ and the kinetic equation would depend on k variables. However, the relevance of such model can be discussed, in particular with comparison to the models of suspension of polymeric fluids where each polymer is represented by beads connected along a chain. The total number of beads k, also called degree of polymerization, can reach for instance 2000 in the case of ductile materials such as plastic films. Hence, in practice, a fine description of a polymer chain is not suitable and it is more convenient to deal with a coarse model where k = 2, see dumbbell model [6,29].
The second result of this paper corresponds to the case where the orientation of the cluster is correlated to its center i.e. ξ = F (t, x). Under additional assumptions, see Theorem 0.2, the derived model is a transport equation for ρ coupled to a Stokes equations for the fluid velocity and pressure (u, p) and a hyperbolic equation for the function describing the evolution of the cluster orientation F . A local existence and uniqueness result for the former system is also presented, see Theorem 0.3.
The starting point is a microscopic model representing sedimentation of N ∈ N * particle pairs in a uniform gravitational field. The pairs are defined as where x i 1 , x i 2 are the centers of the i th pair and R the radius. We define (u N , p N ) as the unique solution to the following Stokes problem : (2) −∆u N + ∇p N = 0, div u N = 0, completed with the no-slip boundary conditions : where (U i 1 , U i 2 ) ∈ R 3 × R 3 , 1 ≤ i ≤ N are the linear velocities. In this model, the angular velocity is neglected and we complete the PDE with the motion equation for each couple of particles : Newton law yields the following equations where inertia is neglected : where m is the mass of the identical particle adjusted for buoyancy, g the gravitational acceleration, F i 1 , F i 2 are the drag forces applied by the fluid on the i th particle : with n the unit outer normal and Σ(u N , p N ) = (∇u N + (∇u N ) ⊤ ) − p N I the stress tensor. In order to formulate our results we introduce the main assumptions on the cloud. 0.1. Assumptions and main results. We assume that the radius is given by R = r 0 2N . In this paper we use the following notations, given a pair of particles B(x 1 , R) and B(x 2 , R): Let T > 0 be fixed. We introduce the empirical density and set ρ N its first marginal: We denote by d min the minimal distance between the centers x i + : We assume that there exists two constants M 1 > M 2 > 1 independent of N such that: We assume that µ N converges in the sense of measures to µ i.e. for all test function We assume that the first marginal of µ denoted by ρ is a probability measure such that ρ ∈ L ∞ (R 3 ) ∩ L 1 (R 3 ). We introduce W ∞ (t) := W ∞ (ρ N (t, ·), ρ(t, ·)) the infinite-Wasserstein distance between ρ N and ρ, see (19) for a definition. We assume that For the first result, we assume that there exists a positive constant E 1 > 0 such that: Regarding the second result, we assume in addition that there exists a positive constant E 2 > 0 such that: Remark 0.1. Since ρ ∈ L ∞ (0; T, L ∞ (R 3 )), this yields a lower bound for the infinite Wasserstein distance for all t ∈ (0, T ) and all N ∈ N * : On the other hand, the definition of the infinite Wasserstein distance ensures that which yields according to (9) (14) sup Assumption (11) is only needed for the second Theorem 0.2. Precisely, under assumption (10), the minimal distance is at least of order C √ N and R ≪ d min . Indeed using (12), (10) Whereas under the additional assumption (11), the threshold for the minimal distance is of order C N 1/3 . Our main results read: (7), (8), (9) and (10) is small enough, µ satisfies the following transport equation : Remark 0.2. The matrix A is defined as A := A 1 + A 2 where A 1 and A 2 are the resistance matrices associated to the sedimentation of a couple of identical spheres, see Section 1.1 for the definition. The term (A) −1 κg represents the mean velocity of a couple of identical particles sedimenting under gravitational field. We assume herein that Remark 0.3. In the case where µ 0 is compactly supported with respect to the second variable ξ uniformly in the first variable x, local existence and uniqueness of the above coupled system can be shown following the result of [21,Chapter 8] for the model of suspension of rod-like particles. In particular, the L 1 norm of the spatial density ρ is conserved in time while the L ∞ norm of ρ(t, ·) is bounded by sup This ensures existence and uniqueness for a small time T and ρ L ∞ (0,T ;L ∞ (R 3 )) is controlled by µ 0 ∞ and sup The second result concerns the case where the vectors along the line of centers ξ i are correlated to the positions of centers x i + . Theorem 0.2. Assume now that ρ ∈ W 1,∞ (R 3 )∩W 1,1 (R 3 ), A −1 ∈ W 2,∞ (R 3 ) and consider the additional assumption (11). Assume that there exists a function There existsT ≤ T independent of N and unique F N ∈ L ∞ (0,T ; W 1,∞ (R 3 )) such that for all t ∈ [0,T ] we have: Moreover, the sequence (F N ) N admits a limit F ∈ L ∞ (0,T ; W 1,∞ (R 3 )). The limit measure µ is of the form µ = ρ ⊗ δ F and the triplet (ρ, F, u) satisfies the following system We finish with a local existence and uniqueness result for the limit model.
As in [32], the idea of proof of Theorem 0.1 and 0.2 is to provide a derivation of the kinetic equation satisfied weakly by µ N . This is done by computing the first order terms of the velocities of each pair: The interaction force Φ is the Oseen tensor, see formula (18). This development is a corollary of the method of reflections which consists in approaching the solution u N of 2N separated particles by the superposition of fields produced by the isolated 2N particle solutions. We refer to [36], [31], [27,Chapter 8] and [11,Section 4], [28] for an introduction to the topic. We also refer to [23] where a converging method of reflections is developed and is used in [22]. In this paper we reproduce the same method of reflections developed in [32,Section 3]. However this method is no longer valid in the case where the minimal distance is comparable to the particle radii. The idea is then to approach the velocity field u N by the superposition of fields produced by the isolated N couple of particles This requires an analysis of the solution of the Stokes equation past a pair of particles. In particular, we need to show that these special solutions have the same decay rate as the Stokeslets, see [32, Section 2.1]. Precisely, in Section 1, we prove that the solution of the Stokes equation past a pair of particles can be approached by the Oseen tensor at first order. The convergence of the method of reflections is ensured under the condition that the minimal distance d min between the centers x i + satisfies and that the distance |x i 1 − x i 2 | for each pair satisfies formula (7). In this paper, we focus only on the derivation of the mesoscopic model. Precisely, we do not tackle the propagation in time of the dilution regime and the mean field approximation. We provide in Propositions B.3 and B.1 some estimates showing that the control on the minimal distance d min depends on the control on the infinite Wasserstein distance W ∞ (ρ N , ρ). However, the gradient of the Oseen tensor appearing in equation (17) leads to a log term in the estimates involving the control of W ∞ (ρ N , ρ), see Proposition B.2. This prevents us from performing a Gronwall argument in order to prove the mean field approximation in the spirit of [14,15]. 0.2. Outline of the paper. The remaining sections of this paper are organized as follows. In section 1 we present an analysis of the particular solution of two translating spheres in a Stokes flow. The main result of this section is the justification of the approximation of this particular solution using the Oseen tensor and proving some decay properties similar to the Stokeslets. In section 2 we present and prove the convergence of the method of reflections using the estimate of Appendix A. We also present two particular cases of the application of this method which are useful later. In section 3 we compute the particle velocities (ẋ i + ,ξ i ) 1≤i≤N using the estimates provided in the previous section. Section 4 is devoted to the proof of the first Theorem 0.1. Precisely, we prove that the discrete density µ N satisfies weakly a transport equation (46) which can be seen as a discrete version of the limit equation (15). In particular, equation (46) is formulated using a discrete convolution operator K N ρ N ∼ Φ * ρ N defined rigorously in Section 4. The convergence proof is obtained by showing that K N ρ N converges to the continuous convolution operator Kρ = Φ * ρ. Convergence estimates of K N ρ N − Kρ are provided in the Appendix B. Section 5 is devoted to the proof of Theorems 0.2 and 0.3. The first step is to prove local existence and uniqueness results for the correlation function F N solution of (48) and also for F the solution of the hyperbolic equation (54). The idea is to apply a fixed-point argument using some stability estimates provided in the last Appendix C. The last part of Section 5 concerns the convergence of the microscopic model to the mesoscopic model. This convergence result is obtained by showing that the sequence F N converges is some sense to F . 0.3. Notations. In this paper, n always refers to the unit outer normal to a surface and dσ denotes the measure integration on the surface of the particles. We recall the definition of the Green's function for the Stokes problem (U, P) where U is also called the Oseen tensor, See Φ(x) = 1 8π Given two probability measures ν 1 , ν 2 , we define the infinite Wasserstein distance as where Π(ν 1 , ν 2 ) is the set of all probability measures on R 3 × R 3 with first marginal ν 1 and second marginal ν 2 . In the case where ν 1 is absolutely continuous with respect to the Lebesgue measure, then according to [4] the following definition holds true In particular, this distance is well adapted to the estimates of the discrete convolution operator K N ρ N defined in (43). Precisely, the infinite Wasserstein distance allows to localise the singularity of the Oseen tensor and is closely related to the minimal distance d min . We refer also to [14,15,4,32] for more details. Given a couple of velocities (U 1 , U 2 ) ∈ R 3 × R 3 we use the following notations 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 , ρ L ∞ (0,T ;L ∞ (R 3 )) , E 1 , E 2 and eventually on κ|g| which is uniformly bounded, see [32].

Two translating spheres in a Stokes flow
In this section, we focus on the analysis of the Stokes problem in R 3 past a pair of particles. Given x 1 , x 2 ∈ R 3 and R 1 , R 2 > 0, such that |x 1 − x 2 | > R 1 + R 2 , we consider two spheres B α := B(x α , R α ) α = 1, 2 and focus on the following Stokes problem: (20) −∆u + ∇p = 0, div u = 0, on R 3 \B 1 ∪B 2 , completed with the no-slip boundary conditions: where U α ∈ R 3 for α = 1, 2. Classical results on the Steady Stokes equations for exterior domains (see [9, Chapter V] for more details) ensure the existence and uniqueness of equations (20) - (21). In this section, we aim to describe the velocity field u in terms of the force applied by the fluid on the particles defined as: We refer to the paper [25] for the following statements. Neglecting angular velocities and torque we emphasize that there exists a linear mapping called resistance matrix satisfying: where A αβ , 1 ≤ α, β ≤ 2, are 3 × 3 matrices depending only on the non-dimensionalized centre-to-centre separation s and the ratio of the spheres' radii λ: each of these matrices is of the form: where I is the 3 × 3 identity matrix and g α,β , h α,β are scalar functions. We refer to the paper of Jeffrey and Onishi [25] where the authors provide a development formulas for g α,β and h α,β given by a convergent power series of |s| −1 . Note that the matrices satisfy Inversly, there exists also a linear mapping called mobility matrix such that The matrices a α,β depend on the same parameters as matrices A α,β and satisfy a formula analogous to (23). They are also symmetric in the sense of formula (24). The resistance and mobility matrices satisfy the following formula: a 11 a 12 a 21 a 22 = I 0 0 I , Again, we refer to [25] for more details.
1.1. Restriction to the case of two identical spheres. We simplify the study by assuming that R 1 = R 2 = R i.e. λ = 1. This means that the resistance matrix depends only on the parameter s which becomes: and we have: Hence we reformulate the resistance matrix as follows: and the mobility matrix: Formula (26) yields the following relations We are interested in providing a formula for the velocity u and showing some decay properties. In this paper we use the notation ( completed with the no-slip boundary conditions: Note that there is no ambiguity regarding the dependence of the solution U[U 1 , U 2 ] with respect to x + and ξ. Indeed, in this paper, since we consider the solutions U[U i 1 , U i 2 ] associated to each cluster B i with some velocities U i 1 , U i 2 , the dependence with respect to the centers x i + and orientations ξ i is implicitly given by the dependence of the velocities with respect to i. The main result of the section is the following Moreover, there exists a positive constant independent of U 1 , U 2 , ξ, x + and depending only on The unique solution (U[U 1 , U 2 ], P [U 1 , U 2 ]) satisfies the following decay property with C(M 1 ) independent of x + , ξ, U 1 and U 2 .
Proof. We consider the case where x + = 0 and R = 1, the generalization to arbitrary x + and R can be obtained by scaling arguments. In what follows we use the short cut (u, p) := (U[U 1 , U 2 ], P [U 1 , U 2 ]) and extend u by U α on B(x α , 1), α = 1, 2 and we have u ∈Ḣ 1 (R 3 ). We consider a regular truncation function χ = 0 on B(0, 2M 1 ) ⊃ B(x 1 , 1) ∪ B(x 2 , 1) and χ = 1 on c B(0, 3M 1 ) and we set Lemma 18] for instance, and satisfies Using Stokes regularity results, see [9, Theorem IV.4.1], combined with some Sobolev embeddings we haveū ∈ C ∞ (R 3 ) and satisfies a Stokes equation on R 3 with a source term f = − div Σ(ū,p) having support in B(0, 3M 1 ) \ B(0, 2M 1 )). Hence we can apply the convolution formula with the Green function Φ and writē Note that u =ū on c B(0, 3M 1 ). We may then apply a Taylor expansion of Φ(· − y) for |x| > 3M 1 and get An integration by parts for the first term yields we recall that in the above computations the unit normal vector n is pointing outward. It remains to estimate the error term, we recall that using the Bogovskii properties and the on the other hand, an integration by parts together with (27) yields For the remaining term we introduce G(x, y) where ψ = 0 on c B(0, 7/2M 1 ) and ψ = 1 on B(0, 3M 1 ). With this construction and since this yields using the decay property of the Oseen tensor for all |x| > 4M 1 and y ∈ we conclude by using the fact that |ξ| ≤ M 1 and the uniform bounds on A 1 , A 2 , see Remark 0.2.

The method of reflections
In this section, we aim to show that the method of reflections holds true in the special case where the minimal distance and the radius R are of the same order. The idea is to approach the velocity field u N by the particular solutions developed in the section above. We recall that u N is the unique solution to the following Stokes problem : completed with the no-slip boundary conditions : Thanks to the superposition principle, the sum of the i , but does not match the boundary conditions. Hence, we define the error term: which satisfies a Stokes equation on We set then for α = 1, 2 and 1 ≤ i ≤ N: , and reproduce the same approximation to obtain: which satisfies a Stokes equation with the following boundary conditions for all 1 ≤ i ≤ N, α = 1, 2 and x ∈ B(x i α , R): By iterating the process, one can show that for all k ≥ 1 we have: where for all α = 1, 2, 1 ≤ i ≤ N and p ≥ 0: The convergence is analogous to the convergence proof in [32, Section 3.1]. We begin by the following estimates that are needed in the computations.  (7), (10) we have for all 1 ≤ i = j ≤ N, 1 ≤ β ≤ 2: The first step is to show that the sequence max |)) converges when p goes to infinity. (7), (8), (10) and the assumption that r 0 ρ

Lemma 2.2. Under assumptions
Proof. According to formulas (32) and Lemma 2.1, we have for all α = 1, 2 and 1 ≤ i ≤ N: where we used Lemma A.1 for k = 1. Hence, the first term in the right-hand side vanishes according to (10) and (14). Finally, if we assume that r 0 ρ is small enough, we obtain the existence of a positive constant K < 1/2 such that: We have the following result.
Proof. The proof is analogous to the convergence proof of [32,Proposition 3.4]. This is due to the fact that the particular solutions have the same decay rate as the Oseen-tensor, see (32).
2.1.1. First case. Given W ∈ R 3 we consider in this part w the unique solution to the Stokes equation (2) completed with the following boundary conditions : We denote by W i,(p) α , α = 1, 2, 1 ≤ i ≤ N, p ∈ N the velocities obtained from the method of reflections applied to the velocity field w. In other words : We aim to show that, in this special case, the sequence of velocities W i,(p) α and the error term U[w (k) * ] are much smaller than before. This is due to the initial vanishing boundary conditions for i = 1. Indeed we have : Proof. We show that the statement holds true for p = 0 then we prove it for all p ≥ 1 by induction. According to formula (34) we have for p = 0: and for i = 1, α = 1, 2, U i,(0) α = 0. Using (30), this yields for i = 1, α = 1, 2: where we used the fact that the radius R is comparable to |x 1 − | thanks to (7). Thus, we denote by C > 0 the maximum between the global constant appearing in (32) and the one in the above estimate.
This shows that the first statement holds true for p = 0. For the second estimate we have |W 1,(1) α | = 0 and for p = 1 we have using the decay rate (32) where we used Lemma A.1 for k = 2 and assumption (10). We define then the constant L > 0 as the constant satisfying: Now for all p ≥ 1, i = 1 we have using again (32) Since Rd 1i d min ≪ r 0 L, the second term can be bounded by (Cr 0 L) p 2 p−2 which yields the expected result because 2 p−1 + 2 p−2 ≤ 2 p . We prove now the second estimate. Let p ≥ 1, using the decay rate (32) : According to these estimates and the definition of L (37), if we assume that r 0 max( ρ L ∞ (0,T ;L ∞ (R 3 )) , ρ ) ) is small enough to have 2LCr 0 < 1 then the following result holds true : This result shows that we can obtain a better estimate for the error term of the method of reflections in this particular case: And for i = 1 we have : , with α = 1, 2 and i = 1, formula (34) yields: We estimate the first term applying Corollary 2.5 and the same arguments as before |W |, We reproduce the same for the second term applying Corollary 2.5: For the last term, according to (30) we have : Thus F 1 2 = −F 1 1 and we obtain using the decay rate of R (31) together with the fact that |x 1 − | is comparable to R thanks to assumption (7) : Gathering all the inequalities we have for i = 1: Analogously for i = 1 we apply Lemma A.1 for k = 3 and obtain up to a constant depending on ρ L ∞ (0,T ;L ∞ (R 3 )) : We have according to formula (34) : where η = 2Cr 0 L < 1 is the constant appearing in Proposition 2.4. Reproducing the same computations as before yields: In the case i = 1 we have: Thanks to these estimates we have the following convergence rate: Proof. Reproducing exactly the same proof as in [32, Proposition 3.4], the main difference appears in the last estimate where we apply Proposition 2.6: Taking the limit when k goes to infinity we get:

AMINA MECHERBET
The term inside brackets is bounded as follows: according to (12).

Second case.
Given W ∈ R 3 we consider in this part w the unique solution to the Stokes equation (2) completed with the following boundary conditions : Denote by W i,(p) α , α = 1, 2, 1 ≤ i ≤ N, p ∈ N the velocities obtained from the method of reflections applied to the velocity field w. In other words : We aim to show that, in this special case, the sequence of velocities W i,(p) α are also smaller than the general case. This is due to the initial boundary conditions which vanish for i = 1. Indeed we have : Proposition 2.8. There exists two positive constants C > 0 and L = L( ρ L ∞ (0,T ;L ∞ (R 3 )) ) such that : Proof. The proof is analogous to the one of Proposition 2.4.

Extraction of the first order terms for the velocities
This section is devoted to the computation of the velocities U i + , U i − for 1 ≤ i ≤ N. The idea of proof is to apply the method of reflections to the velocity field u N as presented above and we set : we also use the following notations for the forces associated to the solutions U 3.1. Preliminary estimates.
Proof. We prove the formula for i = 1 and the same holds true for all 1 ≤ i ≤ N. We set w the unique solution to the Stokes equation (2) completed with the following boundary conditions : with W an arbitrary vector of R 3 . We use the method of reflections to obtain : For the last term we apply again the method of reflections to the velocity field w, see Section 2.1.2. We set: We obtain : Thanks to the method of reflections, the second term on the right hand side can be bounded by . For the first term, direct computations using (27) show that Using Corollary 2.9 we get that This being true for all W ∈ R 3 it yields: Using the definitions of F 1,∞ 1 and F 1,∞ 2 , see (39), this becomes: Recall that A 1 (ξ) and A 2 (ξ) are of the form h 1 (|ξ|)I + h 2 (|ξ|) ξ⊗ξ |ξ| 2 . Moreover, according to formulas (29) A 1 + A 2 (resp. A 1 − A 2 ) is invertible and its inverse is (a 1 + a 2 ) (resp. a 1 − a 2 ). Thus : We use the fact that (A 1 (ξ 1 ) + A 2 (ξ 1 )) −1 is uniformly bounded independently of the particles and N to get On the other hand, as (U we rewrite formula (41) as : Using again formula (34) this yields : ](x 1 2 ), We conclude by recalling that (A 1 + A 2 ) −1 ∞ = A −1 ∞ is uniformly bounded. Applying the same ideas we obtain the following result: Proof. The proof is analogous to the one of Proposition 3.1. The idea is to consider this time w the unique solution to the Stokes equation (2) completed with the following boundary conditions : Proof. First of all, from Propositions 3.1 and 3.2 we can show that the velocities U i α are uniformly bounded with respect to N for all 1 ≤ i ≤ N and α = 1, 2. Indeed, using formula (34) together with the decay properties (32) and Propositions 3.1 and 3.2 we have : This allows us to bound the terms max We recall that using (30) we have , see proof of Propositions 3.2 and 3.1. Hence, we replace F j 1 + F j 2 = 2F j,∞ + by −2mg with 6πRκg = mg and bound the error terms using the decay properties of the Oseen tensor Φ and the field R, see (31) j =i

Now it remains to replace both terms Φ(x j
which is comparable to R according to assumption (7). Gathering all the estimates, the error term is of order √ R which is of order d min according to assumption (10) and Remark 0.1.

Estimates forẋ
Proof. The first formula of Proposition 3.2 together with the uniform bound on the velocities (U i + , U i − ), see proof of Corollary 3.3, yields: We want to estimate the first term, we have using (30) Now recall that, from the proof of Proposition 3.1 we have: Thus, we get the following formula: We recall that mg = 6πRκg = 1 2 6πr 0 N g. Finally we obtain: It remains to bound the error terms. We begin by the first one: We emphasize that for all y ∈ [x i 1 , x i 2 ]: where we used the fact that For the second error term we have: which yields the same estimate as for the first error term. Finally, the last error term gives: where we used the fact that F j,∞

Proof of Theorem 0.1
In order to derive the transport-Stokes equation satisfied at the limit, the idea is to show that the discrete density µ N satisfies weakly a transport equation. We introduce the following notations. Given a density ρ, we define the operator Kρ as: The operator is well defined and is Lipschitz in the case where ρ ∈ L 1 ∩ L ∞ . Moreover, note that Kρ satisfies the Stokes equation −∆K(ρ) + ∇p = 6πr 0 κgρ, on R 3 . Analogously, we define K N ρ N as: where χΦ(·) = χ · d min Φ(·), χ is a truncation function such that χ = 0 on B(0, 1/4) and χ = 1 on c B(0, 1/2).

4.1.
Derivation of the transport-Stokes equation. The transport equation satisfied by µ N is obtained directly using the ODE system derived for each couple (x i + , ξ i ). We recall that: Following the idea of [32, Section 5.2] and using the fact that the centers x i + do not collide thanks to our assumptions, one can show that we can construct two divergence-free velocity fields E N andẼ N such that : and there exists a positive constant independent of N such that This construction yields the following result Proposition 4.1. µ N satisfies weakly the transport equation: We can prove now Theorem 0.1.

4.2.
proof of Theorem 0.1. The proof is a corollary of Proposition 4.1. Indeed, we want to show that for all ψ ∈ C ∞ c (R 3 ) we have: which is obtained directly by passing through the limit in each term of formula (46). Indeed we recall that we have the following estimates:

Proof of theorem 0.2 and 0.3
This section is devoted to the proof of Theorem 0.2 and 0.3. The Lipschitz-like estimates proved in Proposition B.3 suggests a correlation between the vectors along the line of centers ξ i and the centers x i + . In this section, we show in particular that this correlation is well propagated in time.

Derivation of the transport-Stokes equation.
We assume now that there exists a lipschitz function F 0 such that In order to propagate this correlation we search for a function F N (t, ·) ∈ W 1,∞ (R 3 ) such that for all t ∈ [0, T ] we have According to the ODE satisfied by ξ i , see (44), F N must satisfy the following equation The following proposition shows the existence and uniqueness of F N .
Proposition 5.1. There exists T >0 such that for all N ∈ N * , there exists a unique (local) solution F N ∈ L ∞ (0, T ; W 1,∞ (R 3 )) of the following equation Proof. The idea is to apply a fixed-point argument. We define the mapping A which associates to any F ∈ L ∞ (0, T ; W 1,∞ (R 3 )) the unique solution A(F ) =F to the transport equation We define X N as the characteristic flow satisfying : N (s, t, x)).
X N (t, t, x) = x. The Lipschitz property of A −1 , F , K N ρ N and E N ensures the existence, uniqueness and regularity of such a flow, see Proposition B.1 and formula (45). Moreover, direct estimates show that for all 0 ≤ s ≤ t: s)). Hence, we can writê +Ẽ(s, X N (s, t, x))ds.

Direct computations yield
Gathering all the estimates and using Proposition B.1 and the uniform bounds (45), there exists some constants independent of N such that: On the other hand, given F 1 , F 2 ∈ L ∞ (0, T ; W 1,∞ (R 3 )) we set X i the associated characteristic flow and we have . The characteristic flows satisfies This yields We construct the following sequence (F k ) k∈N ⊂ L ∞ (0, T ; W 1,∞ (R 3 )) defined as For T small enough and independent of N, using estimates (51) and (52), the sequence (F k ) k is bounded in L ∞ (0, T ; W 1,∞ (R 3 )) and is a Cauchy sequence in the Banach space ). It remains to show that F = A(F ). The weak formulation of the transport equation writes Uniqueness of the fixed-point is ensured thanks to estimate (52).
Proposition 5.1 and formula (44) yield the following result Corollary 5.2. There exists a unique solution of (48) F N ∈ L ∞ (0, T ; W 1,∞ (R 3 )) such that µ N = (id, F N )#ρ N and ρ N satisfies weakly 5.2. proof of Theorem 0.2 and 0.3. In the previous part we showed the existence of a unique function F N such that: . In order to provide the limit behaviour of the system, we need to extract the limit equation satisfied by F = lim N →∞ F N and to estimate and specify the convergence. It is straightforward that the limit function F should satisfy the following equation: We begin with the proof of local existence and uniqueness of the solution to system (16).
Proof of Theorem 0.3. Let p > 3, F 0 ∈ W 2,p (R 3 ), ρ 0 ∈ W 1,p (R 3 ) having compact support. The idea is to apply a fixed-point argument. We define the operator A which associates to each u ∈ L ∞ (0, T ; W 3,p (R 3 )) the following divergence free velocity where F (u) ∈ L ∞ (0, T ; W 2,p (R 3 )) is the unique solution, see Proposition C.1, to the following equation is the unique solution, see Proposition C.2, to the transport equation On the other hand, since ρ(t, ·) ∈ L p (R 3 ) and is compactly supported, see Remark C.1, we have in particular ρ(t, ·) ∈ L q 1 (R 3 ) ∩ L q 2 (R 3 ) with We apply again [9, Theorem IV.2.1] for q = q 1 (resp. q = q 2 ) to get ∇A(u) ∈ L p (R 3 ) (resp. A(u) ∈ L p (R 3 )) and we have according to [ Hence, since q 1 , q 2 < 3 < p, Holder's inequality yields where sup according to Remark C.1

Finally we have
We recall the following bounds, see Proposition C.2 and Proposition C.1 with C = C( F (u) L ∞ (0,T ;W 2,p (R 3 )) , u L ∞ (0,T ;W 3,p (R 3 )) ). According to Proposition C.1, for a small time interval we have for a fixed λ > 1 On the other hand, gathering the stability estimates of Proposition C.2 and Proposition C.1 and (57) we get for u i ∈ W 3,p (R 3 ), i = 1, 2 We consider the following sequence We set F k := A(u k ), ρ k := ρ(u k ). Previous estimates show that the sequences (u k ) k∈N , (F k ) k∈N , (ρ k ) k∈N are uniformly bounded in L ∞ (0, T ; W 3,p (R 3 )), L ∞ (0, T ; W 2,p (R 3 )), L ∞ (0, T ; W 1,p (R 3 )), respectively, and are Cauchy sequences in L ∞ (0, T ; W 2,p (R 3 )), L ∞ (0, T ; W 1,p (R 3 )), L ∞ (0, T ; L p (R 3 )), respectively for T small enough. Consequently, there exists (u, F, ρ) such that This allows to pass through the limit in the weak formulations of u k and ρ k . In addition, we use the fact that ∇F k converges weakly-* in L ∞ (0, T ; L ∞ (R 3 )) in order to pass through the limit in the weak formulation of F k . Hence, the triplet (u, ρ, F ) satisfies equation (16). We recover the regularity of each term using the a priori bounds. Uniqueness is a consequence of the previous stability estimates.

Proof of Theorem 0.2.
Proof of Theorem 0.2. We recall that W ∞ (ρ N , ρ) → 0 according to (9). We want to show that the triplet (ρ N , F N , K N ρ N ) converges to (ρ, F, Kρ) the unique solution of equation (16). From Proposition B.2 and using the same arguments as in Proposition C.1 we have where W ∞ (s) := W ∞ (ρ N (s, ·), ρ(s, ·)). Hence F N converges to F in L ∞ (0, T ; L ∞ (R 3 )) and K N ρ N converges to Kρ in L ∞ (0, T ; W 1,∞ (R 3 )) if the Wasserstein distance is preserved in finite time. This allows us to pass through the limit in the weak formulation of Appendix A. Some preliminary estimates This section is devoted to the proof of the following lemma which is analogous to [ Proof. We introduce a radial truncation function χ such that χ = 0 on B(0, 1/2) and χ = 1 on c B(0, 3/4). We have for all k ≥ 0: Recall that for all y ∈ B(x, 3W ∞ ) such that |x − T (y)| ≤ d min /2 we have χΦ(x − T (y)) = 0. Hence in all cases we have the following bound for all y ∈ B(x, 3W ∞ ): this yields the following bound Analogously we obtain a similar bound for ∇K N . We focus now on the bound for ∇ 2 K N ρ N . We have We use the same estimates as before to bound the first term by ρ ∞ Using an integration by parts for the first term in the right hand side of (60) we get ∇Φ(x − y)∇ρ(y)dy |∇Φ(x − y)| |ρ(y)|dσ(y) , Finally, for the second term in the right hand side of (60) we have The following convergence estimates are used in the proof of Theorem 0.2.
Proposition B.2 (Convergence estimates). The following estimates hold true: Proof. We use in the proof the shortcut W ∞ := W ∞ (ρ N , ρ). Let x ∈ R 3 , we have We split the integral into two disjoint domains J := {y ∈ supp ρ , |x − y| ≤ 3W ∞ } and its complementary. Note that on J, according to the definition of the truncation function χ, we have χΦ(x − T (y)) = 0 for all y ∈ J such that |x − T (y)| ≤ d min 4 . We can then bound directly the first integral as follows

Direct computations yields
We focus now on the remaining term, note that for all y ∈ c J := c B(x, 3W ∞ ) we have |x − T (y)| ≥ |x − y| − |T (y) − y| ≥ 2W ∞ ≥ d min , which yields that χΦ(x−T (y)) = Φ(x−T (y)) on c J. Moreover, we have |x−T (y)| ≥ 1 2 |x−y| on c J. We have then In the last line we use the fact that 1 |x−y| 2 is integrable on c B(x, 3W ∞ ). The proof for the second estimate is analogous to the first one. The main difference occurs for the last estimate where the log term appears. This is due to the fact that we integrate 1 |x−y| 3 on c B(x, 3W ∞ ).
We present now an estimate for the conservation of the particle configuration. This estimate combined with Proposition B.1 shows that the dilution regime is conserved provided that we have a control on the infinite Wasserstein distance. Proposition B.3. For all 1 ≤ i ≤ N and j = i we have We remark that the conservation of the infinite Wasserstein distance, which is initially of order 1 N 1/3 , ensures the control of the particle distance. Unfortunately, due to the log term appearing in Proposition B.2 we are not able to prove the conservation in time of the infinite Wasserstein distance.
Appendix C. Existence, uniqueness and some stability properties In this section we present some existence, uniqueness and stability estimates.
Multiplying by |D α F | p−1 and integrating by parts the second term using the fact that div(u) = 0, we get Since F 1,∞ F 2,p , u 1,∞ F 2,p , we get up to a constant depending on A −1 2,∞ d dt D α F p p D α F p p ( F 2,p + u 3,p ) + D α F p−1 p F 2,p ( F 2,p + u 3,p ) .
Applying Young's inequality and summing over α = 0, 1, 2 we get F L ∞ (0,T ;W 2,p (R 3 )) F 0 2,p e C(p, F 2,p , u 3,p , A −1 2,∞ )T , which shows that F ∈ L ∞ (0, T ; W 2,p (R 3 )) for a finite time T > 0. Now consider two divergence free velocity fields u 1 , u 2 ∈ L ∞ (0, T ; W 3,p (R 3 )) and denote by F i the solution to (61). We have Multiplying by |F 1 − F 2 | p−1 and integrating by parts the second term in the left hand side using the divergence free property of u, we get For the derivative we have ∂ t (∇F 1 − ∇F 2 ) + ∇(∇F 1 − ∇F 2 )(A −1 (F )κg + u 1 ) Using the same estimates as previously, we obtain d dt where C 1 , C 2 depend on A −1 2,∞ , u i 3,p , F i 2,p . We conclude by integrating with respect to time and apply Gronwall's inequality.