On the correction to Einstein's formula for the effective viscosity

This paper is a follow-up of article [6], on the derivation of accurate effective models for viscous dilute suspensions. The goal is to identify an effective Stokes equation providing a $o(\lambda^2)$ approximation of the exact fluid-particle system, with $\lambda$ the solid volume fraction of the particles. This means that we look for an improvement of Einstein's formula for the effective viscosity in the form $\mu_{eff}(x) = \mu + \frac{5}{2} \mu \rho(x) \lambda + \mu_2(x) \lambda^2$. Under a separation assumption on the particles, we proved in [6] that if a $o(\lambda)^2$ Stokes effective approximation exists, the correction $\mu_2$ is necessarily given by a mean field limit, that can then be studied and computed under further assumptions on the particle configurations. Roughly, we go here from the conditional result of [6] to an unconditional result: we show that such a $o(\lambda^2)$ Stokes approximation indeed exists, as soon as the mean field limit exists. This includes the case of periodic and random stationary particle configurations.


Introduction
We consider a suspension of spherical particles, modelled by a collection of balls B i = B(x i , r n ), 1 ≤ i ≤ n, all included in a fixed compact of R 3 . The centers x i of the balls may (and will) depend on n, but we omit it in the notations. We consider a regime where n is large, and r n ∼ n −1/3 . More precisely, we assume for simplicity that λ := n 4 3 πr 3 n is independent of n. We also assume that the balls occupy a volume of size 1, in the sense that where ρ is a bounded density with support O for a smooth bounded domain O such that |O| = 1. In particular, λ can be interpreted as the solid volume fraction. The suspension is immersed in a viscous fluid. We consider particles light enough so that neglecting inertia of the fluid and the particles is reasonable. Setting Ω n := R 3 \ i B i ,    −div (2µD(u n ) − Ip n ) = g n , on Ω n , div (u n ) = 0, on Ω n , where g n ∈ L 2 (R 3 ) ∩ L 6/5 (R 3 ) models some forcing. The constant vectors u i and ω i are the translation velocity and rotation vector of ball B i . They are unknowns, associated to the newtonian dynamics of the particles: in the absence of inertia, relations are of the form (2) where σ(u, p) = 2D(u) − pI is the newtonian stress tensor, and ν is the unit normal vector pointing outward. These relations correspond to prescribing the force and the torque on each particle. One further assumes decay of u n at infinity, which will be encoded in the functional setting.
As n → +∞, one may expect that some averaging takes place. The hope is to replace the fluid-particle system above by a Stokes equation, with a viscosity coefficient µ ef f = µ ef f (x), different from µ in the domain O, reflecting there the rigidity of the particles. Convergence to such a Stokes equation can indeed be shown through homogenization techniques, if one further assumes periodicity or stationarity assumptions on the distributions of balls: see [3], or [11] for the scalar case. Note that such homogenization results are valid for any λ, but somehow abstract, as the expression of the effective viscosity involves a corrector equation which is not much simpler than the original system. They are moreover restricted to homogeneous distributions. In the present paper, we aim at more explicit formulas for the effective viscosity in the dilute regime, namely when λ is small (but not vanishing as n goes to infinity). We want to show that for n large, the solution u n of (1) has for o(λ 2 ) approximation the solution u of −div (2[µ + µ 1 λ + µ 2 λ 2 ]D(ū) − Ip) = g, on R 3 , div (ū) = 0, on R 3 , for appropriate first and second order corrections µ 1 = µ 1 (x), µ 2 = µ 2 (x). Clearly, µ 1 and µ 2 should be non-zero only in the region O where the suspension is located. Note also that, if the distribution of the particles is anisotropic, µ 1 and µ 2 are not expected to be scalar functions.
In full generality, we look for µ 1 The search for the effective viscosity has a long history, starting from the work of Einstein [4]: he showed that if the suspension is homogeneous, and if the interaction between the particles can be neglected, a o(λ) approximation is given by µ ef f = µ + 5 2 λµ. A rigorous derivation of this formula and of inhomogeneous extensions was later provided under more or less stringent separation assumptions on the particles. We refer to [14,12,8] for periodic distributions of balls, and to [10,13], where the periodicity assumption is relaxed into a lower bound on the minimal distance: See also the most recent paper [7], where formula µ 1 (x) = 5 2 ρ(x)µ is established under mild requirements. Our concern in the present paper is related to a o(λ 2 ) effective approximation, that is beyond Einstein's formula. Such second order effective viscosity has been discussed in several papers, see [2,16,1]. However, one can observe discrepancies between the results, and very different approaches depending on the type of suspensions considered. A more global analysis was initiated by the first author and Matthieu Hillairet in the recent paper [6]. Roughly, this paper shows that under (A1), if lim sup n ||u n − u|| = o(λ 2 ), where || || is a weak norm and where u is a solution of (3) with µ 1 = 5 2 ρ(x)µ, then necessarily: Hence, if a second order effective model exists, the average ν 2 of the second order correction over the domain (that coincides with µ 2 if µ 2 does not depend on x) is given by the mean field limit (4)- (5). The second part of article [6] consists in an analysis of such mean field limit, using ideas introduced by S. Serfaty and her co-authors in the context of Coulomb gases [15]. More explicit formula are provided, notably in the periodic case.
The limitation of the results in [6] is that they hold conditionally to the existence of a second order effective model of type (3). Ideally, one would like to prove the existence of an effective model as soon as the mean field limit in (4)-(5) does exist. This necessary condition is however not enough: indeed ν 2 corresponds to an average over the whole domain O, so that it is unlikely to guarantee the existence of an effective local coefficient µ 2 = µ 2 (x). Nevertheless, as we will show, we can exhibit more local in nature mean field limits, whose existence ensures the existence of an effective model. Moreover, such limits allow to determine µ 2 , and not only its average ν 2 . We introduce It can be seen as a compactly supported distribution on R 3 x × R 3 y , with values in the space Sym Sym 3,σ (R) : We stress that M(x) is a Calderon-Zygmund kernel, hence not integrable. In particular, the last integral must be understood in a weak sense: it can be defined rigorously through the decomposition where the first integral in the decomposition exists in the usual sense, while the second one which allows to give a meaning to µ 2,n , f ⊗ g for much less regular f and g. Our theorem reads as follows.
Assume (A0)-(A1), that g n → g in L 6 5 (R 3 ), and that with µ 2,n defined in (6). Then any accumulation point u λ of u n,λ solves where R λ satisfies for all q ≥ 3 A few remarks are in order.
Remark 1.1. One has directly from estimate (9) that for any p ≤ 3 2 , where u is the solution of the effective model It follows by Sobolev imbedding that u − u λ L r loc = O(λ 7 3 ), for any r = 3p 3−p ≤ 3. Moreover, as will be seen below, (u n,λ ) n∈N is bounded inḢ 1 ∩ L 6 . Combining the last estimate with Rellich's theorem, it follows easily that Remark 1.2. The main assumption of the theorem is the convergence of µ 2,n to µ 2 (x)δ x=y . With regards to the form of µ 2,n , F , cf. (7), this convergence corresponds to the local mean field limits alluded to above. In particular, the necessary condition (4) derived in [6] corresponds to the convergence of µ 2,n , F for the special case F (x, y) = 1.
The outline of the paper is the following. After preliminary results on the Stokes system, gathered in Section 2, we adress the proof of Theorem 1.1 in Section 3. Finally, we turn in Section 4 to the discussion of assumption (A2). Roughly, we show that it is fulfilled by both periodic and random stationary particle distributions that satisfy the separation assumption (A1), and we discuss the corresponding limit µ 2 . We rely there much on article [6]. We notably show that when the particle distributions is given by an isotropic process (plus technical conditions), then µ 2 S : S = 5 2 µ|S| 2 , a result that was not given in [6].

Reminder on the Stokes problem
In this section we recall some properties regarding the Stokes equation on an exterior domain. We denote by (U , P) the Green function of the Stokes equation: U (x) := 1 8π Let A ∈ Sym 3,σ (R). We denote by (V [A] , Q [A]) the solution to the Dirichlet problem given by the explicit formula An important feature of this solution is that, as easily seen from symmetry considerations, it fulfills the extra conditions Moreover we can link V [A] to the Green function through the identity where R[A] is homogeneous of degree −4. Here, ∇U is a third rank tensor defined using Einstein summation convention by Moreover, a simple calculation yields the identity where M was defined in (5), so that We also introduce the extensions Direct computation shows that where s η is the surface measure on the sphere of radius η.
We finish this part by recalling a classical estimate for the Stokes equation. Let w ∈ W 1,2 (∪B i ), divergence-free. We consider the unique solution (v, q) satisfying with the following conditions Using an integration by parts, one can show that v is a minimizer of Combining this minimizing property with [13, Lemma 4.4] we have Proof. The first equality is well-known: it follows from identity ∆v = 2div (D(v)) and integration by parts. For the inequality, by the minimizing property of v, it is enough to construct a divergence-free velocity field u that matches the condition D(u) = D(w) on ∪B i and satisfies the same inequality. Classical considerations about the Bogovskii operator ensure the existence of fields where using Poincaré-Wirtinger inequality we get for all i, with C, C ′ independent of n by scaling considerations. See [9, Lemma 18] for details. We then take u = i u i . Since the balls B(x i , 2r n ) are disjoint by (A1) for λ small enough, u satisfies Moreover, adding a proper rigid vector field to w on each B i , which does not change D(w) on each B i , we can always assume that We conclude by applying [13,Lemma 4.4].
3 Proof of Theorem 1.1 By linearity of the Stokes equation, we can restrict to the case µ = 1. Let q ≥ 3. The goal is to show that any accumulation point u λ of u n,λ satisfies a system of type (8) . For any such φ, we consider φ n satisfying with the following conditions Note that φ n depends implicitly on λ. Similarly, we shall note u n instead of u n,λ for short.
We remind that in the second line of the above computations, the unit normal vector ν is pointing outward the balls. Using equations (24) and (25) we have We get the following relation for all φ using the fact that D(u n ) = 0 on B i : By a simple energy estimate, u n also satisfies where the last inequality comes from the Sobolev embedding and the boundedness of (g n ) n∈N in L 6 5 . Hence, (u n ) n∈N is bounded inḢ 1 ∩ L 6 . Denoting u λ = lim k u n k an accumulation point, we deduce from (26) with n = n k that Note that this limit exists because all other terms in (26) converge when n = n k , k → +∞. Moreover, it is clearly linear in φ as φ n is linear in φ. Furthermore, testing φ n − φ in (24), similar integrations by parts lead to Cauchy Schwartz inequality implies that φ n is bounded uniformly in n inḢ 1 , hence in L 6 . Eventually, as g n → g strongly in L In order to obtain (27), we shall write (φ n , q n ) = (φ 1 n + φ 2 n , q 1 n + q 2 n ), where φ 1 n is a (somehow natural) approximation of φ 1 n and where φ 2 n is a remainder. Namely, we look for an approximation φ 1 n of the form where U , V are defined in (11), (20). The rough idea behind this approximation is that the first term at the right-hand side should take care of the source term in (24), while the second term should take care of the boundary conditions at the balls B i . In particular, the field φ R 3 := ∇U * 2(λ 5 2 |O|ρ + λ 2 µ 2 )D(φ) solves the Stokes equation By looking to (24), it is tempting to take A i = Dφ i , where as φ should be close to this value on the small ball B i . However, this approximation is not accurate enough, and would only allow to recover Einstein's formula. To go beyond, one must account for two extra contributions. The first one is the trace left at the balls by φ R 3 . More precisely, it will be enough to correct the trace of the O(λ) term in φ R 3 . The second contribution is the trace left by all φ j,n , j = i, on ball B i , which corresponds to binary interactions between particles. Again, it will be enough to account for the least decaying term in φ j,n , cf. decomposition (16). This leads to the following definition: where, using relation (18) while, still using (18) Direct computations using formula (21) yield Using the definition of S 1 i and S 2 i we find Moreover, thanks to property (15), it is easily seen that It follows that (φ 2 n , q 2 n ) satisfies div (σ(φ 2 n , q 2 n )) = 0 on Ω n , with boundary condition D(φ 2 n ) = D(ψ 2 n +ψ 2 n ) on ∪ B i , where ψ 2 n andψ 2 n are defined by: We aim now at estimating both terms g φ 1 n and g φ 2 n .

Estimate of φ 2 n
Most of this paragraph is dedicated to the derivation of Before we prove this proposition, we show that it implies Note the extra factor λ 1/2 compared to Proposition 3.1, crucial to obtain a o(λ 2 ) error. Note also that for q ≥ 3, this term is bounded by λ Proof. We introduce the solution u g of the Stokes equation As g ∈ L 3+ε , u g ∈ W 2,3+ε loc , so that D(u g ) is continuous. Integrations by parts yield for any constant vectors u i g , ω i g , 1 ≤ i ≤ n, by the last two relations in (31). As We can apply classical considerations on the Bogovskii operator [5]: for any 1 ≤ i ≤ n, there , by a proper choice of u i g and ω i g , we can ensure the Korn inequality: where the constant C in the last inequality can be taken independent of i and n by translation and scaling arguments. Extending U i g by zero, and denoting U g = U i g , we have for d n > 4r n (which is implied by (A1) for λ small enough): Back to our calculation, we find By using (33) and Cauchy-Schwartz inequality, we end up with so that combining with Proposition 3.1 yields the result.
We now turn to the proof of Proposition 3.1. Proposition 2.1 implies As regards ψ 2 n , we compute As regardsψ 2 n , we use the identities (18) and (19) to write for all For E 4 we have for all q ∈ [2, ∞) because M is a Calderon-Zygmund operator. The estimate of E 5 is more difficult. We shall rely on (A1), and notably make a crucial use of the following lemma, taken from [6] :

Lemma 3.3. [6, Lemma 2.4] *
Under assumption (A1), for all q ∈ (1, ∞), there exists C > 0, such that for all A 1 , . . . , A n in Sym 3,σ (R), In particular, this lemma can be applied to matrices A j = S 2 j , cf. (29). We find where in this computation and all computations below, the constant C may change from line to line. Moreover, by (A0), so that lim sup Besides Lemma 3.3, we shall also make use of the following easy generalization of Young's inequality: We now introduce y i := x i n −1/3 so that |y i − y j | ≥ 1 2 (c + |y i − y j |) with c the constant appearing in (A1). Using decomposition (19) and the homogeneity of each term in this decomposition, we obtain: for all q ≥ 2, q * Only the first inequality is stated in [, Lemma 2.4], but a look at the proof shows that it follows from the second one.
We can then apply Lemma 3.3 to the first term, and apply (39) to the second term (together with the fact that i 1 c+|y i −y j | 4 is uniformly bounded in j). We end up with for any q ≥ 2, where the last bound comes from a Hölder inequality. It remains to bound where we used again the L q continuity of the convolution with M. Combining this inequality with (38), and injecting in the bound for E 5 , we get that for all q ≥ 2, lim sup We now turn to E 3 . We use the fact that D(R[A])(x) = O(|A||x| −5 ). We find where the last line follows from (39). Using again (37), we get eventually We recall now the expressions of E 1 and E 2 on B i : We claim first that lim sup To prove (42), the idea is to regularize D(φ)ρ using a convolution with a mollifier χ η . Denoting || || 0,µ the Hölder semi-norm of exponent µ, we get for all q ≥ 1. Hence we get using the fact that M is a Zygmund-Calderon kernel: for all q ≥ 2, any µ ∈ (0, 1), which vanishes when taking the lim sup in n for fixed η and then taking the limit η → 0.
Eventually, we bound E 2 by: where the last line comes from (39). We end up with lim sup

Estimate of φ 1 n
The goal of this paragraph is to show that The main estimate (27) is then an easy consequence of Propositions 3.2 and 3.4, which, as explained before (27), concludes the proof of Theorem 1.1.
The derivation of Proposition 3.4 will rely strongly on our assumption (A2). At first, let us remind that φ n is bounded uniformly in n inḢ 1 , and so is φ 2 n by Proposition (3.1). It follows that that φ 1 n is bounded uniformly in n inḢ 1 . Hence, by a density argument, it is enough to establish the inequality in Proposition (3.4) for a C ∞ c field g, as long as the constant C at the right-hand side only involves ||g|| L 3+ε . From now on, we assume g to be smooth and compactly supported.
We introduce again the solution (u g , p g ) of (32). As g is smooth, so is u g . Moreover, by elliptic regularity and Sobolev imbedding, we have for a compact K containing all the balls, With (30), we compute: We set µ 2,n := 75 16π See (6) for comparison. We have with the notation that for M a distribution over R 3 × R 3 with values in Sym(Sym 3,σ (R)): From assumption (A0), it easily follows that lim n R 1 = 0 Then, we write The first term at the right-hand side goes to zero as n → +∞, by assumption (A2). For the remainder, we decompose it as The first term is bounded by Note that the constant C may be chosen so that it depends on g only through g L 3+ε , using (44). The second term is bounded by where the fourth, resp. fifth inequality is a consequence of Lemma 3.3, resp. (A0). Hence, lim sup n |R 2,2 | = 0.
As regards the last term, we use the fact that, under (A0),

Discussion of the mean-field limit
We come back in this final section to the main assumption (A2): we discuss examples of particle distributions (x i ) for which µ 2,n converges as in (A2), and discuss how to compute the limit µ 2 , that is the O(λ 2 ) correction to the effective viscosity. This discussion is closely related to the analysis performed in [6], and relies on results established there.
At first, one can notice that the sequence of compactly supported distributions µ 2,n defined in (6) is bounded with respect to n. More precisely: Lemma 4.1. There exists C > 0, such that for any smooth F = F (x, y), for any n, where K is any fixed convex compact containing x 1 , . . . x n for all n.
Proof. We need to bound We introduce again y i := x i n −1/3 so that |y i − y j | ≥ 1 2 (c + |y i − y j |) with c the constant appearing in (A1). We end up with using that for each i, under assumption (A1), To bound I 1,n , we introduce s n := c/(4n 1 3 ), where c is again the constant appearing in (A1). We use the splitting: M(x − y)dxdyF (x j , x j ) =: J 1,n + J 2,n + J 3,n We find Similarly, we find |J 2,n | ≤ C. Eventually, let χ is a truncation which is zero on B(0, 1) and one outside B(0, 2). The last term can be written using the Calderon-Zygmund theorem. This concludes the proof.
Combining the previous lemma with the density of span {ϕ ⊗φ, ϕ,φ ∈ C ∞ (K)} in C ∞ (K × K), we deduce that for µ 2 ∈ L ∞ R 3 , Sym Sym 3,σ (R) , assumption (A2) is satisfied if and only if for all smooth ϕ,φ: or equivalently if and only if for all S ∈ Sym 3,σ (R), for all smooth ϕ,φ: cf. (5). As g S is even, this is the same as : for all S ∈ Sym 3,σ (R), for all smooth ϕ (47) Of course, the existence of a non-trivial limit at the right-hand side of (47) is directly related to the singularity of g s at zero. If g s were smooth or not too singular, one could show that the limit at the right-hand side is zero.
Our first goal is to understand configurations of points (x i ) for which the limit as n → +∞ of 1 exists. This is closely related to the study in [6], in which the special case ϕ = 1 is investigated. This study relies on ideas from homogenization theory, together with the notion of renormalized energy, described in [15] in the context of Coulomb gases. We shall explain how the approach in [6] adapts to a general ϕ. One proceeds in several steps.
Step 1. Expression of the mean field as a (renormalized) energy.
Proceeding as in [6,Lemma 3.1], one can restrict to the case where Under this assumption, proceeding as in [6,Proposition 3.3], one shows that (47), and thus (A2), amounts to: for all smooth ϕ, for all S ∈ Sym 3,σ (R), where To investigate the convergence of W n [ϕ], the main idea is then to associate to this quadratic quantity an energy. This idea is based on the following ) For any f ∈ L 2 (R 3 ), where u S is the solution of the Stokes equation Considering Definition (49), it is tempting to replace f by ϕ(ρ n − ρ) in the lemma, which would yield the formal identity: where h n solves the equation The solution of this equation decaying to zero is explicit, and given by where G S (x) = (S∇) · U (x) = − 3 8π (Sx · x) x |x| 5 , while the corresponding pressure field is |x| 5 . However, this formal identity is not justified, as both sides are infinite: the left-hand side is infinite due to the singularity of the kernel at the diagonal (which is excluded in Definition (49)), and the right-hand side is infinite because the solution h n of (51) is not inḢ 1 (R 3 ), as ρ n is an atomic measure. Still, one can interpret W n [ϕ] in terms of a so-called renormalized energy, intensively studied in the context of Coulomb gases and Ginzburg-Landau systems. See [15] and the references therein for more. One must at first regularize the singular measure. More precisely, following [6], we consider for all η > 0, the field G η S satisfying Note that the condition G η S = G S at |x| = η contained in the first condition can be seen as a Dirichlet condition for the Stokes problem satisfied by G η S inside {|x| < η}. Let us stress that G η S belongs to H 1 loc (R 3 ), and can be seen as a regularization of G S . Accordingly, we set As shown in [6, Lemma 3.5], it solves the modified Stokes equation where Ψ η := 3 A main advantage of this regularization is that h n − h η n is supported in ∪ i B(x i , η). The main result, that is a straightforward adaptation of [6, Propositions 3.7 and 3.8] (dedicated to the special case ϕ = 1) is given by More precisely, for η < c 2 n − 1 3 , where c is the constant appearing in (A1), one has The first formula in this proposition is the rigorous translation of the false identity (50): it connects the mean field W n [ϕ] to an energy. However, two differences are in order: first, the energy of h n being infinite, it is replaced by the energy of the regularization h η n . Second, and more importantly, one must substract the term 1 where F is the mean of F (in the periodic or random sense). Roughly, one has the correspondence The point is that for a model of type (56), the behaviour of the solution h ε is well-understood. For instance, in the periodic case, one has 1. One has and for all L > 0, G S,L is the periodic field with zero mean solution of 2. In the special case of a simple cubic lattice (m = 1), the formula simplifies into , and a ≈ −0, 04655.
In the ergodic stationary case, under the almost sure assumption (A1), a more explicit expression can be obtained in terms of the two-point correlation function ρ 2 (x, y) = r(x − y). We remind that a point process in R 3 has a two-point correlation function ρ 2 ∈ L 1 loc (R 3 × R 3 ) if for all bounded K and smooth function F  One can actually push further this calculation in the case of an isotropic process. The fact that r is radial expresses an isotropy of the point process. The fact that it is constant and equal to m 2 at infinity corresponds to a natural decorrelation at large distances. Note that all assumptions of the proposition are satisfied by the usual hardcore Poisson processes. One could actually relax the hypotheses on r, just assuming fast enough convergence at infinity. (S∇) · G S,1 (z − z ′ ) r(L(z − z ′ ))dzdz ′ = m 2 5 |S| 2 .
We take M large enough so that r(x) = m 2 for |x| ≥ M . We decompose We then perform an integration by parts in variable z. Note that for all z ′ ∈ K 2M/L,1 , the boundary of the ball B(z ′ , M/L) is disjoint from the boundary of K 0,1 . Taking into account that (S∇) · G S,1 (z − z ′ ) is Z d -periodic, the boundary term at ∂K 0,1 vanishes, and eventually Using the explicit expression for G S , we find I L → 3 8π ∂B 1 (Sn · n) 2 m 2 = m 2 5 |S| 2 thanks to the identity ∂B 1 n i n j n k n l = 4π 15 (∂ ij ∂ kl + ∂ ik ∂ jl + ∂ il ∂ jk ). Eventually, for the last term, we find that J L = K 2M/L,1 B(z ′ ,M/L) (S∇) · G S (z − z ′ ) r(L(z − z ′ ))dzdz ′ + o(1) As r is radial and as the mean of (S∇) · G S over spheres is zero, the integral in z at the right-hand side vanishes identically. This concludes the proof.