Analysis of the viscosity of dilute suspensions beyond Einstein's formula

TL;DR

This paper mathematically refines Einstein's formula for suspension viscosity by deriving the second-order correction involving particle interactions, applicable to large systems with specific spatial distributions.

Contribution

It introduces a rigorous derivation of the second-order correction to Einstein's viscosity formula for dilute suspensions, including the limit as the number of particles grows large.

Findings

Derived explicit formulas for the $O(\,\phi^2)$ correction to viscosity.

Extended analysis to infinite particle systems with periodic or ergodic distributions.

Validated the correction's applicability beyond classical Einstein approximation.

Abstract

We provide a mathematical analysis of the effective viscosity of suspensions of spherical particles in a Stokes flow, at low solid volume fraction $\phi$. Our objective is to go beyond the Einstein's approximation $\mu_{eff}=(1+\frac{5}{2}\phi)\mu$. Assuming a lower bound on the minimal distance between the $N$ particles, we are able to identify the $O(\phi^2)$ correction to the effective viscosity, which involves pairwise particle interactions. Applying the methodology developped over the last years on Coulomb gases, we are able to tackle the limit $N \rightarrow +\infty$ of the $O(\phi^2)$-correction, and provide explicit formula for this limit when the particles centers can be described by either periodic or stationary ergodic point processes.

Equations

$$\{\begin{array}{rl}-\mu \Delta u+\nabla p& =0,x\in \mathcal{F},\\ divu& =0,x\in \mathcal{F},\\ u{|}_{B_i}& =u_i+{\omega }_i\times (x-x_i),\end{array}$$

$${\int }_{\partial B_i}{\sigma }_{\mu }(u,p)nds$$

$${\delta }_N=\frac{1}{N}\sum _{i=1}^N{\delta }_{x_i}\underset{N\to +\infty }{\to }f(x)dx\text{weakly,}$$

$$support(f)=\overline{\mathcal{O}},\mathcal{O}\text{ smooth, bounded and open},f{|}_{\mathcal{O}}\in C^1(\overline{\mathcal{O}}).$$

$$\varphi =\frac{4}{3}N\pi a^3/|\mathcal{O}|$$

$$\underset{i\ne j}{\min }|x_i-x_j|\ge c{N}^{-1/3}\text{ for some constant c>0 independent of N}.$$

$$\{\begin{array}{rl}-\mu \Delta u+\nabla p& =0,divu=0,x\in {\mathbb{R}}^3\setminus \mathcal{O},\\ -2div({\mu }^{\prime }D(u^{\prime }))+\nabla p^{\prime }& =0,div{u}^{\prime }=0,x\in \mathcal{O},\end{array}$$

$$u=u^{\prime }\text{at }\partial \mathcal{O},{\sigma }_{\mu }(u,p)n={\sigma }_{{\mu }^{\prime }}(u^{\prime },p^{\prime })n\text{at }\partial \mathcal{O}.$$

$${\mu }_{ijkl}^{\prime }={\mu }_{jikl}^{\prime }={\mu }_{jilk}^{\prime }={\mu }_{lkji}^{\prime },$$

$${\mu }^{\prime }=\mu \mathrm{Id}+\varphi {\mu }_1+{\varphi }^2{\mu }_2+\cdots +{\varphi }^k{\mu }_k+o({\varphi }^k).$$

$${\nu }_i:=\frac{1}{|\mathcal{O}|}{\int }_{\mathcal{O}}{\mu }_i(x)dx$$

$$\underset{|x|\to +\infty }{\lim }u(x)-Sx=0$$

$${\mu }^{\prime }=\mu +\frac{5}{2}\varphi \mu +o(\varphi ).$$

$${\mu }_1=\frac{5}{2}|\mathcal{O}|f(x)\mu .$$

$${\mu }_N=\frac{\mu }{|S{|}^2}{\int }_{\mathcal{F}}|D(u_N){|}^2$$

$${\mu }_N=\mu +\frac{5}{2}\varphi \mu +O({\varphi }^{3/2}).$$

$${\mu }^{\prime }=\mu +\frac{5}{2}\varphi \mu +7.6{\varphi }^2\mu +o({\varphi }^2)$$

$${\nu }_2=\frac{1}{|\mathcal{O}|}{\int }_{\mathcal{O}}{\mu }_2(x)dx$$

$$\underset{N\to +\infty }{\mathrm{limsup}}\Vert u_N-u_{eff}{\Vert }_{H^s(U)}=o({\varphi }^2),\text{as ϕ→0},$$

$${\mathcal{V}}_N:=\frac{75|\mathcal{O}|}{16\pi }\left(\frac{1}{N^2}\sum _{i\ne j}{g}_S(x_i-x_j)-{\int }_{{\mathbb{R}}^3\times {\mathbb{R}}^3}{g}_S(x-y)f(x)f(y)dxdy\right)$$

$$g_S:=-D\left(\frac{S:(x\otimes x)x}{|x{|}^5}\right):S.$$

$$H_N=\frac{1}{N^2}\sum _{i\ne j}g(x_i-x_j)+\frac{1}{N}\sum _{i=1}^{N}V(x_i)$$

$$H(\lambda )={\int }_{{\mathbb{R}}^3\times {\mathbb{R}}^3}g(x-y)d\lambda (x)d\mu (y)+{\int }_{{\mathbb{R}}^3}V(x)d\lambda (x)$$

$${\mathcal{H}}_N={\int }_{{\mathbb{R}}^3\times {\mathbb{R}}^3\setminus \text{Diag}}g(x-y)d({\delta }_N-\lambda )(x)d({\delta }_N-\lambda )(y)$$

$$\text{”}{\int }_{{\mathbb{R}}^3\times {\mathbb{R}}^3}g(x-y)d({\delta }_N-\lambda )(x)d({\delta }_N-\lambda )(y)=\frac{1}{4\pi }{\int }_{{\mathbb{R}}^3}|\nabla h_N{|}^2\text{”}$$

$${\mathcal{H}}_N=\underset{\eta \to 0}{\lim }\left(\frac{1}{4\pi }{\int }_{{\mathbb{R}}^3}|\nabla h_N^{\eta }{|}^2-Ng(\eta )\right)$$

$$\lim_{N\rightarrow+\infty}\mathcal{V}_{N}=\frac{25L^{3}}{2M^{2}}\Bigl{(}\sum_{i\neq j}S{\nabla}\cdot G_{S,L}(a_{i}-a_{j})\>+\>KS{\nabla}\cdot(G_{S,L}-G_{S})(0)\Bigr{)},$$

$$\nu_{2}S:S=\mu\big{(}\alpha\sum_{i=1}^{3}|S_{ii}|^{2}+\beta\sum_{i\neq j}|S_{ij}|^{2}\big{)},$$

$$\underset{N}{\lim }{\mathcal{V}}_N$$

$$=\frac{25}{2}\underset{L\to +\infty }{\lim }\frac{1}{L^3}{\int }_{K_L\times K_L}S\nabla \cdot G_{S,L}(z-z^{\prime })\rho (z-z^{\prime })dzd{z}^{\prime }.$$

$$u_{eff}-u_0=\varphi u_{eff,1}+\cdots +{\varphi }^k{u}_{eff,k}+o({\varphi }^k)\text{in }{\dot{H}}^1({\mathbb{R}}^3)$$

Full text

1 Setting of the problem

Our general concern is the computation of the effective viscosity generated by a suspension of $N$ particles in a fluid flow. We consider spherical particles of small radius $a$, centered at $x_{i,N}$, with $N\geq 1$ and $1\leq i\leq N$. To lighten notations, we write $x_{i}$ instead of $x_{i,N}$, and $B_{i}=B(x_{i},a)$. We assume that the Reynolds number of the fluid flow is small, so that the fluid velocity is governed by the Stokes equation. Moreover, the particles are assumed to be force- and torque-free. If $\mathcal{F}={\mathbb{R}}^{3}\setminus(\cup_{i}B_{i})$ is the fluid domain, governing equations are

[TABLE]

where $\mu$ is the kinematic viscosity, while the constant vectors $u_{i}$ and $\omega_{i}$ are Lagrange multipliers associated to the constraints

[TABLE]

Here, $\sigma_{\mu}(u,p):=2\mu D(u)-pI$ is the usual Cauchy stress tensor. The boundary condition at infinity will be specified later on.

We are interested in a situation where the number of particles is large, $N\gg 1$. We want to understand the additional viscosity created by the particles. Ideally, our goal is to replace the viscosity coefficient $\mu$ in (1.1) by an effective viscosity tensor $\mu^{\prime}$ that would encode the average effect induced by the particles. Note that such replacement can only make sense in the flow region $\mathcal{O}$ in which the particles are distributed in a dense way. For instance, a finite number of isolated particles will not contribute to the effective viscosity, and should not be taken into account in $\mathcal{O}$. The selection of the flow region is formalized through the following hypothesis on the empirical measure:

[TABLE]

Note that we do not ask for regularity of the limit density $f$ over ${\mathbb{R}}^{3}$, but only in restriction to $\mathcal{O}$. Hence, our assumption covers the important case $f=\frac{1}{|\mathcal{O}|}1_{\mathcal{O}}$.

We investigate the classical regime of dilute suspensions, in which the solid volume fraction

[TABLE]

is small, but independent of $N$. Besides (H1), we make the separation hypothesis

[TABLE]

Let us stress that (H2) is compatible with (H1) only if the $L^{\infty}$ norm of $f$ is small enough (roughly less than $1/c^{3}$), which in turn forces $\mathcal{O}$ to be large enough.

Our hope is to replace a model of type (1.1) by a model of the form

[TABLE]

with the usual continuity conditions on the velocity and the stress:

[TABLE]

A priori, $\mu^{\prime}$ could be inhomogeneous (and should be if the density $f$ seen above is itself non-constant over $\mathcal{O}$). It could also be anisotropic, if the cloud of particles favours some direction. With this in mind, it is natural to look for $\mu^{\prime}=\mu^{\prime}(x)$ as a general 4-tensor, with $\sigma^{\prime}=2\mu^{\prime}D(u)$ given in coordinates by $\sigma_{ij}=\mu^{\prime}_{ijkl}D(u)_{kl}$. By standard classical considerations of mechanics, $\mu^{\prime}$ should satisfy the relations

[TABLE]

namely $\mu^{\prime}$ should define a symmetric isomorphism over the space of $3\times 3$ symmetric matrices.

As we consider a situation in which $\phi$ is small, we may expect $\mu^{\prime}$ to be a small perturbation of $\mu$, and hopefully admit an expansion in powers of $\phi$:

[TABLE]

The main mathematical questions are:

• •

Can solutions $u_{N}$ of (1.1)-(1.2) be approximated by solutions $u_{eff}=1_{{\mathbb{R}}^{3}\setminus\mathcal{O}}u+1_{\mathcal{O}}u^{\prime}$ of (1.4)-(1.5), for an appropriate choice of $\mu^{\prime}$ and an appropriate topology ?

• •

If so, does $\mu^{\prime}$ admit an expansion of type (1.6), for some $k$ ?

• •

If so, what are the values of the viscosity coefficients $\mu_{i}$, $1\leq i\leq k$ ?

Let us stress that, in most articles about the effective viscosity of suspensions, it is implicitly assumed that the first two questions have a positive answer, at least for $k=1$ or $2$. In other words, the existence of an effective model is taken for granted, and the point is then to answer the third question, or at least to determine the mean values

[TABLE]

of the viscosity coefficients. As we will see in Section 2, these mean values can be determined from the asymptotic behaviour of some integral quantities $\mathcal{I}_{N}$ as $N\rightarrow+\infty$. These integrals involve the solutions $u_{N}$ of (1.1)-(1.2) with condition at infinity

[TABLE]

where $S$ is an arbitrary symmetric trace-free matrix.

The effective viscosity problem for dilute suspensions of spherical particles has a long history, mostly focused on the first order correction created by the suspension, that is $k=1$ in (1.6). The pioneering work on this problem was due to Einstein [15], not mentioning earlier contributions on the similar conductivity problem by Maxwell [29], Clausius [11], Mossotti [32]. The celebrated Einstein’s formula,

[TABLE]

was derived under the assumption that the particles are homogeneously and isotropically distributed, and neglecting the interactions between particles. In other words, the correction $\mu_{1}=\frac{5}{2}\mu$ is obtained by summing $N$ times the contribution of one spherical particle to the effective stress. The calculation of Einstein will be seen in Section 2. It was later extended to the case of an inhomogeneous suspension by Almog and Brenner [1, page 16], who found

[TABLE]

The mathematical justification of formula (1.9) came much later. As far as we know, the first step in this direction was due to Sanchez-Palencia [38] and Levy and Sanchez-Palencia [28], who recovered Einstein’s formula from homogenization techniques, when the suspension is periodically distributed in a bounded domain. Another justification, based on variational principles, is due to Haines and Mazzucato [19]. They also consider a periodic array of spherical particles in a bounded domain $\Omega$, and define the viscosity coefficient of the suspension in terms of the energy dissipation rate:

[TABLE]

where $u_{N}$ is the solution of (1.1)-(1.2)-(1.8), replacing ${\mathbb{R}}^{3}$ by $\Omega$. Their main result is that

[TABLE]

For preliminary results in the same spirit, see Keller-Rubenfeld [27]. Eventually, a recent work [21] by the second author and Di Wu shows the validity of Einstein’s formula under general assumptions of type (H1)-(H2). See also [33] for a similar recent result.

Our goal in the present paper is to go beyond this famous formula, and to study the second order correction to the effective viscosity, that is $k=2$ in (1.6). Results on this problem have split so far into two settings: periodic distributions, and random distributions of spheres. Many different formula have emerged in the literature, after analytical, numerical and experimental studies. In the periodic case, one can refer to the works [37, 42, 34, 2], or to the more recent work [2], dedicated to the case of spherical inclusions of another Stokes fluid with viscosity $\tilde{\mu}\neq\mu$. Still, in the simple case of a primitive cubic lattice, the expressions for the second order correction differ. In the random case, the most reknowned analysis is due to Batchelor and Green [5], who consider a homogeneous and stationary distribution of spheres, and express the correction $\mu_{2}$ as an ensemble average that involves the $N$-point correlation function of the process. As pointed out by Batchelor and Green, the natural idea when investigating the effective viscosity up to $O(\phi^{2})$ is to replace the $N$-point correlation function by the $2$-point correlation function, but this leads to a divergent integral. To overcome this difficulty, Batchelor and Green develop what they call a renormalization technique, that was developed earlier by Batchelor to determine the sedimentation speed of a dilute suspension. After further analysis of the expression of the two-point correlation function of spheres in a Stokes flow [6], completed by numerical computations, they claim that under a pure strain, the particles induce a viscosity of the form

[TABLE]

Although the result of Batchelor and Green is generally accepted by the fluid mechanics community, the lack of clarity about their renormalization technique has led to debates, see [22, 35, 1].

One main objective in the present paper is to give a rigorous and global mathematical framework for the computation of

[TABLE]

leading to explicit formula in periodic and stationary random settings. We will adopt the point of view of the studies mentioned before: we will assume the validity of an effective model of type (1.4)-(1.5)-(1.6) with $k=2$, and will identify the averaged coefficient $\nu_{2}$.

More precisely, our analysis divides into two parts. The first part, carried in Section 2, has as its main consequence the following

Theorem 1.1**.**

Let $(x_{i})_{1\leq i\leq N}$ a family of points supported in a fixed compact of ${\mathbb{R}}^{3}$, and satisfying (H1)-(H2). For any trace-free symmetric matrix $S$ and any $\phi>0$, let $u_{N}$, resp. $u_{eff}$, the solution of (1.1)-(1.2)-(1.8) with the radius $a$ of the balls defined through (1.3), resp. the solution of (1.4)-(1.5)-(1.8) where $\mu^{\prime}$ obeys (1.6) with $k=2$, $\mu_{1}$ being given in (1.10).

If $u_{N}-u_{eff}=o(\phi^{2})$ in $H^{-\infty}_{loc}({\mathbb{R}}^{3})$, meaning that for all bounded open set $U$, there exists $s\in{\mathbb{R}}$ such that

[TABLE]

then, necessarily, the coefficient $\nu_{2}$ defined in (1.12) satisfies $\>\nu_{2}S:S=\mu\lim_{N\rightarrow+\infty}\mathcal{V}_{N}$ where $\nu_{2}$ was defined in (1.12), and

[TABLE]

with the Calderón-Zygmund kernel

[TABLE]

Roughly, this theorem states that if there is an effective model at order $\phi^{2}$, the mean quadratic correction $\nu_{2}$ is given by the limit of $\mathcal{V}_{N}$, defined in (1.13). Note that the integral at the right-hand side of (1.13) is well-defined: $f\in L^{2}({\mathbb{R}})$ and $f\rightarrow g_{S}\star f$ is a Calderón-Zygmund operator, therefore continuous on $L^{2}({\mathbb{R}}^{3})$. We insist that our result is an if theorem: the limit of (1.13) does not necessarily exist for any configuration of particles $x_{i}=x_{i,N}$ satisfying (H1)-(H2). In particular, it is not clear that an effective model at order $\phi^{2}$ is available for all such configurations.

Still, the second part of our analysis shows that for points associated to stationary random processes (including periodic patterns or Poisson hard core processes), the limit of the functional does exist, and is given by an explicit formula. We shall leave for later investigation the problem of approximating $u_{N}$ by $u_{eff}$ when the limit of $\mathcal{V}_{N}$ exists.

Our study of functional (1.13) is detailed in Sections 3 to 5. It borrows a lot from the mathematical analysis of Coulomb gases, as developped over the last years by Sylvia Serfaty and her coauthors [40, 36, 9]. Although our paper is self-contained, we find useful to give a brief account of this analysis here. As explained in the lecture notes [41], one of its main goals is to understand what configurations of points minimize Coulomb energies of the form

[TABLE]

where $g(x)=\frac{1}{|x|}$ is a repulsive potential of Coulomb type, and $V$ is typically a confining potential. It is well-known, see [41, chapter 2], that under suitable assumptions on $V$, the sequence of functionals $H_{N}$ (seen as a functionals over probability measures by extension by $+\infty$ outside the set of empirical measures) $\Gamma$-converges to the functional

[TABLE]

Hence, the empirical measure $\delta_{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{i}}$ associated to the minimizer $(x_{1},\dots,x_{N})$ of $H_{N}$ converges weakly to the minimizer $\lambda$ of $H$.

In the series of works [40, 36], see also [39] on the Ginzburg-Landau model, Serfaty and her coauthors investigate the next order term in the asymptotic expansion of $\min_{x_{1},\dots,x_{N}}H_{N}$. A keypoint in these works is understanding the behaviour of (the minimum of)

[TABLE]

as $N\rightarrow+\infty$. This is done through the notion of renormalized energy. Roughly, the starting point behind this notion is the (abusive) formal identity

[TABLE]

where $h_{N}$ is the solution of $\Delta h_{N}=4\pi(\delta_{N}-\lambda)$ in ${\mathbb{R}}^{3}$. Of course, this identity does not make sense, as both sides are infinite. On one hand, the left-hand side is not well-defined: the potential $g$ is singular at the diagonal, so that the integral with respect to the product of the empirical measures diverges. On the other hand, the right-hand side is not better defined: as the empirical measure does not belong to $H^{-1}({\mathbb{R}}^{3})$, $h_{N}$ is not in $\dot{H}^{1}({\mathbb{R}}^{3})$.

Still, as explained in [41, chapter 3], one can modify this identity, and show a formula of the form

[TABLE]

where $h_{N}^{\eta}$ is an approximation of $h_{N}$ obtained by regularization of the Dirac masses at the right-hand side of the Laplace equation: $\Delta h_{N}^{\eta}=4\pi(\delta^{\eta}_{N}-\lambda)$ in ${\mathbb{R}}^{3}$. Note the removal of the term $Ng(\eta)$ at the right-hand side of (1.17). This term, which goes to infinity as the parameter $\eta\rightarrow 0$, corresponds to the self-interaction of the Dirac masses: it must be removed, consistently with the fact that the integral defining $\mathcal{H}_{N}$ excludes the diagonal. This explains the term renormalized energy. See [41, chapter 3] for more details.

From there (omitting to discuss the delicate commutation of the limits in $N$ and $\eta$ !), the asymptotics of $\min_{x_{1},\dots,x_{N}}\mathcal{H}_{N}$ can be deduced from the one of $\min_{x_{1},\dots,x_{N}}\int_{{\mathbb{R}}^{3}}|{\nabla}h_{N}^{\eta}|^{2}$, for fixed $\eta$. The next step is to show that such minimum can be expressed as spatial averages of (minimal) microscopic energies, expressed in terms of solutions of the so-called jellium problems: see [41, chapter 4]. These problems, obtained through rescaling and blow-up of the equation on $h_{N}^{\eta}$, are an analogue of cell problems in homogenization. More will be said in Section 4, and we refer to the lecture notes [41] for all necessary complements.

Thus, the main idea in the second part of our paper is to take advantage of the analogy between the functionals $\mathcal{V}_{N}$ and $\mathcal{H}_{N}$ to apply the strategy just described. Doing so, we face specific difficulties: our distribution of points is not minimizing an energy, the potential $g_{S}$ is much more singular than $g$, the reformulation of the functional in terms of an energy is less obvious, etc. Still, we are able to reproduce the same kind of scheme. We introduce in Section 3 an analogue of the renormalized energy. The analogue of the jellium problem is discussed in Section 4. Finally, in Section 5, we are able to tackle the convergence of $\mathcal{V}_{N}$, and give explicit formula for the limit in two cases: the case of a (properly rescaled) $L{\mathbb{Z}}^{3}$-periodic pattern of $M$-spherical particles with centers $a_{1}$, …, $a_{M}$, and the case of a (properly rescaled) hardcore stationary random process with locally integrable two points correlation function $\rho_{2}(y,z)=\rho(y-z)$. In the first case, we show that

[TABLE]

where $G_{S}$ and $G_{S,L}$ are the whole space and $L{\mathbb{Z}}^{3}$-periodic kernels defined respectively in (3.12) and (5.18). See Proposition 5.4. In the special case of a primitive cubic lattice, for which $M=L=1$, we can push further the calculation: we find that

[TABLE]

with $\alpha\approx 9.48$ and $\beta\approx-2,15$, *cf. *Proposition 5.5 for precise expressions. Our result is in agreement with [42]. In the random stationary case, if the process has mean intensity one, we show that

[TABLE]

These formula open the road to numerical computations of the viscosity coefficients of specific processes, and should in particular allow to check the formula found in the literature [5, 35].

Let us conclude this introduction by pointing out that our analysis falls into the general scope of deriving macroscopic properties of dilute suspensions. In this perspective, it can be related to mathematical studies on the drag or sedimentation speed of suspensions, see [25, 13, 23, 24, 30] among many. See also the recent work [14] on the conductivity problem.

2 Expansion of the effective viscosity

The aim of this section is to understand the origin of the functional $\mathcal{V}_{N}$ introduced in (1.13), and to prove Theorem 1.1. The outline is the following. We first consider the effective model (1.4)-(1.5)-(1.6). Given $S$ a symmetric trace-free matrix, and a solution $u_{eff}$ with condition at infinity (1.8), we exhibit an integral quantity $\mathcal{I}_{eff}=\mathcal{I}_{eff}(S)$ that involves $u_{eff}$ and allows to recover (partially) the mean viscosity coefficient $\nu_{2}$. In the next paragraph, we introduce the analogue $\mathcal{I}_{N}$ of $\mathcal{I}_{eff}$, that involves this time the solution $u_{N}$ of (1.1)-(1.2) and (1.8). In brief, we show that if $u_{N}$ is $o(\phi^{2})$ close to $u_{eff}$, then $\mathcal{I}_{N}$ is $o(\phi^{2})$ close to $\mathcal{I}_{eff}$. Finally, we provide an expansion of $\mathcal{I}_{N}$, allowing to express $\nu_{2}$ in terms of $\mathcal{V}_{N}$. Theorem 1.1 follows.

2.1 Recovering the viscosity coefficients in the effective model

Let $k\geq 2$, $\mu^{\prime}$ satisfying (1.6), with viscosity coefficients $\mu_{i}$ that may depend on $x$. Let $S$ symmetric and trace-free. We denote $u_{0}(x)=Sx$. Let $u_{eff}=1_{{\mathbb{R}}^{3}\setminus\mathcal{O}}u+1_{\mathcal{O}}u^{\prime}$ the weak solution in $u_{0}+\dot{H}^{1}({\mathbb{R}}^{3})$ of (1.4)-(1.5)-(1.8). By a standard energy estimate, one can show the expansion

[TABLE]

where the system satisfied by $u_{eff,i}=1_{{\mathbb{R}}^{3}\setminus\mathcal{O}}u_{i}+1_{\mathcal{O}}u^{\prime}_{i}$ is derived by plugging the expansion in (1.4)-(1.5) and keeping terms with power $\phi^{i}$ only. We notably find

[TABLE]

together with the conditions: $u_{1}=0$ at infinity,

[TABLE]

Similarly,

[TABLE]

together with: $u_{2}=0$ at infinity,

[TABLE]

We now define, inspired by formula (4.11.16) in [4],

[TABLE]

where $n$ refers to the outward normal. We will show that

[TABLE]

We first use (1.5) to write

[TABLE]

Using the equations satisfied by $u^{\prime}_{1}$ and $u^{\prime}_{2}$, after integration by parts, we get

[TABLE]

Plugging this last line in the expression for $\mathcal{I}_{eff}$ yields (2.4).

We see through formula (2.4) that the expansion of $\mathcal{I}_{eff}$ in powers of $\phi$ gives access to $\nu_{1}$, and, if $\mu_{1}$ is known, it further gives access to $\nu_{2}$. On the basis of the works [1, 33] and of the recent paper [21] which considers the same setting as ours, it is natural to assume that $\mu_{1}$ is given by (1.10). This implies $\nu_{1}=\frac{5}{2}\mu$. With such expression of $\mu_{1}$, and the form of $f$ specified in (H1), we can check that $u_{S}=(5|\mathcal{O}|)^{-1}u_{eff,1}$ satisfies

[TABLE]

It follows that

[TABLE]

2.2 Recovering the viscosity coefficients in the model with particles

To determine the possible value of the mean viscosity coefficient $\nu_{2}$, we must now relate the functional $\mathcal{I}_{eff}$, based on the effective model, to a functional $\mathcal{I}_{N}$ based on the real model with spherical rigid particles. From now on, we place ourselves under the assumptions of Theorem 1.1. Note that, thanks to hypothesis (H2), the spherical particles do not overlap for $\phi$ small enough, so that a weak solution $u_{N}\in u_{0}+\dot{H}^{1}({\mathbb{R}}^{3})$ of (1.1)-(1.2)-(1.8) exists and is unique.

By integration by parts, for any $R$ such that $\mathcal{O}\Subset B_{R}$, we have

[TABLE]

By analogy with (2.3), and as all particles remain in a fixed compact $K\supset\mathcal{O}$ independent of $N$, we set for any $R$ such that $K\subset B_{R}$:

[TABLE]

which again does not depend on our choice of $R$ by integration by parts. Now, if $u_{eff}$ and $u_{N}$ are $o(\phi^{2})$-close in the sense of Theorem 1.1, then

[TABLE]

Indeed, $u_{N}-u_{eff}$ is a solution of a homogenenous Stokes equation outside $K$. By elliptic regularity, we find that $\limsup_{N\rightarrow+\infty}\|u_{eff}-u_{N}\|_{H^{s}(K^{\prime})}=0$, for any compact $K^{\prime}\subset{\mathbb{R}}^{3}\setminus K$ and any positive $s$. Relation (2.9) follows.

We now turn to the most difficult part of this section, that is expanding $\mathcal{I}_{N}$ in powers of $\phi$. We aim at proving

Proposition 2.1**.**

Let $(x_{i})_{1\leq i\leq N}$, satisfying (H1)-(H2). For $S$ trace-free and symmetric, for $\phi>0$, let $u_{N}$ the solution of (1.1)-(1.2)-(1.8) with the ball radius $a$ defined through (1.3). Let $\mathcal{I}_{N}$ as in (2.8), $\mathcal{V}_{N}$ as in (1.13), and $u_{S}$ the solution of (2.5). One has

[TABLE]

As before, notation $A_{N}=B_{N}+o(\phi^{2})$ means $\limsup_{N}|A_{N}-B_{N}|=o(\phi^{2})$. Obviously, Theorem 1.1 follows directly from (2.6), (2.9) and from the proposition.

To start the proof, we set $v_{N}:=u_{N}-u_{0}$. Note that $v_{N}\in\dot{H}^{1}(\mathcal{F})$ still satisfies the Stokes equation outside the ball, with $v_{N}=0$ at infinity, and $v_{N}=-Sx+u_{i}+\omega_{i}\times(x-x_{i})$ inside $B_{i}$. Moreover, taking into account the identities:

[TABLE]

and

[TABLE]

one has for all $i$:

[TABLE]

From the definition (2.8), we can re-express $\mathcal{I}_{N}$ as

[TABLE]

To obtain an expansion of $\mathcal{I}_{N}$ in powers of $\phi$, we will now approximate $(v_{N},p_{N})$ by some explicit field $(v_{app},p_{app})$, inspired by the method of reflections. This approximation involves the elementary problem:

[TABLE]

The solution of (2.13) is explicit [18], and given by

[TABLE]

with

[TABLE]

The pressure is

[TABLE]

We now introduce

[TABLE]

where

[TABLE]

In short, the first sum at the right-hand side of (2.17) corresponds to a superposition of $N$ elementary solutions, meaning that the interaction between the balls is neglected. This sum satisfies the Stokes equation outside the ball, but creates an error at each ball $B_{i}$, whose leading term is $S_{i}x$. This explains the correction by the second sum at the right-hand side of (2.17). One could of course reiterate the process: as the distance between particles is large compared to their radius, we expect the interactions to be smaller and smaller. This is the principle of the method of reflections that is investigated in [24]. From there, Proposition 2.1 will follow from two facts. Defining

[TABLE]

we will show first that

[TABLE]

and then

[TABLE]

Identity (2.19) follows from a calculation that we now detail. We define

[TABLE]

We have

[TABLE]

To treat $I_{b}$ and $I_{d}$, we rely on the following property, that is checked easily through integration by parts: for any $(v,p)$ solution of Stokes in $B_{i}$, and any trace-free symmetric matrix $S$, $\mathcal{I}_{i}(v,p):S=0$. As for all $i$ and all $j\neq i$, $v^{s}[S](\cdot-x_{j})$ or $v^{s}[S_{j}](\cdot-x_{j})$ is a solution of Stokes inside $B_{i}$, we deduce

[TABLE]

As regards $I_{a}$, we use the following formula which follows from a tedious calculation [18]: for any traceless matrix $S$,

[TABLE]

It follows that

[TABLE]

This term corresponds to the famous Einstein formula for the mean effective viscosity. It is coherent with the expression (1.10) for $\mu_{1}$, which implies $\nu_{1}=\frac{5}{2}\mu$.

Eventually, as regards $I_{c}$, we can use again (2.22), replacing $S$ by $S_{i}$:

[TABLE]

with $g_{S}$ defined in (1.14). In view of (2.21)-(2.23)-(2.24), to conclude that (2.19) holds, it is enough to prove

Lemma 2.2**.**

For any $f\in L^{2}({\mathbb{R}}^{3})$,

[TABLE]

with $g_{S}$ defined in (1.14), and $u_{S}\in\dot{H}^{1}({\mathbb{R}}^{3})$ the solution of (2.5).

Proof. Note that both sides of the identity are continuous over $L^{2}$: the left-hand side is continuous as the Calderón-Zygmund operator $f\rightarrow g_{S}\star f$ is continuous over $L^{2}$, while the right-hand side is continuous by classical elliptic estimates for the Stokes operator. By density, this is therefore enough to assume that $f\in C^{\infty}_{c}({\mathbb{R}}^{3})$. We denote by $U=(U_{ij}),Q=(Q_{j})$ the fondamental solution of the Stokes operator. This means that for all $j$, the vector field $U_{j}=(U_{ij})_{1\leq i\leq N}$ and the scalar field $Q_{j}$ satisfy the Stokes equation

[TABLE]

It is well-known, see [16, page 239], that

[TABLE]

From there, one can deduce the following formula, cf [16, page 290, equation (IV.8.14)]:

[TABLE]

Using the Einstein convention for summation, this implies in turn that

[TABLE]

where we have used that $S$ is trace-free to obtain the third equality. Hence,

[TABLE]

Note that the permutations between the derivatives and the convolution product do not raise any difficulty, as $f\in C_{c}^{\infty}({\mathbb{R}}^{3})$. Now, using $S_{kl}=S_{lk}$, and denoting by ${\rm St}^{-1}$ the convolution with the fundamental solution (inverse of the Stokes operator), we get

[TABLE]

Eventually,

[TABLE]

This concludes the proof of the lemma.

Remark 2.3*.*

By polarization of the previous identity, at least for $f,\tilde{f}$ smooth and decaying enough, one has

[TABLE]

The last step in proving Proposition 2.1, hence Theorem 1.1, is to show the bound (2.20). If $w:=v_{N}-v_{app}$, $q:=p_{N}-p_{app}$,

[TABLE]

Direct verifications show that $v_{app}$, hence $w$, satisfies the same force- and torque-free conditions as $v$. This means that for any family of constant vectors $u_{i}$ and $\omega_{i}$, $1\leq i\leq N$,

[TABLE]

By a proper choice of $u_{i}$ and $\omega_{i}$, we find

[TABLE]

for any family $(\tilde{u}_{i},\tilde{\omega}_{i})$, using this time that $v_{N}$ is force- and torque-free. Let $q\geq 2$. By a proper choice of $(\tilde{u}_{i},\tilde{\omega}_{i})$, by Poincaré and Korn inequalities, one can ensure that for all $i$,

[TABLE]

where

[TABLE]

Note that the factor $\frac{1}{a}$ at the right-hand side is consistent with scaling considerations. Moreover, by standard use of the Bogovskii operator, see [16], there exists a constant $C$ (depending only on the constant $c$ in (H2)) and a field $W\in W^{1,q}(\mathcal{F})$ , zero outside $\cup_{i=1}^{N}B(x_{i},2a)$ satisfying

[TABLE]

We deduce, with $p\leq 2$ the conjugate exponent of $q$:

[TABLE]

By well-known variational properties of the Stokes solution, $\|D(v_{N})\|_{L^{2}}$ minimizes $\|D(v)\|_{L^{2}}$ over the set of all $v$ in $\dot{H}^{1}({\mathbb{R}}^{3})$ satisfying a boundary condition of the form $v|_{B_{i}}=-Sx+u_{i}+\omega_{i}\times(x-x_{i})$ for all $i$. By the same considerations as before, based on the Bogovski operator, we infer that

[TABLE]

so that

[TABLE]

Using this inequality with the first term in (2.31) and applying the Hölder inequality to the second term, we end up with

[TABLE]

To deduce (2.20), it is now enough to prove that for all $q>1$, there exists a constant $C$ independent of $N$ or $\phi$ such that

[TABLE]

Indeed, taking $q>2$, meaning $p<2$, and combining this inequality with (2.32) yields (2.20), more precisely

[TABLE]

In order to show the bound (2.33), we must write down the expression for $w|_{B_{i}}=v_{N}|_{B_{i}}-v_{app}|_{B_{i}}$, where $v_{app}$ was introduced in (2.17). A little calculation, using Taylor’s formula with integral remainder, shows that

[TABLE]

with $w_{i}^{r}$ a rigid vector field (that disappears when taking the symmetric gradient), with

[TABLE]

and with the bilinear application:

[TABLE]

We remind that $v^{s}[S]$ and $v[S]$ were introduced in (2.14) and (2.15), while the matrices $S_{j}$ are defined in (2.18). Note that the matrices $D_{i}$ and $S_{i}$ have the same kind of structure. More precisely, we can define for a collection $(A_{1},\dots,A_{N})$ of $N$ symmetric matrices, an application

[TABLE]

Then, $(S_{1},\dots,S_{N})=\mathcal{A}(S,\dots,S)$ and $(D_{1},\dots,D_{N})=\mathcal{A}(S_{1},\dots,S_{N})=\mathcal{A}^{2}(S,\dots,S)$. Note that for any matrix $A$, the kernel $D(v[A])$, homogenenous of degree $-3$, is of Calderón-Zygmund type. Using this property, we are able to prove in the appendix the following lemma, which is an adaptation of a result by the second author and Di Wu [21]:

Lemma 2.4**.**

For all $1<q<+\infty$, there exists a constant $C$, depending on $q$ and on the constant $c$ in (H2), such that, if $(A^{\prime}_{1},\dots,A^{\prime}_{N})=\mathcal{A}(A_{1},\dots,A_{N})$, then

[TABLE]

We can now proceed to the proof of (2.33). Denoting $w_{i}^{1}:=D_{i}(x-x_{i})$, we find by the lemma:

[TABLE]

Then, we notice that for any matrix $A$, $|D(v^{s}[A]-v[A])(x)|=O(a^{5}|x|^{-5})$. This implies that $w^{2}_{i}:=E_{i}(x-x_{i})$ satisfies

[TABLE]

By assumption (H2), the points $y_{i}:=N^{1/3}x_{i}$ satisfy for all $i\neq j$:

[TABLE]

In particular,

[TABLE]

We then make use of the following easy generalization of Young’s convolution inequality:

[TABLE]

Applied with $a_{ij}=\frac{1}{(c+|y_{i}-y_{j}|)^{5}}$ and $b_{j}=|S|+|S_{j}|$, together with Lemma 2.4, it yields

[TABLE]

It remains to bound the symmetric gradient of $w_{i}^{3}:=\mathbf{F}_{i}|_{x}(x-x_{i},x-x_{i})$. By the expression of $v^{s}$, we get that in $B_{i}$:

[TABLE]

Proceeding as above, we find

[TABLE]

As $D(w)=D(w_{i}^{1})+D(w_{i}^{2})+D(w_{i}^{3})$, cf. (2.34), the previous estimates yield (2.33). This concludes the proof of Proposition 2.1, and therefore the proof of Theorem 1.1.

3 The $\phi^{2}$ correction $\mathcal{V}_{N}$ as a renormalized energy

We start in this section the asymptotic analysis of the viscosity coefficient

[TABLE]

As a preliminary step, we will show that there is no loss of generality in assuming

[TABLE]

We introduce the set

[TABLE]

By (H1)-(H2), it is easily seen that $N_{ext}=o(N)$ as $N\rightarrow+\infty$. We now show

Lemma 3.1**.**

$\mathcal{V}_{N}$* is uniformly bounded in $N$, and*

[TABLE]

goes to zero as $N\rightarrow+\infty$.

Proof. For any open set $U$, we denote \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{U}=\frac{1}{|U|}\int_{U}.

Let $d:=\frac{c}{4}N^{-1/3}\leq\min_{i\neq j}\frac{|x_{i}-x_{j}|}{4}$ by (H2). We write

[TABLE]

For the first term, with $y_{i}:=N^{1/3}x_{i}$ and with (H2) in mind, that is $|y_{i}-y_{j}|\geq c$ for $i\neq j$:

[TABLE]

see (1.14). This yields, by a discrete convolution inequality:

[TABLE]

where we have used that $\sum_{j=1}^{N}\frac{1}{(c+|y_{i}-y_{j}|)^{4}}$ is uniformly bounded in $N$ and in the index $i$ thanks to the separation assumption. By similar arguments, $|II|\leq\mathcal{C}$. As regards the last term, we notice that

[TABLE]

where $F_{N}=\sum_{i=1}^{N}1_{B(x_{i},d)}$. The operator $\mathcal{T}F(x)=\int g_{S}(x-y)F(y)dy$ is a Calderón-Zygmund operator, and therefore continuous over $L^{2}$. As $F_{N}^{2}=F_{N}$ (the balls are disjoint), we find that the $L^{2}$ norm of $F_{N}$ is $(Nd^{3})^{1/2}$ and

[TABLE]

Similarly,

[TABLE]

It follows that $|III|\leq\frac{C}{Nd^{3}}$. With our choice of $d$, the first part of the lemma is proved.

From there, to prove that $\mathcal{V}_{N,ext}$ goes to zero, as $N_{ext}=o(N)$, it is enough to show that

[TABLE]

By symmetry, it is enough that

[TABLE]

This can be shown by a similar decomposition as the previous one. Namely,

[TABLE]

Proceeding as above, we find this time

[TABLE]

which concludes the proof.

Remark 3.2*.*

By Lemma 3.1, there is no restriction assuming (3.1) when studying the asymptotic behaviour of $\mathcal{V}_{N}$. Therefore, we make from now on the assumption (3.1).

As explained in the introduction, the analysis of $\mathcal{V}_{N}$ will rely on the mathematical methods introduced over the last years for Coulomb gases, the core problem being the analysis of a functional of the form (1.15). We shall first reexpress $\mathcal{V}_{N}$ in a similar form. More precisely, we will show

Proposition 3.3**.**

Denoting

[TABLE]

we have $\mathcal{V}_{N}=\mathcal{W}_{N}+\varepsilon(N)$ where $\varepsilon(N)\to 0$ as $N\to\infty.$

Remark 3.4*.*

In the definition of $\mathcal{W}_{N}$, the integrals of the form

[TABLE]

that appear when expanding the product, are understood as

[TABLE]

where ${\rm St}$ is the Stokes operator, see (2.30) and the proof below for an explanation.

Proof. Clearly,

[TABLE]

so that formally

[TABLE]

Note that it is not obvious that this formal decomposition makes sense, because all three quantities at the right-hand side involve integrals of $g_{S}(x-y)$ against product measures of the form $d\delta_{N}(x)f(y)dy$ (or the symmetric one), which may fail to converge because of the singularity of $g_{S}$. To solve this issue, a rigorous path consists in approximating, at fixed $N$, each Dirac mass $\delta_{x_{i}}$ by a (compactly suppported) approximation of unity $\rho_{\eta}(x-x_{i})$, where $\eta>0$ is the approximation parameter and goes to zero. One can then set, for each $\eta$, $\delta_{N}^{\eta}(x):=\frac{1}{N}\sum_{i=1}^{N}\rho_{\eta}(x-x_{i})$, leading to the rigorous decomposition

[TABLE]

where $\mathcal{V}_{N}^{\eta}$, $\mathcal{W}_{N}^{\eta}$ are deduced from $\mathcal{V}_{N}$, $\mathcal{W}_{N}$ replacing the empirical measure by its regularization. It is easy to show that $\lim_{\eta\rightarrow 0}\mathcal{V}_{N}^{\eta}=\mathcal{V}_{N}$. To conclude the proof, we shall establish the following: first,

[TABLE]

the same limit holding for the symmetric term. In particular, (3.2) will show that $\mathcal{W}_{N}=\lim_{\eta\rightarrow 0}\mathcal{W}_{N}^{\eta}$ exists, in the sense given in Remark 3.4. Then, we will prove

[TABLE]

which together with (3.2) will complete the proof of the proposition.

The limit (3.2) follows from identity (2.30). Indeed, for $\eta>0$, this formula yields

[TABLE]

Now, we remark that due to our assumptions on $f$, by elliptic regularity, $h={\rm St}^{-1}(S{\nabla}f)(x)$ is $C^{1}$ inside $\mathcal{O}$. Moreover, in virtue of Remark (3.2), we can assume (3.1). Hence, as $\eta\rightarrow 0$,

[TABLE]

It remains to prove (3.3). In the special case where $f\in C^{r}({\mathbb{R}}^{3})$ for some $r\in(0,1)$ (implying that it vanishes at ${\partial}\mathcal{O}$), classical results on Calderón-Zygmund operators yield that the function $\int_{{\mathbb{R}}^{3}}g_{S}(x-y)f(x)dx=\frac{8\pi}{3}S{\nabla}\cdot h(y)$ is a continuous (even Hölder) bounded function, so (H1) implies straightforwardly

[TABLE]

In the general case where $f$ is discontinuous across ${\partial}\mathcal{O}$, the proof is a bit more involved. The difficulty lies in the fact that some points $x_{i}$ get closer to the boundary as $N\rightarrow+\infty$.

Let ${\varepsilon}>0$. Under (H2), there exists $c^{\prime}>0$ (depending on $c$ only) such that for $N^{-1/3}\leq{\varepsilon}$,

[TABLE]

Let $\chi_{\varepsilon}:{\mathbb{R}}^{3}\rightarrow[0,1]$ be a smooth function such that $\chi_{\varepsilon}=1$ in a $c^{\prime}{\varepsilon}/4$ neighborhood of ${\partial}\mathcal{O}$, $\chi_{\varepsilon}=0$ outside a $c^{\prime}{\varepsilon}/2$ neighborhood of ${\partial}\mathcal{O}$. We write

[TABLE]

By formula (2.30), the second term reads

[TABLE]

with $u_{\varepsilon}={\rm St}^{-1}S{\nabla}((1-\chi_{\varepsilon})f)$. The source term $(1-\chi_{\varepsilon})f$ being $C^{1}$ and compactly supported, $S{\nabla}\cdot u_{\varepsilon}$ is Hölder and bounded, so that, as $N\rightarrow+\infty$, the integral goes to zero by the weak convergence assumption (H1), for any fixed ${\varepsilon}>0$. As regards the first term, we split it again into

[TABLE]

where $v_{\varepsilon}$ is this time the solution of the Stokes equation with source $S{\nabla}(\chi_{\varepsilon}f)$. It is Hölder away from ${\partial}\mathcal{O}$, so that the last term at the right-hand side goes again to zero as $N\rightarrow+\infty$, by assumption (H1).

It remains to handle the first term at the right-hand side. We shall show below that for a proper choice of $\chi_{\varepsilon}$ one has

[TABLE]

Taking advantage of this fact, we write

[TABLE]

where we used property (3.4) to obtain the last inequality. With this bound and the convergence to zero of the other terms for fixed ${\varepsilon}$ and $N\rightarrow+\infty$, the limit (3.3) follows.

We still have to show that ${\nabla}v^{\varepsilon}$ is uniformly bounded in $L^{\infty}$ for a good choice of $\chi_{\varepsilon}$. We borrow here to the analysis of vortex patches in the Euler equation, initiated by Chemin in 2-d [10], extended by Gamblin and Saint-Raymond in 3-d [17]. First, as $\mathcal{O}$ is smooth, one can find a family of five smooth divergence-free vector fields $w_{1},\dots,w_{5}$, tangent at ${\partial}\mathcal{O}$ and non-degenerate in the sense that

[TABLE]

see [17, Proposition 3.2]. We take $\chi_{\varepsilon}$ in the form $\chi(t/{\varepsilon})$, for a coordinate $t$ transverse to the boundary, meaning that ${\partial}_{t}$ is normal at ${\partial}\mathcal{O}$. With this choice and the assumptions on $f$, one checks easily that $\chi_{\varepsilon}f$ is bounded uniformly in ${\varepsilon}$ in $L^{\infty}({\mathbb{R}}^{3})$ and that for all $i$, $w_{i}\cdot{\nabla}(\chi_{\varepsilon}f)$ is bounded uniformly in ${\varepsilon}$ in $C^{0}({\mathbb{R}}^{3})\subset C^{r-1}({\mathbb{R}}^{3})$ for all $r\in(0,1)$. Hence, the norm $\|\chi_{\varepsilon}f\|_{r,W}$ introduced in [17, page 395], where $W=(w_{1},\dots,w_{5})$, is bounded uniformly in ${\varepsilon}$.

We then split the Stokes system

[TABLE]

into the equations:

[TABLE]

and

[TABLE]

Let us show that ${\partial}_{i}{\partial}_{j}\Delta^{-1}(\chi_{\varepsilon}f)$ is bounded uniformly in ${\varepsilon}$ in $L^{\infty}$. Let $\chi\in C^{\infty}_{c}({\mathbb{R}}^{3})$, $\chi\geq 0$, $\chi=1$ near zero. Let for all $m\in{\mathbb{R}}$, $\Lambda^{m}(\xi):=(\chi(\xi)+|\xi|^{2})^{m/2}$. It is easily seen through Fourier transform that for all $s\in{\mathbb{N}}$

[TABLE]

Moreover, by the calculations in [17, page 401], replacing $\omega$ with $\chi_{\varepsilon}f$, we get

[TABLE]

Combining (3.6) and (3.7), we find that

[TABLE]

is bounded uniformly in ${\varepsilon}$ in $L^{\infty}$ and consequently,

[TABLE]

Also, by continuity of Riesz transforms over $L^{p}$, we have

[TABLE]

Now, applying $w_{k}\cdot{\nabla}$ to the equation satisfied by $\Omega_{\varepsilon}$, we obtain for all $1\leq k\leq 5$,

[TABLE]

where $F_{i,j,{\varepsilon}}$, $G_{i,{\varepsilon}}$ and $H_{\varepsilon}$ are combinations of $\Omega_{\varepsilon}$, $\chi_{\varepsilon}f$ and $w_{k}\cdot{\nabla}(\chi_{\varepsilon}f)$. In particular, they are bounded uniformly in ${\varepsilon}$ in $L^{\infty}\cap L^{p}$, for any $1<p<\infty$.

For the first term at the r.h.s., we write with the same cut-off function $\chi$ as before:

[TABLE]

By continuity of $(-\Delta)^{-1}{\partial}_{i}{\partial}_{j}$ over $L^{2}$, the first term, with low frequencies, belongs to $H^{s}$ for any $s$, with uniform bound in ${\varepsilon}$. By continuity of $(1-\chi(D))(-\Delta)^{-1}{\partial}_{i}{\partial}_{j}$ over Hölder spaces (Coifman-Meyer theorem), the second term, with high frequencies, is uniformly bounded in ${\varepsilon}$ in $C^{r-1}({\mathbb{R}}^{3})$, for any $0<r<1$.

For the second and third terms in (3.8), we claim that

[TABLE]

This can be seen easily by expressing these fields as $\sum_{i}{\partial}_{i}\Phi\star G_{i,{\varepsilon}}$ and $\Phi\star H_{\varepsilon}$ with $\Phi$ the fundamental solution, and by using the uniform $L^{p}$ bounds on $G_{i,{\varepsilon}}$ and $H_{\varepsilon}$. Eventually, we find that

[TABLE]

We conclude by [17, Proposition 3.3] that ${\nabla}v_{\varepsilon}$ is bounded in $L^{\infty}({\mathbb{R}}^{3})$ uniformly in ${\varepsilon}$.

3.1 Smoothing

By Proposition 3.3, we are left with understanding the asymptotic behaviour of

[TABLE]

The following field will play a crucial role. For $U,Q$ defined in (2.26), we set

[TABLE]

From (2.27), we have $g_{S}=\frac{8\pi}{3}\,(S{\nabla})\cdot G_{S}$, and $G_{S}$ solves in the sense of distributions

[TABLE]

Moreover, from the explicit expression

[TABLE]

and taking into account the fact that $S$ is symmetric and trace-free, we get

[TABLE]

Let us note that $G_{S}$ is called a point stresslet in the literature, see [18]. It can be interpreted as the velocity field created in a fluid of viscosity $1$ by a point particle whose resistance to a strain is given by $-S$.

We now come back to the analysis of (3.9). Formal replacement of the function $f$ in Lemma 2.2 by $\delta_{N}-f$ yields the formula

[TABLE]

where

[TABLE]

satisfies

[TABLE]

The formula (3.13) is similar to the formula (1.16), and is as much abusive, as both sides are infinite. Still, by an appropriate regularization of the source term $S{\nabla}\sum_{i}\delta_{x_{i}}$, we shall be able in the end to obtain a rigorous formula, convenient for the study of $\mathcal{W}_{N}$. This regularization process is the purpose of the present paragraph.

For any $\eta>0$, we denote $B_{\eta}=B(0,\eta)$, and define $G_{S}^{\eta}$ by:

[TABLE]

[TABLE]

Note that by homogeneity,

[TABLE]

The field $G_{S}^{\eta}$ belongs to $\dot{H}^{1}({\mathbb{R}}^{3})$, and solves

[TABLE]

where $S^{\eta}$ is the measure on the sphere defined by

[TABLE]

with $n=\frac{x}{|x|}$ the unit normal vector pointing outward $B_{\eta}$, $\left[F\right]:=F|_{{\partial}B_{\eta}^{+}}-F|_{{\partial}B_{\eta}^{-}}$ the jump at ${\partial}B_{\eta}$ (with ${\partial}B_{\eta}^{+}$, resp. ${\partial}B_{\eta}^{-}$, the outer, resp. inner boundary of the ball), and $s^{\eta}$ the standard surface measure on $\partial B_{\eta}$. We claim the following

Lemma 3.5**.**

For all $\eta>0$, $\>S^{\eta}={\rm div}\,\Psi^{\eta}\>$ in ${\mathbb{R}}^{3}$, where

[TABLE]

Moreover, $\Psi^{\eta}\rightarrow S\delta$ in the sense of distributions as $\eta\rightarrow 0$, so that $S^{\eta}\rightarrow S{\nabla}\delta$.

Proof of the lemma. From the explicit formula (3.12) for $G_{S}$ and $p_{S}$, we find

[TABLE]

so that

[TABLE]

Using that $S$ is trace-free, one can check from definition (3.21) that $\hbox{div \!}\Psi^{\eta}=0$ in the complement of ${\partial}B_{\eta}$, while:

[TABLE]

where the last equality comes from (3.22). Together with (3.20), it implies the first claim of the lemma.

To compute the limit of $\Psi^{\eta}$ as $\eta\rightarrow 0$, we write $\Psi^{\eta}=\Psi^{\eta}_{1}+\Psi^{\eta}_{2}$, with

[TABLE]

Let $\varphi\in C^{\infty}_{c}({\mathbb{R}}^{3})$ a test function. We can write $\langle\Psi_{1}^{\eta},\varphi\rangle=\langle\Psi_{1}^{\eta},\varphi(0)\rangle+\langle\Psi_{1}^{\eta},\varphi-\varphi(0)\rangle.$ The second term is $O(\eta)$, while the first term can be computed using the elementary formula $\int_{B_{1}}x_{i}x_{j}dx=\frac{4\pi}{15}\delta_{ij}$. We find

[TABLE]

For the second term, using the homogeneity (3.18), we find again that $\lim_{\eta}\langle\Psi_{2}^{\eta},\varphi\rangle=\langle\Psi_{2}^{1},\varphi(0)\rangle$. Note that the pressure $p_{S}^{1}$ is defined up to a constant, so that we can always select the one with zero average. With this choice, we find

[TABLE]

where the sixth equality comes from the elementary formula $\int_{{\partial}B_{1}}n_{i}n_{j}n_{k}n_{l}ds^{1}=\frac{4\pi}{15}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$. The result follows.

For later purpose, we also prove here the

Lemma 3.6**.**

[TABLE]

Proof.

[TABLE]

By (3.18), $\int_{B_{\eta}}|{\nabla}G^{\eta}_{S}|^{2}dx=\frac{1}{\eta^{3}}\int_{B_{1}}|{\nabla}G^{1}_{S}|^{2}dx$. The second term can be computed with (3.12):

[TABLE]

3.2 The renormalized energy

Thanks to the regularization of $S{\nabla}\delta$ introduced in the previous paragraph, cf. Lemma 3.5, we shall be able to set a rigorous alternative to the abusive formula (3.13). Specifically, we shall state an identity involving $\mathcal{W}_{N}$, defined in (3.9), and the energy of the function

[TABLE]

This function solves

[TABLE]

and is a regularization of $h_{N}$, cf. (3.14)-(3.15).

The main result of this section is the

Proposition 3.7**.**

[TABLE]

Proof. We assume that $\eta$ is small enough so that $2\eta<\min_{i\neq j}|x_{i}-x_{j}|$. From the explicit expressions (3.14), (3.25), we find that $h_{N},h_{N}^{\eta}=O(|x|^{-2})$, ${\nabla}(h_{N},h_{N}^{\eta})=O(|x|^{-3})$ and $p_{N},p_{N}^{\eta}=O(|x|^{-3})$ at infinity. As these quantities decay fast enough, we can perform an integration by parts to find

[TABLE]

where we defined $h_{N}^{i}:=h_{N}-G_{S}(x-x_{i})$.

As $h_{N}^{i}$ is smooth over the support of $S^{\eta}(\cdot-x_{i})$, we can apply Lemma 3.5 to obtain

[TABLE]

We can then apply Lemma 3.6 to obtain

[TABLE]

As regards the fourth term, we notice that by our definition (3.16)-(3.17) of $G_{S}^{\eta}$, and the fact that the balls $B(x_{i},\eta)$ are disjoint, the function $h_{N}-h_{N}^{\eta}=\sum_{i}(G_{S}(x-x_{i})-G_{S}^{\eta}(x-x_{i}))$ is zero over $\cup_{i}{\partial}B(x_{i},\eta)$, which is the support of $\sum_{i}S^{\eta}(x-x_{i})$. It follows that

[TABLE]

where we integrated by parts, using that $G_{S}-G_{S}^{\eta}$ is zero outside the balls. Let us notice that the second integral at the right-hand side converges despite the singularity of $S{\nabla}\cdot G_{S}$, using the smoothness of $f$ near $x_{i}$ (by assumption (3.1) and Remark 3.2). Moreover, it goes to zero as $\eta\rightarrow 0$. Using the homogeneity and smoothness properties of $G_{S}^{\eta}$ inside $B^{\eta}$, we also find that the first sum goes to zero with $\eta$, resulting in

[TABLE]

We end up with

[TABLE]

It remains to rewrite properly the right-hand side: we first get

[TABLE]

and integrating by parts

[TABLE]

The last equality was deduced from the identity $g_{S}=\frac{8\pi}{3}\,(S{\nabla})\cdot G_{S}$, see the line after (3.10). The proposition follows.

We can refine the previous proposition by the following

Proposition 3.8**.**

Let $c>0$ the constant in (H2). There exists $C>0$ such that: for all $\alpha<\eta<\frac{c}{2}N^{-1/3}$,

[TABLE]

Proof. One has from (3.25)

[TABLE]

It follows that

[TABLE]

After integration by parts,

[TABLE]

while

[TABLE]

We get

[TABLE]

We note that $G_{S}^{\eta}-G_{S}^{\alpha}$ is zero outside $B_{\eta}$, while $S^{\alpha}+S^{\eta}$ is supported in $B_{\eta}$. Moreover, thanks to (H2), for $\alpha<\eta<\frac{c}{2}$, the balls $B(x_{i},\eta)$ are disjoint. We deduce: $a=0$.

After integration by parts, taking into account that $G_{S}^{\eta}-G_{S}^{\alpha}$ vanishes outside $B_{\eta}$, we can write $b=b_{\eta}-b_{\alpha}$ with

[TABLE]

By assumption (3.1), for $N$ large enough, for all $1\leq i\leq N$ and all $\eta\leq\frac{c}{2}N^{-1/3}$, $B(x_{i},\eta)$ is included in $\mathcal{O}$. Hence, $f$ is $C^{1}$ in $B(x_{i},\eta)$, and

[TABLE]

This results in: $b_{\eta}\leq CN^{2}\eta$.

Similarly, decomposing B(x_{i},\eta)=B(x_{i},\alpha)\cup\Big{(}B(x_{i},\eta)\setminus B(x_{i},\alpha)\Big{)}, we find

[TABLE]

using again that $f$ is Lipschitz over $B(x_{i},\eta)$. We end up with $b_{\alpha}\leq CN^{2}\eta$, and finally $b\leq CN^{2}\eta$.

For the last term $c$ in (3.28), we first notice that as $G_{S}^{\eta}-G_{S}^{\alpha}$ is zero outside $B_{\eta}$:

[TABLE]

where we used Lemma 3.6 in the last line. By the definition of $S^{\alpha}$, the remaining term splits into

[TABLE]

By integration by parts, applied in $B_{\eta}\setminus B_{\alpha}$ for the first term and in $B_{\alpha}$ for the second term, we get

[TABLE]

From there, the conclusion follows easily.

If we let $\alpha\rightarrow 0$ in Proposition 3.8, combining with Propositions 3.27 and 3.3, we find

Corollary 3.9**.**

For all $\eta<\frac{c}{2}N^{-1/3}$,

[TABLE]

where ${\varepsilon}(N)\rightarrow 0$ as $N\rightarrow+\infty$.

This corollary shows that to understand the limit of $\mathcal{V}_{N}$, it is enough to study the limit of

[TABLE]

for $\eta_{N}:=\eta N^{-1/3}$, $\eta<\frac{c}{2}$ fixed. For periodic and more general stationary point processes, this will be possible through an homogenization approach. This homogenization approach involves an analogue of a cell equation, called jellium in the literature on Coulomb gases. We will motivate and introduce this system in the next section.

4 Blown-up system

Formula (3.27) suggests to understand at first the behaviour of $\int_{{\mathbb{R}}^{3}}|{\nabla}h_{N}^{\eta}|^{2}$ at fixed $\eta$, when $N\rightarrow+\infty$. To analyze the system (3.26), a useful intuition can be taken from classical homogenization problems of the form

[TABLE]

with $F(x,y)$ periodic in variable $y$, and $\overline{F}(x):=\int_{{\mathbb{T}}}F(x,y)dy$. Roughly, $\Omega$ would be like $\mathcal{O}$, the small scale ${\varepsilon}$ like $N^{-1/3}$, the term $\frac{1}{{\varepsilon}^{3}}F(x,x/{\varepsilon})$ would correspond to the sum of (regularized) Dirac masses, while the term $\frac{1}{{\varepsilon}^{3}}\overline{F}$ would be an analogue of $Nf$. The factor $\frac{1}{{\varepsilon}^{3}}$ in front of $F$ is put consistently with the fact that $\sum_{i}\delta_{x_{i}}$ has mass $N$. The dependence on $x$ of the source term in (4.1) mimics the possible macroscopic inhomogeneity of the point distribution $\{x_{i}\}$.

In the much simpler model (4.1), standard arguments show that $h_{\varepsilon}$ behaves like

[TABLE]

where $H(x,y)$ satisfies the cell problem

[TABLE]

Let us stress that substracting the term $\frac{1}{{\varepsilon}^{3}}\overline{F}(x)$ in the source term of (4.1) is crucial for the asymptotics (4.2) to hold. It follows that

[TABLE]

Note that the factor ${\varepsilon}^{6}$ in front of the left-hand side is coherent with the factor $\frac{1}{N^{2}}$ at the right-hand side of (3.27). Note also that

[TABLE]

Such average over larger and larger boxes may be still meaningful in more general settings, typically in stochastic homogenization.

Inspired by those remarks, and back to system (3.26), the hope is that some homogenization process may take place, at least locally near each $x\in\mathcal{O}$. More precisely, we hope to recover $\lim_{N}\mathcal{W}_{N}$ by summing over $x\in\mathcal{O}$ some microscopic energy, locally averaged around $x$. This microscopic energy will be deduced from an analogue of the cell problem, called a jellium in the literature on the Ginzburg-Landau model and Coulomb gases.

4.1 Setting of the problem

We will call point distribution a locally finite subset of ${\mathbb{R}}^{3}$. Given a point distribution $\Lambda$, we consider the following problem in ${\mathbb{R}}^{3}$

[TABLE]

Given a solution $H=H(y)$, $P=P(y)$, we introduce for any $\eta>0$

[TABLE]

which satisfies by (3.11), (3.19):

[TABLE]

We remark that, the set $\Lambda$ being locally finite, the sum at the right-hand side of (4.3) or (4.5) is well-defined as a distribution. Also, the sum at the right-hand side of (4.4) is well-defined pointwise, because $G_{S}^{\eta}-G_{S}$ is supported in $B_{\eta}$.

As discussed at the beginning of Section 4, we expect the limit of $\int_{{\mathbb{R}}^{3}}|{\nabla}h_{N}^{\eta}|^{2}$ to be described in terms of quantities of the form

[TABLE]

where $K_{R}:=(-\frac{R}{2},\frac{R}{2})^{3}$, for various $\Lambda$ and solutions $H^{\eta}$ of (4.5). Broadly, the energy concentrated locally around a point $x$ should be understood from a blow-up of the original system (3.26), zooming at scale $N^{-1/3}$ around $x$. Let $x\in\mathcal{O}$ (the center of the blow-up), and $\eta_{N}:=\eta N^{-1/3}$, for a fixed $\eta>0$. If we introduce

[TABLE]

we find that

[TABLE]

System (4.5) corresponds to a formal asymptotics where one replaces $\sum_{i=1}^{N}\delta_{z_{i,N}}$ by $\sum_{i=1}^{\infty}\delta_{z_{i}}$, with $\Lambda=\{z_{i}\}$ a point distribution. Note that, under (H2), we expect this point distribution to be well-separated, meaning that there is $c>0$ such that: for all $z^{\prime}\neq z\in\Lambda$, $|z^{\prime}-z|\geq c$. Still, we insist that this asymptotics is purely formal and requires much more to be made rigorous. Such rigorous asymptotics will be carried in Section 5 for various classes of point configurations.

We now collect several general remarks on the blown-up system (4.3). We start by defining a renormalized energy. For any $L>0$, we denote $K_{L}:=(-\frac{L}{2},\frac{L}{2})^{3}$.

Definition 4.1**.**

Given a point distribution $\Lambda$, we say that a solution $H$ of (4.3) is admissible if for all $\eta>0$, the field $H^{\eta}$ defined by (4.4) satisfies ${\nabla}H^{\eta}\in L^{2}_{loc}({\mathbb{R}}^{3})$.

Given an admissible solution $H$ and $\eta>0$, we say that $H^{\eta}$ is of finite renormalized energy if

[TABLE]

exists in ${\mathbb{R}}$. We say that $H$ is of finite renormalized energy if $H^{\eta}$ is for all $\eta$, and

[TABLE]

exists in ${\mathbb{R}}$.

Remark 4.2*.*

From formula (4.4), it is easily seen that $H$ is admissible if and only if there exists one $\eta>0$ with ${\nabla}H^{\eta}\in L^{2}_{loc}({\mathbb{R}}^{3})$.

Proposition 4.3**.**

If $H_{1}$ and $H_{2}$ are admissible solutions of (4.3) satisfying for some $\eta>0$:

[TABLE]

then ${\nabla}H_{1}$ and ${\nabla}H_{2}$ differ from a constant matrix.

Proof. We set $H:=H_{1}-H_{2}=H_{1}^{\eta}-H_{2}^{\eta}$. It is a solution of the homogeneous Stokes equation with

[TABLE]

By standard elliptic regularity, any solution $v$ of the Stokes equation in the unit ball:

[TABLE]

satisfies for some absolute constant $C$,

[TABLE]

We apply this inequality to $v(x)=H(x_{0}+Rx)$, $x_{0}$ arbitrary. After rescaling, we find that

[TABLE]

As $R\rightarrow+\infty$, the right hand-side goes to zero, which concludes the proof.

Proposition 4.4**.**

Let $\Lambda$ be a well-separated point distribution, meaning there exists $c>0$ such that for all $z^{\prime}\neq z\in\Lambda$, $|z^{\prime}-z|\geq c$. Let $0<\alpha<\eta<\frac{\min(c,1)}{4}$. Let $H$ be an admissible solution of (4.3) such that $H^{\eta}$ is of finite renormalized energy. Then, $H^{\alpha}$ is also of finite renormalized energy, and

[TABLE]

In particular, $H$ is of finite renormalized energy as soon as $H^{\eta}$ is for some $\eta\in(0,\frac{c}{4})$, and $\mathcal{W}({\nabla}H)=\mathcal{W}^{\eta}({\nabla}H)$ for all $\eta<\frac{\min(c,1)}{4}$.

Proof. Let $R>0$. As $\Lambda$ is well-separated,

[TABLE]

From this and the fact that the limit $\mathcal{W}^{\eta}({\nabla}H)$ exists (in ${\mathbb{R}}$), it follows that

[TABLE]

Let $\Omega_{R}$ be an open set such that $K_{R-1}\subset\Omega_{R}\subset K_{R}$ and such that

[TABLE]

where $c^{\prime}$ depends on $c$ only. This implies that $G^{\eta}(\cdot+z)$, $G^{\alpha}(\cdot+z)$ are smooth at ${\partial}\Omega_{R}$ for all $z\in\Lambda$, and that $H^{\eta}$, $H^{\alpha}$ are smooth at ${\partial}\Omega_{R}$.

We now proceed as in the proof of Proposition 3.8. We write

[TABLE]

[TABLE]

After integration by parts, and manipulations similar to those used to show Proposition 3.8, we end up with

[TABLE]

Let us emphasize that the contribution of the boundary terms at ${\partial}\Omega_{R}$ is zero: indeed, thanks to (4.10), $(G_{S}^{\eta}-G_{S}^{\alpha})(\cdot+z)$ is zero at ${\partial}\Omega_{R}$ for any $z\in\Lambda$. Similarly,

[TABLE]

The integral in the right-hand side was computed above, see (3.29) and the lines after:

[TABLE]

Back to (4.11), we find

[TABLE]

We deduce from this identity, (4.8) and (4.9) that

[TABLE]

and replacing $R$ by $R+1$:

[TABLE]

As $\Omega_{R}\subset K_{R}\subset\Omega_{R+1}$, the result follows.

4.2 Resolution of the blown-up system for stationary point processes

As pointed out several times, we follow the strategy described in [41] for the treatment of minimizers and minima of Coulomb energies. But in our effective viscosity problem, the points $x_{i,N}$ do not minimize the analogue $\mathcal{V}_{N}$ of the Coulomb energy $\mathcal{H}_{N}$. Actually, although we consider the steady Stokes equation, our point distribution may be time dependent. More precisely, in many settings, the dynamics of the suspension evolves on a timescale associated with viscous transport (scaling like $a^{2}$, with $a$ the radius of the particle), which is much smaller than the convective time scale (scaling like $a$). This allows to neglect the time derivative in the Stokes equation: system (1.1)-(1.2) corresponds then to a snapshot of the flow at a given time $t$. Even when one is interested in the long time behaviour, the existence of an equilibrium measure for the system of particles is a very difficult problem. To bypass this issue, a usual point of view in the physics literature is to assume that the distribution of points is given by a stationary random process (whose refined description is an issue per se).

We will follow this point of view here, and introduce a class of random point processes for which we can solve (4.3). Let $X={\mathbb{R}}$ or $X={\mathbb{T}}_{L}:={\mathbb{R}}/(L{\mathbb{Z}})$ for some $L>0$. We denote by $Point_{X}$ the set of point distributions in $X^{3}$: an element of $Point_{X}$ is a locally finite subset of $X^{3}$, in particular a finite subset when $X={\mathbb{T}}_{L}$. We endow $Point_{X}$ with the smallest $\sigma$-algebra $\mathcal{P}_{X}$ which makes measurable all the mappings

[TABLE]

Given a probability space $(\Omega,\mathcal{A},P)$, a random point process $\Lambda$ with values in $X^{3}$ is a measurable map from $\Omega$ to $Point_{X}$, see [12]. By pushing forward the probability $P$ with $\Lambda$, we can always assume that the process is in canonical form, that is $\Omega=Point_{X}$, $\mathcal{A}=\mathcal{P}_{X}$, and $\Lambda(\omega)=\omega$.

We shall consider processes that, once in canonical form, are

(P1)

stationary: the probability $P$ on $\Omega$ is invariant by the shifts

[TABLE]

(P2)

ergodic: if $A\in\mathcal{A}$ satisfies $\tau_{y}(A)=A$ for all $y$, then $P(A)=0$ or $P(A)=1$.

(P3)

uniformly well-separated: we mean that there exists $c>0$ such that almost surely, $|z-z^{\prime}|\geq c$ for all $z\neq z^{\prime}$ in $\omega$.

These properties are satisfied in two important contexts:

Example 4.5* (Periodic point distributions).*

Namely, for $L>0$, $a_{1},\dots,a_{M}$ in $K_{L}$, we introduce the set $\Lambda_{0}:=\{a_{1},\dots,a_{M}\}+L{\mathbb{Z}}^{d}$. We can of course identify $\Lambda_{0}$ with a point distribution in $X^{3}$ with $X={\mathbb{T}}_{L}$. We then take $\Omega={\mathbb{T}}_{L}^{3}$, $P$ the normalized Lebesgue measure on ${\mathbb{T}}^{3}_{L}$, and set $\Lambda(\omega):=\Lambda_{0}+\omega$. It is easily checked that this random process satisfies all assumptions. Moreover, a realization of this process is a translate of the initial periodic point distribution $\Lambda_{0}$. By translation, the almost sure results that we will show below (well-posedness of the blown-up system, convergence of $\mathcal{W}_{N}$) will actually yield results for $\Lambda_{0}$ itself.

Example 4.6* (Poisson hard core processes).*

These processes are obtained from Poisson point processes, by removing balls in order to guarantee the hypothesis (P3). For instance, given $c>0$, one can remove from the Poisson process all points $z$ which are not alone in $B(z,c)$. This leads to the so-called Matérn I hard-core process. To increase the density of points while keeping (P3), one can refine the removal process in the following way: for each point $z$ of the Poisson process, one associates an ”age” $u_{z}$, with $(u_{z})$ a family of i.i.d. variables, uniform over $(0,1)$. Then, one retains only the points $z$ that are (strictly) the ”oldest” in $B(z,c)$. This leads to the so-called Matérn II hard-core process. Obviously, these two processes satisfy (P1) by stationarity of the Poisson process, and satisfy (P2) because they have only short range of correlations. For much more on hard core processes, we refer to [8].

The point is now to solve almost surely the blown-up system (4.3) for point processes with properties (P1)-(P2)-(P3). We first state

Proposition 4.7**.**

Let $\Lambda=\Lambda(\omega)$ a random point process with properties (P1)-(P2)-(P3). Let $\eta>0$. For almost every $\omega$, there exists a solution ${\bf H}^{\eta}(\omega,\cdot)$ of (4.5) in $H^{1}_{loc}(X^{3})$ such that

[TABLE]

where $D_{\bf H}^{\eta}\in L^{2}(\Omega)$ is the unique solution of the variational formulation (4.12) below.

Remark 4.8*.*

In the case $X={\mathbb{T}}_{L}$, point distributions and solutions $H^{\eta}$ over $X^{3}$ can be identified with $L{\mathbb{Z}}^{3}$-periodic point distributions and $L{\mathbb{Z}}^{3}$-periodic solutions defined on ${\mathbb{R}}^{3}$. This identification is implicit here and in all that follows.

Proof. We treat the case $X={\mathbb{R}}$, the case $X={\mathbb{T}}_{L}$ follows the same approach. We remind that the process is in canonical form: $\Omega=Point_{\mathbb{R}}$, $\mathcal{A}=\mathcal{P}_{\mathbb{R}}$, $\Lambda(\omega)=\omega$. The idea is to associate to (4.5) a probabilistic variational formulation. This approach is inspired by works of Kozlov [26, 7], see also [3]. Prior to the statement of this variational formulation, we introduce some vocabulary and functional spaces. First, for any ${\mathbb{R}}^{d}$-valued measurable $\phi=\phi(\omega)$, we call a realization of $\phi$ an application

[TABLE]

For $p\in[1,+\infty)$, $\phi\in L^{p}(\Omega)$, as $\tau_{y}$ is measure preserving, we have for all $R>0$ that $\mathbb{E}\int_{K_{R}}|R_{\omega}[\phi]|^{p}=R^{3}\,\mathbb{E}|\phi|^{p}$. Hence, almost surely, $R_{\omega}[\phi]$ is in $L^{p}_{loc}({\mathbb{R}}^{3})$. Also, for $\phi\in L^{\infty}(\Omega)$, one finds that almost surely $R_{\omega}[\phi]\in L^{\infty}_{loc}({\mathbb{R}}^{3})$. It is a consequence of Fatou’s lemma: for all $R>0$,

[TABLE]

We say that $\phi$ is smooth if, almost surely, $R_{\omega}[\phi]$ is. For a smooth function $\phi$, we can define its stochastic gradient ${\nabla}_{\omega}\phi$ by the formula

[TABLE]

where here and below, ${\nabla}={\nabla}_{y}$ refers to the usual gradient (in space). Note that ${\nabla}_{\omega}\phi(\tau_{y}\omega)={\nabla}R_{\omega}[\phi](y)$. One can define similarly the stochastic divergence, curl, etc, and reiterate to define partial stochastic derivatives ${\partial}^{\alpha}_{\omega}$.

Starting from a function $V\in L^{p}(\Omega)$, $p\in[1,+\infty]$ one can build smooth functions through convolution. Namely, for $\rho\in C^{\infty}_{c}({\mathbb{R}}^{3})$, one can define

[TABLE]

which is easily seen to be in $L^{p}(\Omega)$, as

[TABLE]

using that $\tau_{y}$ is measure-preserving. Moreover, it is smooth: we leave to the reader to check

[TABLE]

We are now ready to introduce the functional spaces we need. We set

[TABLE]

We remind that $S^{\eta}=\hbox{div \!}\Psi^{\eta}$, with $\Psi^{\eta}$ defined in (3.21). We introduce

[TABLE]

Note that it is well-defined, as $\Psi^{\eta}$ is supported in $B_{\eta}$ and $\omega$ is a discrete subset. It is measurable: indeed, $\Psi^{\eta}$ is the pointwise limit of a sequence of simple functions of the form $\sum_{i}\alpha_{i}1_{A_{i}}$, where $A_{i}$ are Borel subsets of ${\mathbb{R}}^{3}$. As

[TABLE]

is measurable by definition of the $\sigma$-algebra $\mathcal{A}$, we find that $\mathsf{\Pi}^{\eta}$ is. Moreover, as $\Lambda$ is uniformly well-separated, one has $|\mathsf{\Pi}^{\eta}(\omega)|\leq C\|\Psi^{\eta}\|_{L^{\infty}}$ for a constant $C$ that does not depend on $\omega$, so that $\mathsf{\Pi}^{\eta}$ belongs to $L^{\infty}(\Omega)$.

We now introduce the variational formulation: find $D_{\bf H}^{\eta}\in\mathcal{V}_{\sigma}$ such that for all $D_{\phi}\in\mathcal{V}_{\sigma}$,

[TABLE]

As $\mathcal{V}_{\sigma}$ is a closed subspace of $L^{2}(\Omega)$, existence and uniqueness of a solution comes from the Riesz theorem.

It remains to build a solution of (4.5) almost surely, based on $D_{\bf H}^{\eta}$. Let $\phi_{k}=\phi_{k}(\omega)$ a sequence in $\mathcal{D}_{\sigma}$ such that ${\nabla}_{\omega}\phi_{k}$ converges to $D_{\bf H}^{\eta}$ in $L^{2}(\Omega)$. Let $\rho\in C^{\infty}_{c}({\mathbb{R}}^{3})$. It is easily seen that $\rho\star\phi_{k}$ also belongs to $\mathcal{D}_{\sigma}$ and that ${\partial}^{\alpha}_{\omega}{\nabla}_{\omega}(\rho\star\phi_{k})={\partial}^{\alpha}_{\omega}(\rho\star{\nabla}_{\omega}\phi_{k})$ converges to the smooth function ${\partial}^{\alpha}_{\omega}(\rho\star D_{\bf H}^{\eta})$ in $L^{2}(\Omega)$, for all $\alpha$. In particular, as ${\nabla}_{\omega}\times{\nabla}_{\omega}(\rho\star\phi_{k})=0$, we find that ${\nabla}_{\omega}\times(\rho\star D_{\bf H}^{\eta})=0$. Applying the realization operator $R_{\omega}$, we deduce that

[TABLE]

We recall that $R_{\omega}[D_{\bf H}^{\eta}]$ belongs almost surely to $L^{2}_{loc}({\mathbb{R}}^{3})$, so that ${\nabla}\times R_{\omega}[D_{\bf H}^{\eta}]$ is well-defined in $H^{-1}_{loc}({\mathbb{R}}^{3})$. Taking $\rho=\rho_{n}$ an approximation of the identity, and sending $n$ to infinity, we end up with ${\nabla}\times R_{\omega}[D_{\bf H}^{\eta}]=0$ in ${\mathbb{R}}^{3}$. As curl-free vector fields on ${\mathbb{R}}^{3}$ are gradients, it follows that almost surely, there exists ${\bf H}^{\eta}={\bf H}^{\eta}(\omega,y)$ with

[TABLE]

In the case $X={\mathbb{T}}_{L}$, one can show that the mean of $R_{\omega}[D_{\bf H}^{\eta}]$ is almost surely zero, so that the same result holds. Besides, because the matrices ${\nabla}_{\omega}\phi$, $\phi\in\mathcal{D}_{\sigma}$, have zero trace, the same holds for $D_{\bf H}^{\eta}$. Hence,

[TABLE]

One still has to prove that the first equation of (4.5) is satisfied. Therefore, we use (4.12) with test function $D_{\phi}={\nabla}_{\omega}\phi$, where the smooth function $\phi$ is of the form

[TABLE]

Note that for smooth functions $\varphi,\tilde{\varphi}$, a stochastic integration by parts formula holds:

[TABLE]

Thanks to this formula, we may write

[TABLE]

Similarly, we find

[TABLE]

As this identity is valid for all smooth test fields $\varphi$, we end up with

[TABLE]

Proceeding as above, we find that almost surely

[TABLE]

which can be written

[TABLE]

It follows that there exists ${\bf P}^{\eta}={\bf P}^{\eta}(\omega,y)$, such that

[TABLE]

which concludes the proof of the proposition.

Corollary 4.9**.**

For random point processes with properties (P1)-(P2)-(P3), there exists almost surely a solution $H$ of (4.3) with finite renormalized energy and such that for all $\eta>0$, the gradient field ${\nabla}H^{\eta}$, where $H^{\eta}$ is given by (4.4), coincides with the gradient field ${\nabla}{\bf H}^{\eta}$ of Proposition 4.7. Moreover,

[TABLE]

where $m:={\mathbb{E}}|\Lambda\cap K_{1}|$ is the mean intensity of the point process, the expression at the right-hand side being actually constant for $\eta$ small enough.

Proof. By the definition of the mean intensity and by property (P2), which allows to apply the ergodic theorem (cf. [12, Corollary 12.2.V]), we have almost surely

[TABLE]

Let $\eta_{0}<\frac{\min(c,1)}{4}$ fixed, and ${\bf H}^{\eta_{0}}$ given by the previous proposition. We set

[TABLE]

It is clearly an admissible solution of (4.3). By Proposition 4.4, in order to show that $H$ has almost surely finite renormalized energy, it is enough to show that for one $\eta<\frac{\min(c,1)}{4}$, almost surely, the function $H^{\eta}$ given by (4.4), namely

[TABLE]

has finite renormalized energy. This holds for $\eta=\eta_{0}$, as $H^{\eta_{0}}={\bf H}^{\eta_{0}}$ and the ergodic theorem applies. We then notice that

[TABLE]

We remark that $G^{\eta}_{S}-G_{S}^{\eta_{0}}=0$ outside $B_{\max(\eta,\eta_{0})}$, so that the sum at the r.h.s. has only a finite number of non-zero terms. In the same way as we proved that the function $\mathsf{\Pi}^{\eta}$ belongs to $L^{\infty}(\Omega)$, we get that $\sum_{z\in\omega}{\nabla}(G^{\eta}_{S}-G_{S}^{\eta_{0}})(z)$ defines an element of $L^{\infty}(\Omega)$. Hence, by the ergodic theorem, we have almost surely

[TABLE]

Combining this with (4.13) and Proposition 4.4, we obtain the formula for $\mathcal{W}({\nabla}H)$.

The last step is to prove that for all $\eta>0$, ${\nabla}H^{\eta}={\nabla}{\bf H}^{\eta}$ almost surely. As a consequence of the ergodic theorem, one has almost surely

[TABLE]

Reasoning as in the proof of Proposition 4.3, we find that their gradients differ by a constant:

[TABLE]

Applying again the ergodic theorem, we get that almost surely ${\mathbb{E}}D_{H}^{\eta}={\mathbb{E}}D_{\bf H}^{\eta}+C(\omega)$. As $D_{\bf H}^{\eta}$ belongs to $\mathcal{V}_{\sigma}$, its expectation is easily seen to be zero. To conclude, it remains to prove that ${\mathbb{E}}D_{H}^{\eta}={\mathbb{E}}\sum_{z\in\omega}{\nabla}(G_{S}^{\eta}-G_{S}^{\eta_{0}})(z)$ is zero. Using stationarity, we write, for all $R>0$,

[TABLE]

We remark that for all $z$ outside a $\max(\eta,\eta_{0})$-neighborhood of ${\partial}K_{R}$, $\int_{K_{R}}{\nabla}(G_{S}^{\eta}-G_{S}^{\eta_{0}})(z+\cdot)=\int_{{\partial}K_{R}}n\otimes(G_{S}^{\eta}-G_{S}^{\eta_{0}})(z+\cdot)=0$. It follows from the separation assumption and the $L^{\infty}$ bound on ${\nabla}(G_{S}^{\eta}-G_{S}^{\eta_{0}})$ that

[TABLE]

5 Convergence of $\mathcal{V}_{N}$

This section concludes our analysis of the quadratic correction to the effective viscosity. From Theorem 1.1, we know that this quadratic correction should be given by the limit of $\mathcal{V}_{N}$ as $N$ goes to infinity, where $\mathcal{V}_{N}$ was introduced in (1.13). We show here that the functional $\mathcal{V}_{N}$ has indeed a limit, when the particles are given by the kind of stationary point processes seen in Section 4.

5.1 Proof of convergence

Let ${\varepsilon}>0$ a small parameter, and $\Lambda=\Lambda(\omega)$ a random point process with properties (P1)-(P2)-(P3): stationarity, ergodicity, and uniform separation. As seen in Examples 4.5 and 4.6, this setting covers the case of periodic patterns of points as well as classical hard core processes. We set $N=N({\varepsilon})$ the cardinal of the set

[TABLE]

where $\check{\Lambda}:=-\Lambda$ and where we label the elements arbitrarily. Note that $N$ depends on $\omega$, although it does not appear explicitly. From the fact that $\Lambda$ is uniformly well-separated and from the ergodic theorem (cf. [12, Corollary 12.2.V]), we can deduce that almost surely,

[TABLE]

so that we shall note indifferently $\lim_{{\varepsilon}\rightarrow 0}$ or $\lim_{N\rightarrow+\infty}$. Note that, strictly speaking, $N=N({\varepsilon})$ does not necessarily cover all integer values when ${\varepsilon}\rightarrow 0$, but this is no difficulty.

More generally, for all $\varphi$ smooth and compactly supported in ${\mathbb{R}}^{3}$, ergodicity implies

[TABLE]

which shows that (H1) is satisfied with $f=\frac{1}{|\mathcal{O}|}1_{\mathcal{O}}$. The hypothesis (H2) is also trivially satisfied, as well as (3.1). Our main theorem is

Theorem 5.1**.**

Almost surely,

[TABLE]

with $m$ the mean intensity of the process, and $H$ the solution of (4.3) given in Corollary 4.9.

The rest of the paragraph is dedicated to the proof of this theorem.

Let $\eta$ satisfying $\eta<\frac{\min(c,1)}{4}$ and $\eta<\frac{c}{2}(m|\mathcal{O}|)^{-1/3}$. By (5.1), it follows that almost surely, for ${\varepsilon}$ small enough, ${\varepsilon}\eta<\frac{c}{2}N^{-1/3}$. By Corollary 3.9,

[TABLE]

We denote $h_{\varepsilon}^{\eta}:=h_{N}^{\eta{\varepsilon}}$, see (3.25)-(3.26). Let $H$ be the solution of the blown-up system (4.3) provided by Corollary 4.9, $H^{\eta}$ given in (4.4), and $P^{\eta}$ as in (4.5). We define new fields $\bar{h}^{\eta}_{\varepsilon},\bar{p}^{\eta}_{\varepsilon}$ by the following conditions: $\bar{h}^{\eta}_{\varepsilon}\in\dot{H}^{1}({\mathbb{R}}^{3})$,

[TABLE]

We omit to indicate the dependence in $\omega$ to lighten notations. We claim:

Proposition 5.2**.**

[TABLE]

Proposition 5.3**.**

[TABLE]

Note that, by Proposition 4.4 and our choice of $\eta$, $\mathcal{W}^{\eta}({\nabla}H)$ = $\mathcal{W}({\nabla}H)$. Theorem 5.1 follows directly from this fact, (5.2) and the propositions.

Proof of Proposition 5.2. We know from Corollary 4.9 that

[TABLE]

From this and relation (5.1), we see that the proposition amounts to the statement

[TABLE]

A simple application of the ergodic theorem shows that almost surely

[TABLE]

It remains to show that

[TABLE]

It will be deduced from the well-known fact that the Stokes solution $\bar{h}^{\eta}_{\varepsilon}$ minimizes

[TABLE]

among divergence-free fields $\bar{h}$ in ${\rm ext}\,\mathcal{O}$ satisfying the Dirichlet condition $\bar{h}|_{{\partial}\mathcal{O}}=\bar{h}^{\eta}_{\varepsilon}|_{{\partial}\mathcal{O}}$.

First, we prove that the $H^{1/2}(\partial\mathcal{O})$-norm of ${\varepsilon}^{3}\bar{h}^{\eta}_{\varepsilon}$ goes to zero. In this perspective, we introduce for all $\delta>0$ a function $\chi_{\delta}$ with $\chi_{\delta}=1$ in a $\frac{\delta}{2}$-neighborhood of ${\partial}\mathcal{O}$, $\chi_{\delta}=0$ outside a $\delta$-neighborhood of ${\partial}\mathcal{O}$. We write

[TABLE]

By the ergodic theorem and Corollary 4.9, ${\varepsilon}^{3}{\nabla}\bar{h}^{\eta}_{\varepsilon}={\nabla}_{y}H^{\eta}(\frac{\cdot}{{\varepsilon}})$ converges almost surely weakly in $L^{2}(\mathcal{O})$ to ${\mathbb{E}}D_{\mathbf{H}^{\eta}}=0$. Let $\varphi\in L^{2}(\mathcal{O})$. By standard results on the divergence operator, cf [16], there exists $v\in H^{1}_{0}(\mathcal{O})$ with \hbox{div \!}v=\varphi-\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\mathcal{O}}\varphi, $\|v\|_{H^{1}(\mathcal{O})}\leq C_{\mathcal{O}}\|\varphi\|_{L^{2}(\mathcal{O})}$. As by definition $\bar{h}^{\eta}_{\varepsilon}$ has zero mean over $\mathcal{O}$, it follows that

[TABLE]

Hence, ${\varepsilon}^{3}\bar{h}^{\eta}_{\varepsilon}$ converges weakly to zero in $H^{1}(\mathcal{O})$ and therefore strongly in $L^{2}(\mathcal{O})$. It follows that for any given $\delta$,

[TABLE]

To conclude, it is enough to show that $\limsup_{{\varepsilon}\rightarrow 0}\|{\varepsilon}^{3}{\nabla}\bar{h}^{\eta}_{\varepsilon}\chi_{\delta}\|_{L^{2}(\mathcal{O})}$ goes to zero as $\delta\rightarrow 0$. This comes from

[TABLE]

Finally, $\|{\varepsilon}^{3}\bar{h}^{\eta}_{\varepsilon}\|_{H^{1/2}({\partial}\mathcal{O})}=o(1)$. To conclude that (5.3) holds, we notice that

[TABLE]

By classical results on the right inverse of the divergence operator, see [16], one can find for $R$ such that $\mathcal{O}\Subset B(0,R)$ a solution $\bar{h}$ of the equation

[TABLE]

and such that

[TABLE]

Extending $\bar{h}$ by zero outside $B(0,R)$, we find

[TABLE]

This concludes the proof of the proposition.

Proof of Proposition 5.3. Let $h:=h_{\varepsilon}^{\eta}-\bar{h}_{\varepsilon}^{\eta}$. It satisfies an equation of the form

[TABLE]

where the various source terms will now be defined. First,

[TABLE]

Here, the value of the stress is taken from ${\rm ext}\,\mathcal{O}$, $n$ refers to the normal vector pointing outward $\mathcal{O}$ and $s_{{\partial}}$ refers to the surface measure on ${\partial}\mathcal{O}$. We remind that $\bar{h}^{\eta}_{\varepsilon}\in\dot{H}^{1}({\mathbb{R}}^{3})$ does not jump at the boundary, but its derivatives do, so that one must specify from which side the stress is considered. Then,

[TABLE]

with the value of the stress taken from $\mathcal{O}$, and $n$ as before. Noticing that $S{\nabla}f=-\frac{1}{|\mathcal{O}|}Sn\,s_{\partial}$, we finally set

[TABLE]

where

[TABLE]

Note that the term $R_{3}$ is supported in pieces of spheres. From (3.20), we know that for all $\eta>0$

[TABLE]

This allows to show that the integral of $R_{2}+R_{3}$ is zero. Indeed,

[TABLE]

so that

[TABLE]

The point is now to prove that ${\varepsilon}^{3}\|{\nabla}h\|_{L^{2}({\mathbb{R}}^{3})}\rightarrow 0$ as ${\varepsilon}\rightarrow 0$. From a simple energy estimate, and taking (5.7) into account, we find

[TABLE]

As $(\bar{h}_{\varepsilon}^{\eta},\bar{p}_{\varepsilon}^{\eta})$ is a solution of a homogeneous Stokes equation in ${\rm ext}\,\mathcal{O}$, we get from an integration by parts:

[TABLE]

using the Cauchy-Schwarz inequality and the bound (5.5).

We now wish to show that

[TABLE]

for some $\nu({\varepsilon})$ going to zero with ${\varepsilon}$. More precisely, we will prove that for any divergence-free $\varphi\in\dot{H}^{1}({\mathbb{R}}^{3})$,

[TABLE]

which implies (5.10) by Poincaré inequality. We first notice that

[TABLE]

where

[TABLE]

Then, we use the relation $S^{\eta}=\hbox{div \!}\Psi^{\eta}$, cf. Lemma 3.5 and integrate by parts to get

[TABLE]

For a fixed $\eta$, there is a constant $C$ (depending on $\eta$) such that

[TABLE]

For the last inequality, we have used that all $x_{i}$’s with $i\in I_{\varepsilon}^{\eta}$ belong to an ${\varepsilon}$-neighborhood of ${\partial}\mathcal{O}$, so that $|I_{\varepsilon}^{\eta}|=O({\varepsilon}^{-2})$. Hence,

[TABLE]

Let

[TABLE]

We claim that ${\mathbb{E}}\int_{K_{1}}F_{3}=0$. Indeed, by stationarity, for all $R>0$

[TABLE]

We have used crucially the fact that $\Psi^{\eta}$ is supported in $B(0,\eta)$. The $O(\frac{1}{R})$-term is associated to the points $z\in\Lambda$ which lie in a $\delta$-neighborhood of ${\partial}K_{R}$: see the end of the proof of Corollary 4.9 for similar reasoning. By sending $R$ to infinity, we find that almost surely

[TABLE]

The last step is to compute $\int_{B(0,\eta)}\Psi^{\eta}$, which is independent of $\eta$ by homogeneity. It is in particular equal to $\lim_{\eta\rightarrow 0}\langle\Psi^{\eta},1\rangle$, a limit that was already computed in the proof of Lemma 3.5, cf. (3.23)-(3.24). We get $\int_{B(0,\eta)}\Psi^{\eta}=S$, which shows that ${\mathbb{E}}F_{3}=0$.

By the definition of $F^{\varepsilon}_{3}$, we can write

[TABLE]

where the last inequality follows from (5.1). Plugging this inequality in (5.14), and combining with (5.12), we see that to derive (5.11), it remains to show that almost surely, for all divergence-free fields $\varphi\in H^{1}(\mathcal{O})$,

[TABLE]

where $F^{\varepsilon}:=F_{2}^{\varepsilon}+F_{3}^{\varepsilon}$. Notice that $\hbox{div \!}(F_{2}^{\varepsilon}+F_{3}^{\varepsilon})=0$. We introduce again the functions $\chi_{\delta}$, $\delta>0$, seen above. We get

[TABLE]

For the last term, we take into account that $\varphi$ is divergence-free, so that the pressure disappears: we find

[TABLE]

As seen in (5.4), we have

[TABLE]

and similarly,

[TABLE]

For the first term, we write

[TABLE]

We know that ${\varepsilon}^{3}D(\bar{h}^{\eta}_{\varepsilon})$ goes weakly to zero in $L^{2}(\mathcal{O})$, so that it converges strongly to zero in $H^{-1}(\mathcal{O})$. As ${\nabla}\chi_{\delta}\otimes\varphi$ belongs to $H^{1}_{0}(\mathcal{O})$, we find that for a fixed $\delta$

[TABLE]

Similarly, as ${\mathbb{E}}F_{3}=0$, $F_{3}^{\varepsilon}$ converges weakly to zero in $L^{2}(\mathcal{O})$ and we get

[TABLE]

The last step is to prove that ${\varepsilon}^{3}\overline{p}^{\eta}_{\varepsilon}$ converges weakly to zero in $L^{2}(\mathcal{O})$, which will yield

[TABLE]

As above, for $\phi\in L^{2}(\mathcal{O})$, we introduce $v\in H^{1}_{0}(\mathcal{O})$ such that \hbox{div \!}v=\phi-\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\mathcal{O}}\phi, $\|v\|_{H_{1}(\mathcal{O})}\leq C_{\mathcal{O}}\|\phi\|_{L^{2}(\mathcal{O})}$. Then, using the equation satisfied by $\overline{p}_{\varepsilon}^{\eta}$ in $\mathcal{O}$:

[TABLE]

we find after integration by parts

[TABLE]

This concludes the proof of (5.16), of Proposition 5.3 and of the theorem.

5.2 Formula for periodic point distributions

Theorem 5.1 gives the limit of $\mathcal{V}_{N}$ for properly rescaled stationary and ergodic point processes, under uniform separation of the points. Such setting includes periodic point distributions, as well as Poisson hard core processes. We focus here on the periodic case, for which further explicit formula can be given. For $L>0$, we consider distinct points $a_{1},\dots,a_{M}$ in $K_{L}$, and set $\Lambda_{0}:=\{a_{1},\dots,a_{M}\}+L{\mathbb{Z}}^{d}$, which can be seen as a subset of ${\mathbb{T}}_{L}^{3}$. In Example 4.5, we explained how to build a process on ${\mathbb{T}}_{L}^{3}$ out of $\Lambda_{0}$, with $\Lambda(\omega)=\Lambda_{0}+\omega$, $\omega\in{\mathbb{T}}_{L}^{3}$. By a simple translation, the results above, that are valid for $\Lambda_{0}+\omega$ for a.e. $\omega$, are still valid for $\omega=0$. Thus, for $\Lambda=\Lambda_{0}$, we deduce from Proposition 4.7 the existence of an $L{\mathbb{Z}}^{3}$-periodic solution ${\bf H^{\eta}}$ of (4.5) with ${\nabla}{\bf H^{\eta}}\in L^{2}_{loc}$. If we further assume that ${\bf H^{\eta}}$ is mean-free, it is clearly unique. Then, following Corollary 4.9 and Theorem 5.1, there exists an $L{\mathbb{Z}}^{3}$-periodic solution $H$ of (4.3), such that

[TABLE]

where $H^{\eta}$ is associated to $H$ by (4.4). We have used that in the periodic case, the intensity of the process is $m=\frac{M}{L^{3}}$, while the expectation is simply the average over $K_{L}$.

To make things more explicit, we introduce the periodic Green function $G_{S,L}:{\mathbb{R}}^{3}\rightarrow{\mathbb{R}}^{3}$ satisfying:

[TABLE]

The Green function $G_{S,L}$ is easily expressed in Fourier series. If we write

[TABLE]

a straightforward calculation shows that for all $k\in{\mathbb{Z}}^{3}_{*}$

[TABLE]

where $\pi^{\perp}_{k}$ denotes the projection orthogonally to the line ${\mathbb{R}}k$. Note that the Fourier series for $G_{S,L}$ converges for instance in the quadratic sense.

Proposition 5.4**.**

[TABLE]

Proof. Clearly, the $L{\mathbb{Z}}^{3}$-periodic field defined on $K_{L}$ by $\tilde{H}(y):=\sum_{i=1}^{M}G_{S,L}(y+a_{i})$ is a solution of (4.3), and by Proposition 4.3 ${\nabla}\tilde{H}$ and ${\nabla}H$ differ from a constant matrix. As ${\nabla}(\tilde{H}-H)={\nabla}(\tilde{H}^{\eta}-H^{\eta})$ is the gradient of a periodic function, we have eventually ${\nabla}\tilde{H}={\nabla}H$. Up to adding a constant field to $H$, we can assume that

[TABLE]

Then, if $\eta$ is small enough so that $B(a_{i},\eta)\subset K_{L}$ for all $i$, $H^{\eta}$ is the $L$-periodic field given on $K_{L}$ by

[TABLE]

We integrate by parts to find

[TABLE]

where we have used that the last term of the second line vanishes identically. We then write $G_{S,L}=G_{S}+\phi_{S,L}$ with $\phi_{S,L}$ smooth near [math] to obtain

[TABLE]

Combining with Lemma 3.6 and (LABEL:lim_VN_periodic), we get

[TABLE]

We conclude by the last point of Lemma 3.5 that

[TABLE]

Proposition 5.5** **(Simple cubic lattice).

In the special case where $L=M=1$, we find

[TABLE]

with $\alpha=\frac{5}{2}(1-60a)$, $\beta=\frac{5}{2}(1+40a)$, and $a\approx-0,04655$ is defined in (5.19).

Proof. When $M=L=1$, the formula from the last proposition simplifies into $\lim_{N}\mathcal{V}_{N}=\frac{25}{2}S{\nabla}\cdot\phi_{S,1}(0)$, with $\phi_{S,1}=G_{S,1}-G_{S}$. The periodic Green function $G_{S,1}$ was computed using Fourier series in the last paragraph. We found

[TABLE]

We use formulas from [20], see also [42, equations (64)-(65)]:

[TABLE]

and

[TABLE]

where $c_{1}$ and $c_{2}$ are constants, and

[TABLE]

Note that the formula (5.19) defines implicitly $a$. A numerical computation was carried in [42], see also [34], giving $a\approx-0,04655$.

Inserting in the expression for $G_{S,1}$, we find after a tedious calculation that

[TABLE]

Note that to carry out this calculation, we used the fact that $S$ is trace-free, which leads to the identity

[TABLE]

Moreover, we know from (3.12) that

[TABLE]

We end up with

[TABLE]

and the right formula for $\lim_{N}\mathcal{V}_{N}$.

5.3 Formula in the stationary case with the $2$-point correlation function

We consider here the case of random point processes in ${\mathbb{R}}^{3}$ ($X={\mathbb{R}}$), such that (P1)-(P2)-(P3) hold. We further assume that the mean density is $m=1$. We assume moreover that this point process admits a $2$-point correlation function, that is a function $\rho_{2}=\rho_{2}(x,y)\in L^{1}_{loc}({\mathbb{R}}^{3}\times{\mathbb{R}}^{3})$ such that for all bounded set $K$ and all smooth $F$ in a neighborhood of $K$:

[TABLE]

As the process is stationary, one can write $\rho_{2}(x,y)=\rho(x-y)$. Our goal is to prove the following formula:

Proposition 5.6**.**

Almost surely,

[TABLE]

where $G_{S,L}$ refers to the $L{\mathbb{Z}}^{3}$-periodic Green function introduced in (5.18).

Remark 5.7*.*

We remind that the periodic Green function $G_{S,L}$ has singularities at each point of $L{\mathbb{Z}}^{d}$. But as the sum is restricted to points $z,z^{\prime}$ in $\Lambda\cap K_{L-1}$, $z-z^{\prime}$ is always away from this set of singularities. In the same way, the integral over $K_{L-1}\times K_{L-1}$ in the second equality is well-defined. Under further assumption on the two-point correlation function $\rho$, one could make sense of the integral over $K_{L}\times K_{L}$ and replace the former by the latter.

Proof. Let $\eta$ small enough so that Proposition 4.4 holds. We have

[TABLE]

Let $H_{L}=\sum_{i=1}^{M}G_{S,L}(\cdot+a_{i})$, where $\{a_{1},\dots,a_{M}\}=\Lambda\cap K_{L-1}$. Note that $H_{L}$ is associated to the point process $\Lambda_{L}$ obtained by $L{\mathbb{Z}}^{d}$-periodization of $\Lambda\cap K_{L-1}$. We shall prove below that,

[TABLE]

As $\frac{M}{L^{3}}=\frac{|\Lambda\cap K_{L}|}{L^{3}}\rightarrow 1$ as $L\rightarrow+\infty$, it follows from (5.20) that

[TABLE]

where the last equality comes from Proposition 4.4. One can apply such proposition because the $L{\mathbb{Z}}^{d}$-periodized network $\Lambda_{L}$ has a minimal distance between points that is independent of $L$. This is the reason why we used $K_{L-1}$ instead of $K_{L}$ in the definition of $\Lambda_{L}$. Eventually, by Proposition 5.4,

[TABLE]

Using that

[TABLE]

we get that

[TABLE]

We obtain

[TABLE]

This is the first formula of the proposition. To prove the second one, one can go back to formula (5.21) and take the expectation of both sides. The left-hand side, which is deterministic, is of course unchanged. As regards the r.h.s., one can swap the limit in $L$ and the expectation by invoking the dominated convergence theorem. Indeed, both terms $\frac{1}{L^{3}}\int_{K_{L}}|{\nabla}H^{\eta}_{L}|^{2}$ and \frac{M}{(\eta L)^{3}}\Bigl{(}\int_{B^{1}}|{\nabla}G_{S}^{1}|^{2}+\frac{3}{10\pi}|S|^{2}\Bigr{)} are bounded uniformly in $n$ and in the random parameter $\omega$ (but not uniformly on $\eta$): the first term is bounded through a simple energy estimate, while the second one is bounded thanks to the almost sure separation assumption.

The final step is to prove (5.20) almost surely. We set ${\varepsilon}:=\frac{1}{L}$, and introduce for all $x\in K_{1}$,

[TABLE]

and similarly, for all $x\in K_{1}$,

[TABLE]

where $(H^{\eta},P^{\eta})$ refers to the field built in Proposition 4.7. Clearly,

[TABLE]

while, by the ergodic theorem, one has almost surely:

[TABLE]

It remains to show that

[TABLE]

We notice that the difference $h_{\varepsilon}=\bar{h}^{\eta}_{\varepsilon}-h^{\eta}_{{\varepsilon}}$ satisfies the Stokes equation

[TABLE]

where

[TABLE]

and where we recall that $\Lambda_{L}$ is obtained by $L{\mathbb{Z}}^{3}$-periodization of $\Lambda\cap K_{L-1}$. Testing against ${\varepsilon}^{6}h_{\varepsilon}$, we find

[TABLE]

where

[TABLE]

with

[TABLE]

Note that both $F_{\varepsilon}$ and $G_{\varepsilon}$ are divergence-free.

To handle the first term at the right-hand side of (5.23), we notice that

[TABLE]

resulting in

[TABLE]

As regards the second term, one proceeds exactly as in Paragraph 5.1, replacing $\mathcal{O}$ by $K_{1}$: see the treatment of $F^{2}_{\varepsilon}$ and $F^{3}_{\varepsilon}$, defined in (5.13) and (5.15). One gets in this way that for all divergence-free $\varphi\in H^{1}(K_{1})$,

[TABLE]

As regards the last term, we take into account the periodicity of $H^{\eta}_{L}$ and $\tilde{G}$ to write

[TABLE]

As $\int_{{\partial}K_{1}}\bar{h}^{\eta}_{\varepsilon}\cdot n=0$, we can introduce a solution $\Phi_{\varepsilon}$ of

[TABLE]

Proceeding as in Paragraph 5.1 (replacing $\mathcal{O}$ by $K_{1}$), one can show that $\|{\varepsilon}^{3}\bar{h}^{\eta}_{\varepsilon}\|_{H^{1/2}({\partial}K_{1})}$ goes to zero with ${\varepsilon}$, and so $\|\Phi^{\varepsilon}\|_{H^{1}(K_{1})}$ goes to zero as well. Eventually, we write

[TABLE]

Hence, we find

[TABLE]

which concludes the proof.

Acknowledgements

We express our gratitude to Sylvia Serfaty for explaining to us her work on Coulomb gases and being a source of fruitful suggestions. We acknowledge the support of the SingFlows project, grant ANR-18-CE40-0027 of the French National Research Agency (ANR). D. G.-V. acknowledges the support of the Institut Universitaire de France. M.H. acknowledges the support of Labex Numev Convention grants ANR-10-LABX-20.

Appendix A Proof of Lemma 2.4

For any open set $U$, we denote \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{U}=\frac{1}{|U|}\int_{U}. By (H2), we have

[TABLE]

We write

[TABLE]

with

[TABLE]

Setting $y_{i}=N^{-1/3}x_{i}$, using that for $i\neq j$, $|y_{i}-y_{j}|\geq\frac{1}{2}(c+|y_{i}-y_{j}|)\geq c$,

[TABLE]

From the inequality (2.35), applied with $a_{ij}=\frac{1}{(c+|y_{i}-y_{j}|)^{4}}$ and $b_{j}=A_{j}$, we deduce

[TABLE]

Similarly,

[TABLE]

This leads to

[TABLE]

The last term is the most difficult. We follow [21]. Let us remind that

[TABLE]

Let $\chi_{d}(x)=\chi(x/d)$ a smooth function that is [math] in $B(0,d)$, $1$ outside $B(0,2d)$. Introducing the function $F_{A}=\sum_{j}A_{j}1_{B(x_{j},d)}$, using that $d\leq\min_{i\neq j}\frac{|x_{i}-x_{j}|}{4}$, we can write that

[TABLE]

where $\mathbf{K}(x)$ is an endomorphism of the space of symmetric matrices, defined by

[TABLE]

We then split $A^{\prime}_{i,3}=M_{i}+N_{i}$, with

[TABLE]

By Hölder inequality,

[TABLE]

and so

[TABLE]

The kernel $\chi_{d}\mathbf{K}$ enters the framework of the Calderón-Zygmund theorem, see for instance [31, Chapters 4 and 5]: for all $1<q<+\infty$, the operator \big{(}\chi_{d}\mathbf{K}\bigr{)}\,\star is continuous from $L^{q}({\mathbb{R}}^{3})$ to $L^{q}({\mathbb{R}}^{3})$, with

[TABLE]

We stress that the constant $C_{q}$ depends only on $q$, and not on $d$, as can be seen from the rescaling $x^{\prime}:=x^{\prime}/d$. It follows that

[TABLE]

As the balls $B(x_{j},d)$ are disjoint, $|\sum A_{j}1_{B(x_{j},d)}|^{q}=\sum|A_{j}|^{q}1_{B(x_{j},d)}$, so that $\|F_{A}\|_{L^{q}({\mathbb{R}}^{3})}^{q}=\frac{4\pi}{3}\sum|A_{j}|^{q}d^{3}$, and

[TABLE]

To bound $N_{i}$, we notice that for all $x\in B_{i}$, the support of $x^{\prime}\rightarrow\chi_{d}(x_{i}-x^{\prime})-\chi_{d}(x-x^{\prime})$ is included in

[TABLE]

(remark that by definition of $\phi$, $a$ is less than $d$ for $\phi$ small enough). We get

[TABLE]

so that

[TABLE]

using that for $\phi\ll 1$, $a\ll d$ and \big{|}\ln\big{(}\frac{2d+a}{d-a}\big{)}\big{|} is bounded by an absolute constant.