A local version of Einstein's formula for the effective viscosity of suspensions

TL;DR

This paper rigorously proves a local version of Einstein's effective viscosity formula for dilute suspensions at the level of the Stokes equation, considering a large number of particles with vanishing size and volume fraction.

Contribution

It provides the first proof of Einstein's formula at the Stokes equation level for dilute suspensions with many particles, using a dipole approximation and homogenization techniques.

Findings

Established $L^$ and $L^p$ estimates for microscopic and homogenized solutions.

Validated Einstein's formula in a local setting for suspensions with many particles.

Demonstrated the effectiveness of the dipole approximation and reflection method.

Abstract

We prove a local variant of Einstein's formula for the effective viscosity of dilute suspensions, that is $\mu^\prime=\mu (1+\frac 5 2\phi+o(\phi))$, where $\phi$ is the volume fraction of the suspended particles. Up to now rigorous justifications have only been obtained for dissipation functionals of the flow field. We prove that the formula holds on the level of the Stokes equation (with variable viscosity). We consider a regime where the number $N$ of particles suspended in the fluid goes to infinity while their size $R$ and the volume fraction $\phi=NR^3$ approach zero. We establish $L^\infty$ and $L^p$ estimates for the difference of the microscopic solution to the solution of the homogenized equation. Here we assume that the particles are contained in a bounded region and are well separated in the sense that the minimal distance is comparable to the average one. The main tools for the proof are a dipole approximation of the flow field of the suspension together with the so-called method of reflections and a coarse graining of the volume density.

Equations

$$\mu^{\prime}=\mu\Big{(}1+\frac{5}{2}\phi+o(\phi)\Big{)},$$

$$d_{ij}^N:=\left|X_i^N-X_j^N\right|\text{ for }i\ne j,d^N:=\underset{1\le i,j\le N}{\min }{d}_{ij}^N>2R^N.$$

$${\Omega }^N={\mathbb{R}}^3\setminus \bigcup _{i=1}^N\overline{B_i^N}.$$

$$\Delta u+\nabla p$$

$$\mathrm{div}u$$

$${\int }_{\partial B_i}\sigma ndS$$

$${\int }_{\partial B_i}(x-X_i)\wedge (\sigma n)dS$$

$$u(x)$$

$$\sigma =-p\mathrm{Id}+2eu,eu=\frac{1}{2}(\nabla u+\nabla u^T),$$

$$-\Delta v+\nabla p$$

$$\mathrm{div}v$$

$$v(x)$$

$$v(x)=v(X_i)+\nabla v(X_i)(x-X_i)+o(R).$$

$$d_i(x)=\{\begin{array}{ll}{ϵ}_i(x-X_i)& \text{mbox},forx\in \overline{B_i},\\ \frac{5}{2}{R}^3\frac{(x-X_i)\left((x-X_i)\cdot {ϵ}_i(x-X_i)\right)}{{|x-X_i|}^5}& \\ +R^5\left(\frac{{ϵ}_i(x-X_i)}{{|x-X_i|}^5}-\frac{5}{2}\frac{(x-X_i)\left((x-X_i)\cdot {ϵ}_i(x-X_i)\right)}{{|x-X_i|}^7}\right)& \text{mbox},otherwise.\end{array}$$

$$\tilde{u}=v-\sum _{i=1}^N{d}_i.$$

$$\tilde{u}(x)$$

$$\approx v(x)-{\int }_{{\mathbb{R}}^3}\frac{3}{4\pi }\varphi {\rho }^N(y)\frac{5}{2}\left(\frac{(x-y)\left((x-y)\cdot ev(y)(x-y)\right)}{{\left|x-y\right|}^5}\right)dy.$$

$${\Phi }_{ij}(x)=\frac{1}{8\pi }\left(\frac{{\delta }_{ij}}{\left|x\right|}+\frac{x_i{x}_j}{{\left|x\right|}^3}\right).$$

$${ϵ}_{ki}{\partial }_k{\Phi }_{ij}(x)=-\frac{3}{8\pi }\frac{x_j{x}_k{ϵ}_{ki}{x}_i}{{\left|x\right|}^5}.$$

$${\tilde{u}}_j(x)$$

$$\approx v_j(x)+{\int }_{{\mathbb{R}}^3}5\varphi \rho (y)ev(y{)}_{ki}\left({\partial }_k{\Phi }_{ij}\right)(x-y)dy$$

$$=v_j(x)+{\int }_{{\mathbb{R}}^3}{\Phi }_{ij}(x-y){\mathrm{div}}_y{\left(5\varphi \rho ev(y)\right)}_{i}dy.$$

$$-\Delta \tilde{u}+\nabla p$$

$$-\mathrm{div}\left(\nabla \tilde{u}+5\varphi \rho e\tilde{u}\right)+\nabla \tilde{p}=(1-\varphi \rho )f,\mathrm{div}\tilde{u}=0.$$

$$-\mathrm{div}\left(\left(2+5\varphi \rho \right)e\tilde{u}\right)+\nabla \tilde{p}=(1-\varphi \rho )f,\mathrm{div}\tilde{u}=0.$$

$$-\operatorname{div}\sigma=(1-\phi\rho)f,\quad\text{ where }\quad\sigma=2\Big{(}1+\frac{5}{2}\phi\rho\Big{)}e\tilde{u}-p\operatorname{Id}.$$

$$\mu_{\operatorname{eff}}=\mu\Big{(}1+\frac{5}{2}\phi\rho\Big{)}$$

$$\displaystyle-\Big{(}1+\frac{5}{2}\phi\rho\Big{)}\Delta\tilde{u}+\nabla\tilde{p}$$

$$\mathrm{div}\tilde{u}$$

$$\displaystyle\mbox{There exists $L>0$ such that }\left\lvert X_{i}\right\rvert+R<L\mbox{ for all }i=1,\dots,N.$$

$$\text{mbox}ThereexistsC>0\text{mbox}suchthat{N}^{-\frac{1}{3}}\le Cd.$$

$$\text{mbox}Weassumethat\varphi \log N\to 0\text{mbox}asN\to \infty .$$

$$f\in L^{\frac{6}{5}}({\mathbb{R}}^3)\cap L^{\infty }({\mathbb{R}}^3)$$

$$f\in C^{0,\alpha }({\mathbb{R}}^3)\text{mbox}forsome\alpha >0.$$

$${\rho }^N(x)=\frac{4\pi }{3}\frac{1}{N(s^N{)}^3}n(A_j)\text{mbox}forx\in A_j,$$

$${\rho }^N\rightharpoonup \rho \text{mbox}weakly\ast in{L}^{\infty }({\mathbb{R}}^3)\text{mbox}andinall{L}^p({\mathbb{R}}^3),p\ge 1.$$

Full text

{AMS}

35Q35, 76A99, 76D07, 76D09, 76D09

1 Introduction

In his annus mirabilis, 1905, Einstein published five seminal works contributing to different areas of physics. One of these works was his dissertation "Eine neue Bestimmung der Moleküldimensionen" [Ein06] where he derives a formula for the effective viscosity of a dilute suspension of spheres

[TABLE]

where $\phi$ is the (small) volume fraction of the spheres. He relates it to the formula for the mass diffusivity in order to obtain a formula for the size of the particles in the suspension. Applying this to a solution of sugar in water he is able to estimate the molecular dimensions of sugar, since both viscosity and diffusivity can be measured experimentally.

Today the formula and its validity is still of interest because dilute suspensions appear in a wide scope of applications in physics, chemistry and engineering where the effective properties of such mixtures play a role.

The purpose of this article is to rigorously prove Einstein’s formula on a local level. So instead of considering global energies of the flow field of the suspension we prove that the flow obeys an effective equation incorporating the (possibly non-constant) effective viscosity at every point in space. This is an improvement with regard to prior results because it makes the time evolution of the problem accessible and it provides a proof of Einstein’s formula in a non-periodic setting.

1.1 Setting of the problem

Consider a collection of rigid spherical particles $B_{i}^{N}\coloneqq B_{R^{N}}(X_{i}^{N}),i=1,..,N$ where $X_{i}^{N}\in\mathbb{R}^{3}$ are the particles’ centres and $R^{N}>0$ is the radius of all particles. Let

[TABLE]

This implies that the particles do not intersect nor touch each other. The domain of the suspending material is given by

[TABLE]

For an easier reading we will mostly drop the superscript $N$ in the further discussion.

We assume that $\Omega$ is occupied by a Stokes-fluid with viscosity $\mu=1$, that the particles are inertialess and that the fluid-solid interaction is given by a no-slip boundary condition. This entails the following problem for the fluid velocity $u:\mathbb{R}^{3}\to\mathbb{R}^{3}$:

[TABLE]

where

[TABLE]

and the $V_{i},\omega_{i}\in\mathbb{R}^{3}$ are a priori unknown and must be determined as part of the solution. Here $p$ is the pressure which is a Lagrange multiplier associated to the divergence condition and $f$ is a force density acting on the fluid. In problem (1.2)-(1.7) one can replace $f$ by $f^{N}=\mathds{1}_{\Omega}f$ where $\mathds{1}_{\Omega}$ is the characteristic function of $\Omega$ since the equation holds only in $\Omega$.

1.1.1 Heuristics

We start by a heuristic derivation of the effective flow field. To that aim we consider the solution without particles:

[TABLE]

Since we are in a low volume fraction regime it is reasonable (and will be proven in Theorem 6) to assume that $u$ is already close, in terms of the volume fraction $\phi$, to $v$. In fact it is possible to prove $\left\lVert u-v\right\rVert\lesssim\phi$ for a suitable norm. Then, in order to get a higher-order approximation of $u$, the main point is to satisfy the condition $u=V_{i}+\omega_{i}\wedge(x-X_{i})$ on the particles. On the ball $B_{i}$ the function $v$, up to first order, has the form

[TABLE]

The linear part consists of a skew-symmetric part that induces rotations and that we want to keep, while we need to correct for the symmetric part $\epsilon_{i}=ev(X_{i})$. In order to get closer to a rigid body motion, we subtract from $v$ the (dipole-)function $d_{i}$ that only incorporates the symmetric gradient and is defined by

[TABLE]

Then $v-d_{i}=v(X_{i})+\omega_{i}\wedge(x-X_{i})+o(R)$ in $B_{i}$, where $\omega_{i}$ is determined by the skew-symmetric part of the gradient. Now we want the higher-order approximation $\tilde{u}$ to be close to a rigid body motion on all the particles and set

[TABLE]

Of course for $i\neq j$ the dipole $d_{i}$ will not vanish on $B_{j}$ but since the decay of $d_{i}$ is quadratic we may hope that under some conditions on the particle distribution this effect is comparable to the one coming from higher order terms in the Taylor expansion of $v$ in $B_{i}$. A related approach would be to choose, instead of $\epsilon_{i}$, the coefficients of the dipoles in order to optimize with respect to the distance of $\tilde{u}$ to $u_{N}$, see for example [MM10] for the scalar case, [MMN16] for the vectorial case in the framework of elasticity and the review article [MM18] as well as references therein. We refrain from taking this approach here because it makes the second (homogenization) step harder. Since the $d_{i}$ solve the homogeneous Stokes equation outside $\overline{B_{i}}$ the equation $-\Delta\tilde{u}+\nabla p=f^{N}$ is valid in $\Omega$. Note that $d_{i}$ consists of two parts, one of which decays much more rapidly than the other. Hence we take into account only the first part for the following heuristics. Now let $\phi=RN^{3}$ (this is a slight abuse of notation since $\phi$ denoted the physical volume fraction before). We assume that the rescaled volume density $\rho^{N}=\frac{1}{\phi}\frac{4\pi}{3}R^{3}\sum_{i=1}^{N}\delta_{X_{i}}$, where $\delta_{X_{i}}$ is the dirac measure supported at $X_{i}$, converges in some sense to $\rho$ as $N\to\infty$ so that $\phi\rho$ is the virtual limit volume density. We can write

[TABLE]

Now we introduce the fundamental solution to the Stokes equation

[TABLE]

We will see later (part 3 of the proof of Lemma 5.3) that the following identity holds for symmetric and traceless matrices $\epsilon$, where here, and in the following we use the Einstein convention to always sum over doubly appearing indices:

[TABLE]

Using this we arrive at the following approximation:

[TABLE]

Notice that $\rho^{N}\rightharpoonup\rho$ implies that $\frac{1}{\phi}\mathds{1}_{\Omega}\rightharpoonup\rho$. This yields $f^{N}\approx(1-\phi\rho)f$ for large $N$. Using this and taking $-\Delta$ on both sides we arrive at

[TABLE]

We expect that $\tilde{u}$ is a better approximation of $u$ than $v$, in particular, since $\left\lVert u-v\right\rVert\sim\phi$, we have $\|\tilde{u}-v\|\sim\phi$. Thus, making an error of order $\phi^{2}$, we can replace $v$ by $\tilde{u}$ in the divergence term to obtain the following equation

[TABLE]

We can use the fact that $\operatorname{div}\tilde{u}=0$ (and hence $\operatorname{div}\nabla\tilde{u}^{T}=0$) to write this equation as

[TABLE]

This has the form

[TABLE]

Comparing to the form of the stress tensor for a homogeneous fluid (and keeping in mind that we rescaled the viscosity) this suggests that the effective viscosity of a suspension for small volume fractions of the immersed particles in a material of viscosity $\mu$ is given by

[TABLE]

to first order in $\phi$. Note that, since $\rho$ is typically non-constant, the effective viscosity is a function of the space variable.

In regions where the density $\rho$ is constant, the divergence acting on the part of the transposed gradient vanishes because $\operatorname{div}\tilde{u}=0$. In these regions we recover Einstein’s formula even for the classical form of the Stokes equation:

[TABLE]

The main result of this paper shows, informally stated, that $\frac{u-\tilde{u}}{\phi}$ goes to zero as $\phi\to 0$ in suitable norms. This result is part of the second author’s PhD thesis [Sch19]. Some standard results that are not mentioned in this article and extended proofs can be found in [Sch19].

The same computation as above can be done for the analogous but simpler electrostatical problem of perfectly conducting spheres suspended in a dielectric medium. The corresponding result is that the effective permittivity is given by $\eta_{\operatorname{eff}}=\eta(1+3\phi\rho)$. We also refer to [Sch19] for details.

2 The main result

In this section we state our assumptions on the particle configuration and give a precise formulation of our main result.

2.1 Assumptions

We set $\phi=NR^{3}$ and assume that the following requirements are met by the sequence of particle configurations:

[TABLE]

Note that (2.2) and (2.3) imply that the particles are well separated in the sense that $d\geq 4R$ for large $N$. We need assumption (2.3) which is a bit stronger than the minimal assumption $\phi\to 0$ as $N\to\infty$ to account for the logarithmic divergence of certain sums and integrals at the origin.

We will make the following assumptions on $f$:

[TABLE]

In order to state the result that compares the microscopic solution of problem (1.2)-(1.7) to the solutions of certain homogenized problems, it is necessary to define a limit volume density. It will prove useful to introduce a coarse grained density as in [NV06].

Definition 2.1**.**

Let $s^{N}>0$ be a sequence such that $s^{N}\log N\to 0$ as $N\to\infty$. Let $\mathbb{R}^{3}$ be decomposed into half-open disjoint cubes $A_{j}$ of side length $s^{N}$ where $j\in\mathbb{Z}^{3}$. We define the rescaled averaged particle volume density $\rho^{N}$ by

[TABLE]

*where $n(A_{j})$ is the number of particle centers $X_{i}$ in $A_{j}$. *

Note that $\phi\rho^{N}$ is the local volume density of the particles in each cube. By assumption (2.2) this vanishes in the limit $N\to\infty$ which forces us to rescale by the volume fraction, in order to obtain a quantity proportional to the number density, that does not necessarily converge to zero. Since all particles are contained in a big ball (assumption (2.1)), $\rho^{N}$ will, for large $N$ be compactly supported in $B_{L+1}(0)$. By assumption (2.2) we have that $\rho^{N}$ is uniformly bounded in $L^{\infty}$. For the following we will assume without loss of generality that

[TABLE]

In addition we need the regularity assumption

[TABLE]

which is needed, because in equation (1.13) the derivative coming from the divergence might fall onto the (non-constant) density. This produces a lower order term that can only be controlled with assumption (2.8).

As the dipoles used for the approximation of the microscopic problem are singular, we are forced to define a domain that cuts out a boundary layer around the particles. Let $\delta^{N}>0$ such that $\frac{1}{\left({\delta^{N}}\right)^{2}N}\to 0$ and $\frac{\delta}{d}\to 0$ as $N\to\infty$. In particular $N^{-\frac{1}{2}}\leq C\delta^{N}\leq C^{\prime}N^{-\frac{1}{3}}$. Then we define $r=\max(2R,\delta^{N})$ and introduce

[TABLE]

2.2 Statement of main result

For the following let $\dot{H}^{1}$ be the closure of functions in $C^{\infty}_{c}(\mathbb{R}^{3},\mathbb{R}^{3})$ with respect to the $L^{2}$ norm of the gradient and let $\dot{H}^{-1}$ be its dual. In order to incorporate the incompressibility condition we define $\dot{H}_{\sigma}^{1}=\{w\in\dot{H}^{1}:\operatorname{div}w=0\}$ and denote its dual by $\dot{H}_{\sigma}^{-1}$.

Our main result is the following local version of Einstein’s formula.

Theorem 1**.**

The solution $\bar{u}\in\dot{H}_{\sigma}^{1}$ to the equation

[TABLE]

is close to $u$, the solution of (1.2)-(1.7), in the following sense:

[TABLE]

Furthermore, if $U\subset\mathbb{R}^{3}$ is of finite measure and $p\in[1,\frac{3}{2}]$, then

[TABLE]

Remark 2**.**

Note that

[TABLE]

Therefore, even in $L^{\infty}$, the solution of the homogenized equation is close to $u$ on scale $\phi$ in a volume that is asymptotically the whole $\mathbb{R}^{3}$.

*Since Einstein’s formula is only asymptotic for $\phi\to 0$ and the effective viscosity at $\phi=0$ is $\mu$ there is no hope to obtain convergence to the solution of a single non-trivial limit problem. Therefore the comparison to a family of limiting problems is the right approach. Note that already the difference of $u$ to the homogeneous solution with visocsity $\mu$ is of order $\phi$. So only after dividing the difference by $\phi$ the result is meaningful. *

2.3 Prior and related results

While the number of rigorous mathematical papers regarding Einstein’s formula is limited, there is a multitude of physics literature. As a first generalization ellipsoidal particles have been considered by Jeffery [Jef22] and later by Hinch and Leal [LH71, HL72], while drops of another fluid (with finite viscosity) suspended in a surrounding fluid are considered for the first time in [Tay32].

In [KRM67] the authors establish several extremum principles for the Stokes flow including fairly general boundary conditions and rigid particles. They use those principles to obtain, among other results, bounds and an asymptotic formula for the effective viscosity in the low concentration regime and for high concentrations when the particles are situated on a lattice. In the same year in [FA67] another result is given for high concentrations. Numerical research can be found in [NK84], in which arbitrary concentrations are considered and also asymptotic formulas for high concentrations are obtained. [BBP05] considers the case of highly concentrated suspensions and uses a so-called network approximation.

In [KRM67] the effective viscosity is obtained by comparing the dissipation rate of the suspension and of a homogenous fluid with different viscosity subject to pure strain boundary conditions. The pure strain boundary condition is imposed for a domain that becomes infinite in the limit in order to circumvent boundary effects. This disadvantage is overcome only in 2012 by Haines and Mazzucato [HM12] when they rigorously prove, simultaneously bounding the power of the next order term, that

[TABLE]

They consider a fixed domain with pure strain boundary conditions with particle positions fixed to a lattice and also compare dissipation rates.

In [LSP85] the periodic homogenization of the Navier-Stokes equation is discussed. For the first time, the effective viscosity is not determined by an asymptotic or a dissipation functional method, but as a prefactor of the strain in the homogenized equation. In their paper the authors derive a homogenized Navier-Stokes equation up to terms of order $\phi$ that includes the term ([LSP85, p. 13])

[TABLE]

Almog and Brenner [AB98] consider non-constant volume fraction and ensemble averages and obtain an effective viscosity field $\mu(x)$ which confirms Einstein’s formula. Also here the effective viscosity appears inside the Stokes equation. They also recover the results up to $\phi^{2}$ with a second factor $6.95$. Although [LSP85] and [AB98] derive equations with an effective viscosity and are in that respect similar to our approach, both results are not completely rigorous.

A second order correction to the viscosity is first considered by Batchelor and Green. In [BG72] they calculate the second order correction to the viscosity for a random distribution of spheres to be $7.6\phi^{2}$ with an estimated error of the numerical factor of $10\%$ which comes from numerical and asymptotic evaluation of an, in principle, known function. [AGKL12] obtain the term $\frac{5}{2}\phi^{2}$ (the numerical value is computed wrongly in the paper despite a correct formula) and the recent [GH19] confirm this result (see also below).

Shortly after the first submission of our article, several parallel contributions appeared, that together yield a substantially clearer picture of the matter than the one that existed during the writing of this article. In [HW19] Einstein’s formula is justified in a quantitative way and for general shapes with controlled diameter. There one can also find a hint on how to generalize our method to general shapes. The idea is to compare the dissipation energy of the shape to the one of an enclosing sphere and derive estimates from there. [GH19] contains a conditional result on the second order correction and in particular they recover the term $\frac{5}{2}\phi^{2}$ obtained in [AGKL12] (apparently the paper contains a miscalculation and $\frac{5}{2}$ is the corrected result) in the case of periodically arranged spheres. [DG19] considers the more abstract problem of a general effective viscosity tensor for finite volume fraction, i.e. they consider the limit $N\to\infty$ without simultaneously letting $\phi\to 0$. They derive, in the stationary and ergodic random setting a formula similar to the one obtained in [MK06] for the corresponding electrostatic problem, at the same time providing corrector results. Finally, [Ger19] gives a short and clean proof for a result that justifies the analogue of Einstein’s formula in the electrostatic setting (and seems to work similarly for the Stokes problem). Let us stress here that although quantifying convergences and considering more general cases in comparison to this article, all results in [GH19, HW19, Ger19, DG19] provide closeness results in low $L^{p}$ spaces $p\leq 2$. However, in order to apply the results to the dynamic problem and to consider the evolution of the particle density it is crucial to have estimates in $L^{\infty}$. We provide those in the whole space for the relevant approximation (Theorem 6) and for a domain that converges to $\mathbb{R}^{3}$ in case of the homogenized velocity field (Remark 2).

3 Strategy of the proof and preliminaries

3.1 Strategy of proof

The proof is divided into two parts. In the first part we justify the dipole approximation while in the second part we make the heuristic computations of Subsection 1.1 rigorous. To that end we use successive approximations $u\to v_{1}\to\tilde{u}\to\hat{u}\to\bar{u}$.

In Section 4 we prove that the explicit dipole approximation $\tilde{u}$ from (1.12) is actually close to the microscopic solution $u$. For this we first define a related but abstract dipole approximation $v_{1}$ defined via projections to subspaces of $\dot{H}_{\sigma}^{1}$ incorporating the rigid body boundary conditions on the particles (Subsection 4.1). The so-called method of reflections then gives closeness of $v_{1}$ to $u$ (see Theorem 6). The method of reflections for several particles was first introduced by Smoluchowski in [Smo11] and is used extensively in the physics literature. Rigorous proofs for the convergence of the method suited for the treatment of sedimenting particles are given in [Luk89, JO04, HV18, Hoe18]. See also [MNP00, MMP02] where a reflection method is used for a single particle in a domain with boundary. We adapt here the method developed in [Hoe18] for dipole approximations. Since we take into account rotations for the particles in the Stokes case, some adjustments have to be made and it is necessary to establish a Korn and a Korn-Poincaré inequality for balls with integrated boundary conditions. Using carefully obtained characterizations of the projections (Subsection 4.2) we then show that inside the particles $v_{1}$ and $\tilde{u}$ are already close and that, again, using the decay of the dipoles, this can be extended to the domain outside the particles (see Subsection 4.3).

Closeness of $v_{1}$ to $\tilde{u}$ (see Lemma 4.24 and Lemma 4.26) is a consequence of the fact that both the abstract and the explicit dipoles and hence also their difference are dipoles and thus have good decay properties. This is shown in Subsection 4.4.

The second part of the proof consists in proving the closeness of $\tilde{u}$ to the solution $\bar{u}$ of the Stokes equation with Einstein viscosity (2.10). This is done in Section 5. We use an intermediate approximation $\hat{u}$ which is the solution to the equation

[TABLE]

To prove that $\tilde{u}$ is close to $\hat{u}$ (Lemma 5.3) we use the fact that for every point in space the contributions of the particles in a moderately large region around this point are negligible. But further away the number density $\rho^{N}$ looks approximately like $\rho$ which allows passage from sum to integral. The proof relies heavily on the representation of solutions as convolutions with the fundamental solutions and involves various estimates regarding these convolution integrals.

In order to replace $v$ by $\bar{u}$ in equation (3.1) we first prove that $v$ is already close to $\bar{u}$ namely that $\left\lVert v-\bar{u}\right\rVert\leq C\phi$. This is achieved by standard regularity arguments and estimates of the solutions of the homogenized equation in terms of the right hand side. With the same methods it is then possible to prove that the solution to (3.1), $\hat{u}$, is close to the solution of the final equation (2.10), $\bar{u}$, since their difference satisfies an equation with a right hand side that is already small (Lemma 5.9).

The strategy of first providing a suitable approximation to the microscopic solution from which the homogenized problem can be derived was already executed in several works. See for example [Hoe18] in the case of sedimentation or [AZ17] in the case of acoustic waves.

3.2 Notation

In addition to $\dot{H}^{1}$ and $\dot{H}_{\sigma}^{1}$ introduced above in Section 2.2 we introduce some further notation. Here for all spaces we will drop the target space $\mathbb{R}^{3}$ in the notation. We denote the $\dot{H}^{1}$ pairing with $\left\langle{\cdot,\cdot}\right\rangle$ while we write $\left({\cdot,\cdot}\right)$ for the $L^{p}-L^{q}$ pairing where $\frac{1}{p}+\frac{1}{q}=1$. For two functions $u,v\in\dot{H}^{1}$ this means $\left\langle{u,v}\right\rangle=\left({\nabla u,\nabla v}\right)$. For $f\in\dot{H}^{-1}$ we write write $(f,\varphi)$ for $f[\varphi]$ which coincides with the classical notation if $f\in L^{\frac{6}{5}}(\mathbb{R}^{3})$.

We denote the symmetric gradient of $w$ by $ew=\frac{1}{2}\left({\nabla w+\nabla w^{T}}\right)$. Note that for functions in $\dot{H}_{\sigma}^{1}$ the $L^{2}$ pairing of the gradients is the same up to a factor of $2$ as the $L^{2}$ pairing of the symmetric gradients.

Also note that there is a unique vector $\omega u(x)\in\mathbb{R}^{3}$ such that $\nabla u(x)\,y=eu(x)\,y+\omega u(x)\wedge y$ for all $y\in\mathbb{R}^{3}$. This is just a consequence of the fact that for every skewsymmetric matrix $S\in\mathbb{R}^{3\times 3},S=-S^{T}$ there exists $\omega\in\mathbb{R}^{3}$ such that $Sy=\omega\wedge y$ for all $y\in\mathbb{R}^{3}$.

By the Riesz theorem, for any $f\in\dot{H}_{\sigma}^{-1}$ there is a $w\in\dot{H}_{\sigma}^{1}$ such that

[TABLE]

By [Gal94, Lemma V.1.1] we have that if equation (3.2) holds for all $\varphi\in\dot{H}_{\sigma}^{1}$, then, there exists $p\in L^{2}(\mathbb{R}^{3})$ such that

[TABLE]

We then say that

[TABLE]

in the weak sense. The solution operator $S^{-1}:\dot{H}_{\sigma}^{-1}\to\dot{H}_{\sigma}^{1}$ that maps $f$ to $w$ is an isometric isomorphism and its inverse $S$ is the so-called Stokes operator.

The solution operator $S^{-1}$ is given by $S^{-1}f=\Phi\ast f$ where $\Phi$ is the fundamental solution of the Stokes equation, the so-called Oseen tensor

[TABLE]

The corresponding pressure such that $-\Delta S^{-1}f+\nabla p=f$ is given by

[TABLE]

Since the pressure $p$ is merely a Lagrange multiplier ensuring that the velocity field is solenoidal we will write $p$ for every appearing pressure, so that it may change between different equations but also from line to line in one computation.

In the following universal constants $C>0$ will often appear in statements. They never depend on $N,R,d$ and $X_{1},\dots,X_{N}$ and other $N$-dependent quantities but possibly on $f$ unless otherwise stated. When constants appear they might change their value from line to line without indication.

For all spaces we will use a [math] as subscript to indicate that the boundary values of that function vanish. Also, for any classical Sobolev spaces, the subscript $\sigma$ indicates that the weak divergence vanishes.

3.3 Weak formulation of the problem

Let $f\in L^{\frac{6}{5}}(\mathbb{R}^{3})\cap L^{2}(\mathbb{R}^{3})$. A function $u$ is a weak solution of problem (1.2)-(1.7) if $u\in\dot{H}_{\sigma}^{1}$ (which implies (1.3),(1.7)), if $u$ is a rigid body motion on all $\overline{B_{i}}$ for $i=1,..,N$ (this is (1.6)), if for all $\varphi\in\dot{H}^{1}_{\sigma,0}(\Omega)$

[TABLE]

and if (1.4), (1.5) are satisfied. Here $\nabla u$ and $p$ are a priori only in $L^{2}(\Omega)$ and (1.4), (1.5) may seem a bit ambiguous at first glance since $\sigma$ is only in $L^{2}(\Omega)$ a priori. However, we have that $\operatorname{div}\sigma=f^{N}\in L^{2}$ and hence $\sigma n\in H^{-\frac{1}{2}}(\partial B_{i})$. That $\operatorname{div}\sigma\in L^{2}(\Omega)$ is a simple consequence of (3.5) and the so-called reciprocal principle (see, e.g. [HB65]): For any $p\in L^{2}(\mathbb{R}^{3})$ and $w\in\dot{H}_{\sigma}^{1}$ we write $\sigma=2ew-p\operatorname{Id}$. Then for $v,w\in\dot{H}_{\sigma}^{1}$ we have

[TABLE]

In the second step we used that the scalar product of a symmetric and a skew-symmetric matrix is zero ($\nabla v=ev+\left({\nabla v}\right)^{\operatorname{skew}}$), while in the second, we used, that $v$ is divergence-free whence $\operatorname{Id}\cdot\nabla v=\operatorname{div}v=0$. The name reciprocal principle comes from the fact that the same equality holds for interchanged $w,v$. Note that, if $w$ satisfies (1.6), then, because $ew=0$ in $B_{i}$ for all $i=1,\dots,N$ we can write

[TABLE]

Existence of a weak solution to (1.2)-(1.7) follows by minimizing the energy

[TABLE]

in the space of functions that are rigid body motions inside the particles:

[TABLE]

4 The dipole approximation

It is useful not to start with the explicit dipole approximation from (1.12), but to consider dipoles which can be characterized using a variational formulation which makes the comparison to the microscopic solution simpler. In order to do so we adapt the theory developed in [Hoe18]. With the exception of the Korn and Korn-Poincaré inequality in Lemma 4.7 and Corollary 2, which are not needed in [Hoe18], all statements in Subsections 4.1, 4.2 and 4.3 concerning the abstract dipoles are statements/ideas from [Hoe18] adapted to the situation of rigid body motions. Since the proofs are very similar we just state deviations and refer the reader to [Hoe18] and [Sch19] for the detailed proofs.

4.1 Approximation by abstract dipoles

The solution $v$ of the particle-free problem (1.8)-(1.10) is given by $v=\Phi\ast f^{N}$ and is the minimizer of the energy $E$ from (3.7) in the space $\dot{H}_{\sigma}^{1}$.

On the other hand, the solution $u$ to problem (1.2)-(1.7) is the minimizer of $E$ in the space $W$ (see equation (3.8)). The fact that $u$ minimizes $E$ in $W$ means that $u$ is the orthogonal projection of $v$ from $\dot{H}_{\sigma}^{1}$ to the subspace $W$. We call $P:\dot{H}_{\sigma}^{1}\to W$ the orthogonal projection so that $u=Pv$. This implies that

[TABLE]

By choosing a suitable function $w$ one can thus get an estimate for $\left\lVert u-v\right\rVert_{\dot{H}^{1}}$. This is part of the proof of Theorem 6.

In order to control the velocity field locally we need to consider the $L^{\infty}$ norm, though. To get $L^{\infty}$ estimates it is useful to work with the method of reflections. This is due to the fact that the projection onto $W$ is not so easy to characterize. The method of reflections works with solutions of single particle problems. The single particle spaces involved, are much easier to characterize than $W$. For this we first define the particle wise version of $W$:

[TABLE]

Since $W_{i}$ is a closed subspace of $\dot{H}_{\sigma}^{1}$ there is an orthogonal projection $P_{i}:\dot{H}_{\sigma}^{1}\to W_{i}$. Notice that $W=\cap_{i=1}^{N}W_{i}$. The orthogonal complement of $W_{i}$ has a useful characterization:

Lemma 4.1**.**

[TABLE]

Proof 4.2**.**

Take any element $s$ of the right hand side. Then for any $w\in W_{i}$, using (3.6), we obtain (let $w=V+\omega\wedge(x-X_{i})$ on $\overline{B_{i}}$):

[TABLE]

*On the other hand, since the first line is zero for any $s\in\left({W_{i}}\right)^{\perp}$ we have $0=-\operatorname{div}\sigma[s]=-\Delta s+\nabla p$ in $\mathbb{R}^{3}\setminus\overline{B_{i}}$ by considering $w\in\dot{H}^{1}_{\sigma,0}(\mathbb{R}^{3}\setminus\overline{B_{i}})$. Then we can conclude by using $V,\omega=e_{1},e_{2},e_{3}$ in the last line. *

A function with the property that

[TABLE]

is usually called a dipole since the first moment of the force distribution vanishes inside the ball.

We come back to our goal to approximate $u$ by $v$. We already know that, following the idea from Subsection 1.1, it makes sense to subtract from $v$ at every ball the dipole preventing $v$ from being a rigid body motion . Let $Q_{i}=\operatorname{Id}-P_{i}$ be the orthogonal projection onto $W_{i}^{\perp}$. We know that $v-Q_{i}v=P_{i}v\in W_{i}$ is a rigid body motion on the ball $B_{i}$, hence $Q_{i}v$ is the dipole we are looking for. As explained in Subsection 1.1, subtracting $Q_{i}v$ only helps with the boundary condition on $B_{i}$, so we have to subtract the dipole for all balls which gives rise to the first approximation v_{1}=\big{(}\operatorname{Id}-\sum_{i=1}^{N}Q_{i}\big{)}v Of course this approximation will not be constant on all balls, since the additional $Q_{j}v$ for $j\neq i$ will have non-vanishing contributions on $B_{i}$. But since the dipoles solve the homogeneous Stokes equation outside $B_{i}$, the function $v_{1}$ still satisfies the same equation as $v$ and $u$ in $\Omega$. In order to get even closer to $u$ let us repeat the process of subtracting dipoles in order to make the functions closer to rigid body motions on the balls. This leads to approximations $v_{k}$ given by

[TABLE]

The idea is that taking $k\to\infty$ one should have P=\lim_{k\to\infty}\Big{(}\operatorname{Id}-\sum_{i=1}^{N}Q_{i}\Big{)}^{k}. This would imply that $v_{k}\to Pv=u$ as $k\to\infty$.

Before attempting to prove that the $v_{k}$ converge to $u$ we need a better understanding of the projections $P_{i}$ and $Q_{i}$ respectively.

4.2 Characterization of $W_{i}^{\perp}$ and $P_{i}$

Lemma 4.3**.**

For $w\in W_{i}^{\perp}$ we have

[TABLE]

Hence, for $w\in\dot{H}_{\sigma}^{1}$ the projection to $W_{i}$ satisfies $P_{i}w(x)=V+\omega\wedge(x-X_{i})$ for all $x\in\overline{B_{i}}$ with

[TABLE]

Remark 1**.**

Here and in the following we set

[TABLE]

Proof 4.4**.**

The idea is to use explicit test functions, for spatial motion and rotation separately, for which we know the drag on the sphere explicitly. See [Hoe18], Lemma 3.10, for the first parts of (4.2) and (4.3).

For the second part of (4.2) take $\varphi$ such that $-\Delta\varphi+\nabla p=0$, $\operatorname{div}\varphi=0$ in $\mathbb{R}^{3}\setminus\overline{B_{i}}$ and $\varphi=\omega\wedge(x-X_{i})$, $\omega\in\mathbb{R}^{3}$ on $\overline{B_{i}}$. Then $\varphi\in W_{i}$ and

[TABLE]

for $\left\lvert x-X_{i}\right\rvert>R$. The drag on the sphere is then given by $(\sigma[\varphi]n)=-\frac{3}{R}\omega\wedge(x-X_{i})$. Then

[TABLE]

Setting $\omega=e_{1},e_{2},e_{3}$ we arrive at the second part of (4.2).

Now for $w\in\dot{H}_{\sigma}^{1}$ we know that $w-P_{i}w=Q_{i}w\in W_{i}^{\perp}$. On the other hand there are $V,\omega\in\mathbb{R}^{3}$ such that $P_{i}w(x)=V+\omega\wedge(x-X_{i})$ for all $x\in\overline{B_{i}}$.

Also, using the vector rule $A\wedge(B\wedge C)=(A\cdot C)B-(A\cdot B)C$ we have:

[TABLE]

and therefore

[TABLE]

As a consequence of this characterization we obtain a Poincaré inequality and a (first) Korn inequality and consequently a Korn-Poincaré inequality on $W_{i}^{\perp}$.

Lemma 4.5**.**

Let $r>0$ and $X\in\mathbb{R}^{3}$. Let $p\in[1,\infty]$ and let

[TABLE]

For $p>1$ there is a constant $C>0$ that does not depend on $X$ or $r$ such that for all $w\in H^{1,p}_{X,r}$:

[TABLE]

Proof 4.6** (Proof of Lemma 4.5).**

*By scaling and translation it is enough to prove the inequality for $X=0$ and $r=1$. It is well known, that for closed cones in which $\nabla w=0$ implies that $w=0$, such a Poincaré inequality holds. *

Lemma 4.7**.**

Let $r>0$ and $X\in\mathbb{R}^{3}$. Let $p\in[1,\infty]$ and let

[TABLE]

For $p\in(1,\infty)$ there is a constant $C>0$ that does not depend on $X$ or $r$, such that for all $w\in H^{1,p}_{\sigma,X,r}$:

[TABLE]

Proof 4.8** (Proof of Lemma 4.7).**

Again, we only need to prove the result for $X=0,r=1$, the case of general $X$ and $r$ follows by translation and rescaling. In the following we will omit the domain $B_{1}(0)$ in the spaces.

$H^{1,p}_{\sigma,0,1}$* is a closed cone. The proof of this Korn inequality is indirect and uses the same idea as the proof of general Poincaré inequalities on closed cones. For $w\in H^{1,p}_{\sigma,0,1}$ we have that $ew=0$ implies $w=0$. This follows from a combination of the well-known fact that $ew=0$ implies that $w$ is skew-symmetric affine and the boundary conditions in $H^{1,p}_{\sigma,0,1}$.*

Now, for the sake of contradiction, assume that there is no $C>0$ such that (4.4) holds for all $w\in H^{1,p}_{\sigma,0,1}$. Then there is a sequence $w_{k}\in H^{1,p}_{\sigma,0,1}$ such that $\left\lVert\nabla w_{k}\right\rVert_{L^{p}}\geq k\left\lVert ew_{k}\right\rVert_{L^{p}}$. By rescaling we can arrange that $\left\lVert w_{k}\right\rVert_{W^{1,p}}=1$ for all $k\in\mathbb{N}$. But then there is a subsequence (again denoted by $w_{k}$) and $w_{\ast}\in W^{1,p}$ such that $w_{k}\rightharpoonup w_{\ast}$ in $W^{1,p}$. Note that since $H^{1,p}_{\sigma,0,1}$ is convex we have that $w_{\ast}\in H^{1,p}_{\sigma,0,1}$. On the other hand we know that

[TABLE]

*Hence $ew_{k}\to 0$ in $L^{p}$. Since at the same time $ew_{k}\rightharpoonup ew_{\ast}$ we get $ew_{\ast}=0$ . By our foregoing considerations this implies that $w_{\ast}=0$.

*By compact embedding we know that $w_{k}\to 0$ strongly in $L^{p}$. But to reach a contradiction we also need that the full gradient $\nabla w_{k}\to 0$ strongly, not only the symmetrized part. The key idea is, that by Korn’s second inequality the gradients are already close to some constant skew-symmetric matrices for which we have strong compactness.

Indeed, ([KO88], §2, Theorem 8) there exist matrices $A_{k}\in\mathbb{R}^{3\times 3}_{\operatorname{skew}}$ such that

[TABLE]

Since $\nabla w_{k}$ is bounded in $L^{p}$ this implies that the sequence $(A_{k})_{k}$ must be bounded in $\mathbb{R}^{3\times 3}_{\operatorname{skew}}$. But then there is a subsequence (again denoted $A_{k}$) such that $A_{k}\to A_{\ast}$. Furthermore we can pick the subsequence in such a way that

[TABLE]

Then we have

[TABLE]

*Therefore, $\nabla w_{k}$ converges strongly to $A_{\ast}$ in $L^{p}$. But at the same time $\nabla w_{k}\rightharpoonup 0$ weakly in $L^{p}$. This yields $A_{\ast}=0$ and $\nabla w_{k}\to 0$ strongly in $L^{p}$ and $w_{k}\to 0$ strongly in $W^{1,p}$. This is a contradiction. *

Combining Lemma 4.5 and Lemma 4.7 as well as the Sobolev-embedding $W^{1,q}\hookrightarrow L^{\infty}$ for $q>3$ we obtain

Corollary 2**.**

Let $p\in(1,\infty]$. There is a constant $C>0$ that does not depend on $X$ or $r$ such that for all $w\in H^{1,p}_{\sigma,X,r}$:

[TABLE]

We know that elements of $w\in W_{i}^{\perp}$ solve the homogeneous Stokes equation outside $\overline{B_{i}}$. It is well-known that this coincides with minimizing the respective norm.

Lemma 4.9**.**

*Let $s\in\dot{H}_{\sigma}^{1}$ and $-\Delta s+\nabla p=0$ on $\mathbb{R}^{3}\setminus V$ for some closed $V$ with Lipschitz boundary. Then $s$ minimizes $\left\lVert ew\right\rVert_{L^{2}(\mathbb{R}^{3}\setminus V)}$ among all $w\in\dot{H}_{\sigma}^{1}$ with $w=s$ on $V$. *

By choosing a suitable competitor we can therefore get estimates of the norm of $w\in W_{i}^{\perp}$ in terms of its values in $\overline{B_{i}}$. This competitor can be constructed using extension operators.

Lemma 4.10**.**

There is a constant $C>0$ that does not depend on $X$ or $r$, and an extension operator $E_{X,r}:H^{1,2}_{\sigma,X,r}\to H_{\sigma,0}^{1}(B_{2r}(X))$ such that

[TABLE]

Proof 4.11**.**

We can get $E_{X,r}$ from $E_{0,1}$ by translation and scaling without changing $C$. Let $E_{0,1}:H^{1}_{\sigma}(B_{1}(0))\to H_{\sigma,0}^{1}(B_{2}(0))$ be a continuous extension operator. This can be constructed from the usual extension operators by using the Bogowskii-operator (see [Bog80]. Then we have

[TABLE]

*where the second to last inequality came from Lemma 4.5 while the last is the Korn inequality from Lemma 4.7. *

Using $E_{X_{i},R}w$ as a competitor in Lemma 4.9 and employing Lemma 4.10 we obtain

Corollary 3**.**

There is a constant $C>0$ such that for all $w\in W_{i}^{\perp}$

[TABLE]

We can furthermore prove the following decay properties of the dipoles

Lemma 4.12**.**

There is a constant $C>0$ such that for all $w\in W_{i}^{\perp}$ and for all $x\in\mathbb{R}^{3}\setminus B_{2R}(X_{i})$ we have

[TABLE]

Proof 4.13**.**

*The main idea of the proof is to use that the convolution of the fundamental solution is the inverse of the application of the Stokes operator. Therefore $w$ can be represented as the convolution of $Sw$ and the fundamental solution. For any function $w$ with compactly supported $Sw$ this gives a $\frac{1}{\left\lvert x\right\rvert}$-like decay. The fact that $Sw$ integrates to zero due to the dipole property gives rise to the additional power in the decay. For details see Lemma 3.11 in [Hoe18]. *

4.3 Convergence of the method of reflections

The ultimate goal of the method of reflections is, to prove that the symmetrized gradients of the approximations $v_{k}$ approach [math] in $L^{\infty}$ inside the particles. If $w\in\dot{H}^{1}$ has an uniformly bounded gradient in $B_{i}$ we can get the following estimates

Corollary 4**.**

There is a constant $C>0$ such that for all $w\in\dot{H}_{\sigma}^{1}\cap W^{1,\infty}(B_{i})$ and for all $x\in\mathbb{R}^{3}\setminus B_{2R}(X_{i})$ we have

[TABLE]

Proof 4.14**.**

*We know that $Q_{i}w=w-P_{i}w$ and since $P_{i}w$ is a rigid body motion on $B_{i}$ we have $eP_{i}w=0$ in $B_{i}$ and hence $eQ_{i}w=ew$ in $B_{i}$. Using that $Q_{i}w\in W_{i}^{\perp}$, and combining Corollary 3 with Lemma 4.12 yields (4.7) and (4.8). *

In order to prove the main approximation statement of this section we need some estimates for recurring sums (also see Lemma 2.1 of [JO04] for the first two inequalities):

Lemma 4.15**.**

There exists $C>0$ such that for any $x\in\mathbb{R}^{3}$ with $1\leq i\leq N$ such that $\left\lvert x-X_{i}\right\rvert\leq\left\lvert x-X_{j}\right\rvert$ for all $1\leq j\leq N$ it holds that

[TABLE]

Proof 4.16**.**

Without loss of generality we assume $i=1$ and $X_{1}=0$. We order the balls in such a way that $\left\lvert X_{1}\right\rvert\leq\dots\leq\left\lvert X_{N}\right\rvert$. Since $d_{ij}\geq d$ the balls $B\left({X_{i},\frac{d}{2}}\right)$ and $B\left({X_{j},\frac{d}{2}}\right)$ do not intersect. Moreover for any $2\leq i\leq N$ we have

[TABLE]

We compare the left and the right volume to obtain $i\left({\frac{d}{2}}\right)^{3}\leq(2\left\lvert X_{i}\right\rvert)^{3}$ and $\left\lvert X_{i}\right\rvert\geq\frac{1}{4}d\;i^{\frac{1}{3}}$. Now

[TABLE]

*We deduce now (4.9) from this and the fact that $\left\lvert x-X_{j}\right\rvert\geq\frac{1}{2}\left\lvert X_{i}-X_{j}\right\rvert$. The other estimates follow similarly. *

We will also use a maximum modulus theorem for the Stokes equation:

Lemma 4.17** ([MRS99]).**

Let $\Omega\subset\mathbb{R}^{3}$ be an exterior domain (a domain with bounded complement) and assume that $g\in C^{0}(\Omega^{C})$ satisfies $\int_{\partial\Omega}g\cdot n\;\operatorname{d}\!S=0.$

Then, there is a constant $C>0$ that depends only on $\Omega$, such that the unique solution $u\in\dot{H}_{\sigma}^{1}(\mathbb{R}^{3})$ of the Dirichlet problem

[TABLE]

*satisfies $\left\lVert u\right\rVert_{L^{\infty}}\leq C\left\lVert g\right\rVert_{L^{\infty}}$. *

Remark 5**.**

We only have to use this statement for $\Omega$ being the exterior of a ball. The constant is invariant under translation and scaling of the ball.

We are now able to prove the convergence of the method of reflections.

Theorem 6**.**

There is $\varepsilon>0$ such that if $\phi\log N<\varepsilon$, we have $v_{k}\to u$ in $\dot{H}^{1}$ and $L^{\infty}(\mathbb{R}^{3})$ for $k\to\infty$ and in particular

[TABLE]

Proof 4.18**.**

The proof is in principle similar to the proof of Proposition 3.12 in [Hoe18] and we refer to that article for details. We explain briefly the strategy.

We first have to establish that for all $v_{k}$ we have $Pv_{k}=u$ where $P$ is the projection onto $W$. Using this we know that $\left\lVert u-v_{k}\right\rVert_{\dot{H}^{1}}\leq\left\lVert w-v_{k}\right\rVert_{\dot{H}^{1}}$ for all $w\in W$. With the use of the extension operator $E$ we then construct suitable competitors in $W$ to obtain

[TABLE]

Using that $ev_{k}-eQ_{i}v_{k}=0$ on $B_{i}$ and the decay of $\nabla Q_{j}v_{k}$ for $j\neq i$ (Corollary 4) together with Lemma 4.15 gives an estimate of the form $\left\lVert ev_{k+1}\right\rVert_{L^{\infty}(\cup B_{i})}\leq C\phi\log N\left\lVert ev_{k}\right\rVert_{L^{\infty}(\cup B_{i})}$. By the smallness of $\phi\log N$ we obtain $v_{k}\to u$ in $\dot{H}^{1}$. For the convergence in $L^{\infty}$ we can use the same decay techniques together with Lemma 4.17.

We just show that $ev\in L^{\infty}(\cup B_{j})$ which is needed for the proof. For this we estimate

[TABLE]

*Therefore actually $\nabla v\in L^{\infty}(\mathbb{R}^{3})$. *

4.4 The explicit dipole approximation

The dipole approximation $v_{1}$ is close enough to $u$. The next step is to show that $v_{1}$ is close to the explicit dipole approximation $\tilde{u}$ from (1.12). For this we need a technical lemma that gives us additional decay for differences of terms from the fundamental solution.

Lemma 4.19**.**

There is a constant $C>0$ such that for $x,z\in\mathbb{R}^{3}$ with $\left\lvert x-X_{i}\right\rvert\leq\frac{1}{2}\left\lvert X_{i}-z\right\rvert$ and all $n>l\geq 0$ and $k\neq i$ we have the following estimate:

[TABLE]

*Here, for $y\in\mathbb{R}^{3}$, the tensor $y^{\otimes l}\in\mathbb{R}^{3^{l}}$ is given by $y^{\otimes l}_{i_{1}\dots i_{l}}=y_{i_{1}}\dots y_{i_{l}}$. *

Proof 4.20**.**

We have $\left\lvert x-z\right\rvert\geq\frac{1}{2}\left\lvert X_{i}-z\right\rvert$. First we set $g(y)=\frac{(y-z)^{\otimes l}}{\left\lvert y-z\right\rvert^{n}}$. Observe that for $y\in B_{\left\lvert x-X_{i}\right\rvert}(X_{i})$:

[TABLE]

since $\left\lvert y-z\right\rvert\geq\frac{1}{2}\left\lvert X_{i}-z\right\rvert$. Then

[TABLE]

The symmetrized gradient of the dipoles $d_{i}$ from (1.11) and $Q^{P}_{i}v$ inside the balls is $ev(X_{i})$ and $ev(x)$ respectively. This is why we need to control the oscillation of $\nabla v$ inside the balls.

Lemma 4.21**.**

There is constant $C>0$ such that for all $x\in\mathbb{R}^{3}$ and $a=\left\lvert x-X_{i}\right\rvert$

[TABLE]

*Recall that $\alpha$ is the Hölder exponent of $f$. *

Remark 7**.**

In particular, taking $a<R$ we have

[TABLE]

*Recall that $\alpha$ is the Hölder exponent of $f$. The exponent $\frac{1}{4}$ is chosen for ease of notation. Any exponent smaller than $\frac{1}{3}$ works. *

Proof 4.22** (Proof of Lemma 4.21).**

If $v$ was in $C^{2}$ the statement would be an easy consequence of the Taylor expansion. But we cannot expect $v$ to be two times differentiable since $f^{N}$ is zero on the particles and hence in general not continuous. Still, the function that is subtracted from $f$ to obtain $f^{N}$ is supported only on the particles and hence should have a contribution vanishing with $\phi$. This can be made clear by writing $\nabla v=\nabla\Phi\ast f^{N}$. For the part that remains we expect Hölder continuity, since $f$ is Hölder continuous. We know:

[TABLE]

We split the term into a part with the pure $f$ and a part that incorporates the characteristic functions of the particles. Let us first estimate the latter for $\left\lvert x\right\rvert<2L$:

[TABLE]

Here we used that $\left\lvert\nabla\Phi(x)\right\rvert\leq C\frac{1}{\left\lvert x\right\rvert^{2}}\in L^{\frac{4}{3}}(B_{3L}(0))$. If $\left\lvert x\right\rvert\geq 2L$ then $\left\lvert x-y\right\rvert\geq L$ and we can estimate $\nabla\Phi$ in $L^{\infty}$ obtaining the $\phi$ term. For the other term and $\left\lvert x\right\rvert<L$ we compute:

[TABLE]

*Here we used Lemma 4.19 in the first step and Hölder continuity of $f$ as well as decay properties of $\nabla\Phi$. If $\left\lvert x\right\rvert\geq L$ we can just use the $L^{\infty}$ bound on $\nabla v$ to estimate the difference by a constant. Combining all estimates we arrive at (4.13). *

Corollary 8**.**

For $x\in\mathbb{R}^{3}\setminus B_{2R}(X_{i})$and $d_{i}$ from (1.11) we have

[TABLE]

Proof 4.23**.**

By computation we see that $d_{i}\in W_{i}^{\perp}$ and hence $d_{i}(x)-Q_{i}v\in W_{i}^{\perp}$. Applying Corollary 4 we obtain

[TABLE]

And for $y\in B_{i}$, with Lemma 4.21, we have:

[TABLE]

It remains to use Corollary 8 to prove that the difference between the whole dipole approximations $v_{1}$ and $\tilde{u}$ is small. Because of the singularities of the dipoles we are now forced to consider the domain $\Omega^{N}_{\delta}$ from (2.9).

Lemma 4.24**.**

[TABLE]

Proof 4.25**.**

We have $v_{1}-\tilde{u}=\sum_{i=1}^{N}\left({Q_{i}v-d_{i}}\right)$. Take $x\in\Omega^{N}_{\delta}$ and let $X_{i}$ be the closest centre point. Then by Corollary 8

[TABLE]

For the dipoles, that are further away, we use Corollary 8 and Lemma 4.15 to get

[TABLE]

We now want to get a similar estimate in $\mathbb{R}^{3}\setminus\Omega_{\delta}^{N}$, the region around the particles. Remark 7 gives us a hint what kind of closeness we can hope for. Let us imagine for a moment that both the first and the second derivative of $v$ are Hölder continuous (i.e. $v$ solves the problem with $f$ and not with $f^{N}$ as a right hand side) and hence bounded. This improves the (optimal) estimate from 4.21 only slightly to contain $R$ instead of $R^{\alpha}$ as one of the smallest terms. Even when multiplied by the $R$ from the Poincaré inequality this is, in general, not of the type $\phi o(1)$. So instead of an estimate in $L^{\infty}$ we aim for an estimate in $L^{p}$ for some $p\geq 1$. In order to estimate the difference $v_{1}-\tilde{u}$ in this space we only need to consider the difference of the dipoles originated at the closest particle because we have $L^{\infty}$ control over the remaining terms. So we find approximately:

[TABLE]

where we used Lemma 4.5 and Remark 7 in the third line. For this to be of type $\phi o(1)$ we need $p\leq\frac{3}{2}$. Notice that on the other hand the dipoles decay like $\frac{1}{\left\lvert x\right\rvert^{2}}$. This is only in $L^{p}(\mathbb{R}^{3}\setminus B_{L}(0))$ for $p>\frac{3}{2}$. Using the explicit form of the $o(1)$ term, it is possible to find a common $p=\frac{3}{2}+\varepsilon$ for which we have closeness on the whole space.

Lemma 4.26**.**

Let $U\subset\mathbb{R}^{3}$ be of finite measure. For $p\in[1,\frac{3}{2}]$ it holds:

[TABLE]

Proof 4.27**.**

First of all note that by Lemma 4.24 we only need to prove the statement for $U=\cup_{i=1}^{N}B_{r}(X_{i})$. Let $x\in U$, then $x\in B_{r}(X_{i})$ for one and only one $i$ because of Assumption (2.2). By the proof of Theorem 4.24 we only need to consider the dipole $d_{i}-Q_{i}v$ because the other dipoles behave exactly the same as outside $B_{r}(X_{i})$, giving $L^{\infty}$ estimates. If $x\in B_{2R}(X_{i})$, by the maximum modulus theorem from [MRS99] plus corollary 2 as well as Lemma 4.21:

[TABLE]

By the computation done before the lemma this bounds $\left\lVert v_{1}-\tilde{u}\right\rVert_{L^{p}(\cup_{i=1}^{N}B_{2R}(X_{i}))}$ by $\phi o(1)$. If $r=2R$ we are done. Otherwise observe that for $\left\lvert x-X_{i}\right\rvert\in(2R,\delta)$ we can use Corollary 8 and get

[TABLE]

This was the computation for $p\neq\frac{3}{2}$. If $p=\frac{3}{2}$ we get

[TABLE]

*Here we used $-N^{-\frac{1}{3}}\phi^{\frac{1}{4}}\log R\leq CR^{-\frac{1}{2}}N^{-\frac{1}{3}}\phi^{\frac{1}{4}}\leq C\phi^{-\frac{1}{6}}N^{-\frac{1}{6}}\phi^{\frac{1}{4}}=o(1)$. *

5 Homogenization

5.1 From the microscopic approximation to a homogenized equation

The solution $v$ to the reference problem without particles involves the individual particles on the right hand side and is therefore not in a good form for the treatment of the limit problem. Therefore we first prove that $v$ is close to the solution $\hat{v}$ of the following problem:

[TABLE]

Lemma 5.1**.**

For the solution $v$ of problem (1.8)-(1.10) and the solution $\hat{v}$ of problem (5.1)-(5.3) it holds that

[TABLE]

Let $U\subset\mathbb{R}^{3}$ be of finite measure and $p\in[1,3]$. Then

[TABLE]

Proof 5.2**.**

The idea of the proof is to represent $v,\hat{v}$ in terms of the fundamental solution. Since there is a $\phi$ in front of $\rho$, and $\rho^{N}$ (see 2.6) is close to $\rho$ (see assumption (2.4)) in a weak sense (when convoluted with the fundamental solution), we can interchange $\rho$ by $\rho^{N}$. For given $x\in\mathbb{R}^{3}$ we can ignore terms at regions that are close to $x$ since they are small anyway. For regions further away from $x$, the number density $\rho^{N}$ looks approximately like the rescaled sum of the characteristic functions of the particles.

We write $v,\hat{v}$ by means of the fundamental solution to see that

[TABLE]

The proof comes in four parts. The first part shows that we can replace $\rho$ by $\rho^{N}$ and that particles close to $x$ can be ignored, the second part establishes the closeness of the functions in $L^{\infty}$ while the third part is concerned with the closeness of the gradients in $L^{\infty}$. In the last part the necessary $L^{p}$ estimates are shown.

Part 1: We can replace $\rho$ by $\rho^{N}$ since

[TABLE]

because $\rho^{N}\rightharpoonup\rho$ in $L^{p}(B_{L+1}(0))$ and $f\Phi(x-\cdot)\in L^{q}(B_{L+1}(0))$ where $q$ is the Hölder dual of $p$ and hence $q<\frac{3}{2}$.

Let $X_{i}$ be the closest centre point to $x$. Then we can ignore the $i$th term in the sum:

[TABLE]

and we can replace $\mathds{1}_{B_{k}}$ by $\frac{4\pi}{3}R^{3}\delta_{X_{k}}$ using Lemmas 4.19 and 4.15 as well as the Hölder continuity of $f$:

[TABLE]

Therefore, what is left to show is

[TABLE]

We can ignore the contributions by particles in the range $s$ :

[TABLE]

Here we used that $\rho^{N}$ is uniformly bounded and in the range $s$ there can only be a number of particles $\leq C\frac{s^{3}}{d^{3}}$ and then applied Lemma 4.15 with $N=\frac{s^{3}}{d^{3}}$.

The above reasoning applies to particles in the range of $3s$ in the same way. This means we can ignore all cubes $A_{j}$ that intersect the boundary $\partial B_{s}(x)$ since they will be included in $B_{3s}(x)$ anyway.

Part 2: Therefore estimating the difference above reduces to estimating the difference of appropriately grouped terms in the sum to its corresponding parts (the cube $A_{j}$) of the integral. This means we want to estimate

[TABLE]

Using Hölder-continuity of $f$ and Lemma 4.19 we have

[TABLE]

Hence, using Lemma 4.15:

[TABLE]

Part 3: In order to understand that the estimate holds for the gradient note that

[TABLE]

*We can now reproduce steps 1 and 2 from the proof above using the following facts i)-vi) in that order:

1. i)

$\rho_{N}\rightharpoonup\rho$* in $L^{p}(B_{L+1}(0))$ and $f\nabla\Phi(x-\cdot)\in L^{q}(B_{L+1}(0))$, since $\frac{1}{\left\lvert x\right\rvert^{2}}$ is $q$-integrable for *$q<\frac{3}{2}$ 2. ii)

$\frac{1}{\delta^{2}N}\to 0$, and therefore

[TABLE] 3. iii)

The appearing sums over third powers are well behaved:

[TABLE] 4. iv)

Terms is range $s$ can be disregarded:

[TABLE] 5. v)

The sum over squares is a good as the sum over first powers:

[TABLE] 6. vi)

The sum over the third powers can be controlled since $s$ approaches zero fast enough:

[TABLE]

This gives

[TABLE]

Part 4: In order to get the $L^{p}$ result simply notice that we can use the $L^{\infty}$ results everywhere even where $x\in B_{r}(X_{i})$ as long as we did not use that $\left\lvert x-X_{i}\right\rvert>r$. In fact this was used only once so that we have to look at the following term again when $x\in B_{r}(X_{i})$:

[TABLE]

If $\left\lvert x-X_{i}\right\rvert>2R$ then this is smaller than $CR^{3}\frac{1}{\left\lvert x-X_{i}\right\rvert}$. If $\left\lvert x-X_{i}\right\rvert\leq 2R$ it scales like $R^{2}$. Integrating the $p$th power of the left hand side over the union of the $B_{r}(X_{i})$ gives

[TABLE]

Now we can establish the first closeness result for the solution of the homogenized equation.

Lemma 5.3**.**

Let $\hat{u}$ be the solution to

[TABLE]

and let $\tilde{u}$ be the explicit Stokes dipole approximation. Then we have

[TABLE]

Let $U\subset\mathbb{R}^{3}$ be of finite measure and $p\in[1,\frac{3}{2}]$. Then

[TABLE]

Proof 5.4**.**

In principle we employ the same strategy as in the proof of Lemma 5.1. We represent $\hat{u}$ in terms of the fundamental solution and then we use Lemma 5.1 as well as $\rho^{N}\rightharpoonup\rho$ to show, that the difference of the sum and the integral is small.

We write $\hat{u}$ (componentwise) in terms of the fundamental solution:

[TABLE]

In order to show the closeness of $\tilde{u}$ to the representation from (5.6) let $x\in\Omega^{N}_{\delta}$ be given. We have

[TABLE]

Where ${d_{k}}_{j}$ is the $j$-th component of the explicit dipole at particle $k$ (see (1.11)) and not the distance. Taking into account Lemma 5.1 it remains to prove that

[TABLE]

The proof is divided in several parts. The first part shows that we can replace $\hat{v}$ by $v$. In the second part we show that the closest particle as well as the fast decaying parts of the dipoles can be ignored. The next part determines the explicit form of the gradient of the fundamental solution when applied to a symmetric, trace-free matrix. In the fourth part it is shown that we can replace $\rho$ by $\rho^{N}$ and ignore close particles. The fifth part establishes the closeness of the functions in $L^{\infty}$ while the last part is concerned with the $L^{p}$ result.

Part 1: Using Lemma 5.1 we have

[TABLE]

Therefore it suffices to prove

[TABLE]

Part 2: Let $X_{i}$ be the closest centre point to $x$. Then we can ignore the $i$th term in the sum:

[TABLE]

Next we look at the fast decaying terms of $d_{k}$:

[TABLE]

Thus we can ignore these terms and only consider the slowly decaying terms of the form

[TABLE]

Part 3: We now derive an expression for $ev(y)\nabla\Phi(x-y)$. We compute the derivative of the Oseen-Tensor to be

[TABLE]

Take any symmetric, trace free matrix $\epsilon$. Then

[TABLE]

Let us replace $R^{3}=\phi\frac{1}{N}$. Then we are left to show that

[TABLE]

Part 4: With the same type of argument as in Part 1 of the proof of 5.1 we can replace $\rho$ by $\rho^{N}$, we can leave out a ball of size $s$ around $x$ in the integral and we can ignore the parts of the sum where $\left\lvert X_{k}-x\right\rvert\leq s$. It remains to show:

[TABLE]

We can employ the same reasoning as above to exclude all particles in the range of $3s$. This means we can ignore all cubes $A_{j}$ that intersect the boundary $\partial B_{s}(x)$ since they will eventually be included in $B_{3s}(x)$ anyway.

Part 5: Therefore estimating the difference above reduces to estimating the difference of appropriately grouped terms in the sum to its corresponding parts (the cube $A_{j}$) of the integral. I.e. we need to estimate

[TABLE]

Notice that the cubes with $A_{j}\cap B_{2L}(0)=\emptyset$ have no contribution since there $\rho^{N}=0$.

Looking at one term of the sum we are left to estimate

[TABLE]

We now use the definition of $\rho^{N}$ to write this as

[TABLE]

We can replace $ev(y)$ by $ev(X_{k})$ in the integral since for the difference, by Lemma 4.21, we have:

[TABLE]

and hence

[TABLE]

Since the number of particles in one $A_{j}$ is bounded by $Ns^{3}$, adding this up we obtain

[TABLE]

Using Lemma 4.19 for $y,X_{k}\in A_{j}$, we obtain

[TABLE]

Summing up over $j$ gives

[TABLE]

Part 6: As in Part 4 of the proof of Lemma 5.1 we can use the $L^{\infty}$ results everywhere even where $x\in B_{r}(X_{i})$ as long as we do not use that $\left\lvert x-X_{i}\right\rvert>r$. This was only used once so that we have to look at $d_{i}(x)$ again when $x\in B_{r}(X_{i})$: If $\left\lvert x-X_{i}\right\rvert>R$, this is smaller than $CR^{3}\frac{1}{\left\lvert x-X_{i}\right\rvert^{2}}$. If $\left\lvert x-X_{i}\right\rvert\leq R$ it scales like $R$. Integrating the $p$th power of this over the union of the $B_{r}(X_{i})$ gives

[TABLE]

5.2 Passage to the Stokes equation with variable viscosity

In order to obtain the final result we want to replace the $\hat{v}$ in equation (5.4) by $\hat{u}$. First we establish a regularity lemma:

Lemma 5.5**.**

There is a constant $C>0$ such that the following holds. Let $g\in L^{2}(\mathbb{R}^{3})$ be compactly supported in $B_{2L}(0)$. Let $w\in\dot{H}_{\sigma}^{1}$ solve

[TABLE]

*Then, $w\in L^{\infty}(\mathbb{R}^{3})$ and $\left\lVert w\right\rVert_{L^{\infty}(\mathbb{R}^{3})}\leq C\left\lVert g\right\rVert_{L^{2}(\mathbb{R}^{3})}$. *

Proof 5.6**.**

We apply the fundamental solution to write

[TABLE]

Now we establish existence and estimates for the final equation:

Lemma 5.7**.**

There is a constant $C>0$ such that the following holds. The equation

[TABLE]

has a solution in $\dot{H}_{\sigma}^{1}$ and for small $\phi$ we have $\left\lVert\nabla\bar{u}\right\rVert_{L^{2}}\leq\left\lVert f\right\rVert_{L^{\frac{6}{5}}}$. Moreover, the gradient of the solution satisfies $\nabla\bar{u}\in H^{1}(\mathbb{R}^{3})$. The estimate for $\nabla^{2}\bar{u}$ is given by

[TABLE]

Proof 5.8**.**

Consider the weak formulation

[TABLE]

where $\varphi\in\dot{H}_{\sigma}^{1}$. Existence follows from Lax-Milgram theorem. For $\phi\leq\left\lVert\rho\right\rVert_{L^{\infty}}^{-1}$ we get the estimate for the gradient by setting $\varphi=\bar{u}$ and estimating the right hand side like

[TABLE]

The estimate on the second derivative is obtained by rewriting the weak formulation as

[TABLE]

*and then testing with difference quotients $\varphi=-D^{-h}_{k}D^{h}_{k}\bar{u}$, where for any function $g$ $D^{h}_{k}g(x)=\frac{1}{h}\left({g(x+he_{k})-g(x)}\right)$. *

Lemma 5.9**.**

There is a constant $C>0$ such that the following holds. The weak solutions in $\dot{H}_{\sigma}^{1}$ to the equations

[TABLE]

*differ on scale $\phi^{2}$, i.e. $\left\lVert\hat{u}-\bar{u}\right\rVert_{L^{\infty}(\mathbb{R}^{3})}\leq C\phi^{2}$. *

Proof 5.10**.**

By subtracting equation (5.1) from equation (5.9) we obtain:

[TABLE]

Hence, for the difference $\bar{u}-\hat{v}$, we get:

[TABLE]

Testing with $\bar{u}-\hat{v}$ gives

[TABLE]

On the other hand we know that $\nabla\bar{u}\in H^{1}$ and by the same argument $\nabla\hat{v}\in H^{1}$ so that we can test equation (5.10) by $-\Delta\left({\bar{u}-\hat{v}}\right)$ in order to obtain

[TABLE]

This proves that $\left\lVert\nabla\left({\bar{u}-\hat{v}}\right)\right\rVert_{H^{1}}\leq C\phi$. Now we subtract the equations for $\bar{u}$ and $\hat{u}$ to obtain for the difference $w=\bar{u}-\hat{u}$:

[TABLE]

This means that

[TABLE]

The right hand side is compactly supported in $B_{2L}(0)$ and in $L^{2}$. By Lemma 5.5 this means that

[TABLE]

Proof 5.11** (Proof of Theorem 1).**

*The statement follows by combining Theorem 6, Lemma 4.24, Lemma 4.26, Lemma 5.3 and Lemma 5.9. Note that we do not need a separate $L^{p}$ statement in Theorem 6 and in 5.9 since we have control over the $L^{\infty}$ norm of the difference on the whole space. *

Aknowledgements

We thank Juan Velázquez for suggesting to look at the dipole approximations as well as for many fruitful discussions and Richard Höfer for his hint on how to use the method of reflections to its full potential.

The authors acknowledge support through the CRC 1060 The mathematics of emergent effects at the University of Bonn that is funded through the German Science Foundation (DFG).