# Effective viscosity of random suspensions without uniform separation

Mitia Duerinckx

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mitia Duerinckx. Effective viscosity of random suspensions without uniform separation. Annales de l'Institut Henri Poincaré C, Analyse non linéaire, 2021, 39 (5), pp.1009-1052. 10.4171/AIHPC/25 . hal-03456132

HAL Id: hal-03456132

https://hal.science/hal-03456132

Submitted on 23 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Effective viscosity of random suspensions
without uniform separation

Mitia Duerinckx

**Abstract.** This work is devoted to the definition and the analysis of the effective viscosity associated with a random suspension of small rigid particles in a steady Stokes fluid. While previous works on the topic have conveniently assumed that particles are uniformly separated, we relax this restrictive assumption in the form of mild moment bounds on interparticle distances.

## 1. Introduction

Consider a colloidal suspension of small rigid particles in a Stokes fluid. Suspended particles act as obstacles, hindering the fluid flow and thus increasing the viscosity. In a recent contribution ([7]) with Gloria, we show in terms of homogenization theory that the suspension behaves at leading order like a Stokes fluid with some *effective* viscosity, and in [8] we establish optimal error estimates. In [6] we analyze the value of this effective viscosity in the low-density regime, in particular establishing the so-called Einstein formula and improving on several recent works on the topic ([14, 17–19, 22]). In [10] we further investigate the collective sedimentation of suspended particles under gravity. In all those contributions, a crucial technical assumption is that particles are uniformly separated, which is necessary in various arguments, for instance when appealing to trace estimates and regularity theory at particle boundaries. This separation assumption is however unsatisfactory from the physical viewpoint, as it is incompatible with the steady-state behavior, e.g. [1, 2], and the present contribution aims at relaxing it as much as possible in the form of mild inverse moment bounds on interparticle distances. We focus on the definition of the effective viscosity and on the qualitative homogenization result, and we further provide general tools that can be used to adapt some more advanced results; see e.g. [6, Section 2] and [18, Section 5] on the validity of Einstein's formula in the low-density regime without uniform separation.

In the case of smooth particles with some nondegeneracy condition, we essentially show in three dimensions that the effective viscosity is well defined provided that $E[\rho^{-1}] < \infty$, where $\rho$ stands for the distance between two neighboring particles, and we prove qualitative homogenization under the stronger condition $E[\rho^{-3/2}] < \infty$. Although

<footnote n="2020">2020 Mathematics Subject Classification. 35R60, 76M50, 35Q35, 76D07.
Keywords. Random suspensions.</footnote>
# Effective viscosity of random suspensions
without uniform separation

Mitia Duerinckx

**Abstract.** This work is devoted to the definition and the analysis of the effective viscosity associated with a random suspension of small rigid particles in a steady Stokes fluid. While previous works on the topic have conveniently assumed that particles are uniformly separated, we relax this restrictive assumption in the form of mild moment bounds on interparticle distances.

## 1. Introduction

Consider a colloidal suspension of small rigid particles in a Stokes fluid. Suspended particles act as obstacles, hindering the fluid flow and thus increasing the viscosity. In a recent contribution ([7]) with Gloria, we show in terms of homogenization theory that the suspension behaves at leading order like a Stokes fluid with some *effective* viscosity, and in [8] we establish optimal error estimates. In [6] we analyze the value of this effective viscosity in the low-density regime, in particular establishing the so-called Einstein formula and improving on several recent works on the topic ([14, 17–19, 22]). In [10] we further investigate the collective sedimentation of suspended particles under gravity. In all those contributions, a crucial technical assumption is that particles are uniformly separated, which is necessary in various arguments, for instance when appealing to trace estimates and regularity theory at particle boundaries. This separation assumption is however unsatisfactory from the physical viewpoint, as it is incompatible with the steady-state behavior, e.g. [1, 2], and the present contribution aims at relaxing it as much as possible in the form of mild inverse moment bounds on interparticle distances. We focus on the definition of the effective viscosity and on the qualitative homogenization result, and we further provide general tools that can be used to adapt some more advanced results; see e.g. [6, Section 2] and [18, Section 5] on the validity of Einstein's formula in the low-density regime without uniform separation.

In the case of smooth particles with some nondegeneracy condition, we essentially show in three dimensions that the effective viscosity is well defined provided that $E[\rho^{-1}] < \infty$, where $\rho$ stands for the distance between two neighboring particles, and we prove qualitative homogenization under the stronger condition $E[\rho^{-3/2}] < \infty$. Although

<footnote n="2020">2020 Mathematics Subject Classification. 35R60, 76M50, 35Q35, 76D07.
Keywords. Random suspensions.</footnote>
likely optimal in a general stationary ergodic setting, these moment bounds on interparticle distances are still quite restrictive and unphysical; cf. [1, 2]. We may draw the link with the well-known paradox of absence of solid-solid contacts in a three-dimensional Stokes flow, which is related to flaws in the modeling: real-life solid particles are slightly elastic, their boundaries display some roughness, and no-slip boundary conditions are not exactly valid; see e.g. [16] and references therein. Such corrections are not considered in the present contribution and we rather provide a detailed analysis of the ideal Stokes model. In [9], with Gloria, we investigate another line of research: under suitable mixing conditions, large clusters of close particles are unlikely in view of subcritical percolation, which can be exploited to prove homogenization without any conditions on interparticle distances. Finer geometric information might also be used in the spirit of [15].

Our approach in this contribution is mainly inspired by the work of Jikov ([24, 25]) on the homogenization problem for scalar elliptic equations with stiff inclusions; see also [20, Section 8.6]. In that scalar setting, however, required moment bounds on interparticle distances are much milder and only logarithmic moments are required in three dimensions. We emphasize two main differences:

* First, and most importantly, the incompressibility constraint in the present Stokes problem brings important rigidity and leads to completely different scalings. This is easily understood by noting that the incompressibility constraint can be eliminated by writing the Stokes equations as fourth-order elliptic equations on the vector potential; see e.g. [12]. As in [16], spatial cutoffs in this situation are then naturally to be performed on the vector potential, so that one derivative of cutoff functions is lost with respect to scalar and compressible settings, which explains the different scalings; see the proof of Proposition 3.1.

* Second, the vectorial character of the Stokes problem prohibits the use of scalar truncations: in contrast with e.g. [20, Section 8.6], this forces us to appeal to the Sobolev embedding and further deteriorates the required moment conditions.

We note some similarities with the homogenization problem for elliptic systems with degenerate random coefficients, e.g. [3–5, 11], where similar inverse moment conditions are required on coefficients.

Before stating our main results, we close this introduction by recalling the formulation of the Stokes model for a viscous fluid in the presence of a random suspension of small rigid particles, e.g. [7]. We denote by $d \ge 2$ the space dimension, and we consider a random ensemble of particles $I = \bigcup_n I_n \subset \mathbb{R}^d$. Stationarity, ergodicity, and regularity assumptions are postponed to Section 2. In order to model a dense suspension of small particles, we rescale the random set $I$ by a small parameter $\varepsilon > 0$ and consider $\varepsilon I = \bigcup_n \varepsilon I_n$. We then view these small particles $\{\varepsilon I_n\}_n$ as suspended in a solvent described by the steady Stokes equation: in a reference domain $U \subset \mathbb{R}^d$, given an internal force $f \in L^2(U)^d$, the fluid velocity $u_\varepsilon \in H^1(U \setminus \varepsilon I)^d$ satisfies

$$ - \Delta u_{\varepsilon} + \nabla S_{\varepsilon} = f, \quad \text{div}(u_{\varepsilon}) = 0, \quad \text{in } U \setminus \varepsilon I, \qquad (1.1) $$
likely optimal in a general stationary ergodic setting, these moment bounds on interparticle distances are still quite restrictive and unphysical; cf. [1, 2]. We may draw the link with the well-known paradox of absence of solid-solid contacts in a three-dimensional Stokes flow, which is related to flaws in the modeling: real-life solid particles are slightly elastic, their boundaries display some roughness, and no-slip boundary conditions are not exactly valid; see e.g. [16] and references therein. Such corrections are not considered in the present contribution and we rather provide a detailed analysis of the ideal Stokes model. In [9], with Gloria, we investigate another line of research: under suitable mixing conditions, large clusters of close particles are unlikely in view of subcritical percolation, which can be exploited to prove homogenization without any conditions on interparticle distances. Finer geometric information might also be used in the spirit of [15].

Our approach in this contribution is mainly inspired by the work of Jikov ([24, 25]) on the homogenization problem for scalar elliptic equations with stiff inclusions; see also [20, Section 8.6]. In that scalar setting, however, required moment bounds on interparticle distances are much milder and only logarithmic moments are required in three dimensions. We emphasize two main differences:

* First, and most importantly, the incompressibility constraint in the present Stokes problem brings important rigidity and leads to completely different scalings. This is easily understood by noting that the incompressibility constraint can be eliminated by writing the Stokes equations as fourth-order elliptic equations on the vector potential; see e.g. [12]. As in [16], spatial cutoffs in this situation are then naturally to be performed on the vector potential, so that one derivative of cutoff functions is lost with respect to scalar and compressible settings, which explains the different scalings; see the proof of Proposition 3.1.

* Second, the vectorial character of the Stokes problem prohibits the use of scalar truncations: in contrast with e.g. [20, Section 8.6], this forces us to appeal to the Sobolev embedding and further deteriorates the required moment conditions.

We note some similarities with the homogenization problem for elliptic systems with degenerate random coefficients, e.g. [3–5, 11], where similar inverse moment conditions are required on coefficients.

Before stating our main results, we close this introduction by recalling the formulation of the Stokes model for a viscous fluid in the presence of a random suspension of small rigid particles, e.g. [7]. We denote by $d \ge 2$ the space dimension, and we consider a random ensemble of particles $I = \bigcup_n I_n \subset \mathbb{R}^d$. Stationarity, ergodicity, and regularity assumptions are postponed to Section 2. In order to model a dense suspension of small particles, we rescale the random set $I$ by a small parameter $\varepsilon > 0$ and consider $\varepsilon I = \bigcup_n \varepsilon I_n$. We then view these small particles $\{\varepsilon I_n\}_n$ as suspended in a solvent described by the steady Stokes equation: in a reference domain $U \subset \mathbb{R}^d$, given an internal force $f \in L^2(U)^d$, the fluid velocity $u_\varepsilon \in H^1(U \setminus \varepsilon I)^d$ satisfies

$$ - \Delta u_{\varepsilon} + \nabla S_{\varepsilon} = f, \quad \text{div}(u_{\varepsilon}) = 0, \quad \text{in } U \setminus \varepsilon I, \qquad (1.1) $$
with $u_ε = 0$ on $\partial U$. (We implicitly assume here that no particle intersects the boundary.) The pressure field is only defined up to an additive constant and we choose $S_ε \in L^1(U \setminus \varepsilon I)$ with $\int_{U \setminus \varepsilon I} S_ε = 0$. Next, no-slip boundary conditions are imposed at particle boundaries: since particles are constrained to have rigid motions, this amounts to letting the velocity field $u_ε$ be extended inside particles, $u_ε \in H^1(U)^d$, with the rigidity constraint

$$D(u_\varepsilon) = 0 \quad \text{in } \varepsilon I, \qquad (1.2)$$

where $D(u_ε)$ stands for the symmetric gradient of $u_ε$. In other words, this condition means that the velocity field $u_ε$ coincides with a rigid motion $x \mapsto V_{ε,n} + \Theta_{ε,n}x$ inside each particle $εI_n$, for some $V_{ε,n} \in \mathbb{R}^d$ and some skew-symmetric matrix $\Theta_{ε,n} \in \mathbb{R}^{d \times d}$. Finally, assuming that the particles have the same mass density as the fluid, or in the absence of gravity, buoyancy forces vanish, and the force and torque balances on each particle take the form

$$\int_{\varepsilon \partial I_n} \sigma(u_\varepsilon, S_\varepsilon) \nu = 0, \qquad (1.3)$$

$$\int_{\varepsilon \partial I_n} \Theta x \cdot \sigma(u_\varepsilon, S_\varepsilon) \nu = 0 \quad \text{for all skew-symmetric } \Theta \in \mathbb{R}^{d \times d}, \qquad (1.4)$$

where $\sigma(u_ε, S_ε)$ is the Cauchy stress tensor

$$\sigma(u_\varepsilon, S_\varepsilon) = 2 D(u_\varepsilon) - S_\varepsilon \mathrm{Id}, \qquad (1.5)$$

and where $\nu$ stands for the outward unit normal vector at the particle boundaries. These equations (1.1)–(1.5) have the following weak formulation:

$$2 \int_{U} D(g) : D(u_{\varepsilon}) = \int_{U} g \cdot f, \quad \forall g \in C_{c}^{1}(U)^{d} : \operatorname{div}(g) = 0, D(g)|_{\varepsilon I} = 0.$$

This Stokes problem can also be viewed as a model for incompressible linear elasticity with stiff inclusions.

## Notation

* For vector fields $u, u'$ and matrix fields $T, T'$, we set $(\nabla u)_{ij} = \nabla_j u_i$, $\operatorname{div}(T) = \nabla_j T_{ij}$, $T : T' = T_{ij} T'_{ij}$, $(u \otimes u')_{ij} = u_i u'_j$, where we systematically use Einstein's summation convention on repeated indices. For a matrix $E$, we write $\nabla_E u = E : \nabla u$.

* For a velocity field $u$ and pressure field $S$, we denote by $(D(u))_{ij} = \frac{1}{2}(\nabla_j u_i + \nabla_i u_j)$ the symmetric gradient and by $\sigma(u, S) = 2D(u) - S \mathrm{Id}$ the Cauchy stress tensor. At particle boundaries, we let $\nu$ denote the outward unit normal vector.

* We denote by $M^{\text{sym}} \subset \mathbb{R}^{d \times d}$ the subset of symmetric matrices, by $M_0^{\text{sym}}$ the subset of symmetric trace-free matrices, and by $M^{\text{skew}}$ the subset of skew-symmetric matrices. We also write $L^p(\mathbb{R}^d)^{d \times d}_{\text{sym}} = L^p(\mathbb{R}^d; M^{\text{sym}})$.
with $u_{\varepsilon} = 0$ on $\partial U$. (We implicitly assume here that no particle intersects the boundary.) The pressure field is only defined up to an additive constant and we choose $S_{\varepsilon} \in L^1(U \setminus \varepsilon I)$ with $\int_{U \setminus \varepsilon I} S_{\varepsilon} = 0$. Next, no-slip boundary conditions are imposed at particle boundaries: since particles are constrained to have rigid motions, this amounts to letting the velocity field $u_{\varepsilon}$ be extended inside particles, $u_{\varepsilon} \in H^1(U)^d$, with the rigidity constraint

$$D(u_{\varepsilon}) = 0 \quad \text{in } \varepsilon I, \qquad (1.2)$$

where $D(u_{\varepsilon})$ stands for the symmetric gradient of $u_{\varepsilon}$. In other words, this condition means that the velocity field $u_{\varepsilon}$ coincides with a rigid motion $x \mapsto V_{\varepsilon,n} + \Theta_{\varepsilon,n}x$ inside each particle $\varepsilon I_n$, for some $V_{\varepsilon,n} \in \mathbb{R}^d$ and some skew-symmetric matrix $\Theta_{\varepsilon,n} \in \mathbb{R}^{d \times d}$. Finally, assuming that the particles have the same mass density as the fluid, or in the absence of gravity, buoyancy forces vanish, and the force and torque balances on each particle take the form

$$\int_{\varepsilon \partial I_n} \sigma(u_{\varepsilon}, S_{\varepsilon}) \nu = 0, \qquad (1.3)$$

$$\int_{\varepsilon \partial I_n} \Theta x \cdot \sigma(u_{\varepsilon}, S_{\varepsilon}) \nu = 0 \quad \text{for all skew-symmetric } \Theta \in \mathbb{R}^{d \times d}, \qquad (1.4)$$

where $\sigma(u_{\varepsilon}, S_{\varepsilon})$ is the Cauchy stress tensor

$$\sigma(u_{\varepsilon}, S_{\varepsilon}) = 2 D(u_{\varepsilon}) - S_{\varepsilon} \mathrm{Id}, \qquad (1.5)$$

and where $\nu$ stands for the outward unit normal vector at the particle boundaries. These equations (1.1)–(1.5) have the following weak formulation:

$$2 \int_{U} D(g) : D(u_{\varepsilon}) = \int_{U} g \cdot f, \quad \forall g \in C_{c}^{1}(U)^{d} : \operatorname{div}(g) = 0, D(g)|_{\varepsilon I} = 0.$$

This Stokes problem can also be viewed as a model for incompressible linear elasticity with stiff inclusions.

## Notation

* For vector fields $u, u'$ and matrix fields $T, T'$, we set $(\nabla u)_{ij} = \nabla_j u_i$, $\operatorname{div}(T) = \nabla_j T_{ij}$, $T : T' = T_{ij} T'_{ij}$, $(u \otimes u')_{ij} = u_i u'_j$, where we systematically use Einstein's summation convention on repeated indices. For a matrix $E$, we write $\nabla_E u = E : \nabla u$.

* For a velocity field $u$ and pressure field $S$, we denote by $(D(u))_{ij} = \frac{1}{2}(\nabla_j u_i + \nabla_i u_j)$ the symmetric gradient and by $\sigma(u, S) = 2 D(u) - S \mathrm{Id}$ the Cauchy stress tensor. At particle boundaries, we let $\nu$ denote the outward unit normal vector.

* We denote by $M^{\text{sym}} \subset \mathbb{R}^{d \times d}$ the subset of symmetric matrices, by $M_0^{\text{sym}}$ the subset of symmetric trace-free matrices, and by $M^{\text{skew}}$ the subset of skew-symmetric matrices. We also write $L^p(\mathbb{R}^d)^{d \times d}_{\text{sym}} = L^p(\mathbb{R}^d; M^{\text{sym}})$.
* We denote by $C \ge 1$ any constant that only depends on the dimension $d$, on the reference domain $U$, and on the parameters appearing in the different assumptions (in particular on $\delta$ in $(H_\delta^0)-(H_\delta^t)$ below). The value of the constant $C$ is allowed to change from one line to another. We use the notation $\lesssim$ (resp. $\gtrsim$) for $\le C \times$ (resp. $\ge \frac{1}{C} \times$) up to such a multiplicative constant $C$. We add subscripts to $C$, $\lesssim$, $\gtrsim$ to indicate dependence on other parameters.

* The ball centered at $x$ of radius $r$ in $\mathbb{R}^d$ is denoted by $B_r(x)$, and we simply write $B(x) = B_1(x)$, $B_r = B_r(0)$, and $B = B_1(0)$.

## 2. Main results

We focus on the case $d > 2$ for the statement of the main results, while the two-dimensional case has some important differences and is briefly discussed in Remark 3.4.

### 2.1. Assumptions

We start with the construction and suitable assumptions on the random ensemble of particles. Given an underlying probability space $(\Omega, \mathcal{A}, \mathbb{P})$, let $\mathcal{P} = \{x_n\}_n$ be a random point process on $\mathbb{R}^d$, with a given enumeration, consider a collection of random shapes $\{I_n^\circ\}_n$, where each $I_n^\circ$ is a connected random Borel subset of the unit ball $B$,<fnref n="1" /> and define the corresponding random inclusions $I_n := x_n + I_n^\circ$. We then consider the random set $\mathcal{I} := \bigcup_n I_n$, which is assumed to satisfy the following general conditions, for some deterministic constant $\delta > 0$.

**Assumption (H$_\delta^\circ$) – General conditions.**

* Stationarity and ergodicity: The point process $\mathcal{P} = \{x_n\}_n$ and the associated random set $\mathcal{I}$ are stationary and ergodic.<fnref n="2" />

* Uniform $C^2$ regularity: Random shapes $\{I_n^\circ\}_n$ almost surely satisfy interior and exterior ball conditions with radius $\delta$.

* Hardcore condition: There holds $\bar{I}_n \cap \bar{I}_m = \emptyset$ almost surely for all $n \ne m$.

When particles are close, not only does their distance matter, but also the order of their quasi-contact. We therefore need to refine the above hardcore condition, and we focus on the case of smooth particles with uniformly nonosculating boundaries. This is expressed

<footnote n="1">¹Letting $\mathcal{B}(\mathbb{R}^d)$ denote the Borel $\sigma$-algebra on $\mathbb{R}^d$, we recall that a map $I^\circ: \Omega \to \mathcal{B}(\mathbb{R}^d): \omega \mapsto I^\circ(\omega)$ is a random Borel subset of $\mathbb{R}^d$ if the set $\{(\omega, x) : \omega \in \Omega, x \in I^\circ(\omega)\}$ belongs to the product $\sigma$-algebra $\mathcal{A} \times \mathcal{B}(\mathbb{R}^d)$, or alternatively if the indicator function $1_{I^\circ}$ is $(\mathcal{A} \times \mathcal{B}(\mathbb{R}^d))$-measurable on $\Omega \times \mathbb{R}^d$.</footnote>

<footnote n="2">²Stationarity means that the laws of the translated point process $x + \mathcal{P}$ and of the translated random Borel set $x + \mathcal{I}$ do not depend on the shift $x \in \mathbb{R}^d$. Ergodicity then means that, if a measurable function of $\mathcal{P}$ or $\mathcal{I}$ is almost surely unchanged when $\mathcal{P}$ or $\mathcal{I}$ is replaced by $x + \mathcal{P}$ or $x + \mathcal{I}$ for any $x \in \mathbb{R}^d$, then the function is almost surely constant.</footnote>
* We denote by $C \ge 1$ any constant that only depends on the dimension $d$, on the reference domain $U$, and on the parameters appearing in the different assumptions (in particular on $\delta$ in $(H_\delta^0)-(H_\delta^t)$ below). The value of the constant $C$ is allowed to change from one line to another. We use the notation $\lesssim$ (resp. $\gtrsim$) for $\le C \times$ (resp. $\ge \frac{1}{C} \times$) up to such a multiplicative constant $C$. We add subscripts to $C$, $\lesssim$, $\gtrsim$ to indicate dependence on other parameters.

* The ball centered at $x$ of radius $r$ in $\mathbb{R}^d$ is denoted by $B_r(x)$, and we simply write $B(x) = B_1(x)$, $B_r = B_r(0)$, and $B = B_1(0)$.

## 2. Main results

We focus on the case $d > 2$ for the statement of the main results, while the two-dimensional case has some important differences and is briefly discussed in Remark 3.4.

### 2.1. Assumptions

We start with the construction and suitable assumptions on the random ensemble of particles. Given an underlying probability space $(\Omega, \mathcal{A}, \mathbb{P})$, let $\mathcal{P} = \{x_n\}_n$ be a random point process on $\mathbb{R}^d$, with a given enumeration, consider a collection of random shapes $\{I_n^\circ\}_n$, where each $I_n^\circ$ is a connected random Borel subset of the unit ball $B$,<fnref n="1" /> and define the corresponding random inclusions $I_n := x_n + I_n^\circ$. We then consider the random set $\mathcal{I} := \bigcup_n I_n$, which is assumed to satisfy the following general conditions, for some deterministic constant $\delta > 0$.

**Assumption (H$_\delta^\circ$) – General conditions.**

* Stationarity and ergodicity: The point process $\mathcal{P} = \{x_n\}_n$ and the associated random set $\mathcal{I}$ are stationary and ergodic.<fnref n="2" />

* Uniform $C^2$ regularity: Random shapes $\{I_n^\circ\}_n$ almost surely satisfy interior and exterior ball conditions with radius $\delta$.

* Hardcore condition: There holds $\bar{I}_n \cap \bar{I}_m = \emptyset$ almost surely for all $n \ne m$.

When particles are close, not only does their distance matter, but also the order of their quasi-contact. We therefore need to refine the above hardcore condition, and we focus on the case of smooth particles with uniformly nonosculating boundaries. This is expressed

<footnote n="1">¹Letting $\mathcal{B}(\mathbb{R}^d)$ denote the Borel $\sigma$-algebra on $\mathbb{R}^d$, we recall that a map $I^\circ: \Omega \to \mathcal{B}(\mathbb{R}^d): \omega \mapsto I^\circ(\omega)$ is a random Borel subset of $\mathbb{R}^d$ if the set $\{(\omega, x) : \omega \in \Omega, x \in I^\circ(\omega)\}$ belongs to the product $\sigma$-algebra $\mathcal{A} \times \mathcal{B}(\mathbb{R}^d)$, or alternatively if the indicator function $1_{I^\circ}$ is $(\mathcal{A} \times \mathcal{B}(\mathbb{R}^d))$-measurable on $\Omega \times \mathbb{R}^d$.</footnote>

<footnote n="2">²Stationarity means that the laws of the translated point process $x + \mathcal{P}$ and of the translated random Borel set $x + \mathcal{I}$ do not depend on the shift $x \in \mathbb{R}^d$. Ergodicity then means that, if a measurable function of $\mathcal{P}$ or $\mathcal{I}$ is almost surely unchanged when $\mathcal{P}$ or $\mathcal{I}$ is replaced by $x + \mathcal{P}$ or $x + \mathcal{I}$ for any $x \in \mathbb{R}^d$, then the function is almost surely constant.</footnote>
below in the form of some “parabolic” version of a cone condition. While always satisfied
in the case of spherical particles, this excludes for instance the case of particles that would
almost touch on flat components, as it would correspond to a contact of infinite order;
see Figures 1–2. Note that our analysis is easily adapted to intermediate situations with
contacts of any fixed order: this would lead to stronger moment conditions on interparticle
distances and is not pursued here.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** The figure displays a configuration with close particles satisfying Assumption (H'_δ). Disjoint neighborhoods {I_n^+}n are represented as light gray areas around the particles. The zooms on the neighborhoods of quasi-contact points show that particle boundaries are not osculating, as prescribed by Assumption (H'_δ), with parabolic domains delimited by dotted lines.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Figure 2.** The figure displays examples of configurations of close particles that are forbidden by Assumption (H'_δ) as their boundaries are osculating to infinite order.</figcaption></figure>

Before we actually state relevant geometric conditions, we need to introduce some
further notation. First, we construct neighborhoods {$I_n^+$}$_n$ of the particles {$I_n$}$_n$ in the
form of truncated Voronoi cells,

$$
I_n^+ := (I_n + B_\delta) \cap \{x \in \mathbb{R}^d : \operatorname{dist}(x, I_n) < \inf_{m:m \neq n} \operatorname{dist}(x, I_m)\}. \quad (2.1)
$$

In view of the uniform $C^2$ regularity of the particles, cf. (H°_δ), it is easily checked that these
neighborhoods {$I_n^+$}n are uniformly Lipschitz (with Lipschitz constant bounded by C/δ).
Next, we define “model” parabolic domains that are enclosed by close paraboloids with
different radii: given a distance ρ ≥ 0 and radii a₂ > a₁ > 0, we set

$$
\begin{equation}
\begin{aligned}
\Gamma_{a_1, a_2}^{+}(\rho) &:= B_{\delta} \cap \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : -\rho - \frac{1}{a_2} |x'|^2 < x_1 < \frac{1}{a_1} |x'|^2\}, \\
\Gamma_{a_1, a_2}^{-}(\rho) &:= B_{\delta} \cap \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : -\rho - \frac{1}{a_1} |x'|^2 < x_1 < -\frac{1}{a_2} |x'|^2\},
\end{aligned}
\tag{2.2}
\end{equation}
$$
below in the form of some “parabolic” version of a cone condition. While always satisfied
in the case of spherical particles, this excludes for instance the case of particles that would
almost touch on flat components, as it would correspond to a contact of infinite order;
see Figures 1–2. Note that our analysis is easily adapted to intermediate situations with
contacts of any fixed order: this would lead to stronger moment conditions on interparticle
distances and is not pursued here.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The figure displays a configuration with close particles satisfying Assumption (H'_δ). Disjoint neighborhoods {I_n^+}n are represented as light gray areas around the particles. The zooms on the neighborhoods of quasi-contact points show that particle boundaries are not osculating, as prescribed by Assumption (H'_δ), with parabolic domains delimited by dotted lines.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The figure displays examples of configurations of close particles that are forbidden by Assumption (H'_δ) as their boundaries are osculating to infinite order.</figcaption></figure>

Before we actually state relevant geometric conditions, we need to introduce some
further notation. First, we construct neighborhoods {$I_n^+$}$_n$ of the particles {$I_n$}$_n$ in the
form of truncated Voronoi cells,

$$
I_n^+ := (I_n + B_\delta) \cap \{x \in \mathbb{R}^d : \operatorname{dist}(x, I_n) < \inf_{m:m \neq n} \operatorname{dist}(x, I_m)\}. \quad (2.1)
$$

In view of the uniform $C^2$ regularity of the particles, cf. $(H_\delta^\circ)$, it is easily checked that these
neighborhoods $\{I_n^+\}_n$ are uniformly Lipschitz (with Lipschitz constant bounded by $C/\delta$).
Next, we define “model” parabolic domains that are enclosed by close paraboloids with
different radii: given a distance $\rho \ge 0$ and radii $a_2 > a_1 > 0$, we set

$$
\begin{equation}
\begin{aligned}
\Gamma_{a_1, a_2}^{+}(\rho) &:= B_{\delta} \cap \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : -\rho - \frac{1}{a_2} |x'|^2 < x_1 < \frac{1}{a_1} |x'|^2\}, \\
\Gamma_{a_1, a_2}^{-}(\rho) &:= B_{\delta} \cap \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : -\rho - \frac{1}{a_1} |x'|^2 < x_1 < -\frac{1}{a_2} |x'|^2\},
\end{aligned}
\tag{2.2}
\end{equation}
$$
In these terms, we formulate the following geometric condition, for some deterministic constant $\delta > 0$. It states that such parabolic domains can be included in the interparticle spacing $I_n^+ \setminus I_n$ in the neighborhood of quasi-contact points, and the condition $\frac{1}{a_1} - \frac{1}{a_2} \ge \delta$ means that paraboloids can be chosen to be $\delta$-uniformly not osculating; see Figures 1–2.

**Assumption (H'_g)** – Uniform nondegeneracy of contact points. For all $n$, for all $x \in \partial I_n$, there exists $0 \le \rho \le \delta$, there exist radii $a_2 > a_1 \ge \delta$ with $\frac{1}{a_1} - \frac{1}{a_2} \ge \delta$, and there exists a rotation $Q \in O(d)$, such that the rotated parabolic domain $x + Q\Gamma_{a_1,a_2}^+(\rho)$ or $x + Q\Gamma_{a_1,a_2}^-(\rho)$ is contained in $I_n^+ \setminus I_n$.

Finally, we turn to assumptions on interparticle distances. For all $n$, the (half) interparticle distance from $I_n$ is given by

$$ \rho_n^\circ := \min_{m:m \ne n} \frac{1}{2} \mathrm{dist}(I_n, I_m). \qquad (2.3) $$

While previous works on the Stokes model (1.1)–(1.5) have focused on the convenient case of uniformly separated particles, that is, $\inf_n \rho_n^\circ > 0$, the present contribution aims at showing that this can be substantially weakened in the form of mild inverse moment bounds. For that purpose, under Assumption (H'_g), we first need to introduce a better suited notion of interparticle distance $\rho_n \le \rho_n^\circ$: for all $x \in \partial I_n$, we let $\rho_n(x)$ denote the supremum of the admissible choices of $\rho$ in (H'_g), and we then define

$$ \rho_n := \inf_{x \in \partial I_n} \rho_n(x). \qquad (2.4) $$

## 2.2. Construction of correctors

We start with the definition of correctors for the Stokes problem (1.1)–(1.5), thus adapting [7, Proposition 2.1] to the present setting without uniform particle separation. The proof relies on the construction of a suitable admissible test function for the variational problem (2.6) below, and we believe that the moment condition (2.5) is optimal in general. As is shown in the proof, existence and uniqueness of the corrector $\psi_E$ also hold under (2.5) with $\eta = 0$, but existence of a stationary pressure field is based on a weak compactness argument in $L^1+(\Omega)$ and therefore requires $\eta > 0$. Contacts between particles are allowed in dimension $d > 5$ as no moment condition is required in that case.

**Theorem 1** (Correctors). Let $d > 2$. On top of Assumptions (H'_g) and (H'_g), assume that interparticle distances $\{\rho_n\}_n$, cf. (2.4), satisfy the following moment condition, for some $\eta > 0$:

$$ \begin{align} \text{for } d < 5: \quad & \sum_n \mathbb{E}[\rho_n^{-\frac{5-d}{2}-\eta} \mathbb{1}_{0 \in I_n}] < \infty, \tag{2.5} \\ \text{for } d = 5: \quad & \sum_n \mathbb{E}[|\log \rho_n|^{1+\eta} \mathbb{1}_{0 \in I_n}] < \infty, \notag \end{align} $$

while no moment condition is required in dimension $d > 5$. Then, for all $E \in \mathbb{M}_0^{\mathrm{sym}}$, there
In these terms, we formulate the following geometric condition, for some deterministic constant $\delta > 0$. It states that such parabolic domains can be included in the interparticle spacing $I_n^+ \setminus I_n$ in the neighborhood of quasi-contact points, and the condition $\frac{1}{a_1} - \frac{1}{a_2} \ge \delta$ means that paraboloids can be chosen to be $\delta$-uniformly not osculating; see Figures 1–2.

**Assumption (H'_g)** – Uniform nondegeneracy of contact points. For all $n$, for all $x \in \partial I_n$, there exists $0 \le \rho \le \delta$, there exist radii $a_2 > a_1 \ge \delta$ with $\frac{1}{a_1} - \frac{1}{a_2} \ge \delta$, and there exists a rotation $Q \in O(d)$, such that the rotated parabolic domain $x + Q\Gamma_{a_1,a_2}^+(\rho)$ or $x + Q\Gamma_{a_1,a_2}^-(\rho)$ is contained in $I_n^+ \setminus I_n$.

Finally, we turn to assumptions on interparticle distances. For all $n$, the (half) interparticle distance from $I_n$ is given by

$$ \rho_n^\circ := \min_{m:m \ne n} \frac{1}{2} \mathrm{dist}(I_n, I_m). \qquad (2.3) $$

While previous works on the Stokes model (1.1)–(1.5) have focused on the convenient case of uniformly separated particles, that is, $\inf_n \rho_n^\circ > 0$, the present contribution aims at showing that this can be substantially weakened in the form of mild inverse moment bounds. For that purpose, under Assumption (H'_g), we first need to introduce a better suited notion of interparticle distance $\rho_n \le \rho_n^\circ$: for all $x \in \partial I_n$, we let $\rho_n(x)$ denote the supremum of the admissible choices of $\rho$ in (H'_g), and we then define

$$ \rho_n := \inf_{x \in \partial I_n} \rho_n(x). \qquad (2.4) $$

## 2.2. Construction of correctors

We start with the definition of correctors for the Stokes problem (1.1)–(1.5), thus adapting [7, Proposition 2.1] to the present setting without uniform particle separation. The proof relies on the construction of a suitable admissible test function for the variational problem (2.6) below, and we believe that the moment condition (2.5) is optimal in general. As is shown in the proof, existence and uniqueness of the corrector $\psi_E$ also hold under (2.5) with $\eta = 0$, but existence of a stationary pressure field is based on a weak compactness argument in $L^1(\Omega)$ and therefore requires $\eta > 0$. Contacts between particles are allowed in dimension $d > 5$ as no moment condition is required in that case.

**Theorem 1** (Correctors). Let $d > 2$. On top of Assumptions (H'_g) and (H'_g), assume that interparticle distances $\{\rho_n\}_n$, cf. (2.4), satisfy the following moment condition, for some $\eta > 0$:

$$ \begin{align} \text{for } d < 5: \quad & \sum_n \mathbb{E}[\rho_n^{-\frac{5-d}{2}-\eta} \mathbb{1}_{0 \in I_n}] < \infty, \tag{2.5} \\ \text{for } d = 5: \quad & \sum_n \mathbb{E}[|\log \rho_n|^{1+\eta} \mathbb{1}_{0 \in I_n}] < \infty, \notag \end{align} $$

while no moment condition is required in dimension $d > 5$. Then, for all $E \in \mathbb{M}_0^{\mathrm{sym}}$, there
exists a unique minimizer $D(\psi_E)$ of the variational problem

$$ \inf\{\mathbb{E}[|D(\psi) + E|^2] : \psi \in L^2(\Omega; H_{\text{loc}}^1(\mathbb{R}^d)^d), \nabla\psi \text{ stationary,} \\ \qquad \text{div}(\psi) = 0, (D(\psi) + E)|_I = 0, \mathbb{E}[D(\psi)] = 0\}, \quad (2.6) $$

and the minimum value defines a positive-definite symmetric linear map $\bar{\mathbf{B}}$ on $\mathbb{M}_0^{\text{sym}}$, which is the so-called effective viscosity,

$$ E: \bar{\mathbf{B}} E := \mathbb{E}[|D(\psi_E) + E|^2]. \quad (2.7) $$

Moreover, the minimizer $D(\psi_E)$ can be characterized by the following PDE: there exist a unique random vector field $\psi_E \in L^2(\Omega; H_{\text{loc}}^1(\mathbb{R}^d)^d)$, with anchoring $\int_B \psi_E = 0$, and a unique associated pressure field $\Sigma_E \in L^1(\Omega; L_{\text{loc}}^1(\mathbb{R}^d \setminus I))$, such that

* the following equations are almost surely satisfied in the strong sense:

$$ \left\{ \begin{array}{ll} -\Delta\psi_E + \nabla\Sigma_E = 0 & \text{in } \mathbb{R}^d \setminus I, \\ \text{div}(\psi_E) = 0 & \text{in } \mathbb{R}^d, \\ D(\psi_E + Ex) = 0 & \text{in } I, \\ \int_{\partial I_n} \sigma(\psi_E + Ex, \Sigma_E)\nu = 0 & \forall n, \\ \int_{\partial I_n} \Theta(x-x_n) \cdot \sigma(\psi_E + Ex, \Sigma_E)\nu = 0 & \forall n, \forall \Theta \in \mathbb{M}^{\text{skew}}; \end{array} \right. \quad (2.8) $$

* $\nabla\psi_E$ and $\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I}$ are stationary, with the following estimates, for some $\eta > 0$:

$$ \begin{aligned} \mathbb{E}[|\nabla\psi_E|^2] &\lesssim |E|^2, & \mathbb{E}[|\nabla\psi_E|] &= 0, \\ \mathbb{E}[|\Sigma_E|^{1+\eta} \mathbb{1}_{\mathbb{R}^d \setminus I}] &\lesssim |E|^{1+\eta}, & \mathbb{E}[|\Sigma_E| \mathbb{1}_{\mathbb{R}^d \setminus I}] &= 0. \end{aligned} $$

In particular, the following convergences hold almost surely as $\varepsilon \downarrow 0$:

$$ \begin{aligned} (\nabla \psi_E)(\cdot_e) &\to 0 & \text{weakly in } L_{\text{loc}}^2(\mathbb{R}^d), \\ (\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I})(\cdot_e) &\to 0 & \text{weakly in } L_{\text{loc}}^{1+\eta}(\mathbb{R}^d), \\ \varepsilon \psi_E(\cdot_e) &\to 0 & \text{strongly in } L_{\text{loc}}^q(\mathbb{R}^d), & \text{for all } q < \frac{2d}{d-2}. \end{aligned} \quad (2.9) $$

In contrast with the case of uniformly separated particles, cf. [7, Proposition 2.1], we emphasize that under the moment condition (2.5) the pressure field $\Sigma_E$ above is only defined in $L^{1+\eta}(\Omega)$ for some $\eta > 0$, and not in $L^2(\Omega)$. Improving on this integrability naturally requires a stronger moment condition, as shown in the following.

**Proposition 2** (Integrability of the pressure). Let $d > 2$ and let $\gamma := \frac{2d(d+1)}{d^2+5d-2}$ for abbreviation. On top of Assumptions $(H_\delta^\circ)$ and $(H'_\delta)$, given $1 < \alpha < 2$, let one of the following conditions hold for interparticle distances $\{\rho_n\}_n$:
exists a unique minimizer $D(\psi_E)$ of the variational problem

$$ \inf\{\mathbb{E}[|D(\psi) + E|^2] : \psi \in L^2(\Omega; H_{\text{loc}}^1(\mathbb{R}^d)^d), \nabla\psi \text{ stationary,} \\ \qquad \text{div}(\psi) = 0, (D(\psi) + E)|_I = 0, \mathbb{E}[D(\psi)] = 0\}, \quad (2.6) $$

and the minimum value defines a positive-definite symmetric linear map $\bar{\mathbf{B}}$ on $\mathbb{M}_0^{\text{sym}}$, which is the so-called effective viscosity,

$$ E: \bar{\mathbf{B}} E := \mathbb{E}[|D(\psi_E) + E|^2]. \quad (2.7) $$

Moreover, the minimizer $D(\psi_E)$ can be characterized by the following PDE: there exist a unique random vector field $\psi_E \in L^2(\Omega; H_{\text{loc}}^1(\mathbb{R}^d)^d)$, with anchoring $\int_B \psi_E = 0$, and a unique associated pressure field $\Sigma_E \in L^1(\Omega; L_{\text{loc}}^1(\mathbb{R}^d \setminus I))$, such that

* the following equations are almost surely satisfied in the strong sense:

$$ \left\{
\begin{array}{ll}
-\Delta\psi_E + \nabla\Sigma_E = 0 & \text{in } \mathbb{R}^d \setminus I, \\
\text{div}(\psi_E) = 0 & \text{in } \mathbb{R}^d, \\
D(\psi_E + Ex) = 0 & \text{in } I, \\
\displaystyle\int_{\partial I_n} \sigma(\psi_E + Ex, \Sigma_E)\nu = 0 & \forall n, \\
\displaystyle\int_{\partial I_n} \Theta(x-x_n) \cdot \sigma(\psi_E + Ex, \Sigma_E)\nu = 0 & \forall n, \forall \Theta \in \mathbb{M}^{\text{skew}};
\end{array}
\right.
\quad (2.8) $$

* $\nabla\psi_E$ and $\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I}$ are stationary, with the following estimates, for some $\eta > 0$:

$$ \begin{aligned}
\mathbb{E}[|\nabla\psi_E|^2] &\lesssim |E|^2, & \mathbb{E}[|\nabla\psi_E|] &= 0, \\
\mathbb{E}[|\Sigma_E|^{1+\eta} \mathbb{1}_{\mathbb{R}^d \setminus I}] &\lesssim |E|^{1+\eta}, & \mathbb{E}[|\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I}|] &= 0.
\end{aligned} $$

In particular, the following convergences hold almost surely as $\varepsilon \downarrow 0$:

$$ \begin{aligned}
(\nabla \psi_E)(\dot{\varepsilon}) &\to 0 \quad \text{weakly in } L_{\text{loc}}^2(\mathbb{R}^d), \\
(\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I})(\dot{\varepsilon}) &\to 0 \quad \text{weakly in } L_{\text{loc}}^{1+\eta}(\mathbb{R}^d), \\
\varepsilon \psi_E(\dot{\varepsilon}) &\to 0 \quad \text{strongly in } L_{\text{loc}}^q(\mathbb{R}^d), \text{ for all } q < \frac{2d}{d-2}.
\end{aligned} \quad (2.9) $$

In contrast with the case of uniformly separated particles, cf. [7, Proposition 2.1], we emphasize that under the moment condition (2.5) the pressure field $\Sigma_E$ above is only defined in $L^{1+\eta}(\Omega)$ for some $\eta > 0$, and not in $L^2(\Omega)$. Improving on this integrability naturally requires a stronger moment condition, as shown in the following.

**Proposition 2** (Integrability of the pressure). Let $d > 2$ and let $\gamma := \frac{2d(d+1)}{d^2+5d-2}$ for abbreviation. On top of Assumptions $(H_\delta^\circ)$ and $(H'_\delta)$, given $1 < \alpha < 2$, let one of the following conditions hold for interparticle distances $\{\rho_n\}_n$:
• in the case $\alpha \le \frac{d}{d-1}$ with $d \le 5$, assume that

$$ \text{for } d < 5: \quad \sum_n \mathbb{E}[\rho_n^{-\frac{\alpha}{2-\alpha} \frac{5-d}{2}} 1_{0 \in I_n}] < \infty, $$

$$ \text{for } d = 5: \quad \sum_n \mathbb{E}[|\log \rho_n|^{\frac{\alpha}{2-\alpha}} 1_{0 \in I_n}] < \infty, $$

• in the case $\alpha < \gamma$ with $d > 5$, no moment condition is required;

• in the case $\frac{d}{d-1} < \gamma < \alpha < 2$, with $\alpha \ne \frac{d}{d-2}$, assume that

$$ \sum_n \mathbb{E}[\rho_n^{-\frac{\alpha}{2-\alpha}(\frac{1}{\gamma}-\frac{1}{\alpha})(d+1)} 1_{0 \in I_n}] < \infty. $$

Then for all $E \in \mathbb{M}_0^{\text{sym}}$ the pressure field $\Sigma_E$ constructed in Theorem 1 satisfies

$$ \mathbb{E}[|\Sigma_E|^\alpha 1_{\mathbb{R}^d \setminus I}] \lesssim |E|^\alpha, $$

and there holds almost surely $(\Sigma_E 1_{\mathbb{R}^d \setminus I})(\frac{1}{\varepsilon}) \to 0$ weakly in $L_{\text{loc}}^\alpha(\mathbb{R}^d)$ as $\varepsilon \downarrow 0$.

## 2.3. Homogenization result

We turn to the homogenization result for the Stokes problem (1.1)–(1.5). For that purpose, we first define admissible random ensembles of particles in a given bounded Lipschitz domain $U \subset \mathbb{R}^d$: the proof indeed requires controlling the distances of particles to the boundary $\partial U$ similarly to interparticle distances. We let $\mathcal{N}_\varepsilon(U) \subset \mathbb{N}$ denote a random subset of indices such that

$$ \{n : I_n \subset \frac{1}{\varepsilon}U, \text{dist}(I_n, \partial \frac{1}{\varepsilon}U) \ge \delta\} \subset \mathcal{N}_\varepsilon(U) \subset \{n : I_n \subset \frac{1}{\varepsilon}U\}, $$

and we define the associated random ensemble of particles in $U$,

$$ I_{\varepsilon}(U) := \bigcup_{n \in \mathcal{N}_{\varepsilon}(U)} \varepsilon I_n. \qquad (2.10) $$

In this setting we consider corresponding neighborhoods $\{I_{n;U,\varepsilon}^+\}_n$ of the particles $\{I_n\}_n$,

$$ I_{n;U,\varepsilon}^{+} := I_{n}^{+} \cap \frac{1}{\varepsilon}U, $$

we assume that Assumption (H'_δ) holds with neighborhoods $\{I_n^+\}_n$ replaced by $\{I_{n;U,\varepsilon}^+\}_n$, and we define the corresponding distances $\{\rho_{n;U,\varepsilon}\}_n$ as in (2.4).

With this notation we may now formulate the homogenization result for (1.1)–(1.5). The proof is based on a div-curl argument together with an extension result for fluxes as inspired by the work of Jikov ([24, 25]). Due to nonuniform particle separation, extended fluxes are only controlled in $L^\alpha$ for some integrability $\alpha < 2$ depending on the moment condition on interparticle distances; see Theorem 4. In view of the Sobolev embedding, Jikov’s div-curl argument can then be performed provided $\alpha \ge \frac{2d}{d+2}$. This restriction leads
* in the case $\alpha \le \frac{d}{d-1}$ with $d \le 5$, assume that

$$ \text{for } d < 5: \quad \sum_n \mathbb{E}[\rho_n^{-\frac{\alpha}{2-\alpha} \frac{5-d}{2}} 1_{0 \in I_n}] < \infty, $$

$$ \text{for } d = 5: \quad \sum_n \mathbb{E}[|\log \rho_n|^{\frac{\alpha}{2-\alpha}} 1_{0 \in I_n}] < \infty, $$

* in the case $\alpha < \gamma$ with $d > 5$, no moment condition is required;

* in the case $\frac{d}{d-1} < \gamma < \alpha < 2$, with $\alpha \ne \frac{d}{d-2}$, assume that

$$ \sum_n \mathbb{E}[\rho_n^{-\frac{\alpha}{2-\alpha}(\frac{1}{\gamma}-\frac{1}{\alpha})(d+1)} 1_{0 \in I_n}] < \infty. $$

Then for all $E \in \mathbb{M}_0^{\text{sym}}$ the pressure field $\Sigma_E$ constructed in Theorem 1 satisfies

$$ \mathbb{E}[|\Sigma_E|^\alpha 1_{\mathbb{R}^d \setminus I}] \lesssim |E|^\alpha, $$

and there holds almost surely $(\Sigma_E 1_{\mathbb{R}^d \setminus I})(\frac{1}{\epsilon}) \to 0$ weakly in $L_{\text{loc}}^\alpha(\mathbb{R}^d)$ as $\epsilon \downarrow 0$.

## 2.3. Homogenization result

We turn to the homogenization result for the Stokes problem (1.1)–(1.5). For that purpose, we first define admissible random ensembles of particles in a given bounded Lipschitz domain $U \subset \mathbb{R}^d$: the proof indeed requires controlling the distances of particles to the boundary $\partial U$ similarly to interparticle distances. We let $\mathcal{N}_\epsilon(U) \subset \mathbb{N}$ denote a random subset of indices such that

$$ \{n : I_n \subset \frac{1}{\epsilon}U, \text{dist}(I_n, \partial \frac{1}{\epsilon}U) \ge \delta\} \subset \mathcal{N}_\epsilon(U) \subset \{n : I_n \subset \frac{1}{\epsilon}U\}, $$

and we define the associated random ensemble of particles in $U$,

$$ I_{\epsilon}(U) := \bigcup_{n \in \mathcal{N}_{\epsilon}(U)} \epsilon I_n. \qquad (2.10) $$

In this setting we consider corresponding neighborhoods $\{I_{n;U,\epsilon}^+\}_n$ of the particles $\{I_n\}_n$,

$$ I_{n;U,\epsilon}^{+} := I_{n}^{+} \cap \frac{1}{\epsilon}U, $$

we assume that Assumption (H'_δ) holds with neighborhoods $\{I_n^+\}_n$ replaced by $\{I_{n;U,\epsilon}^+\}_n$, and we define the corresponding distances $\{\rho_{n;U,\epsilon}\}_n$ as in (2.4).

With this notation we may now formulate the homogenization result for (1.1)–(1.5). The proof is based on a div-curl argument together with an extension result for fluxes as inspired by the work of Jikov ([24, 25]). Due to nonuniform particle separation, extended fluxes are only controlled in $L^\alpha$ for some integrability $\alpha < 2$ depending on the moment condition on interparticle distances; see Theorem 4. In view of the Sobolev embedding, Jikov’s div-curl argument can then be performed provided $\alpha \ge \frac{2d}{d+2}$. This restriction leads
to the moment condition (2.11) below, which is expected to be optimal in general and coincides with the one in Proposition 2 with $\alpha = \frac{2d}{d+2}$. We emphasize that this condition becomes more stringent in large dimensions as the Sobolev exponent $\frac{2d}{d+2}$ increases to 2. Not surprisingly, the condition is stronger than the one for the existence of the corrector in Theorem 1 since defining correctors only requires constructing an admissible test function for the variational problem (2.6).

**Theorem 3 (Homogenization result).** Let $d > 2$. On top of Assumption (H°), given a bounded Lipschitz domain $U \subset \mathbb{R}^d$, let Assumption (H') hold for $\{I_{n;U,\varepsilon}^+\}_n$, and assume that interparticle distances $\{\rho_{n;U,\varepsilon}\}_n$ satisfy, almost surely,

$$ 
\begin{aligned}
\text{for } d=3: \quad & \limsup_{\varepsilon \downarrow 0} \varepsilon^d \sum_{n \in \mathcal{N}_\varepsilon(U)} (\rho_{n;U,\varepsilon})^{-\frac{3}{2}} < \infty, \\
\text{for } d \ge 4: \quad & \limsup_{\varepsilon \downarrow 0} \varepsilon^d \sum_{n \in \mathcal{N}_\varepsilon(U)} (\rho_{n;U,\varepsilon})^{-\frac{d}{2}-1} < \infty,
\end{aligned} 
\quad (2.11) $$

where in the case $d=6$ the exponent $\frac{d}{2}-1=2$ must be replaced by some exponent $> 2$. Denote by $\lambda := \mathbb{E}[1_I]$ the volume fraction of the suspension, let $\psi, \Sigma, \bar{\mathbf{B}}$ be defined as in Theorem 1, and define the following effective constant $\bar{\mathbf{b}} \in \mathbb{M}_0^{\text{sym}}$: for all $E \in \mathbb{M}_0^{\text{sym}}$,

$$ \bar{\mathbf{b}}: E := \frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1_{I_n}}{|I_n|} \int_{\partial I_n} (x-x_n) \cdot \sigma(\psi_E + Ex, \Sigma_E) \nu \right]. \quad (2.12) $$

Given an internal force $f \in L^2(U)^d$, let the velocity field $u_\varepsilon \in L^2(\Omega; H_0^1(U)^d)$ and the associated pressure field $S_\varepsilon \in L^1(\Omega; L^1(U \setminus I_\varepsilon(U)))$, with anchoring $\int_{U \setminus I_\varepsilon(U)} S_\varepsilon = 0$, be almost surely the unique solutions of the Stokes problem (1.1)–(1.5), that is,

$$ 
\left\{
\begin{array}{ll}
-\Delta u_{\varepsilon} + \nabla S_{\varepsilon} = f & \text{in } U \setminus I_{\varepsilon}(U), \\
\text{div}(u_{\varepsilon}) = 0 & \text{in } U, \\
D(u_{\varepsilon}) = 0 & \text{in } I_{\varepsilon}(U), \\
\displaystyle\int_{\varepsilon \partial I_n} \sigma(u_{\varepsilon}, S_{\varepsilon}) \nu = 0 & \forall n, \\
\displaystyle\int_{\varepsilon \partial I_n} \Theta(x - \varepsilon x_n) \cdot \sigma(u_{\varepsilon}, S_{\varepsilon}) \nu = 0 & \forall n, \forall \Theta \in \mathbb{M}^{\text{skew}}.
\end{array}
\right.
\quad (2.13) $$

Then we have almost surely, as $\varepsilon \downarrow 0$,

$$ 
\begin{array}{l}
u_{\varepsilon} - \bar{u} \to 0 \quad \text{weakly in } H_0^1(U), \\
(S_{\varepsilon} - \bar{S} - \bar{\mathbf{b}} : D(\bar{u})) 1_{U \setminus I_{\varepsilon}(U)} \to 0 \quad \text{weakly in } L^{\frac{2d}{d+2}}(U),
\end{array}
$$

where the limiting velocity field $\bar{u} \in H_0^1(U)^d$ and the associated pressure field $\bar{S} \in L^2(U)$, with anchoring $\int_U \bar{S} = 0$, are the unique solutions of the homogenized equation

$$ 
\left\{
\begin{array}{ll}
-\text{div}(2\bar{\mathbf{B}} \, \text{D}(\bar{u})) + \nabla \bar{S} = (1-\lambda)f & \text{in } U, \\
\text{div}(\bar{u}) = 0 & \text{in } U.
\end{array}
\right.
\quad (2.14) $$
to the moment condition (2.11) below, which is expected to be optimal in general and coincides with the one in Proposition 2 with $\alpha = \frac{2d}{d+2}$. We emphasize that this condition becomes more stringent in large dimensions as the Sobolev exponent $\frac{2d}{d+2}$ increases to 2. Not surprisingly, the condition is stronger than the one for the existence of the corrector in Theorem 1 since defining correctors only requires constructing an admissible test function for the variational problem (2.6).

**Theorem 3 (Homogenization result).** Let $d > 2$. On top of Assumption (H°), given a bounded Lipschitz domain $U \subset \mathbb{R}^d$, let Assumption (H') hold for $\{I_{n;U,\varepsilon}^+\}_n$, and assume that interparticle distances $\{\rho_{n;U,\varepsilon}\}_n$ satisfy, almost surely,

$$ 
\begin{aligned}
& \text{for } d=3: \quad \limsup_{\varepsilon \downarrow 0} \varepsilon^d \sum_{n \in \mathcal{N}_\varepsilon(U)} (\rho_{n;U,\varepsilon})^{-\frac{3}{2}} < \infty, \\
& \text{for } d \ge 4: \quad \limsup_{\varepsilon \downarrow 0} \varepsilon^d \sum_{n \in \mathcal{N}_\varepsilon(U)} (\rho_{n;U,\varepsilon})^{-\frac{d}{2}-1} < \infty,
\end{aligned} 
\quad (2.11) $$

where in the case $d=6$ the exponent $\frac{d}{2}-1=2$ must be replaced by some exponent $> 2$. Denote by $\lambda := \mathbb{E}[1_I]$ the volume fraction of the suspension, let $\psi, \Sigma, \bar{\mathbf{B}}$ be defined as in Theorem 1, and define the following effective constant $\bar{\mathbf{b}} \in \mathbb{M}_0^{\text{sym}}$: for all $E \in \mathbb{M}_0^{\text{sym}}$,

$$ \bar{\mathbf{b}}: E := \frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1_{I_n}}{|I_n|} \int_{\partial I_n} (x-x_n) \cdot \sigma(\psi_E + Ex, \Sigma_E) \nu \right]. \quad (2.12) $$

Given an internal force $f \in L^2(U)^d$, let the velocity field $u_\varepsilon \in L^2(\Omega; H_0^1(U)^d)$ and the associated pressure field $S_\varepsilon \in L^1(\Omega; L^1(U \setminus I_\varepsilon(U)))$, with anchoring $\int_{U \setminus I_\varepsilon(U)} S_\varepsilon = 0$, be almost surely the unique solutions of the Stokes problem (1.1)–(1.5), that is,

$$ 
\left\{
\begin{array}{ll}
-\Delta u_{\varepsilon} + \nabla S_{\varepsilon} = f & \text{in } U \setminus I_{\varepsilon}(U), \\
\text{div}(u_{\varepsilon}) = 0 & \text{in } U, \\
D(u_{\varepsilon}) = 0 & \text{in } I_{\varepsilon}(U), \\
\displaystyle\int_{\varepsilon \partial I_n} \sigma(u_{\varepsilon}, S_{\varepsilon}) \nu = 0 & \forall n, \\
\displaystyle\int_{\varepsilon \partial I_n} \Theta(x - \varepsilon x_n) \cdot \sigma(u_{\varepsilon}, S_{\varepsilon}) \nu = 0 & \forall n, \forall \Theta \in \mathbb{M}^{\text{skew}}.
\end{array}
\right.
\quad (2.13) $$

Then we have almost surely, as $\varepsilon \downarrow 0$,

$$ 
\begin{array}{l}
u_{\varepsilon} - \bar{u} \to 0 \quad \text{weakly in } H_0^1(U), \\
(S_{\varepsilon} - \bar{S} - \bar{\mathbf{b}} : D(\bar{u})) 1_{U \setminus I_{\varepsilon}(U)} \to 0 \quad \text{weakly in } L^{\frac{2d}{d+2}}(U),
\end{array}
 $$

where the limiting velocity field $\bar{u} \in H_0^1(U)^d$ and the associated pressure field $\bar{S} \in L^2(U)$, with anchoring $\int_U \bar{S} = 0$, are the unique solutions of the homogenized equation

$$ 
\left\{
\begin{array}{ll}
-\text{div}(2\bar{\mathbf{B}} \, \text{D}(\bar{u})) + \nabla \bar{S} = (1-\lambda)f & \text{in } U, \\
\text{div}(\bar{u}) = 0 & \text{in } U.
\end{array}
\right.
\quad (2.14) $$
In addition, provided that $f \in L^p(U)^d$ for some $p > d$, the following corrector results hold almost surely, as $\varepsilon \downarrow 0$:

$$ 
\begin{align} 
& \left\| u_{\varepsilon} - \bar{u} - \sum_{E \in \mathcal{E}} \varepsilon \psi_E \left( \frac{\cdot}{\varepsilon} \right) \nabla_E \bar{u} \right\|_{H^1(U)} \to 0, \tag{2.15} \\ 
& \inf_{c \in \mathbb{R}} \left\| S_{\varepsilon} - \bar{S} - \bar{b} : D(\bar{u}) - \sum_{E \in \mathcal{E}} \Sigma_E \left( \frac{\cdot}{\varepsilon} \right) \nabla_E \bar{u} - c \right\|_{L^{\frac{2d}{d+2}}(U \setminus \mathcal{I})} \to 0. \nonumber
\end{align} 
$$

where $\mathcal{E}$ stands for an orthonormal basis of $M_0^{\text{sym}}$.

## 2.4. Further technical tools

On top of the definition of the effective viscosity and the above qualitative homogenization result, we wish to further extend more advanced results such as the validity of Einstein's formula for the effective viscosity at low density ([6, 18]), optimal quantitative error estimates for homogenization ([8]), and the analysis of sedimentation ([10]). For these aims, we provide a couple of technical tools for the analysis of suspensions without uniform separation. These tools are used in [6, Section 2] and [18, Section 5] for the validity of Einstein's formula.

We start with the following extension result for fluxes in the presence of rigid particles, which constitutes the main technical tool in our proof of Theorem 3. Starting from a notion of flux $q$ that accounts for the behavior outside rigid particles, we construct an extension $\tilde{q}$ that is defined nontrivially inside the particles in such a way that the continuity equation holds globally; cf. (2.18). For that purpose, one views the suspension of rigid particles as the limit of a suspension of droplets with diverging shear viscosity, and extended fluxes are then naturally defined as limits of corresponding fluxes; see Remark 4.2. This construction is inspired by a corresponding scalar result by Jikov ([24, 25]) in the context of scalar elliptic equations with stiff inclusions (see also [20, Section 3.5]), but additional care is needed here to deal with the incompressibility constraint.

**Theorem 4** (Extension of fluxes). Let $d > 2$, let Assumptions (H$_{\delta}^{\circ}$) and (H$_{\delta}^{\prime}$) hold, and let a realization of the random set $\mathcal{I}$ be fixed. Given $\beta \in (1, \infty)$ and $f \in L_{\text{loc}}^{\beta}(\mathbb{R}^d)^d$, let $q \in L_{\text{loc}}^{\beta}(\mathbb{R}^d)^{d \times d}_{\text{sym}}$ with $\text{tr}(q) = 0$ satisfy

$$ 2 \int_{\mathbb{R}^d} D(g) : q = \int_{\mathbb{R}^d} g \cdot f \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \text{div}(g) = 0, \ D(g)|_{\mathcal{I}} = 0. \quad (2.16) $$

Then, for all $\alpha$, $r$ chosen as

$$ r \ge \frac{\beta}{\beta - 1}, \quad 1 < \alpha \le \beta \wedge \frac{dr\beta}{r(d-\beta) + d\beta}, \quad \text{with} \quad \begin{cases} r < \frac{d\beta}{\beta - d} & \text{if } \beta > d, \\ r \ne \frac{d\beta}{d\beta - d - \beta} & \text{if } \beta > \frac{d}{d-1}, \\ \alpha < \frac{d}{d-1} & \text{if } r = \frac{\beta}{\beta-1}, \end{cases} \quad (2.17) $$
In addition, provided that $f \in L^p(U)^d$ for some $p > d$, the following corrector results hold almost surely, as $\varepsilon \downarrow 0$:

$$ \begin{gathered} \left\| u_{\varepsilon} - \bar{u} - \sum_{E \in \mathcal{E}} \varepsilon \psi_E \left(\frac{\cdot}{\varepsilon}\right) \nabla_E \bar{u} \right\|_{H^1(U)} \to 0, \quad (2.15) \\ \inf_{c \in \mathbb{R}} \left\| S_{\varepsilon} - \bar{S} - \bar{b} : D(\bar{u}) - \sum_{E \in \mathcal{E}} \Sigma_E \left(\frac{\cdot}{\varepsilon}\right) \nabla_E \bar{u} - c \right\|_{L^{\frac{2d}{d+2}}(U \setminus \mathcal{I})} \to 0. \end{gathered} $$

where $\mathcal{E}$ stands for an orthonormal basis of $M_0^{\text{sym}}$.

## 2.4. Further technical tools

On top of the definition of the effective viscosity and the above qualitative homogenization result, we wish to further extend more advanced results such as the validity of Einstein's formula for the effective viscosity at low density ([6, 18]), optimal quantitative error estimates for homogenization ([8]), and the analysis of sedimentation ([10]). For these aims, we provide a couple of technical tools for the analysis of suspensions without uniform separation. These tools are used in [6, Section 2] and [18, Section 5] for the validity of Einstein's formula.

We start with the following extension result for fluxes in the presence of rigid particles, which constitutes the main technical tool in our proof of Theorem 3. Starting from a notion of flux $q$ that accounts for the behavior outside rigid particles, we construct an extension $\tilde{q}$ that is defined nontrivially inside the particles in such a way that the continuity equation holds globally; cf. (2.18). For that purpose, one views the suspension of rigid particles as the limit of a suspension of droplets with diverging shear viscosity, and extended fluxes are then naturally defined as limits of corresponding fluxes; see Remark 4.2. This construction is inspired by a corresponding scalar result by Jikov ([24, 25]) in the context of scalar elliptic equations with stiff inclusions (see also [20, Section 3.5]), but additional care is needed here to deal with the incompressibility constraint.

**Theorem 4** (Extension of fluxes). Let $d > 2$, let Assumptions (H$_{\delta}^{\circ}$) and (H$_{\delta}^{\prime}$) hold, and let a realization of the random set $\mathcal{I}$ be fixed. Given $\beta \in (1, \infty)$ and $f \in L_{\text{loc}}^{\beta}(\mathbb{R}^d)^d$, let $q \in L_{\text{loc}}^{\beta}(\mathbb{R}^d)^{d \times d}_{\text{sym}}$ with $\text{tr}(q) = 0$ satisfy

$$ 2 \int_{\mathbb{R}^d} D(g) : q = \int_{\mathbb{R}^d} g \cdot f \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \text{div}(g) = 0, \ D(g)|_{\mathcal{I}} = 0. \quad (2.16) $$

Then, for all $\alpha$, $r$ chosen as

$$ r \ge \frac{\beta}{\beta - 1}, \quad 1 < \alpha \le \beta \wedge \frac{dr\beta}{r(d-\beta) + d\beta}, \quad \text{with} \quad \begin{cases} r < \frac{d\beta}{\beta - d} & \text{if } \beta > d, \\ r \ne \frac{d\beta}{d\beta - d - \beta} & \text{if } \beta > \frac{d}{d-1}, \\ \alpha < \frac{d}{d-1} & \text{if } r = \frac{\beta}{\beta-1}, \end{cases} \quad (2.17) $$
there exists an extension $\tilde{q} \in L_{\text{loc}}^{\alpha}(\mathbb{R}^d)^{d \times d}$ with $\text{tr}(\tilde{q}) = 0$, as well as an associated pressure field $\tilde{S} \in L_{\text{loc}}^{\alpha}(\mathbb{R}^d)$, such that

$$ \tilde{q}|_{\mathbb{R}^d \setminus \mathcal{I}} = q|_{\mathbb{R}^d \setminus \mathcal{I}}, \quad \text{and} \quad -\text{div}(2\tilde{q} - \tilde{S} \text{Id}) = f \quad \text{in } \mathbb{R}^d, \qquad (2.18) $$

and such that, for all $R \ge 1$, the estimate

$$ \begin{aligned} & \|\tilde{q}\|_{L^\alpha(B_R)} + \left\| \tilde{S} - \int_{B_R} \tilde{S} \right\|_{L^\alpha(B_R)} \\ & \lesssim_{\alpha, \beta, r} \Lambda(B_R; r, \frac{\beta\alpha}{\beta-\alpha}) (\|f\|_{L^{d+\beta}(\hat{B}_R)} + \|q\|_{L^\beta(\hat{B}_R \setminus \mathcal{I})}) \end{aligned} \qquad (2.19) $$

holds, where we have set $\hat{B}_R := B_R \cup \bigcup_{n: I_n \cap B_R \neq \emptyset} I_n^+$ and

$$ \Lambda(D; r, p) := \left( |D| + \sum_{n: I_n \cap D \neq \emptyset} \mu_r(\rho_n)^p \right)^{\frac{1}{p}}, \qquad (2.20) $$

in terms of

$$ \mu_r(\rho) := \begin{cases} \rho^{\frac{d+1}{2r} - \frac{3}{2}} & : r > \frac{d+1}{3}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{3}, \\ 1 & : r < \frac{d+1}{3}. \end{cases} \qquad (2.21) $$

As applications of this extension result, we establish a trace estimate at particle boundaries and a version of Caccioppoli's inequality.

**Corollary 5** (Trace estimate). Let $d > 2$, let Assumptions (H) and (H') hold, and let a realization of the random set $\mathcal{I}$ be fixed. Let the velocity field $u \in H_{\text{loc}}^{1}(\mathbb{R}^d)^d$ and the associated pressure field $S \in L_{\text{loc}}^{1}(\mathbb{R}^d \setminus \mathcal{I})$ satisfy the homogeneous Stokes problem

$$ \left\{
\begin{array}{ll}
-\Delta u + \nabla S = 0 & \text{in } \mathbb{R}^d \setminus \mathcal{I}, \\
\text{div}(u) = 0 & \text{in } \mathbb{R}^d, \\
D(u) = 0 & \text{in } \mathcal{I}, \\
\displaystyle\int_{\partial I_n} \sigma(u, S)v = 0 & \forall n, \\
\displaystyle\int_{\partial I_n} \Theta(x-x_n) \cdot \sigma(u, S)v = 0 & \forall n, \forall \Theta \in \mathbb{M}^{\text{skew}}.
\end{array}
\right.
\qquad (2.22) $$

Then for all $n$ and $g \in W^{1,\infty}(I_n^+)^d$ we have for all $\eta > 0$,

$$ \left| \int_{\partial I_n} g \cdot \sigma(u, S) v \right| \lesssim_{\eta} \|g\|_{W^{1,\infty}(I_n^+ \setminus I_n)} \left( \int_{I_n^+ \setminus I_n} |D(u)|^2 \right)^{\frac{1}{2}} \\ \times \begin{cases} \rho_n^{\frac{1}{4d}(d+1)(d+2) - \frac{5}{2} - \eta} & : d \le 6, \\ 1 & : d > 6. \end{cases} $$
there exists an extension $\tilde{q} \in L_{\text{loc}}^{\alpha}(\mathbb{R}^d)^{d \times d}$ with $\text{tr}(\tilde{q}) = 0$, as well as an associated pressure field $\tilde{S} \in L_{\text{loc}}^{\alpha}(\mathbb{R}^d)$, such that

$$ \tilde{q}|_{\mathbb{R}^d \setminus \mathcal{I}} = q|_{\mathbb{R}^d \setminus \mathcal{I}}, \quad \text{and} \quad -\text{div}(2\tilde{q} - \tilde{S} \text{Id}) = f \quad \text{in } \mathbb{R}^d, \qquad (2.18) $$

and such that, for all $R \ge 1$, the estimate

$$ \begin{aligned} & \|\tilde{q}\|_{L^\alpha(B_R)} + \left\| \tilde{S} - \int_{B_R} \tilde{S} \right\|_{L^\alpha(B_R)} \\ & \lesssim_{\alpha, \beta, r} \Lambda(B_R; r, \frac{\beta\alpha}{\beta-\alpha}) (\|f\|_{L^{d+\beta}(\hat{B}_R)} + \|q\|_{L^\beta(\hat{B}_R \setminus \mathcal{I})}) \end{aligned} \qquad (2.19) $$

holds, where we have set $\hat{B}_R := B_R \cup \bigcup_{n: I_n \cap B_R \neq \emptyset} I_n^+$ and

$$ \Lambda(D; r, p) := \left( |D| + \sum_{n: I_n \cap D \neq \emptyset} \mu_r(\rho_n)^p \right)^{\frac{1}{p}}, \qquad (2.20) $$

in terms of

$$ \mu_r(\rho) := \begin{cases} \rho^{\frac{d+1}{2r} - \frac{3}{2}} & : r > \frac{d+1}{3}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{3}, \\ 1 & : r < \frac{d+1}{3}. \end{cases} \qquad (2.21) $$

As applications of this extension result, we establish a trace estimate at particle boundaries and a version of Caccioppoli's inequality.

**Corollary 5** (Trace estimate). Let $d > 2$, let Assumptions (H) and (H') hold, and let a realization of the random set $\mathcal{I}$ be fixed. Let the velocity field $u \in H_{\text{loc}}^{1}(\mathbb{R}^d)^d$ and the associated pressure field $S \in L_{\text{loc}}^{1}(\mathbb{R}^d \setminus \mathcal{I})$ satisfy the homogeneous Stokes problem

$$ \left\{
\begin{array}{ll}
-\Delta u + \nabla S = 0 & \text{in } \mathbb{R}^d \setminus \mathcal{I}, \\
\text{div}(u) = 0 & \text{in } \mathbb{R}^d, \\
D(u) = 0 & \text{in } \mathcal{I}, \\
\displaystyle\int_{\partial I_n} \sigma(u, S)v = 0 & \forall n, \\
\displaystyle\int_{\partial I_n} \Theta(x-x_n) \cdot \sigma(u, S)v = 0 & \forall n, \forall \Theta \in \mathbb{M}^{\text{skew}}.
\end{array}
\right.
\qquad (2.22) $$

Then for all $n$ and $g \in W^{1,\infty}(I_n^+)^d$ we have for all $\eta > 0$,

$$ \left| \int_{\partial I_n} g \cdot \sigma(u, S) v \right| \lesssim_{\eta} \|g\|_{W^{1,\infty}(I_n^+ \setminus I_n)} \left( \int_{I_n^+ \setminus I_n} |D(u)|^2 \right)^{\frac{1}{2}} \\ \times \begin{cases} \rho_n^{\frac{1}{4d}(d+1)(d+2) - \frac{5}{2} - \eta} & : d \le 6, \\ 1 & : d > 6. \end{cases} $$
**Corollary 6** (Caccioppmoli's inequality). Let $d > 2$, let Assumptions $(H_\delta^\circ)$ and $(H_\delta')$ hold, and let a realization of the random set $\mathcal{I}$ be fixed. Then, for all $\eta > 0$, there exists $s < \frac{2d}{d-2}$ such that any solution $(u, S)$ of the homogeneous Stokes problem (2.22) satisfies, for all $R \ge 5$ and $K \ge 1$,

$$ \left( \fint_{B_R} |\nabla u|^2 \right)^{\frac{1}{2}} \\ \lesssim_{s,\eta} \left( KR^{-1} \left( \fint_{B_{2R}} |u - \fint_{B_{2R}} u|^s \right)^{\frac{1}{s}} + (K^{-1} + R^{-\frac{d}{2}(\frac{1}{s} - \frac{d-2}{2d})}) \left( \fint_{B_{2R}} |\nabla u|^2 \right)^{\frac{1}{2}} \right) \\ \times \begin{cases} 1 + R^{-d} \sum_{n: I_n \cap B_{2R} \neq \emptyset} \rho_n^{\frac{1}{4}(d+1)(d+2) - \frac{5}{2}d - \eta} & : d \le 5, \\ 1 + R^{-d} \sum_{n: I_n \cap B_{2R} \neq \emptyset} \rho_n^{1 - \frac{d}{2} - \eta} & : d > 5. \end{cases} $$

### 3. Extension of fluxes

This section is devoted to the proof of Theorem 4. The argument relies on the following local extension result for incompressible fields, which is of independent interest.

**Proposition 3.1.** Let $d > 2$, let Assumptions $(H_\delta^\circ)$ and $(H_\delta')$ hold, and let a realization of the random set $\mathcal{I}$ be fixed. Let $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$. Then, for all $n$, there exists an extension operator $P_n$ such that for all $g \in C_b^1(I_n)^d$ with $\text{div}(g) = 0$ the extension $P_n g \in W_0^{1,s}(I_n^+)^d$ satisfies

$$ \mathrm{D}(P_n g)|_{I_n} = \mathrm{D}(g), \quad \text{and} \quad \mathrm{div}(P_n g) = 0 \text{ in } I_n^+, \qquad (3.1) $$

and for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \| \nabla P_n g \|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \| \mathrm{D}(g) \|_{L^p(I_n)}, \qquad (3.2) $$

where we recall the notation (2.21) for $\mu_r$.

For future reference, we also highlight the following key tool for pressure estimates. It follows from the above local extension result combined with a standard use of the Bogovskii operator. Note that the restriction on the geometry of the domain $D$ and the associated constant $K(D)$ can be refined as e.g. in [13, Lemma III.3.2 and Theorem III.3.1].

**Lemma 3.2.** Let $d > 2$, let Assumptions $(H_\delta^\circ)$ and $(H_\delta')$ hold, and let a realization of the random set $\mathcal{I}$ be fixed. Let $D \subset \mathbb{R}^d$ be a bounded Lipschitz domain that is star-like with respect to every point in some ball of radius $R_0$, and set $K(D) := \frac{1}{R_0} \text{diam}(D)$. Let $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$. Then, for all $h \in C_b(D)$ with $\int_{D \setminus I} h = 0$, there exists $z \in W_0^{1,s}(D)^d$ such that

$$ \mathrm{D}(z)|_I = 0, \quad \text{and} \quad \mathrm{div}(z) = h \mathbb{1}_{D \setminus I} \text{ in } D, $$
**Corollary 6** (Caccioppoli's inequality). Let $d > 2$, let Assumptions $(H_\delta^\circ)$ and $(H_\delta')$ hold, and let a realization of the random set $\mathcal{I}$ be fixed. Then, for all $\eta > 0$, there exists $s < \frac{2d}{d-2}$ such that any solution $(u, S)$ of the homogeneous Stokes problem (2.22) satisfies, for all $R \ge 5$ and $K \ge 1$,

$$ \left( \int_{B_R} |\nabla u|^2 \right)^{\frac{1}{2}} \\ \lesssim_{s,\eta} \left( KR^{-1} \left( \int_{B_{2R}} |u - f| |u|^s \right)^{\frac{1}{s}} + (K^{-1} + R^{-\frac{d}{2}(\frac{1}{s} - \frac{d}{2d})}) \left( \int_{B_{2R}} |\nabla u|^2 \right)^{\frac{1}{2}} \right) \\ \times \begin{cases} 1 + R^{-d} \sum_{n: I_n \cap B_{2R} \neq \emptyset} \rho_n^{\frac{1}{4}(d+1)(d+2) - \frac{5}{2}d - \eta} & : d \le 5, \\ 1 + R^{-d} \sum_{n: I_n \cap B_{2R} \neq \emptyset} \rho_n^{1 - \frac{d}{2} - \eta} & : d > 5. \end{cases} $$

### 3. Extension of fluxes

This section is devoted to the proof of Theorem 4. The argument relies on the following local extension result for incompressible fields, which is of independent interest.

**Proposition 3.1.** Let $d > 2$, let Assumptions $(H_\delta^\circ)$ and $(H_\delta')$ hold, and let a realization of the random set $\mathcal{I}$ be fixed. Let $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$. Then, for all $n$, there exists an extension operator $P_n$ such that for all $g \in C_b^1(I_n)^d$ with $\mathrm{div}(g) = 0$ the extension $P_n g \in W_0^{1,s}(I_n^+)^d$ satisfies

$$ \mathrm{D}(P_n g)|_{I_n} = \mathrm{D}(g), \quad \text{and} \quad \mathrm{div}(P_n g) = 0 \text{ in } I_n^+, \qquad (3.1) $$

and for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla P_n g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \|\mathrm{D}(g)\|_{L^p(I_n)}, \qquad (3.2) $$

where we recall the notation (2.21) for $\mu_r$.

For future reference, we also highlight the following key tool for pressure estimates. It follows from the above local extension result combined with a standard use of the Bogovskii operator. Note that the restriction on the geometry of the domain $D$ and the associated constant $K(D)$ can be refined as e.g. in [13, Lemma III.3.2 and Theorem III.3.1].

**Lemma 3.2.** Let $d > 2$, let Assumptions $(H_\delta^\circ)$ and $(H_\delta')$ hold, and let a realization of the random set $\mathcal{I}$ be fixed. Let $D \subset \mathbb{R}^d$ be a bounded Lipschitz domain that is star-like with respect to every point in some ball of radius $R_0$, and set $K(D) := \frac{1}{R_0} \mathrm{diam}(D)$. Let $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$. Then, for all $h \in C_b(D)$ with $\int_{D\setminus I} h = 0$, there exists $z \in W_0^{1,s}(D)^d$ such that

$$ \mathrm{D}(z)|_I = 0, \quad \text{and} \quad \mathrm{div}(z) = h \mathbb{1}_{D\setminus I} \text{ in } D, $$
and for all $p \ge s \lor \frac{drs}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla z\|_{L^s(D)} \lesssim_{p,r,s} K(D)^{d+1} \Lambda(D; r, \frac{ps}{p-s}) \|h\|_{L^p(D\setminus I)} \quad (3.3) $$

where we recall the definition (2.20)–(2.21) of $\Lambda$.

### 3.1. Cutoff functions

We start with the construction of suitable cutoff functions for the inclusions $\{I_n\}_n$ in their neighborhoods $\{I_n^+\}_n$. The open subsets $\{J_n^j\}_j$ in the statement below are neighborhoods of quasi-contact points, that is, neighborhoods where $\partial I_n$ and $\partial I_n^+$ are very close; see Figure 3. The proof is inspired by the work of Jikov on homogenization problems with stiff inclusions, e.g. [20, Section 3.2], and is also analogous to computations by Gérard-Varet and Hillairet in [16] for the drag force on a sphere close to a wall. This result is easily adapted beyond Assumption (H$_\delta^+$) to cover higher-order quasi-contacts between the particles, then leading to a worse dependence on the distance $\rho_n$.

<figure><img src="image_3.png" /><figcaption>**Figure 3.** This displays a configuration of close particles. Disjoint neighborhoods $\{I_n^+\}_n$ are represented around the particles, and suitable neighborhoods $\{J_n^j\}_j$ of quasi-contact points are drawn in light gray.</figcaption></figure>

**Lemma 3.3 (Cutoff functions).** Let Assumptions (H$_\delta^+$) and (H$_\delta^+$) hold, and let a realization of the random set $\mathcal{I}$ be fixed. For all $n$, there exists a function $w_n \in W_0^{1,\infty}(I_n^+; [0, 1])$ with $w_n|_{I_n} = 1$ such that for all $r \ge 1$,

$$ \|\nabla w_n\|_{L^r(I_n^+)} \lesssim_r \begin{cases} \rho_n^{\frac{d+1}{2r}-1} & : r > \frac{d+1}{2}, \\ |\log \rho_n|^r & : r = \frac{d+1}{2}, \\ 1 & : r < \frac{d+1}{2}, \end{cases} \quad (3.4) $$

and

$$ \|\nabla^2 w_n\|_{L^r(I_n^+)} \lesssim_r \begin{cases} \rho_n^{\frac{d+1}{2r}-2} & : r > \frac{d+1}{4}, \\ |\log \rho_n|^r & : r = \frac{d+1}{4}, \\ 1 & : r < \frac{d+1}{4}. \end{cases} \quad (3.5) $$
and for all $p \ge s \lor \frac{drs}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla z\|_{L^s(D)} \lesssim_{p,r,s} K(D)^{d+1} \Lambda(D; r, \frac{ps}{p-s}) \|h\|_{L^p(D\setminus I)} \quad (3.3) $$

where we recall the definition (2.20)–(2.21) of $\Lambda$.

### 3.1. Cutoff functions

We start with the construction of suitable cutoff functions for the inclusions $\{I_n\}_n$ in their neighborhoods $\{I_n^+\}_n$. The open subsets $\{J_n^j\}_j$ in the statement below are neighborhoods of quasi-contact points, that is, neighborhoods where $\partial I_n$ and $\partial I_n^+$ are very close; see Figure 3. The proof is inspired by the work of Jikov on homogenization problems with stiff inclusions, e.g. [20, Section 3.2], and is also analogous to computations by Gérard-Varet and Hillairet in [16] for the drag force on a sphere close to a wall. This result is easily adapted beyond Assumption (H$^{\circ}_{\delta}$) to cover higher-order quasi-contacts between the particles, then leading to a worse dependence on the distance $\rho_n$.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** This displays a configuration of close particles. Disjoint neighborhoods $\{I_n^+\}_n$ are represented around the particles, and suitable neighborhoods $\{J_n^j\}_j$ of quasi-contact points are drawn in light gray.</figcaption></figure>

**Lemma 3.3 (Cutoff functions).** Let Assumptions (H$^{\circ}_{\delta}$) and (H$^{\prime}_{\delta}$) hold, and let a realization of the random set $\mathcal{I}$ be fixed. For all $n$, there exists a function $w_n \in W_0^{1,\infty}(I_n^+; [0, 1])$ with $w_n|_{I_n} = 1$ such that for all $r \ge 1$,

$$ \|\nabla w_n\|_{L^r(I_n^+)} \lesssim_r \begin{cases} \rho_n^{\frac{d+1}{2r}-1} & : r > \frac{d+1}{2}, \\ |\log \rho_n|^r & : r = \frac{d+1}{2}, \\ 1 & : r < \frac{d+1}{2}, \end{cases} \quad (3.4) $$

and

$$ \|\nabla^2 w_n\|_{L^r(I_n^+)} \lesssim_r \begin{cases} \rho_n^{\frac{d+1}{2r}-2} & : r > \frac{d+1}{4}, \\ |\log \rho_n|^r & : r = \frac{d+1}{4}, \\ 1 & : r < \frac{d+1}{4}. \end{cases} \quad (3.5) $$
In addition, there is a collection ${J_n^j}_{j=1}^{M_n}$ of open subsets of the form $J_n^j = B(x_n^j, \frac{1}{C}\delta) \cap I_n^+$, with $M_n \lesssim 1$ and $\text{dist}(J_n^j, J_n^k) \ge \frac{1}{C}\delta$ for all $j \ne k$, such that

$$ \|w_n\|_{W^{2,\infty}(I_n^+ \cup \bigcup_{j=1}^{M_n} J_n^j)} \lesssim 1, $$

and for all $r \ge 1$,

$$ \max_{1 \le j \le M_n} \|\cdot - x_n^j |\nabla^2 w_n\|_{L^r(J_n^j)} \lesssim_r \begin{cases} \rho_n^{\frac{d+1}{2r}-\frac{3}{2}} & : r > \frac{d+1}{3}, \\ |\log \rho_n|^{\frac{1}{r}} & : r = \frac{d+1}{3}, \\ 1 & : r < \frac{d+1}{3}. \end{cases} \quad (3.6) $$

*Proof.* Under Assumption (H'_δ), the construction of the neighborhoods ${J_n^j}_{j=1}^{M_n}$ is transparent, cf. Figure 3, and we only need to construct $w_n$ in one of those sets. In view of the definition of the parabolic domains $\Gamma_{a_1,a_2}^{\pm}(\rho)$, cf. (2.2), it suffices to construct a cutoff function $w_{a_1,a_2}^{\rho}$ in $B_{\delta} \subset \mathbb{R} \times \mathbb{R}^{d-1}$ such that $w_{a_1,a_2}^{\rho} = 0$ for $x_1 < -\rho + \frac{1}{a_2}|x'|^2$ and $w_{a_1,a_2}^{\rho} = 1$ for $x_1 > \frac{1}{a_1}|x'|^2$. By assumption we consider $a_2 > a_1 \ge \delta$ with $\frac{1}{a_1} - \frac{1}{a_2} \ge \delta$, and by scaling it suffices to consider $a_1 = 1$. More precisely, we consider the set

$$ E = \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : x_1 \ge -1, |x'| \le \frac{1}{2}\}, $$

and, given $\rho > 0$ and $a > 1$ with $1 - \frac{1}{a} \ge \delta$, we construct a cutoff function $w_a^\rho \in C_b^{1,1}(E)$ such that $w_a^\rho = 0$ for $x_1 < -\rho + \frac{1}{a}|x'|^2$ and $w_a^\rho = 1$ for $x_1 > |x'|^2$, and such that for all $r \ge 1$,

$$ \|\nabla w_a^\rho\|_{L^r(E)} \lesssim_r \begin{cases} \rho^{\frac{d+1}{2r}-1} & : r > \frac{d+1}{2}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{2}, \\ 1 & : r < \frac{d+1}{2}, \end{cases} \quad (3.7) $$

$$ \|\nabla^2 w_a^\rho\|_{L^r(E)} \lesssim_r \begin{cases} \rho^{\frac{d+1}{2r}-2} & : r > \frac{d+1}{4}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{4}, \\ 1 & : r < \frac{d+1}{4}, \end{cases} \quad (3.8) $$

$$ \| \cdot \|_{\nabla^2 w_a^\rho} \|_{L^r(E)} \lesssim_r \begin{cases} \rho^{\frac{d+1}{2r}-\frac{3}{2}} & : r > \frac{d+1}{3}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{3}, \\ 1 & : r < \frac{d+1}{3}. \end{cases} \quad (3.9) $$

In other words, we need to construct a suitable interpolation between 1 and 0 in the domain enclosed by the two parabolas,

$$ \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : -\rho + \frac{1}{a}|x'|^2 < x_1 < |x'|^2, |x'| \le \frac{1}{2}\}. $$

As we aim to construct a $C^{1,1}$ test function, we cannot use linear interpolation: instead of the linear function $h^0(t) = t$ with $h^0(0) = 0$ and $h^0(1) = 1$, we rather consider as in [16] the cubic function

$$ h(t) := t^2(3 - 2t), $$
In addition, there is a collection ${\{J_n^j\}_{j=1}^{M_n}}$ of open subsets of the form $J_n^j = B(x_n^j, \frac{1}{C}\delta) \cap I_n^+$, with $M_n \lesssim 1$ and $\text{dist}(J_n^j, J_n^k) \ge \frac{1}{C}\delta$ for all $j \ne k$, such that

$$ \|w_n\|_{W^{2,\infty}(I_n^+ \cup \bigcup_{j=1}^{M_n} J_n^j)} \lesssim 1, $$

and for all $r \ge 1$,

$$ \max_{1 \le j \le M_n} \|\cdot - x_n^j |\nabla^2 w_n\|_{L^r(J_n^j)} \lesssim_r \begin{cases} \rho_n^{\frac{d+1}{2r}-\frac{3}{2}} & : r > \frac{d+1}{3}, \\ |\log \rho_n|^{\frac{1}{r}} & : r = \frac{d+1}{3}, \\ 1 & : r < \frac{d+1}{3}. \end{cases} \quad (3.6) $$

*Proof.* Under Assumption (H'_δ), the construction of the neighborhoods ${\{J_n^j\}_{j=1}^{M_n}}$ is transparent, cf. Figure 3, and we only need to construct $w_n$ in one of those sets. In view of the definition of the parabolic domains $\Gamma_{a_1,a_2}^{\pm}(\rho)$, cf. (2.2), it suffices to construct a cutoff function $w_{a_1,a_2}^{\rho}$ in $B_{\delta} \subset \mathbb{R} \times \mathbb{R}^{d-1}$ such that $w_{a_1,a_2}^{\rho} = 0$ for $x_1 < -\rho + \frac{1}{a_2}|x'|^2$ and $w_{a_1,a_2}^{\rho} = 1$ for $x_1 > \frac{1}{a_1}|x'|^2$. By assumption we consider $a_2 > a_1 \ge \delta$ with $\frac{1}{a_1} - \frac{1}{a_2} \ge \delta$, and by scaling it suffices to consider $a_1 = 1$. More precisely, we consider the set

$$ E = \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : x_1 \ge -1, |x'| \le \frac{1}{2}\}, $$

and, given $\rho > 0$ and $a > 1$ with $1 - \frac{1}{a} \ge \delta$, we construct a cutoff function $w_a^\rho \in C_b^{1,1}(E)$ such that $w_a^\rho = 0$ for $x_1 < -\rho + \frac{1}{a}|x'|^2$ and $w_a^\rho = 1$ for $x_1 > |x'|^2$, and such that for all $r \ge 1$,

$$ \|\nabla w_a^\rho\|_{L^r(E)} \lesssim_r \begin{cases} \rho^{\frac{d+1}{2r}-1} & : r > \frac{d+1}{2}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{2}, \\ 1 & : r < \frac{d+1}{2}, \end{cases} \quad (3.7) $$

$$ \|\nabla^2 w_a^\rho\|_{L^r(E)} \lesssim_r \begin{cases} \rho^{\frac{d+1}{2r}-2} & : r > \frac{d+1}{4}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{4}, \\ 1 & : r < \frac{d+1}{4}, \end{cases} \quad (3.8) $$

$$ \| |\cdot| \nabla^2 w_a^\rho \|_{L^r(E)} \lesssim_r \begin{cases} \rho^{\frac{d+1}{2r}-\frac{3}{2}} & : r > \frac{d+1}{3}, \\ |\log \rho|^{\frac{1}{r}} & : r = \frac{d+1}{3}, \\ 1 & : r < \frac{d+1}{3}. \end{cases} \quad (3.9) $$

In other words, we need to construct a suitable interpolation between 1 and 0 in the domain enclosed by the two parabolas,

$$ \{(x_1, x') \in \mathbb{R} \times \mathbb{R}^{d-1} : -\rho + \frac{1}{a}|x'|^2 < x_1 < |x'|^2, |x'| \le \frac{1}{2}\}. $$

As we aim to construct a $C^{1,1}$ test function, we cannot use linear interpolation: instead of the linear function $h^0(t) = t$ with $h^0(0) = 0$ and $h^0(1) = 1$, we rather consider as in [16] the cubic function

$$ h(t) := t^2(3 - 2t), $$
with $h(0) = 0$, $h(1) = 1$, and $h'(0) = h'(1) = 0$. We then define

$$w_a^\rho(x) := w_a^\rho(x_1, x') := \begin{cases} 0 & : x_1 \le -\rho + \frac{1}{a}|x'|^2, \\ h\left(\frac{1}{\theta_a^\rho(x')}(\rho + x_1 - \frac{1}{a}|x'|^2)\right) & : -\rho + \frac{1}{a}|x'|^2 \le x_1 \le |x'|^2, \\ 1 & : x_1 \ge |x'|^2, \end{cases}$$

where for abbreviation we denote by $\theta_a^\rho$ the horizontal distance between the two parabolas,

$$\theta_a^\rho(x') := \rho + (1 - \frac{1}{a})|x'|^2.$$

We check that $w_a^\rho$ belongs to $C^{1,1}(E)$ and it remains to establish the estimates (3.7)–(3.9). Recalling the assumption $1 - \frac{1}{a} \ge \delta$, a direct computation shows that there holds for $-\rho + \frac{1}{a}|x'|^2 \le x_1 \le |x'|^2$,

$$|\nabla w_a^\rho(x)| \lesssim (\rho + |x'|^2)^{-1}, \quad |\nabla^2 w_a^\rho(x)| \lesssim (\rho + |x'|^2)^{-2}. \quad (3.10)$$

We start with the proof of (3.7). Using (3.10), evaluating the integral over $x_1$, and using radial coordinates, we find

$$\int_E |\nabla w_a^\rho|^r \lesssim_r \int_{|x'| \le \frac{1}{2}} (\rho + |x'|^2)^{1-r} dx' \lesssim \int_0^{\frac{1}{2}} \frac{s^{d-2}}{(\rho + s^2)^{r-1}} ds,$$

which proves (3.7) after evaluating the integral. The proof of (3.8) follows the same lines and is skipped.

We turn to the proof of (3.9). For $-\rho + \frac{1}{a}|x'|^2 \le x_1 \le |x'|^2$ and $|x'| \le \frac{1}{2}$, we find

$$|x| \le |x_1| + |x'| \lesssim \rho + |x'|.$$

Combining this with (3.10), evaluating the integral over $x_1$, and using radial coordinates, we find

$$\int_E |\cdot|^r |\nabla^2 w_a^\rho|^r \lesssim_r \rho^r \int_{|x'| \le \frac{1}{2}} (\rho + |x'|^2)^{1-2r} dx' + \int_{|x'| \le \frac{1}{2}} |x'|^r (\rho + |x'|^2)^{1-2r} dx' \\ \lesssim \rho^r \int_0^{\frac{1}{2}} \frac{s^{d-2}}{(\rho + s^2)^{2r-1}} ds + \int_0^{\frac{1}{2}} \frac{s^{d+r-2}}{(\rho + s^2)^{2r-1}} ds,$$

which proves (3.9) after evaluating the integrals.

## 3.2. Proof of Proposition 3.1

Starting from a Stein extension of $g$, the argument relies on the cutoff function constructed in Lemma 3.3 in order to make this extension vanish on the boundary $\partial I_n^+$. A naive cutoff would however break the incompressibility property and cause serious trouble, especially close to quasi-contact points. Instead, in the spirit of [12], taking inspiration from calculations by Gérard-Varet and Hillairet in [16], we take cutoffs at the level of the vector potential. We split the proof into three main steps.
with $h(0) = 0$, $h(1) = 1$, and $h'(0) = h'(1) = 0$. We then define

$$ w_a^\rho(x) := w_a^\rho(x_1, x') := \begin{cases} 0 & : x_1 \le -\rho + \frac{1}{a}|x'|^2, \\ h\left(\frac{1}{\theta_a^\rho(x')}(\rho + x_1 - \frac{1}{a}|x'|^2)\right) & : -\rho + \frac{1}{a}|x'|^2 \le x_1 \le |x'|^2, \\ 1 & : x_1 \ge |x'|^2, \end{cases} $$

where for abbreviation we denote by $\theta_a^\rho$ the horizontal distance between the two parabolas,

$$ \theta_a^\rho(x') := \rho + (1 - \frac{1}{a})|x'|^2. $$

We check that $w_a^\rho$ belongs to $C^{1,1}(E)$ and it remains to establish the estimates (3.7)–(3.9). Recalling the assumption $1 - \frac{1}{a} \ge \delta$, a direct computation shows that there holds for $-\rho + \frac{1}{a}|x'|^2 \le x_1 \le |x'|^2$,

$$ |\nabla w_a^\rho(x)| \lesssim (\rho + |x'|^2)^{-1}, \quad |\nabla^2 w_a^\rho(x)| \lesssim (\rho + |x'|^2)^{-2}. \quad (3.10) $$

We start with the proof of (3.7). Using (3.10), evaluating the integral over $x_1$, and using radial coordinates, we find

$$ \int_E |\nabla w_a^\rho|^r \lesssim_r \int_{|x'| \le \frac{1}{2}} (\rho + |x'|^2)^{1-r} dx' \lesssim \int_0^{\frac{1}{2}} \frac{s^{d-2}}{(\rho + s^2)^{r-1}} ds, $$

which proves (3.7) after evaluating the integral. The proof of (3.8) follows the same lines and is skipped.

We turn to the proof of (3.9). For $-\rho + \frac{1}{a}|x'|^2 \le x_1 \le |x'|^2$ and $|x'| \le \frac{1}{2}$, we find

$$ |x| \le |x_1| + |x'| \lesssim \rho + |x'|. $$

Combining this with (3.10), evaluating the integral over $x_1$, and using radial coordinates, we find

$$ \begin{aligned} \int_E |\cdot|^r |\nabla^2 w_a^\rho|^r & \lesssim_r \rho^r \int_{|x'| \le \frac{1}{2}} (\rho + |x'|^2)^{1-2r} dx' + \int_{|x'| \le \frac{1}{2}} |x'|^r (\rho + |x'|^2)^{1-2r} dx' \\ & \lesssim \rho^r \int_0^{\frac{1}{2}} \frac{s^{d-2}}{(\rho + s^2)^{2r-1}} ds + \int_0^{\frac{1}{2}} \frac{s^{d+r-2}}{(\rho + s^2)^{2r-1}} ds, \end{aligned} $$

which proves (3.9) after evaluating the integrals.

## 3.2. Proof of Proposition 3.1

Starting from a Stein extension of $g$, the argument relies on the cutoff function constructed in Lemma 3.3 in order to make this extension vanish on the boundary $\partial I_n^+$. A naive cutoff would however break the incompressibility property and cause serious trouble, especially close to quasi-contact points. Instead, in the spirit of [12], taking inspiration from calculations by Gérard-Varet and Hillairet in [16], we take cutoffs at the level of the vector potential. We split the proof into three main steps.
Step 1. Extension to $I_n^+$.

Given $g \in C_b^1(I_n)^d$ with $\text{div}(g) = 0$, we construct an extension $P_n^1 g \in H_0^1(I_n + B)^d$ such that

$$ \mathrm{D}(P_n^1 g)|_{I_n} = \mathrm{D}(g), \quad \text{and} \quad \text{div}(P_n^1 g) = 0 \text{ in } I_n + B, \qquad (3.11) $$

and for all $1 < s < \infty$,

$$ \|\nabla P_n^1 g\|_{L^s(I_n+B)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.12) $$

For that purpose, let us first choose $V_g \in \mathbb{R}^d$ and $\Theta_g \in \mathbb{M}^{\text{skew}}$ such that Korn's inequality yields for all $1 < s < \infty$,

$$ \|g - V_g - \Theta_g x\|_{W^{1,s}(I_n)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.13) $$

Next, in view of the $C^2$ regularity of $I_n$, cf. Assumption (H$_6^0$), we can choose a Stein extension $P_n^0 g \in C_b^1(I_n + B)^d$ with $P_n^0 g|_{I_n} = g - V_g - \Theta_g x$, such that for all $s \ge 1$,

$$ \|P_n^0 g\|_{W^{1,s}(I_n+B)} \lesssim_s \|g - V_g - \Theta_g x\|_{W^{1,s}(I_n)}, $$

and thus, by (3.13), for all $1 < s < \infty$,

$$ \|P_n^0 g\|_{W^{1,s}(I_n+B)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.14) $$

It remains to apply a cutoff to $P_n^0 g$ to make it vanish on the boundary $\partial(I_n + B)$ while keeping the properties in (3.11). For that purpose, choose a cutoff function $\chi \in C_c^\infty(I_n + B)$ with $\chi|_{I_n} = 1$ and $\|\chi\|_{W^{1,\infty}(I_n+B)} \lesssim 1$. By a standard construction based on the Bogovskii operator, e.g. [13, Theorem III.3.1], since the compatibility relation

$$ \int_{(I_n+B)\setminus I_n} \mathrm{div}(\chi P_n^0 g) = -\int_{\partial I_n} g \cdot \nu = -\int_{I_n} \mathrm{div}(g) = 0 $$

holds, there exists $z_n(g) \in H_0^1((I_n + B) \setminus I_n)^d$ such that

$$ \mathrm{div}(z_n(g)) = \mathrm{div}(\chi P_n^0 g) \quad \text{in} \ (I_n + B) \setminus I_n, $$

and for all $1 < s < \infty$,

$$ \|\nabla z_n(g)\|_{L^s((I_n+B)\setminus I_n)} \lesssim_s \|\mathrm{div}(\chi P_n^0 g)\|_{L^s((I_n+B)\setminus I_n)}. $$

Expanding the divergence in the right-hand side of this estimate, and using (3.14), we find for all $1 < s < \infty$,

$$ \|\nabla z_n(g)\|_{L^s((I_n+B)\setminus I_n)} \lesssim_s \|P_n^0 g\|_{W^{1,s}(I_n+B)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.15) $$

Now we define

$$ P_n^1 g := \chi P_n^0 g - z_n(g) \quad \in H_0^1(I_n + B)^d, $$
Step 1. Extension to $I_n^+$.

Given $g \in C_b^1(I_n)^d$ with $\text{div}(g) = 0$, we construct an extension $P_n^1 g \in H_0^1(I_n + B)^d$ such that

$$ \mathrm{D}(P_n^1 g)|_{I_n} = \mathrm{D}(g), \quad \text{and} \quad \mathrm{div}(P_n^1 g) = 0 \text{ in } I_n + B, \qquad (3.11) $$

and for all $1 < s < \infty$,

$$ \|\nabla P_n^1 g\|_{L^s(I_n+B)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.12) $$

For that purpose, let us first choose $V_g \in \mathbb{R}^d$ and $\Theta_g \in \mathbb{M}^{\text{skew}}$ such that Korn's inequality yields for all $1 < s < \infty$,

$$ \|g - V_g - \Theta_g x\|_{W^{1,s}(I_n)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.13) $$

Next, in view of the $C^2$ regularity of $I_n$, cf. Assumption (H$_6^0$), we can choose a Stein extension $P_n^0 g \in C_b^1(I_n + B)^d$ with $P_n^0 g|_{I_n} = g - V_g - \Theta_g x$, such that for all $s \ge 1$,

$$ \|P_n^0 g\|_{W^{1,s}(I_n+B)} \lesssim_s \|g - V_g - \Theta_g x\|_{W^{1,s}(I_n)}, $$

and thus, by (3.13), for all $1 < s < \infty$,

$$ \|P_n^0 g\|_{W^{1,s}(I_n+B)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.14) $$

It remains to apply a cutoff to $P_n^0 g$ to make it vanish on the boundary $\partial(I_n + B)$ while keeping the properties in (3.11). For that purpose, choose a cutoff function $\chi \in C_c^\infty(I_n + B)$ with $\chi|_{I_n} = 1$ and $\|\chi\|_{W^{1,\infty}(I_n+B)} \lesssim 1$. By a standard construction based on the Bogovskii operator, e.g. [13, Theorem III.3.1], since the compatibility relation

$$ \int_{(I_n+B)\setminus I_n} \mathrm{div}(\chi P_n^0 g) = -\int_{\partial I_n} g \cdot \nu = -\int_{I_n} \mathrm{div}(g) = 0 $$

holds, there exists $z_n(g) \in H_0^1((I_n + B) \setminus I_n)^d$ such that

$$ \mathrm{div}(z_n(g)) = \mathrm{div}(\chi P_n^0 g) \quad \text{in} \ (I_n + B) \setminus I_n, $$

and for all $1 < s < \infty$,

$$ \|\nabla z_n(g)\|_{L^s((I_n+B)\setminus I_n)} \lesssim_s \|\mathrm{div}(\chi P_n^0 g)\|_{L^s((I_n+B)\setminus I_n)}. $$

Expanding the divergence in the right-hand side of this estimate, and using (3.14), we find for all $1 < s < \infty$,

$$ \|\nabla z_n(g)\|_{L^s((I_n+B)\setminus I_n)} \lesssim_s \|P_n^0 g\|_{W^{1,s}(I_n+B)} \lesssim_s \|\mathrm{D}(g)\|_{L^s(I_n)}. \qquad (3.15) $$

Now we define

$$ P_n^1 g := \chi P_n^0 g - z_n(g) \quad \in H_0^1(I_n + B)^d, $$
which indeed satisfies $P_n^1 g|_{I_n} = P_n^0 g|_{I_n} = g - V_g - \Theta_g x$ and $\text{div}(P_n^1 g) = 0$, hence (3.11). In addition, combining (3.14) and (3.15) yields for all $1 < s < \infty$,

$$ \|\nabla P_n^1 g\|_{L^s(I_n+B)} \lesssim \|P_n^0 g\|_{W^{1,s}(I_n+B)} + \|\nabla z_n(g)\|_{L^s((I_n+B)\setminus I_n)} \lesssim_s \|D(g)\|_{L^s(I_n)}, $$

that is, (3.12).

Step 2. Matrix potential for $P_n^1 g$.

We construct a matrix field $\sigma[P_n^1 g] \in C^1(\mathbb{R}^d)_{\text{skew}}^{d \times d}$ that decays at infinity such that

$$ \text{div}(\sigma[P_n^1 g]|_{I_n^+}) = P_n^1 g|_{I_n^+}, \quad (3.16) $$

and such that for all $\frac{d}{d-1} < s < \infty$ and $p > d$,

$$ \|\nabla\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|D(g)\|_{L^{\frac{ds}{d+s}}(I_n)}, \quad (3.17) $$

$$ \|\nabla\sigma[P_n^1 g]\|_{L^\infty(\mathbb{R}^d)} \lesssim_p \|D(g)\|_{L^p(I_n)}. $$

For that purpose, we extend $P_n^1 g$ by 0 to $\mathbb{R}^d$, viewing it as a compactly supported element of $H^1(\mathbb{R}^d)^d$, and for all $i, j$ we define $\nabla\sigma_{ij}[P_n^1 g] \in L^2(\mathbb{R}^d)^d$ as the unique solution of

$$ -\Delta\sigma_{ij}[P_n^1 g] = \partial_i(P_n^1 g)_j - \partial_j(P_n^1 g)_i \quad \text{in } \mathbb{R}^d. \quad (3.18) $$

In view of (3.12), Calderón–Zygmund potential theory yields $\nabla\sigma_{ij}[P_n^1 g] \in W^{1,s}(\mathbb{R}^d)^d$ for all $1 < s < \infty$. Moreover, as $P_n^1 g$ is compactly supported, Riesz potential theory ensures that $\sigma_{ij}[P_n^1 g]$ can itself be uniquely chosen as a decaying element in $C^1(\mathbb{R}^d)$. Uniqueness and the form of the right-hand side in (3.18) ensure that $\sigma[P_n^1 g]$ is skew-symmetric. Taking the divergence in (3.18), and using that $\text{div}(P_n^1 g) = 0$, we find

$$ -\Delta \text{div}(\sigma[P_n^1 g]) = -\Delta P_n^1 g \quad \text{in } \mathbb{R}^d, $$

which entails

$$ \text{div}(\sigma[P_n^1 g]) = P_n^1 g, $$

that is, (3.16). It remains to check (3.17). First, for all $\frac{d}{d-1} \le s < \infty$ and $p > d$, the Sobolev embedding gives

$$ \|\nabla\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|\nabla^2\sigma[P_n^1 g]\|_{L^{\frac{ds}{d+s}}(\mathbb{R}^d)}, $$

$$ \|\nabla\sigma[P_n^1 g]\|_{L^\infty(\mathbb{R}^d)} \lesssim_p \|\nabla\sigma[P_n^1 g]\|_{W^{1,p}(\mathbb{R}^d)}. $$

Second, for all $1 < s < \infty$, Calderón–Zygmund potential theory for (3.18) gives

$$ \|\nabla^2\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|\nabla P_n^1 g\|_{L^s(\mathbb{R}^d)}, \quad (3.19) $$

$$ \|\nabla\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|P_n^1 g\|_{L^s(\mathbb{R}^d)}. $$

Combining these two ingredients, appealing to Poincaré's inequality for $P_n^1 g$ supported in $I_n + B$, and using (3.12), claim (3.17) follows. For future reference, we note that a
which indeed satisfies $P_n^1 g|_{I_n} = P_n^0 g|_{I_n} = g - V_g - \Theta_g x$ and $\text{div}(P_n^1 g) = 0$, hence (3.11). In addition, combining (3.14) and (3.15) yields for all $1 < s < \infty$,

$$ \|\nabla P_n^1 g\|_{L^s(I_n+B)} \lesssim \|P_n^0 g\|_{W^{1,s}(I_n+B)} + \|\nabla z_n(g)\|_{L^s((I_n+B)\setminus I_n)} \lesssim_s \|D(g)\|_{L^s(I_n)}, $$

that is, (3.12).

**Step 2.** Matrix potential for $P_n^1 g$.

We construct a matrix field $\sigma[P_n^1 g] \in C^1(\mathbb{R}^d)_{\text{skew}}^{d \times d}$ that decays at infinity such that

$$ \text{div}(\sigma[P_n^1 g]|_{I_n^+}) = P_n^1 g|_{I_n^+}, \quad (3.16) $$

and such that for all $\frac{d}{d-1} < s < \infty$ and $p > d$,

$$ \|\nabla\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|D(g)\|_{L^{\frac{ds}{d+s}}(I_n)}, \quad (3.17) $$

$$ \|\nabla\sigma[P_n^1 g]\|_{L^\infty(\mathbb{R}^d)} \lesssim_p \|D(g)\|_{L^p(I_n)}. $$

For that purpose, we extend $P_n^1 g$ by 0 to $\mathbb{R}^d$, viewing it as a compactly supported element of $H^1(\mathbb{R}^d)^d$, and for all $i, j$ we define $\nabla\sigma_{ij}[P_n^1 g] \in L^2(\mathbb{R}^d)^d$ as the unique solution of

$$ -\Delta\sigma_{ij}[P_n^1 g] = \partial_i(P_n^1 g)_j - \partial_j(P_n^1 g)_i \quad \text{in } \mathbb{R}^d. \quad (3.18) $$

In view of (3.12), Calderón–Zygmund potential theory yields $\nabla\sigma_{ij}[P_n^1 g] \in W^{1,s}(\mathbb{R}^d)^d$ for all $1 < s < \infty$. Moreover, as $P_n^1 g$ is compactly supported, Riesz potential theory ensures that $\sigma_{ij}[P_n^1 g]$ can itself be uniquely chosen as a decaying element in $C^1(\mathbb{R}^d)$. Uniqueness and the form of the right-hand side in (3.18) ensure that $\sigma[P_n^1 g]$ is skew-symmetric. Taking the divergence in (3.18), and using that $\text{div}(P_n^1 g) = 0$, we find

$$ -\Delta \text{div}(\sigma[P_n^1 g]) = -\Delta P_n^1 g \quad \text{in } \mathbb{R}^d, $$

which entails

$$ \text{div}(\sigma[P_n^1 g]) = P_n^1 g, $$

that is, (3.16). It remains to check (3.17). First, for all $\frac{d}{d-1} \le s < \infty$ and $p > d$, the Sobolev embedding gives

$$ \|\nabla\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|\nabla^2\sigma[P_n^1 g]\|_{L^{\frac{ds}{d+s}}(\mathbb{R}^d)}, $$

$$ \|\nabla\sigma[P_n^1 g]\|_{L^\infty(\mathbb{R}^d)} \lesssim_p \|\nabla\sigma[P_n^1 g]\|_{W^{1,p}(\mathbb{R}^d)}. $$

Second, for all $1 < s < \infty$, Calderón–Zygmund potential theory for (3.18) gives

$$ \|\nabla^2\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|\nabla P_n^1 g\|_{L^s(\mathbb{R}^d)}, \quad (3.19) $$

$$ \|\nabla\sigma[P_n^1 g]\|_{L^s(\mathbb{R}^d)} \lesssim_s \|P_n^1 g\|_{L^s(\mathbb{R}^d)}. $$

Combining these two ingredients, appealing to Poincaré's inequality for $P_n^1 g$ supported in $I_n + B$, and using (3.12), claim (3.17) follows. For future reference, we note that a
similar argument also gives, for all $\frac{d}{d-1} < s < \infty$ and $p > d$,

$$ \begin{aligned} \| P_n^1 g \|_{L^s(\mathbb{R}^d)} &\lesssim_s \| \mathbf{D}(g) \|_{L^{\frac{ds}{d+s}}(I_n)}, \\ \| P_n^1 g \|_{L^\infty(\mathbb{R}^d)} &\lesssim_p \| \mathbf{D}(g) \|_{L^p(I_n)}. \end{aligned} \quad (3.20) $$

**Step 3. Conclusion.**

Recall the cutoff function $w_n \in H_0^1(I_n^+)$ that we constructed in Lemma 3.3, as well as the collection of neighborhoods of quasi-contact points $J_n^j = B(x_n^j, \frac{1}{C}\delta) \cap I_n^+$, $1 \le j \le M_n$. Recalling that $\text{dist}(J_n^j, J_n^k) \ge \frac{1}{C}\delta$ for $j \ne k$, we further define enlarged neighborhoods

$$ J_n^j \subset J_n^{j,+} := B(x_n^j, \frac{6}{5C}\delta) \cap I_n^+ \subset J_n^{j,++} := B(x_n^j, \frac{7}{5C}\delta) \cap I_n^+, $$

which then satisfy $\text{dist}(J_n^{j,++}, J_n^{k,++}) \ge \frac{1}{5C}\delta$ for $j \ne k$, and we abbreviate these as

$$ J_n := \bigcup_{j=1}^{M_n} J_n^j, \quad J_n^+ := \bigcup_{j=1}^{M_n} J_n^{j,+}, \quad J_n^{++} := \bigcup_{j=1}^{M_n} J_n^{j,++}. $$

We split the proof into two further substeps, first constructing the extension $P_n g$ close to quasi-contact points in $J_n^+$, and then completing the construction globally.

**Substep 3.1. Construction of $P_n g$ close to quasi-contact points.**

Given $1 < s \le r < \infty$ with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we construct a vector field $P_n^2 g \in H_0^1(I_n^+)^d$ such that

$$ P_n^2 g|_{I_n \cap J_n^+} = P_n^1 g|_{I_n \cap J_n^+}, \quad \text{and} \quad \text{div}(P_n^2 g) = 0 \text{ in } I_n^+, \quad (3.21) $$

and for all $p \ge s \vee \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \| \nabla P_n^2 g \|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \| \mathbf{D}(g) \|_{L^p(I_n)}, \quad (3.22) $$

where we recall the notation (2.21) for $\mu_r$. For all $j$ we first choose a smooth cutoff function $\chi_n^j \in C^\infty(I_n^+)$ such that

$$ \chi_n^j|_{J_n^{j,+}} = 1, \quad \chi_n^j|_{I_n^+ \setminus J_n^{j,++}} = 0, \quad \|\chi_n^j\|_{W^{2,\infty}(I_n^+)} \lesssim 1. $$

Given a collection of matrices $\{\Theta_n^j\}_{j=1}^{M_n} \subset \mathbb{M}^{\text{skew}}$ to be fixed later, we then define

$$ P_n^2 g := \sum_{j=1}^{M_n} \text{div}(w_n \chi_n^j (\sigma[P_n^1 g] - \Theta_n^j)). \quad (3.23) $$

By definition of $w_n$, this is supported in $I_n^+$ and satisfies, in view of (3.16),

$$ P_n^2 g|_{I_n \cap J_n^+} = \text{div}(\sigma[P_n^1 g])|_{I_n \cap J_n^+} = P_n^1 g|_{I_n \cap J_n^+}. $$
similar argument also gives, for all $\frac{d}{d-1} < s < \infty$ and $p > d$,

$$ \begin{aligned} \| P_n^1 g \|_{L^s(\mathbb{R}^d)} &\lesssim_s \| \mathbf{D}(g) \|_{L^{\frac{ds}{d+s}}(I_n)}, \\ \| P_n^1 g \|_{L^\infty(\mathbb{R}^d)} &\lesssim_p \| \mathbf{D}(g) \|_{L^p(I_n)}. \end{aligned} \quad (3.20) $$

### Step 3. Conclusion.

Recall the cutoff function $w_n \in H_0^1(I_n^+)$ that we constructed in Lemma 3.3, as well as the collection of neighborhoods of quasi-contact points $J_n^j = B(x_n^j, \frac{1}{C}\delta) \cap I_n^+, 1 \le j \le M_n$. Recalling that $\text{dist}(J_n^j, J_n^k) \ge \frac{1}{C}\delta$ for $j \ne k$, we further define enlarged neighborhoods

$$ J_n^j \subset J_n^{j,+} := B(x_n^j, \frac{6}{5C}\delta) \cap I_n^+ \subset J_n^{j,++} := B(x_n^j, \frac{7}{5C}\delta) \cap I_n^+, $$

which then satisfy $\text{dist}(J_n^{j,++}, J_n^{k,++}) \ge \frac{1}{5C}\delta$ for $j \ne k$, and we abbreviate these as

$$ J_n := \bigcup_{j=1}^{M_n} J_n^j, \quad J_n^+ := \bigcup_{j=1}^{M_n} J_n^{j,+}, \quad J_n^{++} := \bigcup_{j=1}^{M_n} J_n^{j,++}. $$

We split the proof into two further substeps, first constructing the extension $P_n g$ close to quasi-contact points in $J_n^+$, and then completing the construction globally.

#### Substep 3.1. Construction of $P_n g$ close to quasi-contact points.

Given $1 < s \le r < \infty$ with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we construct a vector field $P_n^2 g \in H_0^1(I_n^+)^d$ such that

$$ P_n^2 g|_{I_n \cap J_n^+} = P_n^1 g|_{I_n \cap J_n^+}, \quad \text{and} \quad \text{div}(P_n^2 g) = 0 \text{ in } I_n^+, \quad (3.21) $$

and for all $p \ge s \vee \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \| \nabla P_n^2 g \|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \| \mathbf{D}(g) \|_{L^p(I_n)}, \quad (3.22) $$

where we recall the notation (2.21) for $\mu_r$. For all $j$ we first choose a smooth cutoff function $\chi_n^j \in C^\infty(I_n^+)$ such that

$$ \chi_n^j|_{J_n^{j,+}} = 1, \quad \chi_n^j|_{I_n^+ \setminus J_n^{j,++}} = 0, \quad \|\chi_n^j\|_{W^{2,\infty}(I_n^+)} \lesssim 1. $$

Given a collection of matrices $\{\Theta_n^j\}_{j=1}^{M_n} \subset \mathbb{M}^{\text{skew}}$ to be fixed later, we then define

$$ P_n^2 g := \sum_{j=1}^{M_n} \text{div}(w_n \chi_n^j (\sigma[P_n^1 g] - \Theta_n^j)). \quad (3.23) $$

By definition of $w_n$, this is supported in $I_n^+$ and satisfies, in view of (3.16),

$$ P_n^2 g|_{I_n \cap J_n^+} = \text{div}(\sigma[P_n^1 g])|_{I_n \cap J_n^+} = P_n^1 g|_{I_n \cap J_n^+}. $$
Moreover, since $\sigma[P_n^1 g] - \Theta_n^j$ is skew-symmetric, we obviously have $\text{div}(P_n^2 g) = 0$. It remains to estimate the norm of $\nabla P_n^2 g$. To this aim, using (3.16) again, we compute

$$ \nabla P_n^2 g = \sum_{j=1}^{M_n} \nabla(w_n \chi_n^j P_n^1 g) + \sum_{j=1}^{M_n} \nabla((\sigma[P_n^1 g] - \Theta_n^j)\nabla(w_n \chi_n^j)). $$

Expanding the gradients, smuggling in the weights $x \mapsto |x - x_n^j|$, and using Hölder's inequality, we find for all $r \ge s \ge 1$,

$$ \begin{aligned} \|\nabla P_n^2 g\|_{L^s(I_n^+)} &\lesssim \|\nabla P_n^1 g\|_{L^s(I_n^+)} + \|\nabla(w_n \chi_n^j)\|_{L^r(I_n^+)} \|(P_n^1 g, \nabla\sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(I_n^+)} \\ &\quad + \sum_{j=1}^{M_n} \|\cdot - x_n^j|\nabla^2(w_n \chi_n^j)\|_{L^r(J_n^{j,++})} \\ &\quad \times \|\cdot - x_n^j|^{-1}(\sigma[P_n^1 g] - \Theta_n^j)\|_{L^{\frac{rs}{r-s}}(J_n^{j,++})}, \end{aligned} \quad (3.24) $$

and thus, inserting the estimates of Lemma 3.3 for norms of the cutoff function $w_n$, and recalling the definition (2.21) of $\mu_r$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_r \mu_r(\rho_n) \left( \|\nabla P_n^1 g\|_{L^s(I_n^+)} + \|(P_n^1 g, \nabla\sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(I_n^+)} \right. \\ \left. \qquad + \sum_{j=1}^{M_n} \|\cdot - x_n^j|^{-1}(\sigma[P_n^1 g] - \Theta_n^j)\|_{L^{\frac{rs}{r-s}}(J_n^{j,++})} \right). \quad (3.25) $$

We estimate the right-hand side in two different ways, corresponding to two different choices of the constants $\{\Theta_n^j\}_{j=1}^{M_n}$ and allowing for complementary ranges of exponents.

* Case 1. Choosing $\Theta_n^j = 0$ for all $j$, given $1 < s < d$ and $r > \frac{ds}{d-s}$, with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we obtain for all $p \ge s \vee \frac{dr s}{d(r-s)+rs}$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{r,s} \mu_r(\rho_n) \|D(g)\|_{L^p(I_n)}. \quad (3.26) $$

For that purpose, we appeal to Hardy's inequality in the following form (see e.g. [23, Sections 1.3 and 12.8]): for all $x_0 \in \mathbb{R}^d$ and $1 \le p < d$,

$$ \|\cdot - x_0|^{-1}\sigma[P_n^1 g]\|_{L^p(\mathbb{R}^d)} \lesssim_p \|\nabla\sigma[P_n^1 g]\|_{L^p(\mathbb{R}^d)}. $$

Choosing $\Theta_n^j = 0$, and inserting this estimate into (3.25), we find for all $r \ge s \ge 1$ with $\frac{rs}{r-s} < d$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{r,s} \mu_r(\rho_n) \left( \|\nabla P_n^1 g\|_{L^s(\mathbb{R}^d)} + \|(P_n^1 g, \nabla\sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(\mathbb{R}^d)} \right), $$

and claim (3.26) follows from (3.12), (3.17), and (3.20).

* Case 2. Choosing $\Theta_n^j = \sigma[P_n^1 g](x_n^j)$ for all $j$, given $1 < s \le r < \infty$, with $r < \frac{ds}{d-s}$ if $s < d$, we obtain for all $p \ge s \vee \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \|D(g)\|_{L^p(I_n)}. \quad (3.27) $$
Moreover, since $\sigma[P_n^1 g] - \Theta_n^j$ is skew-symmetric, we obviously have $\text{div}(P_n^2 g) = 0$. It remains to estimate the norm of $\nabla P_n^2 g$. To this aim, using (3.16) again, we compute

$$ \nabla P_n^2 g = \sum_{j=1}^{M_n} \nabla(w_n \chi_n^j P_n^1 g) + \sum_{j=1}^{M_n} \nabla((\sigma[P_n^1 g] - \Theta_n^j)\nabla(w_n \chi_n^j)). $$

Expanding the gradients, smuggling in the weights $x \mapsto |x - x_n^j|$, and using Hölder's inequality, we find for all $r \ge s \ge 1$,

$$ \begin{aligned} \|\nabla P_n^2 g\|_{L^s(I_n^+)} &\lesssim \|\nabla P_n^1 g\|_{L^s(I_n^+)} + \|\nabla(w_n \chi_n^j)\|_{L^r(I_n^+)} \|(P_n^1 g, \nabla\sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(I_n^+)} \\ &\quad + \sum_{j=1}^{M_n} \|\cdot - x_n^j|\nabla^2(w_n \chi_n^j)\|_{L^r(J_n^{j,++})} \\ &\quad \times \|\cdot - x_n^j|^{-1}(\sigma[P_n^1 g] - \Theta_n^j)\|_{L^{\frac{rs}{r-s}}(J_n^{j,++})}, \end{aligned} \quad (3.24) $$

and thus, inserting the estimates of Lemma 3.3 for norms of the cutoff function $w_n$, and recalling the definition (2.21) of $\mu_r$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_r \mu_r(\rho_n) \left( \|\nabla P_n^1 g\|_{L^s(I_n^+)} + \|(P_n^1 g, \nabla\sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(I_n^+)} \right. \\ \left. \qquad + \sum_{j=1}^{M_n} \|\cdot - x_n^j|^{-1}(\sigma[P_n^1 g] - \Theta_n^j)\|_{L^{\frac{rs}{r-s}}(J_n^{j,++})} \right). \quad (3.25) $$

We estimate the right-hand side in two different ways, corresponding to two different choices of the constants $\{\Theta_n^j\}_{j=1}^{M_n}$ and allowing for complementary ranges of exponents.

* Case 1. Choosing $\Theta_n^j = 0$ for all $j$, given $1 < s < d$ and $r > \frac{ds}{d-s}$, with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we obtain for all $p \ge s \vee \frac{dr s}{d(r-s)+rs}$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{r,s} \mu_r(\rho_n) \|D(g)\|_{L^p(I_n)}. \quad (3.26) $$

For that purpose, we appeal to Hardy's inequality in the following form (see e.g. [23, Sections 1.3 and 12.8]): for all $x_0 \in \mathbb{R}^d$ and $1 \le p < d$,

$$ \|\cdot - x_0|^{-1}\sigma[P_n^1 g]\|_{L^p(\mathbb{R}^d)} \lesssim_p \|\nabla\sigma[P_n^1 g]\|_{L^p(\mathbb{R}^d)}. $$

Choosing $\Theta_n^j = 0$, and inserting this estimate into (3.25), we find for all $r \ge s \ge 1$ with $\frac{rs}{r-s} < d$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{r,s} \mu_r(\rho_n) \left( \|\nabla P_n^1 g\|_{L^s(\mathbb{R}^d)} + \|(P_n^1 g, \nabla\sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(\mathbb{R}^d)} \right), $$

and claim (3.26) follows from (3.12), (3.17), and (3.20).

* Case 2. Choosing $\Theta_n^j = \sigma[P_n^1 g](x_n^j)$ for all $j$, given $1 < s \le r < \infty$, with $r < \frac{ds}{d-s}$ if $s < d$, we obtain for all $p \ge s \vee \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \|D(g)\|_{L^p(I_n)}. \quad (3.27) $$
For that purpose, we appeal to Hardy's inequality in the following form (see e.g. [23, Sections 1.3 and 12.8]): for all $x_0 \in I_n^+$ and $d < p \le \infty$,

$$ \| |\cdot - x_0|^{-1} (\sigma[P_n^1 g] - \sigma[P_n^1 g](x_0)) \|_{L^p(I_n^+)} \lesssim_p \|\nabla \sigma[P_n^1 g]\|_{L^p(\mathbb{R}^d)}. $$

Choosing $\Theta_n^j = \sigma[P_n^1 g](x_n^j)$, and inserting this estimate into (3.25), we find for all $r \ge s \ge 1$ with $\frac{rs}{r-s} > d$,

$$ \|\nabla P_n^2 g\|_{L^s(I_n^+)} \lesssim_{r,s} \mu_r(\rho_n)(\|\nabla P_n^1 g\|_{L^s(\mathbb{R}^d)} + \|(P_n^1 g, \nabla \sigma[P_n^1 g])\|_{L^{\frac{rs}{r-s}}(\mathbb{R}^d)}), $$

and claim (3.27) follows from (3.12), (3.17), and (3.20).

Combining (3.26) and (3.27), and choosing the constants $\{\Theta_n^j\}_{j=1}^{M_n}$ accordingly in the definition (3.23) of $P_n^2 g$, claim (3.22) follows.

Substep 3.2. Construction of $P_n g$ away from contact points.
Let $1 < s \le r < \infty$ be fixed, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$.
Choosing a cutoff function $\zeta \in C_c^\infty(I_n^+ \setminus J_n)$ with $\zeta|_{I_n \setminus J_n^+} = 1$ and $\|\zeta\|_{W^{1,\infty}(I_n^+ \setminus J_n)} \lesssim 1$, we consider the vector field

$$ Q_n g := \zeta(P_n^1 g - P_n^2 g), $$

and we note that in view of (3.21) it satisfies

$$ Q_n g|_{I_n} = (P_n^1 g - P_n^2 g)|_{I_n}, \qquad (3.28) $$

hence in particular $\text{div}(Q_n g)|_{I_n} = 0$. As this yields the relation

$$ \int_{I_n^+ \setminus (I_n \cup J_n)} \text{div}(Q_n g) = - \int_{\partial(I_n \setminus J_n)} (Q_n g) \cdot \nu = - \int_{I_n \setminus J_n} \text{div}(Q_n g) = 0, $$

we can appeal to the same construction based on the Bogovskii operator as in Step 1: there exists $t_n(g) \in H_0^1(I_n^+ \setminus (I_n \cup J_n))^d$ such that

$$ \text{div}(t_n(g)) = \text{div}(Q_n g) \quad \text{in } I_n^+ \setminus (I_n \cup J_n), $$

and

$$ \|\nabla t_n(g)\|_{L^s(I_n^+ \setminus (I_n \cup J_n))} \lesssim_s \|\text{div}(Q_n g)\|_{L^s(I_n^+ \setminus J_n)}. \qquad (3.29) $$

Here comes the restriction to $d > 2$ as the set $I_n^+ \setminus (I_n \cup J_n)$ is typically not connected in dimension $d = 2$; see Remark 3.4 below. In these terms, we finally define

$$ P_n g := P_n^2 g + Q_n g - t_n(g) \quad \in H_0^1(I_n^+), $$

which satisfies, in view of (3.28),

$$ P_n g|_{I_n} = P_n^2 g|_{I_n} + (P_n^1 g - P_n^2 g)|_{I_n} = P_n^1 g|_{I_n}, $$
For that purpose, we appeal to Hardy's inequality in the following form (see e.g. [23, Sections 1.3 and 12.8]): for all $x_0 \in I_n^+$ and $d < p \le \infty$,

$$ \| | \cdot - x_0|^{-1} (\sigma[P_n^1 g] - \sigma[P_n^1 g](x_0)) \|_{L^p(I_n^+)} \lesssim_p \| \nabla \sigma [P_n^1 g] \|_{L^p(\mathbb{R}^d)}. $$

Choosing $\Theta_n^j = \sigma[P_n^1 g](x_n^j)$, and inserting this estimate into (3.25), we find for all $r \ge s \ge 1$ with $\frac{rs}{r-s} > d$,

$$ \| \nabla P_n^2 g \|_{L^s(I_n^+)} \lesssim_{r,s} \mu_r(\rho_n) (\| \nabla P_n^1 g \|_{L^s(\mathbb{R}^d)} + \|(P_n^1 g, \nabla \sigma [P_n^1 g])\|_{L^{\frac{rs}{r-s}}(\mathbb{R}^d)}), $$

and claim (3.27) follows from (3.12), (3.17), and (3.20).

Combining (3.26) and (3.27), and choosing the constants $\{\Theta_n^j\}_{j=1}^{M_n}$ accordingly in the definition (3.23) of $P_n^2 g$, claim (3.22) follows.

Substep 3.2. Construction of $P_n g$ away from contact points.
Let $1 < s \le r < \infty$ be fixed, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$.
Choosing a cutoff function $\zeta \in C_c^\infty(I_n^+ \setminus J_n)$ with $\zeta|_{I_n \setminus J_n^+} = 1$ and $\|\zeta\|_{W^{1,\infty}(I_n^+ \setminus J_n)} \lesssim 1$, we consider the vector field

$$ Q_n g := \zeta(P_n^1 g - P_n^2 g), $$

and we note that in view of (3.21) it satisfies

$$ Q_n g|_{I_n} = (P_n^1 g - P_n^2 g)|_{I_n}, \quad (3.28) $$

hence in particular $\mathrm{div}(Q_n g)|_{I_n} = 0$. As this yields the relation

$$ \int_{I_n^+ \setminus (I_n \cup J_n)} \mathrm{div}(Q_n g) = - \int_{\partial(I_n \setminus J_n)} (Q_n g) \cdot \nu = - \int_{I_n \setminus J_n} \mathrm{div}(Q_n g) = 0, $$

we can appeal to the same construction based on the Bogovskii operator as in Step 1: there exists $t_n(g) \in H_0^1(I_n^+ \setminus (I_n \cup J_n))^d$ such that

$$ \mathrm{div}(t_n(g)) = \mathrm{div}(Q_n g) \quad \text{in } I_n^+ \setminus (I_n \cup J_n), $$

and

$$ \| \nabla t_n(g) \|_{L^s(I_n^+ \setminus (I_n \cup J_n))} \lesssim_s \| \mathrm{div}(Q_n g) \|_{L^s(I_n^+ \setminus J_n)}. \quad (3.29) $$

Here comes the restriction to $d > 2$ as the set $I_n^+ \setminus (I_n \cup J_n)$ is typically not connected in dimension $d = 2$; see Remark 3.4 below. In these terms, we finally define

$$ P_n g := P_n^2 g + Q_n g - t_n(g) \quad \in H_0^1(I_n^+), $$

which satisfies, in view of (3.28),

$$ P_n g|_{I_n} = P_n^2 g|_{I_n} + (P_n^1 g - P_n^2 g)|_{I_n} = P_n^1 g|_{I_n}, $$
and also $\mathrm{div}(P_n g) = 0$ by definition of $t_n(g)$. In addition, combining (3.29) with the definition of $Q_n g$, we find

$$ \|\nabla P_n g\|_{L^s(I_n^+)} \lesssim_s \|\nabla P_n^2 g\|_{L^s(I_n^+)} + \|\nabla Q_n g\|_{L^s(I_n^+)} \lesssim \|(P_n^1 g, P_n^2 g)\|_{W^{1,s}(I_n^+)} $$

hence, using Poincaré's inequality and inserting (3.12) and (3.22), for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla P_n g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \|D(g)\|_{L^p(I_n)}. $$

This concludes the proof.

<figure><img src="image_4.png" /></figure>

**Remark 3.4** (Two-dimensional case). The restriction to $d > 2$ is due to the impossibility of fixing a stream function $\sigma[P_n^1 g]$ that would vanish at all quasi-contact points. More precisely, in Case 2 of the above proof, we adapt the stream function $\sigma[P_n^1 g]$ locally by making it vanish at each quasi-contact point (cf. choice of $\Theta_n^j$ in (3.23)), and modifications are then glued together in $I_n^+ \setminus J_n$ while the field must remain divergence-free and keep the same symmetric gradient in $I_n$. In two dimensions this is not possible since $I_n^+ \setminus (I_n \cup J_n)$ is not connected whenever $I_n$ has multiple quasi-contact points. Due to this geometric rigidity in two dimensions, the above proof is no longer valid: we must abandon the cancellation of the stream function at quasi-contact points and rather consider the extension operator

$$ \tilde{P}g := \mathrm{div}(w_n \sigma[P_n^1 g]). $$

The bound (3.24) then becomes, for all $r \ge s \ge 1$, with $r < \frac{2s}{2-s}$ if $s < 2$,

$$ \begin{aligned} \|\nabla \tilde{P} g\|_{L^s(I_n^+)} &\lesssim \|\nabla P_n^1 g\|_{L^s(I_n^+)} \\ &\quad + \|\nabla w_n\|_{W^{1,r}(I_n^+)} (\|P_n^1 g\|_{L^{\frac{rs}{r-s}}(\mathbb{R}^2)} + \|\sigma[P_n^1 g]\|_{W^{1,\frac{rs}{r-s}}(\mathbb{R}^2)}) \\ &\lesssim_{r,s} \|\nabla P_n^1 g\|_{L^s(I_n^+)} + \|\nabla w_n\|_{W^{1,r}(I_n^+)} \|P_n^1 g\|_{L^{\frac{rs}{r-s}}(I_n+B)}, \end{aligned} $$

where we used the Sobolev embedding, the bound (3.19) on $\sigma[P_n^1 g]$, and Jensen's inequality. Combining this with (3.5), (3.12), and (3.20), we deduce for all $r \ge s \ge 1$, with $r < \frac{2s}{2-s}$ if $s < 2$, and for all $p \ge s \lor \frac{2rs}{2(r-s)+rs}$, with $p > 2$ if $r = s$,

$$ \|\nabla \tilde{P} g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \rho_n^{\frac{3}{2r}-2} \|D(g)\|_{L^p(I_n)}. $$

Replacing Proposition 3.1 by this extension result would lead to corresponding two-dimensional versions of our main results; we skip the detail for brevity.

### 3.3. Proof of Lemma 3.2

The starting point is the following standard construction based on the Bogovskii operator, e.g. [13, Theorem III.3.1]: given a domain $D$ as in the statement, and given $h \in C_b(D)$ with $\int_{D \setminus I} h = 0$, there exists $z^0 \in H_0^1(D)^d$ such that

$$ \mathrm{div}(z^0) = h \mathbf{1}_{D \setminus I} \quad \text{in } D, $$
and also $\mathrm{div}(P_n g) = 0$ by definition of $t_n(g)$. In addition, combining (3.29) with the definition of $Q_n g$, we find

$$ \|\nabla P_n g\|_{L^s(I_n^+)} \lesssim_s \|\nabla P_n^2 g\|_{L^s(I_n^+)} + \|\nabla Q_n g\|_{L^s(I_n^+)} \lesssim \|(P_n^1 g, P_n^2 g)\|_{W^{1,s}(I_n^+)} $$

hence, using Poincaré's inequality and inserting (3.12) and (3.22), for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \|\nabla P_n g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \mu_r(\rho_n) \|D(g)\|_{L^p(I_n)} $$

This concludes the proof.

<figure><img src="image_6.png" /></figure>

**Remark 3.4** (Two-dimensional case). The restriction to $d > 2$ is due to the impossibility of fixing a stream function $\sigma[P_n^1 g]$ that would vanish at all quasi-contact points. More precisely, in Case 2 of the above proof, we adapt the stream function $\sigma[P_n^1 g]$ locally by making it vanish at each quasi-contact point (cf. choice of $\Theta_n^j$ in (3.23)), and modifications are then glued together in $I_n^+ \setminus J_n$ while the field must remain divergence-free and keep the same symmetric gradient in $I_n$. In two dimensions this is not possible since $I_n^+ \setminus (I_n \cup J_n)$ is not connected whenever $I_n$ has multiple quasi-contact points. Due to this geometric rigidity in two dimensions, the above proof is no longer valid: we must abandon the cancellation of the stream function at quasi-contact points and rather consider the extension operator

$$ \tilde{P}g := \mathrm{div}(w_n \sigma[P_n^1 g]). $$

The bound (3.24) then becomes, for all $r \ge s \ge 1$, with $r < \frac{2s}{2-s}$ if $s < 2$,

$$ \begin{aligned} \|\nabla \tilde{P} g\|_{L^s(I_n^+)} &\lesssim \|\nabla P_n^1 g\|_{L^s(I_n^+)} \\ &\quad + \|\nabla w_n\|_{W^{1,r}(I_n^+)} (\|P_n^1 g\|_{L^{\frac{rs}{r-s}}(\mathbb{R}^2)} + \|\sigma[P_n^1 g]\|_{W^{1,\frac{rs}{r-s}}(\mathbb{R}^2)}) \\ &\lesssim_{r,s} \|\nabla P_n^1 g\|_{L^s(I_n^+)} + \|\nabla w_n\|_{W^{1,r}(I_n^+)} \|P_n^1 g\|_{L^{\frac{rs}{r-s}}(I_n+B)}, \end{aligned} $$

where we used the Sobolev embedding, the bound (3.19) on $\sigma[P_n^1 g]$, and Jensen's inequality. Combining this with (3.5), (3.12), and (3.20), we deduce for all $r \ge s \ge 1$, with $r < \frac{2s}{2-s}$ if $s < 2$, and for all $p \ge s \lor \frac{2rs}{2(r-s)+rs}$, with $p > 2$ if $r = s$,

$$ \|\nabla \tilde{P} g\|_{L^s(I_n^+)} \lesssim_{p,r,s} \rho_n^{\frac{3}{2r}-2} \|D(g)\|_{L^p(I_n)}. $$

Replacing Proposition 3.1 by this extension result would lead to corresponding two-dimensional versions of our main results; we skip the detail for brevity.

### 3.3. Proof of Lemma 3.2

The starting point is the following standard construction based on the Bogovskii operator, e.g. [13, Theorem III.3.1]: given a domain $D$ as in the statement, and given $h \in C_b(D)$ with $\int_{D \setminus I} h = 0$, there exists $z^0 \in H_0^1(D)^d$ such that

$$ \mathrm{div}(z^0) = h \mathbf{1}_{D \setminus I} \quad \text{in } D, $$
and for all $1 < s < \infty$,

$$ \|\nabla z^0\|_{L^s(D)} \lesssim_s K(D)^{d+1} \|h\|_{L^s(D\setminus I)}. \quad (3.30) $$

Next, given $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we appeal to the extension operator $P_n$ that we constructed in Proposition 3.1, and we define

$$ z := z^0 - \sum_n P_n(z^0|_{I_n}). $$

By the properties of $P_n$, we find

$$ D(z)|_I = 0, \quad \text{and} \quad \text{div}(z) = \text{div}(z^0) = h \mathbb{1}_{D \setminus h} \quad \text{in } D, $$

and for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \begin{aligned} \|\nabla z\|_{L^s(D)}^s &\lesssim_s \|\nabla z^0\|_{L^s(D)}^s + \sum_{n: I_n \cap D \ne \emptyset} \|\nabla P_n(z^0|_{I_n})\|_{L^s(I_n^+)}^s \\ &\lesssim_{p,r,s} \|\nabla z^0\|_{L^s(D)}^s + \sum_{n: I_n \cap D \ne \emptyset} \mu_r(\rho_n)^s \|D(z^0)\|_{L^p(I_n)}^s \\ &\lesssim \left( |D| + \sum_{n: I_n \cap D \ne \emptyset} \mu_r(\rho_n)^{\frac{ps}{p-s}} \right)^{1-\frac{s}{p}} \|\nabla z^0\|_{L^p(D)}^s, \end{aligned} $$

where the last bound follows from Hölder's inequality. Combined with (3.30), this yields the conclusion. ■

### 3.4. Proof of Theorem 4

We split the proof into three steps.

*Step 1.* Given $q, S, f$ as in (2.16), and given $1 < \beta < \infty$ and $\alpha, r$ as in (2.17), we show that for all $n$ there exists $z_n \in W^{1,\alpha}(I_n)^d$ such that

$$ 2 \int_{I_n} D(g) : D(z_n) = \int_{I_n^+} g \cdot f - 2 \int_{I_n^+ \setminus I_n} D(g) : q \quad \forall g \in C_c^1(I_n^+)^d : \operatorname{div}(g) = 0 \quad (3.31) $$

and

$$ \|D(z_n)\|_{L^\alpha(I_n)} \lesssim_{\alpha, \beta, r} \mu_r(\rho_n)(\|f\|_{W^{-1,\beta}(I_n^+)} + \|q\|_{L^\beta(I_n^+ \setminus I_n)}). \quad (3.32) $$

While the left-hand side in (3.31) only involves the restriction $g|_{I_n} \in C_b^1(I_n)^d$ of the test function $g$, the right-hand side involves its extension $g \in C_c^1(I_n^+)^d$. In view of condition (2.16), the choice of the extension does not matter. Given $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we recall the extension operator $P_n$ that we constructed in Proposition 3.1, and problem (3.31) then reads

$$ 2 \int_{I_n} D(g) : D(z_n) = \mathcal{F}_n(g) \quad \forall g \in C_b^1(I_n)^d : \operatorname{div}(g) = 0, \quad (3.33) $$
and for all $1 < s < \infty$,

$$ \|\nabla z^0\|_{L^s(D)} \lesssim_s K(D)^{d+1} \|h\|_{L^s(D\setminus I)}. \quad (3.30) $$

Next, given $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we appeal to the extension operator $P_n$ that we constructed in Proposition 3.1, and we define

$$ z := z^0 - \sum_n P_n(z^0|_{I_n}). $$

By the properties of $P_n$, we find

$$ D(z)|_I = 0, \quad \text{and} \quad \operatorname{div}(z) = \operatorname{div}(z^0) = h \mathbb{1}_{D \setminus h} \quad \text{in } D, $$

and for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \begin{aligned} \|\nabla z\|_{L^s(D)}^s &\lesssim_s \|\nabla z^0\|_{L^s(D)}^s + \sum_{n: I_n \cap D \ne \emptyset} \|\nabla P_n(z^0|_{I_n})\|_{L^s(I_n^+)}^s \\ &\lesssim_{p,r,s} \|\nabla z^0\|_{L^s(D)}^s + \sum_{n: I_n \cap D \ne \emptyset} \mu_r(\rho_n)^s \|D(z^0)\|_{L^p(I_n)}^s \\ &\lesssim \left( |D| + \sum_{n: I_n \cap D \ne \emptyset} \mu_r(\rho_n)^{\frac{ps}{p-s}} \right)^{1-\frac{s}{p}} \|\nabla z^0\|_{L^p(D)}^s, \end{aligned} $$

where the last bound follows from Hölder's inequality. Combined with (3.30), this yields the conclusion. ■

### 3.4. Proof of Theorem 4

We split the proof into three steps.

**Step 1.** Given $q, S, f$ as in (2.16), and given $1 < \beta < \infty$ and $\alpha, r$ as in (2.17), we show that for all $n$ there exists $z_n \in W^{1,\alpha}(I_n)^d$ such that

$$ 2 \int_{I_n} D(g) : D(z_n) = \int_{I_n^+} g \cdot f - 2 \int_{I_n^+ \setminus I_n} D(g) : q \quad \forall g \in C_c^1(I_n^+)^d : \operatorname{div}(g) = 0 \quad (3.31) $$

and

$$ \|D(z_n)\|_{L^\alpha(I_n)} \lesssim_{\alpha, \beta, r} \mu_r(\rho_n)(\|f\|_{W^{-1,\beta}(I_n^+)} + \|q\|_{L^\beta(I_n^+ \setminus I_n)}). \quad (3.32) $$

While the left-hand side in (3.31) only involves the restriction $g|_{I_n} \in C_b^1(I_n)^d$ of the test function $g$, the right-hand side involves its extension $g \in C_c^1(I_n^+)^d$. In view of condition (2.16), the choice of the extension does not matter. Given $1 < s \le r < \infty$, with $r \ne \frac{ds}{d-s}$ if $s < d$, and with $r < \frac{ds}{d+s-ds}$ if $s < \frac{d}{d-1}$, we recall the extension operator $P_n$ that we constructed in Proposition 3.1, and problem (3.31) then reads

$$ 2 \int_{I_n} D(g) : D(z_n) = \mathcal{F}_n(g) \quad \forall g \in C_b^1(I_n)^d : \operatorname{div}(g) = 0, \quad (3.33) $$
where we have abbreviated

$$ \mathcal{F}_n(g) := \int_{I_n^+} (P_n g) \cdot f - 2 \int_{I_n^+ \setminus I_n} \mathrm{D}(P_n g) : q. $$

By Proposition 3.1, we find for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \begin{aligned} |\mathcal{F}_n(g)| &\lesssim (\|f\|_{W^{-1,s'}(I_n^+)} + \|q\|_{L^{s'}(I_n^+ \setminus I_n)}) \|\nabla P_n g\|_{L^s(I_n^+)} \\ &\lesssim_{p,r,s} \mu_r(\rho_n)(\|f\|_{W^{-1,s'}(I_n^+)} + \|q\|_{L^{s'}(I_n^+ \setminus I_n)}) \|\mathrm{D}(g)\|_{L^p(I_n)}. \end{aligned} $$

Appealing to the $L^{p'}$ theory for the Stokes equation, e.g. [13, Section IV.6], we deduce that there exists a solution $z_n \in W^{1,p'}(I_n)^d$ of problem (3.33) (unique up to a rigid motion), and that it satisfies

$$ \|\mathrm{D}(z_n)\|_{L^{p'}(I_n)} \lesssim_{p,r,s} \mu_r(\rho_n)(\|f\|_{W^{-1,s'}(I_n^+)} + \|q\|_{L^{s'}(I_n^+ \setminus I_n)}). $$

Setting $\alpha := p'$ and $\beta := s'$, this yields claim (3.31)–(3.32).

Step 2. Construction of extended flux.

Given $1 < \beta < \infty$ and $\alpha, r$ as in (2.17), define

$$ \tilde{q} := q \mathbb{1}_{\mathbb{R}^d \setminus I} + \sum_n \mathrm{D}(z_n) \mathbb{1}_{I_n}, $$

with $z_n$ as constructed in Step 1. Given $g \in C_c^1(\mathbb{R}^d)^d$ with $\mathrm{div}(g) = 0$, we may decompose

$$ g = g_{\circ} + \sum_{n} g_{n}, \quad g_{\circ} := g - \sum_{n} P_{n}g, \quad g_{n} := P_{n}g. $$

Using (2.16) with test function $g_\circ$, and using (3.33) with test function $g_n$, we are led to the following integral identity:

$$ 2 \int_{\mathbb{R}^d} \mathrm{D}(g) : \tilde{q} = \int_{\mathbb{R}^d} g \cdot f \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \mathrm{div}(g) = 0. \quad (3.34) $$

Next we prove bound (2.19) for $\tilde{q}$. Given a bounded domain $D \subset \mathbb{R}^d$, summing (3.32) over all particles, appealing to Hölder's inequality, and using the Sobolev embedding $L^{d\beta/(d+\beta)} \hookrightarrow W^{-1,\beta}$, we find

$$ \|\tilde{q}\|_{L^\alpha(D)}^\alpha \le \|q\|_{L^\alpha(D \setminus I)}^\alpha + \sum_{n: I_n \cap D \ne \emptyset} \|\mathrm{D}(z_n)\|_{L^\alpha(I_n)}^\alpha \quad (3.35) $$

$$ \lesssim_{\alpha, \beta, r} \|q\|_{L^\alpha(D \setminus I)}^\alpha + \sum_{n: I_n \cap D \neq \emptyset} \mu_r(\rho_n)^\alpha (\|f\|_{W^{-1,\beta}(I_n^+)}^\beta + \|q\|_{L^\beta(I_n^+ \setminus I_n)}^\beta)^{\frac{\alpha}{\beta}} $$

$$ \lesssim \left( |D| + \sum_{n: I_n \cap D \neq \emptyset} \mu_r(\rho_n)^{\frac{\beta\alpha}{\beta-\alpha}} \right)^{1-\frac{\alpha}{\beta}} (\|f\|_{L^{\frac{d\beta}{d+\beta}}(\hat{D})}^\alpha + \|q\|_{L^\beta(\hat{D} \setminus I)}^\alpha) $$

where we recall the notation $\hat{D} = D \cup \bigcup_{n: I_n \cap D \neq \emptyset} I_n^+$.
where we have abbreviated

$$ \mathcal{F}_n(g) := \int_{I_n^+} (P_n g) \cdot f - 2 \int_{I_n^+ \setminus I_n} \mathrm{D}(P_n g) : q. $$

By Proposition 3.1, we find for all $p \ge s \lor \frac{dr s}{d(r-s)+rs}$, with $p > d$ if $r = s$,

$$ \begin{aligned} |\mathcal{F}_n(g)| &\lesssim (\|f\|_{W^{-1,s'}(I_n^+)} + \|q\|_{L^{s'}(I_n^+ \setminus I_n)}) \|\nabla P_n g\|_{L^s(I_n^+)} \\ &\lesssim_{p,r,s} \mu_r(\rho_n)(\|f\|_{W^{-1,s'}(I_n^+)} + \|q\|_{L^{s'}(I_n^+ \setminus I_n)}) \|\mathrm{D}(g)\|_{L^p(I_n)}. \end{aligned} $$

Appealing to the $L^{p'}$ theory for the Stokes equation, e.g. [13, Section IV.6], we deduce that there exists a solution $z_n \in W^{1,p'}(I_n)^d$ of problem (3.33) (unique up to a rigid motion), and that it satisfies

$$ \|\mathrm{D}(z_n)\|_{L^{p'}(I_n)} \lesssim_{p,r,s} \mu_r(\rho_n)(\|f\|_{W^{-1,s'}(I_n^+)} + \|q\|_{L^{s'}(I_n^+ \setminus I_n)}). $$

Setting $\alpha := p'$ and $\beta := s'$, this yields claim (3.31)–(3.32).

Step 2. Construction of extended flux.

Given $1 < \beta < \infty$ and $\alpha, r$ as in (2.17), define

$$ \tilde{q} := q \mathbb{1}_{\mathbb{R}^d \setminus I} + \sum_n \mathrm{D}(z_n) \mathbb{1}_{I_n}, $$

with $z_n$ as constructed in Step 1. Given $g \in C_c^1(\mathbb{R}^d)^d$ with $\mathrm{div}(g) = 0$, we may decompose

$$ g = g_{\circ} + \sum_{n} g_n, \quad g_{\circ} := g - \sum_{n} P_n g, \quad g_n := P_n g. $$

Using (2.16) with test function $g_{\circ}$, and using (3.33) with test function $g_n$, we are led to the following integral identity:

$$ 2 \int_{\mathbb{R}^d} \mathrm{D}(g) : \tilde{q} = \int_{\mathbb{R}^d} g \cdot f \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \mathrm{div}(g) = 0. \quad (3.34) $$

Next we prove bound (2.19) for $\tilde{q}$. Given a bounded domain $D \subset \mathbb{R}^d$, summing (3.32) over all particles, appealing to Hölder's inequality, and using the Sobolev embedding $L^{d\beta/(d+\beta)} \hookrightarrow W^{-1,\beta}$, we find

$$ \|\tilde{q}\|_{L^\alpha(D)}^\alpha \le \|q\|_{L^\alpha(D \setminus I)}^\alpha + \sum_{n: I_n \cap D \ne \emptyset} \|\mathrm{D}(z_n)\|_{L^\alpha(I_n)}^\alpha \quad (3.35) $$

$$ \lesssim_{\alpha, \beta, r} \|q\|_{L^\alpha(D \setminus I)}^\alpha + \sum_{n: I_n \cap D \neq \emptyset} \mu_r(\rho_n)^\alpha (\|f\|_{W^{-1,\beta}(I_n^+)}^\beta + \|q\|_{L^\beta(I_n^+ \setminus I_n)}^\beta)^{\frac{\alpha}{\beta}} $$

$$ \lesssim \left( |D| + \sum_{n: I_n \cap D \neq \emptyset} \mu_r(\rho_n)^{\frac{\beta\alpha}{\beta-\alpha}} \right)^{1-\frac{\alpha}{\beta}} (\|f\|_{L^{\frac{d\beta}{d+\beta}}(\hat{D})}^\alpha + \|q\|_{L^\beta(\hat{D} \setminus I)}^\alpha) $$

where we recall the notation $\hat{D} = D \cup \bigcup_{n: I_n \cap D \neq \emptyset} I_n^+$.
Step 3. Construction of extended pressure.

In view of e.g. [20, Proposition 12.10], the relation (3.34) for the extension $\tilde{q}$ ensures the existence of an associated pressure field $\tilde{S} \in L_{\text{loc}}^1(\mathbb{R}^d)$, uniquely defined up to a global additive constant, such that

$$ \int_{\mathbb{R}^d} D(g) : (2\tilde{q} - \tilde{S} \text{ Id}) = \int_{\mathbb{R}^d} g \cdot f \quad \forall g \in C_c^1(\mathbb{R}^d)^d, \qquad (3.36) $$

that is, $-\text{div}(2\tilde{q} - \tilde{S} \text{ Id}) = f$ in $\mathbb{R}^d$. It remains to prove bound (2.19) for $\tilde{S}$. For all $R \ge 1$, by a standard use of the Bogovskii operator, e.g. [13, Theorem III.3.1], we can construct $z_R \in W_0^{1,\alpha'}(B_R)^d$ such that

$$ \text{div}(z_R) = \left( T_R |T_R|^{\alpha-2} - \int_{B_R} T_R |T_R|^{\alpha-2} \right) \mathbb{1}_{B_R}, \quad T_R := \tilde{S} - \int_{B_R} \tilde{S}, $$

and

$$ \|\nabla z_R\|_{L^{\alpha'}(B_R)} \lesssim_{\alpha} \|T_R|T_R|^{\alpha-2}\|_{L^{\alpha'}(B_R)} \lesssim \left\| \tilde{S} - \int_{B_R} \tilde{S} \right\|_{L^{\alpha}(B_R)}^{\alpha-1}. $$

Testing (3.36) with $g = z_R$, we find

$$ \int_{\mathbb{R}^d} \tilde{S} \text{ div}(z_R) = 2 \int_{\mathbb{R}^d} D(z_R) : \tilde{q} - \int_{\mathbb{R}^d} z_R \cdot f, $$

and thus, using the properties of $z_R$,

$$ \left\| \tilde{S} - \int_{B_R} \tilde{S} \right\|_{L^{\alpha}(B_R)} \lesssim \|f\|_{W^{-1,\alpha}(B_R)} + \|\tilde{q}\|_{L^{\alpha}(B_R)}. $$

Combined with (3.35), this yields the conclusion (2.19).

<figure><img src="image_5.png" /></figure>

# 4. Homogenization

This section is devoted to the proofs of Theorems 1 and 3. While Tartar's oscillating test function method as used in [7] is not quite appropriate to the present setting without uniform separation, we provide an alternative argument based on div-curl ideas combined with the extension result in Theorem 4, as inspired by the work of Jikov ([24, 25]) on homogenization problems with stiff inclusions (see also [20, Section 3.2]).

## 4.1. Construction of correctors

We start with the proof of Theorem 1, which we shall deduce from our results in [7] for uniformly separated particles, via an approximation argument together with suitable a priori estimates. The improved pressure estimates in Proposition 2 are deduced simultaneously.
Step 3. Construction of extended pressure.

In view of e.g. [20, Proposition 12.10], the relation (3.34) for the extension $\tilde{q}$ ensures the existence of an associated pressure field $\tilde{S} \in L_{\text{loc}}^1(\mathbb{R}^d)$, uniquely defined up to a global additive constant, such that

$$ \int_{\mathbb{R}^d} D(g) : (2\tilde{q} - \tilde{S} \text{ Id}) = \int_{\mathbb{R}^d} g \cdot f \quad \forall g \in C_c^1(\mathbb{R}^d)^d, \qquad (3.36) $$

that is, $-\text{div}(2\tilde{q} - \tilde{S} \text{ Id}) = f$ in $\mathbb{R}^d$. It remains to prove bound (2.19) for $\tilde{S}$. For all $R \ge 1$, by a standard use of the Bogovskii operator, e.g. [13, Theorem III.3.1], we can construct $z_R \in W_0^{1,\alpha'}(B_R)^d$ such that

$$ \text{div}(z_R) = \left( T_R |T_R|^{\alpha-2} - \int_{B_R} T_R |T_R|^{\alpha-2} \right) \mathbb{1}_{B_R}, \quad T_R := \tilde{S} - \int_{B_R} \tilde{S}, $$

and

$$ \|\nabla z_R\|_{L^{\alpha'}(B_R)} \lesssim_{\alpha} \|T_R|T_R|^{\alpha-2}\|_{L^{\alpha'}(B_R)} \lesssim \left\| \tilde{S} - \int_{B_R} \tilde{S} \right\|_{L^{\alpha}(B_R)}^{\alpha-1}. $$

Testing (3.36) with $g = z_R$, we find

$$ \int_{\mathbb{R}^d} \tilde{S} \text{ div}(z_R) = 2 \int_{\mathbb{R}^d} D(z_R) : \tilde{q} - \int_{\mathbb{R}^d} z_R \cdot f, $$

and thus, using the properties of $z_R$,

$$ \left\| \tilde{S} - \int_{B_R} \tilde{S} \right\|_{L^{\alpha}(B_R)} \lesssim \|f\|_{W^{-1,\alpha}(B_R)} + \|\tilde{q}\|_{L^{\alpha}(B_R)}. $$

Combined with (3.35), this yields the conclusion (2.19).

<figure><img src="image_7.png" /></figure>

# 4. Homogenization

This section is devoted to the proofs of Theorems 1 and 3. While Tartar's oscillating test function method as used in [7] is not quite appropriate to the present setting without uniform separation, we provide an alternative argument based on div-curl ideas combined with the extension result in Theorem 4, as inspired by the work of Jikov ([24, 25]) on homogenization problems with stiff inclusions (see also [20, Section 3.2]).

## 4.1. Construction of correctors

We start with the proof of Theorem 1, which we shall deduce from our results in [7] for uniformly separated particles, via an approximation argument together with suitable a priori estimates. The improved pressure estimates in Proposition 2 are deduced simultaneously.
*Proofs of Theorem 1 and Proposition 2.* Given $2 \le r \ne \frac{2d}{d-2}$ and $1 \le \alpha \le 2 \wedge \frac{2dr}{r(d-2)+2d}$, with $\alpha < \frac{d}{d-1}$ if $r=2$, we assume that interparticle distances satisfy

$$ \sum_n \mathbb{E}[\mu_r(\rho_n)^{\frac{2\alpha}{2-\alpha}} 1_{0 \in I_n}] < \infty, \quad (4.1) $$

and we shall then prove Theorem 1, with pressure $\Sigma_E$ in $L^\alpha(\Omega)$ provided $\alpha > 1$. Optimizing in $r$ further yields Proposition 2. We split the proof into two main steps.

*Step 1.* Approximations with uniformly separated particles.

For $0 < \kappa \le \frac{\delta}{2}$, we consider the restricted inclusions

$$ I_n^\kappa := \{x \in I_n : \mathrm{dist}(x, \partial I_n) > \kappa\}, \quad \mathcal{I}^\kappa := \bigcup_n I_n^\kappa, $$

which still satisfy Assumptions (H) and (H') with $\delta$ replaced by $\frac{\delta}{2}$ and with minimal interparticle distance $\rho_n^\circ \ge \kappa$ (cf. (2.3)), that is, $(I_n^\kappa + \kappa B) \cap (I_m^\kappa + \kappa B) = \emptyset$ for all $n \ne m$. In this context with uniformly separated particles, we may apply [7, Proposition 2.1], which ensures the existence and uniqueness of a corrector $\psi_E^\kappa$ and of an associated pressure $\Sigma_E^\kappa$ that satisfy the different properties stated in Theorem 1 with $\mathcal{I}$ replaced by $\mathcal{I}^\kappa$. In addition, we show that the following moment bounds hold uniformly with respect to the parameter $\kappa > 0$: for $\alpha$ as in the moment condition (4.1),

$$ \mathbb{E}[|\nabla \psi_E^\kappa|^2] \lesssim |E|^2, \quad (4.2) $$

$$ \mathbb{E}[|\Sigma_E^\kappa|^\alpha 1_{\mathbb{R}^d \setminus \mathcal{I}^\kappa}] \lesssim |E|^\alpha. \quad (4.3) $$

These two estimates are established in the following two substeps.

*Substep 1.1.* Proof of (4.2).

In terms of the extension operator $P_n$ that we constructed in Proposition 3.1, we consider the stationary random vector field

$$ \phi_E^\circ := -\sum_n P_n(E(x - x_n)) $$

and we show that it satisfies

$$ (D(\phi_E^\circ) + E)|_{\mathcal{I}} = 0, \quad \mathrm{div}(\phi_E^\circ) = 0, \quad \mathbb{E}[|D(\phi_E^\circ)|^2] \lesssim |E|^2, \quad \mathbb{E}[D(\phi_E^\circ)] = 0. \quad (4.4) $$

The first two properties follow from the construction of $P_n$ with $D(P_n(E(x - x_n)))|_{I_n} = E$ and $\mathrm{div}(P_n(E(x - x_n))) = 0$. Next, stationarity allows us to write

$$ \mathbb{E}[|D(\phi_E^\circ)|^2] = \mathbb{E}\left[\int_B |D(\phi_E^\circ)|^2 \right] = \mathbb{E}\left[\frac{1}{|B|} \sum_n \int_{I_n^+ \cap B} |D(P_n(E(x - x_n)))|^2 \right], $$

hence, by Proposition 3.1,

$$ \mathbb{E}[|D(\phi_E^\circ)|^2] \lesssim |E|^2 \mathbb{E}\left[\frac{1}{|B|} \sum_{n: I_n^+ \cap B \neq \emptyset} \mu_2(\rho_n)^2\right]. $$
*Proofs of Theorem 1 and Proposition 2.* Given $2 \le r \ne \frac{2d}{d-2}$ and $1 \le \alpha \le 2 \wedge \frac{2dr}{r(d-2)+2d}$, with $\alpha < \frac{d}{d-1}$ if $r=2$, we assume that interparticle distances satisfy

$$ \sum_n \mathbb{E}[\mu_r(\rho_n)^{\frac{2\alpha}{2-\alpha}} 1_{0 \in I_n}] < \infty, \quad (4.1) $$

and we shall then prove Theorem 1, with pressure $\Sigma_E$ in $L^\alpha(\Omega)$ provided $\alpha > 1$. Optimizing in $r$ further yields Proposition 2. We split the proof into two main steps.

*Step 1.* Approximations with uniformly separated particles.

For $0 < \kappa \le \frac{\delta}{2}$, we consider the restricted inclusions

$$ I_n^\kappa := \{x \in I_n : \mathrm{dist}(x, \partial I_n) > \kappa\}, \quad \mathcal{I}^\kappa := \bigcup_n I_n^\kappa, $$

which still satisfy Assumptions (H) and (H') with $\delta$ replaced by $\frac{\delta}{2}$ and with minimal interparticle distance $\rho_n^\circ \ge \kappa$ (cf. (2.3)), that is, $(I_n^\kappa + \kappa B) \cap (I_m^\kappa + \kappa B) = \emptyset$ for all $n \ne m$. In this context with uniformly separated particles, we may apply [7, Proposition 2.1], which ensures the existence and uniqueness of a corrector $\psi_E^\kappa$ and of an associated pressure $\Sigma_E^\kappa$ that satisfy the different properties stated in Theorem 1 with $\mathcal{I}$ replaced by $\mathcal{I}^\kappa$. In addition, we show that the following moment bounds hold uniformly with respect to the parameter $\kappa > 0$: for $\alpha$ as in the moment condition (4.1),

$$ \mathbb{E}[|\nabla \psi_E^\kappa|^2] \lesssim |E|^2, \quad (4.2) $$

$$ \mathbb{E}[|\Sigma_E^\kappa|^\alpha 1_{\mathbb{R}^d \setminus \mathcal{I}^\kappa}] \lesssim |E|^\alpha. \quad (4.3) $$

These two estimates are established in the following two substeps.

*Substep 1.1.* Proof of (4.2).

In terms of the extension operator $P_n$ that we constructed in Proposition 3.1, we consider the stationary random vector field

$$ \phi_E^\circ := -\sum_n P_n(E(x - x_n)) $$

and we show that it satisfies

$$ (D(\phi_E^\circ) + E)|_{\mathcal{I}} = 0, \quad \mathrm{div}(\phi_E^\circ) = 0, \quad \mathbb{E}[|D(\phi_E^\circ)|^2] \lesssim |E|^2, \quad \mathbb{E}[D(\phi_E^\circ)] = 0. \quad (4.4) $$

The first two properties follow from the construction of $P_n$ with $D(P_n(E(x - x_n)))|_{I_n} = E$ and $\mathrm{div}(P_n(E(x - x_n))) = 0$. Next, stationarity allows us to write

$$ \mathbb{E}[|D(\phi_E^\circ)|^2] = \mathbb{E}\left[\int_B |D(\phi_E^\circ)|^2 \right] = \mathbb{E}\left[\frac{1}{|B|} \sum_n \int_{I_n^+ \cap B} |D(P_n(E(x - x_n)))|^2 \right], $$

hence, by Proposition 3.1,

$$ \mathbb{E}[|D(\phi_E^\circ)|^2] \lesssim |E|^2 \mathbb{E}\left[\frac{1}{|B|} \sum_{n: I_n^+ \cap B \neq \emptyset} \mu_2(\rho_n)^2\right]. $$
This can then be estimated, by stationarity, as

$$ \mathbb{E}[|D(\phi_E^\circ)|^2] \lesssim |E|^2 \mathbb{E}\left[\frac{1}{|B|} \int_{B_3} \sum_n \mu_2(\rho_n)^2 \mathbb{1}_{I_n}\right] \lesssim |E|^2 \mathbb{E}\left[\sum_n \mu_2(\rho_n)^2 \mathbb{1}_{0 \in I_n}\right], $$

so that the third property in (4.4) follows from the moment assumption (4.1) with $r = 2$ and $\alpha = 1$. Finally, stationarity allows us to write, for all $R > 0$,

$$ \mathbb{E}[D(\phi_E^\circ)] = \mathbb{E}\left[\int_{B_R} D(\phi_E^\circ)\right], $$

hence, inserting the definition of $D(\phi_E^\circ) \in L^2(\Omega)$, and letting $R \uparrow \infty$ to neglect boundary terms,

$$ \mathbb{E}[D(\phi_E^\circ)] = - \lim_{R \uparrow \infty} \mathbb{E}\left[|B_R|^{-1} \sum_{n: I_n^+ \subset B_R} \int_{I_n^+} D(P_n(E(x-x_n)))\right]. $$

Combined with the observation that $\int_{I_n^+} D(P_n(E(x-x_n))) = 0$, this concludes the proof of the last property in (4.4).

With this construction at hand, and noting that $I^\kappa \subset I$, testing the variational problem (2.6) for $\psi_E^\kappa$ with the test function $\phi_E^\circ$ yields

$$ \mathbb{E}[|D(\psi_E^\kappa) + E|^2] \le \mathbb{E}[|D(\phi_E^\circ) + E|^2] \lesssim |E|^2. \quad (4.5) $$

It remains to turn this into an a priori estimate on the full gradient $\nabla \psi_E^\kappa$. For that purpose, we decompose

$$ |\nabla \psi_E^\kappa|^2 = 2|D(\psi_E^\kappa)|^2 - \nabla_j(\psi_E^\kappa)_i \nabla_i(\psi_E^\kappa)_j. \quad (4.6) $$

For all $R \ge 1$, choose a smooth averaging function $\chi_R \in C_c^\infty(\mathbb{R}^d; \mathbb{R}^+)$ such that $\chi_R$ is constant in $B_R$, vanishes outside $B_{2R}$, and satisfies $\int_{\mathbb{R}^d} \chi_R = 1$ and $|\nabla \chi_R| \lesssim R^{-d-1}$. An integration by parts together with the constraint $\text{div}(\psi_E^\kappa) = 0$ yields

$$ \int_{\mathbb{R}^d} \chi_R \nabla_j (\psi_E^\kappa)_i \nabla_i (\psi_E^\kappa)_j = - \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_E^\kappa) : \nabla \psi_E^\kappa, $$

and thus, by definition of $\chi_R$ and by scaling,

$$ \left| \int_{\mathbb{R}^d} \chi_R \nabla_j (\psi_E^\kappa)_i \nabla_i (\psi_E^\kappa)_j \right| \lesssim \|R^{-1} \psi_E^\kappa(R \cdot)\|_{L^2(B_2)} \|\nabla \psi_E^\kappa(R \cdot)\|_{L^2(B_2)}. $$

Passing to the limit $R \uparrow \infty$ and appealing to the ergodic theorem, in view of the stationarity of $\nabla \psi_E^\kappa$ and the sublinearity of $\psi_E^\kappa$, cf. (2.9), we deduce $\mathbb{E}[\nabla_j(\psi_E^\kappa)_i \nabla_i(\psi_E^\kappa)_j] = 0$, so that the decomposition (4.6) entails $\mathbb{E}[|\nabla \psi_E^\kappa|^2] = 2 \mathbb{E}[|D(\psi_E^\kappa)|^2]$ and bound (4.5) yields claim (4.2).

Substep 1.2. Proof of (4.3).

We appeal to Lemma 3.2 in the following form (with $s = 2$ and $p = \alpha'$, with $\alpha, r$ as in the moment condition (4.1)): there exists $z_R^\kappa \in H_0^1(B_R)^d$ such that $D(z_R^\kappa)|_{I^\kappa} = 0$ and

$$ \text{div}(z_R^\kappa) = \left( T_R^\kappa |T_R^\kappa|^{\alpha-2} - \int_{B_R \setminus I^\kappa} T_R^\kappa |T_R^\kappa|^{\alpha-2} \right) \mathbb{1}_{B_R \setminus I^\kappa}, \quad T_R^\kappa := \Sigma_E^\kappa - \int_{B_R \setminus I^\kappa} \Sigma_E^\kappa, $$
This can then be estimated, by stationarity, as

$$ \mathbb{E}[|D(\phi_E^\circ)|^2] \lesssim |E|^2 \mathbb{E}\left[\frac{1}{|B|} \int_{B_3} \sum_n \mu_2(\rho_n)^2 \mathbb{1}_{I_n}\right] \lesssim |E|^2 \mathbb{E}\left[\sum_n \mu_2(\rho_n)^2 \mathbb{1}_{0 \in I_n}\right], $$

so that the third property in (4.4) follows from the moment assumption (4.1) with $r = 2$ and $\alpha = 1$. Finally, stationarity allows us to write, for all $R > 0$,

$$ \mathbb{E}[D(\phi_E^\circ)] = \mathbb{E}\left[\int_{B_R} D(\phi_E^\circ)\right], $$

hence, inserting the definition of $D(\phi_E^\circ) \in L^2(\Omega)$, and letting $R \uparrow \infty$ to neglect boundary terms,

$$ \mathbb{E}[D(\phi_E^\circ)] = - \lim_{R \uparrow \infty} \mathbb{E}\left[ |B_R|^{-1} \sum_{n: I_n^+ \subset B_R} \int_{I_n^+} D(P_n(E(x-x_n))) \right]. $$

Combined with the observation that $\int_{I_n^+} D(P_n(E(x-x_n))) = 0$, this concludes the proof of the last property in (4.4).

With this construction at hand, and noting that $I^\kappa \subset I$, testing the variational problem (2.6) for $\psi_E^\kappa$ with the test function $\phi_E^\circ$ yields

$$ \mathbb{E}[|D(\psi_E^\kappa) + E|^2] \le \mathbb{E}[|D(\phi_E^\circ) + E|^2] \lesssim |E|^2. \quad (4.5) $$

It remains to turn this into an a priori estimate on the full gradient $\nabla \psi_E^\kappa$. For that purpose, we decompose

$$ |\nabla \psi_E^\kappa|^2 = 2|D(\psi_E^\kappa)|^2 - \nabla_j(\psi_E^\kappa)_i \nabla_i(\psi_E^\kappa)_j. \quad (4.6) $$

For all $R \ge 1$, choose a smooth averaging function $\chi_R \in C_c^\infty(\mathbb{R}^d; \mathbb{R}^+)$ such that $\chi_R$ is constant in $B_R$, vanishes outside $B_{2R}$, and satisfies $\int_{\mathbb{R}^d} \chi_R = 1$ and $|\nabla \chi_R| \lesssim R^{-d-1}$. An integration by parts together with the constraint $\text{div}(\psi_E^\kappa) = 0$ yields

$$ \int_{\mathbb{R}^d} \chi_R \nabla_j (\psi_E^\kappa)_i \nabla_i (\psi_E^\kappa)_j = - \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_E^\kappa) : \nabla \psi_E^\kappa, $$

and thus, by definition of $\chi_R$ and by scaling,

$$ \left| \int_{\mathbb{R}^d} \chi_R \nabla_j (\psi_E^\kappa)_i \nabla_i (\psi_E^\kappa)_j \right| \lesssim \|R^{-1} \psi_E^\kappa(R \cdot)\|_{L^2(B_2)} \|\nabla \psi_E^\kappa(R \cdot)\|_{L^2(B_2)}. $$

Passing to the limit $R \uparrow \infty$ and appealing to the ergodic theorem, in view of the stationarity of $\nabla \psi_E^\kappa$ and the sublinearity of $\psi_E^\kappa$, cf. (2.9), we deduce $\mathbb{E}[\nabla_j(\psi_E^\kappa)_i \nabla_i(\psi_E^\kappa)_j] = 0$, so that the decomposition (4.6) entails $\mathbb{E}[|\nabla \psi_E^\kappa|^2] = 2 \mathbb{E}[|D(\psi_E^\kappa)|^2]$ and bound (4.5) yields claim (4.2).

Substep 1.2. Proof of (4.3).

We appeal to Lemma 3.2 in the following form (with $s = 2$ and $p = \alpha'$, with $\alpha, r$ as in the moment condition (4.1)): there exists $z_R^\kappa \in H_0^1(B_R)^d$ such that $D(z_R^\kappa)|_{I^\kappa} = 0$ and

$$ \text{div}(z_R^\kappa) = \left( T_R^\kappa |T_R^\kappa|^{\alpha-2} - \int_{B_R \setminus I^\kappa} T_R^\kappa |T_R^\kappa|^{\alpha-2} \right) \mathbb{1}_{B_R \setminus I^\kappa}, \quad T_R^\kappa := \Sigma_E^\kappa - \int_{B_R \setminus I^\kappa} \Sigma_E^\kappa. $$
and such that

$$
\begin{align*}
\| \nabla z_R^\kappa \|_{L^2(B_R)} &\lesssim_{\alpha,r} \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \| T_R^\kappa |T_R^\kappa|^{\alpha-2} \|_{L^{\alpha'}(B_R \setminus I^\kappa)} \\
&\lesssim \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \left\| \Sigma_E^\kappa - \int_{B_R \setminus I^\kappa} \Sigma_E^\kappa \right\|_{L^\alpha(B_R \setminus I^\kappa)}^{\alpha-1}.
\end{align*}
$$

Testing the corrector equation (2.8) for $(\psi_E^\kappa, \Sigma_E^\kappa)$ with this test function $z_R^\kappa$, we find

$$
\int_{B_R} \Sigma_E^\kappa \operatorname{div}(z_R^\kappa) = 2 \int_{B_R} D(z_R^\kappa) : D(\psi_E^\kappa),
$$

and thus, using the above properties of $z_R^\kappa$,

$$
\left\| \Sigma_E^\kappa - \int_{B_R \setminus I^\kappa} \Sigma_E^\kappa \right\|_{L^\alpha(B_R \setminus I^\kappa)} \lesssim_{\alpha,r} \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \|D(\psi_E^\kappa)\|_{L^2(B_R)}. \quad (4.7)
$$

Dividing both sides by $R^{d/\alpha}$, recalling the definition (2.20)–(2.21) of $\Lambda$, passing to the limit $R \uparrow \infty$, and appealing to the ergodic theorem, recalling that $\nabla\psi_E^\kappa$ and $\Sigma_E^\kappa$ are stationary with vanishing expectation, and using (4.2), we deduce

$$
\mathbb{E}[|\Sigma_E^\kappa|^\alpha 1_{\mathbb{R}^d \setminus I^\kappa}]^{\frac{1}{\alpha}} \lesssim |E| \left( 1 + \sum_n \mathbb{E}[\mu_r(\rho_n)^{\frac{2\alpha}{2-\alpha}} 1_{0 \in I_n}] \right)^{\frac{2-\alpha}{2\alpha}},
$$

and claim (4.3) follows from the moment assumption (4.1).

Step 2. Conclusion.

In view of the uniform bounds (4.2) and (4.3), provided $\alpha > 1$, we may consider some
weak limit point $(\nabla\psi_E, \Sigma_E 1_{\mathbb{R}^d\setminus I})$ of $\{((\nabla\psi_E^\kappa, \Sigma_E^\kappa)_{\kappa>0}) \in L^2(\Omega)^{d\times d} \times L^\alpha(\Omega)$
as $\kappa \downarrow 0$. It follows that $\nabla\psi_E$ is stationary with vanishing expectation and finite second
moments, that it satisfies $\mathrm{div}(\psi_E) = 0$ and $(\mathrm{D}(\psi_E) + E)|_I = 0$, and that $\mathrm{D}(\psi_E)$ is the
unique solution of the limiting variational problem (2.6). Moreover, passing to the limit in
the weak formulation of (2.8), we find

$$
2 \int_{\mathbb{R}^d} \mathrm{D}(g) : \mathrm{D}(\psi_E) = \int_{\mathbb{R}^d} \Sigma_E \mathrm{div}(g) \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \mathrm{D}(g)|_I = 0, \quad (4.8)
$$

hence, in particular,

$$
- \Delta \psi_E + \nabla \Sigma_E = 0 \quad \text{in } \mathbb{R}^d \setminus \mathcal{I}. \tag{4.9}
$$

The pressure field $\Sigma_E$ in this equation is uniquely defined up to a global constant in view
of the almost sure connectedness of $\mathbb{R}^d \setminus \mathcal{I}$, and is thus fully determined by the condi-
tion $\mathbb{E}[\Sigma_E 1_{\mathbb{R}^d \setminus \mathcal{I}}] = 0$. In addition, in view of the regularity of the particle boundaries,
cf. Assumption (H)\textsubscript{0}, the regularity theory for the Stokes equation (e.g. [13, Section IV])
entails that $(\psi_E, \Sigma_E)$ is $\mathbb{C}^2$ smooth in $\mathbb{R}^d \setminus \mathcal{I}$ up to the boundary, and equation (4.9) is
thus satisfied in the strong sense. Next, for all $n$, for all $V \in \mathbb{R}^d$ and $\Theta \in \mathbb{M}^{\text{skew}}$, in terms
of the cutoff function $w_n$ that we constructed in Lemma 3.3, we may test equation (4.8)
and such that

$$
\begin{align*}
\| \nabla z_R^\kappa \|_{L^2(B_R)} &\lesssim_{\alpha,r} \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \| T_R^\kappa |T_R^\kappa|^{\alpha-2} \|_{L^{\alpha'}(B_R \setminus I^\kappa)} \\
&\lesssim \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \left\| \Sigma_E^\kappa - \int_{B_R \setminus I^\kappa} \Sigma_E^\kappa \right\|_{L^\alpha(B_R \setminus I^\kappa)}^{\alpha-1}.
\end{align*}
$$

Testing the corrector equation (2.8) for $(\psi_E^\kappa, \Sigma_E^\kappa)$ with this test function $z_R^\kappa$, we find

$$
\int_{B_R} \Sigma_E^\kappa \operatorname{div}(z_R^\kappa) = 2 \int_{B_R} \mathrm{D}(z_R^\kappa) : \mathrm{D}(\psi_E^\kappa),
$$

and thus, using the above properties of $z_R^\kappa$,

$$
\left\| \Sigma_E^\kappa - \int_{B_R \setminus I^\kappa} \Sigma_E^\kappa \right\|_{L^\alpha(B_R \setminus I^\kappa)} \lesssim_{\alpha,r} \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \|D(\psi_E^\kappa)\|_{L^2(B_R)}. \quad (4.7)
$$

Dividing both sides by $R^{d/\alpha}$, recalling the definition (2.20)-(2.21) of $\Lambda$, passing to the limit $R \uparrow \infty$, and appealing to the ergodic theorem, recalling that $\nabla\psi_E^\kappa$ and $\Sigma_E^\kappa$ are stationary with vanishing expectation, and using (4.2), we deduce

$$
\mathbb{E}[|\Sigma_E^\kappa|^\alpha 1_{\mathbb{R}^d \setminus I^\kappa}]^{\frac{1}{\alpha}} \lesssim |E| \left( 1 + \sum_n \mathbb{E}[\mu_r(\rho_n)^{\frac{2\alpha}{2-\alpha}} 1_{0 \in I_n}] \right)^{\frac{2-\alpha}{2\alpha}},
$$

and claim (4.3) follows from the moment assumption (4.1).

Step 2. Conclusion.

In view of the uniform bounds (4.2) and (4.3), provided $\alpha > 1$, we may consider some
weak limit point $(\nabla\psi_E, \Sigma_E 1_{\mathbb{R}^d\setminus I})$ of $\{((\nabla\psi_E^\kappa, \Sigma_E^\kappa)_{\kappa>0}) \in L^2(\Omega)^{d\times d} \times L^\alpha(\Omega)$
as $\kappa \downarrow 0$. It follows that $\nabla\psi_E$ is stationary with vanishing expectation and finite second
moments, that it satisfies $\mathrm{div}(\psi_E) = 0$ and $(\mathrm{D}(\psi_E) + E)|_I = 0$, and that $\mathrm{D}(\psi_E)$ is the
unique solution of the limiting variational problem (2.6). Moreover, passing to the limit in
the weak formulation of (2.8), we find

$$
2 \int_{\mathbb{R}^d} \mathrm{D}(g) : \mathrm{D}(\psi_E) = \int_{\mathbb{R}^d} \Sigma_E \operatorname{div}(g) \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \mathrm{D}(g)|_I = 0, \quad (4.8)
$$

hence, in particular,

$$
- \Delta \psi_E + \nabla \Sigma_E = 0 \quad \text{in } \mathbb{R}^d \setminus \mathcal{I}. \tag{4.9}
$$

The pressure field $\Sigma_E$ in this equation is uniquely defined up to a global constant in view
of the almost sure connectedness of $\mathbb{R}^d \setminus \mathcal{I}$, and is thus fully determined by the condi-
tion $\mathbb{E}[\Sigma_E 1_{\mathbb{R}^d \setminus \mathcal{I}}] = 0$. In addition, in view of the regularity of the particle boundaries,
cf. Assumption (H)\textsubscript{0}, the regularity theory for the Stokes equation (e.g. [13, Section IV])
entails that $(\psi_E, \Sigma_E)$ is $\mathbb{C}^2$ smooth in $\mathbb{R}^d \setminus \mathcal{I}$ up to the boundary, and equation (4.9) is
thus satisfied in the strong sense. Next, for all $n$, for all $V \in \mathbb{R}^d$ and $\Theta \in \mathbb{M}^{\text{skew}}$, in terms
of the cutoff function $w_n$ that we constructed in Lemma 3.3, we may test equation (4.8)
with $g = w_n(V + \Theta(x - x_n)) \in W_0^{1,\alpha'}(I_n^+)$, which indeed satisfies $D(g)|_I = 0$, and an integration by parts then yields

$$
\begin{align*}
0 &= \int_{\mathbb{R}^d} D(w_n(V + \Theta(x - x_n))) : \sigma(\psi_E + Ex, \Sigma_E) \\
&= - \int_{\partial I_n} (V + \Theta(x - x_n)) \cdot \sigma(\psi_E + Ex, \Sigma_E) \nu,
\end{align*}
$$

showing that the boundary conditions in (2.8) are almost surely satisfied in a pointwise
sense. Finally, the weak convergence of $(\nabla \psi_E, \Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I})(\frac{1}{\varepsilon})$ to 0 in (2.9) follows from
$E[(\nabla \psi_E, \Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I})] = 0$ by the ergodic theorem, while the sublinearity of $\psi_E$ in form
of the strong convergence of $\varepsilon \psi_E(\frac{1}{\varepsilon})$ to 0 is a standard result for random fields with
stationary gradient having vanishing expectation, e.g. [20, Section 7].
■

**4.2. Extension of fluxes**

Applying Theorem 4 to the corrector $\psi_E$, cf. (2.8), and to the solution $u_\varepsilon$ of the Stokes problem (2.13), we obtain the following useful extension result for the fluxes:

$$
q_E := D(\psi_E) + E, \quad p_\varepsilon := D(u_\varepsilon).
$$

**Corollary 4.1.** On top of Assumptions (H<sub>δ</sub><sup>o</sup>) and (H<sub>δ</sub><sup>r</sup>), given a bounded Lipschitz domain U ⊂ ℝ<sup>d</sup>, and given 2 ≤ r ≠ 2>d−2 and 1 < α ≤ 2 ∧ 2dr>r(d−2)+2d, with α < d>d−1 if r = 2, assume that interparticle distances satisfy the following moment condition, almost surely:

$$
\limsup_{\epsilon \downarrow 0} \epsilon^d \sum_{n \in \mathcal{N}_\epsilon(U)} \mu_r(\rho_n; U, \epsilon)^{\frac{2\alpha}{2-\alpha}} < \infty. \quad (4.10)
$$

Then the following properties hold:

(i) For all $E \in \mathbb{M}_0^{\text{sym}}$, there exist a stationary element $\tilde{q}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d)^{d \times d})$ with $\text{tr}(\tilde{q}_E) = 0$, and associated stationary pressure field $\tilde{\Sigma}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d))$, such that almost surely,

$$
\begin{equation}
\begin{aligned}
(\tilde{q}_E, \tilde{\Sigma}_E)|_{\mathbb{R}^d \setminus I} &= (q_E, \Sigma_E)|_{\mathbb{R}^d \setminus I}, \\
&\quad - \text{div}(2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}) = 0 \quad \text{in } \mathbb{R}^d,
\end{aligned}
\tag{4.11}
\end{equation}
$$

and

$$
\| \tilde{q}_E \|_{L^\alpha(\Omega)} + \| \tilde{\Sigma}_E - \mathbb{E}[\tilde{\Sigma}_E] \|_{L^\alpha(\Omega)} \lesssim_{\alpha,r} |E|.
$$

(ii) There exists $\tilde{p}_\epsilon \in L^\alpha(\Omega; L^\alpha(U))^{d \times d}_{\text{sym}}$ with $\text{tr}(\tilde{p}_\epsilon) = 0$, and an associated pressure field $\tilde{S}_\epsilon \in L^\alpha(\Omega; L^\alpha(U))$, such that almost surely,

$$
\begin{equation}
\begin{split}
(\tilde{p}_\epsilon, \tilde{S}_\epsilon)|_{U \setminus I_\epsilon(U)} &= (p_\epsilon, S_\epsilon)|_{U \setminus I_\epsilon(U)}, \\
&\quad - \operatorname{div}(2\tilde{p}_\epsilon - \tilde{S}_\epsilon \operatorname{Id}) = f \mathbb{1}_{U \setminus I_\epsilon(U)} \quad \text{in } U,
\end{split}
\tag{4.12}
\end{equation}
$$

and

$$
\limsup_{\epsilon \downarrow 0} \left( \| \tilde{p}_\epsilon \|_{L^\alpha(U)} + \left\| \tilde{S}_\epsilon - \int_U \tilde{S}_\epsilon \right\|_{L^\alpha(U)} \right) \lesssim_{U, \alpha, r} \| f \|_{L^{\frac{2d}{d+2}}(U)}.
$$
with $g = w_n(V + \Theta(x - x_n)) \in W_0^{1,\alpha'}(I_n^+)$, which indeed satisfies $D(g)|_I = 0$, and an integration by parts then yields

$$ 
\begin{aligned}
0 &= \int_{\mathbb{R}^d} D(w_n(V + \Theta(x - x_n))) : \sigma(\psi_E + Ex, \Sigma_E) \\
&= -\int_{\partial I_n} (V + \Theta(x - x_n)) \cdot \sigma(\psi_E + Ex, \Sigma_E) \nu,
\end{aligned}
$$

showing that the boundary conditions in (2.8) are almost surely satisfied in a pointwise sense. Finally, the weak convergence of $(\nabla \psi_E, \Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I})(\frac{\cdot}{\varepsilon})$ to 0 in (2.9) follows from $\mathbb{E}[(\nabla \psi_E, \Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I})] = 0$ by the ergodic theorem, while the sublinearity of $\psi_E$ in form of the strong convergence of $\varepsilon\psi_E(\frac{\cdot}{\varepsilon})$ to 0 is a standard result for random fields with stationary gradient having vanishing expectation, e.g. [20, Section 7]. ■

## 4.2. Extension of fluxes

Applying Theorem 4 to the corrector $\psi_E$, cf. (2.8), and to the solution $u_\varepsilon$ of the Stokes problem (2.13), we obtain the following useful extension result for the fluxes:

$$ q_E := D(\psi_E) + E, \quad p_\varepsilon := D(u_\varepsilon). $$

**Corollary 4.1.** On top of Assumptions (H) and (H'), given a bounded Lipschitz domain $U \subset \mathbb{R}^d$, and given $2 \le r \ne \frac{2d}{d-2}$ and $1 < \alpha \le 2 \wedge \frac{2dr}{r(d-2)+2d}$, with $\alpha < \frac{d}{d-1}$ if $r = 2$, assume that interparticle distances satisfy the following moment condition, almost surely:

$$ \limsup_{\varepsilon \downarrow 0} \varepsilon^d \sum_{n \in \mathcal{N}_\varepsilon(U)} \mu_r(\rho_{n;U,\varepsilon})^{\frac{2\alpha}{2-\alpha}} < \infty. \quad (4.10) $$

Then the following properties hold:

(i) For all $E \in \mathbb{M}_0^{\text{sym}}$, there exist a stationary element $\tilde{q}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d)^{d \times d})$ with $\text{tr}(\tilde{q}_E) = 0$, and associated stationary pressure field $\tilde{\Sigma}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d))$, such that almost surely,

$$ (\tilde{q}_E, \tilde{\Sigma}_E)|_{\mathbb{R}^d \setminus I} = (q_E, \Sigma_E)|_{\mathbb{R}^d \setminus I}, \\ - \text{div}(2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}) = 0 \quad \text{in } \mathbb{R}^d, \quad (4.11) $$

and

$$ \|\tilde{q}_E\|_{L^\alpha(\Omega)} + \|\tilde{\Sigma}_E - \mathbb{E}[\tilde{\Sigma}_E]\|_{L^\alpha(\Omega)} \lesssim_{\alpha,r} |E|. $$

(ii) There exists $\tilde{p}_\varepsilon \in L^\alpha(\Omega; L^\alpha(U))^{d \times d}_{\text{sym}}$ with $\text{tr}(\tilde{p}_\varepsilon) = 0$, and an associated pressure field $\tilde{S}_\varepsilon \in L^\alpha(\Omega; L^\alpha(U))$, such that almost surely,

$$ (\tilde{p}_\varepsilon, \tilde{S}_\varepsilon)|_{U \setminus I_\varepsilon(U)} = (p_\varepsilon, S_\varepsilon)|_{U \setminus I_\varepsilon(U)}, \\ - \text{div}(2\tilde{p}_\varepsilon - \tilde{S}_\varepsilon \text{Id}) = f \mathbb{1}_{U \setminus I_\varepsilon(U)} \quad \text{in } U, \quad (4.12) $$

and

$$ \limsup_{\varepsilon \downarrow 0} \left( \| \tilde{p}_\varepsilon \|_{L^\alpha(U)} + \left\| \tilde{S}_\varepsilon - \int_U \tilde{S}_\varepsilon \right\|_{L^\alpha(U)} \right) \lesssim_{U, \alpha, r} \| f \|_{L^{\frac{2d}{d+2}}(U)}. $$
*Proof.* We split the proof into two steps.

Step 1. Proof of (i)

The corrector equation (2.8) ensures that the flux $q_E = D(\psi_E) + E \in L^2(\Omega; L_{\text{loc}}^2(\mathbb{R}^d)^{d \times d})$ satisfies $\text{tr}(q_E) = 0$ and

$$ \int_{\mathbb{R}^d} D(g) : q_E = 0 \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \text{div}(g) = 0, D(g)|_{I} = 0. $$

Given $\alpha, r$ as in (4.10), Theorem 4 provides an extension $\tilde{q}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d)^{d \times d})$ with $\text{tr}(\tilde{q}_E) = 0$, and an associated pressure field $\tilde{\Sigma}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d))$, such that

$$ \tilde{q}_E|_{\mathbb{R}^d \setminus I} = q_E|_{\mathbb{R}^d \setminus I}, \quad \text{and} \quad -\text{div}(2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}) = 0 \text{ in } \mathbb{R}^d, \quad (4.13) $$

and such that the following estimate holds, for all $R \ge 1$:

$$ \|\tilde{q}_E\|_{L^\alpha(B_R)} + \left\|\tilde{\Sigma}_E - \int_{B_R} \tilde{\Sigma}_E\right\|_{L^\alpha(B_R)} \lesssim_{\alpha,r} \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \|q_E\|_{L^2(\hat{B}_R \setminus I)}. \quad (4.14) $$

In addition, the construction in the proof of Theorem 4 ensures that $\tilde{q}_E$ can be chosen stationary. Since $\tilde{q}_E$ coincides with $q_E = D(\psi_E) + E$ on $\mathbb{R}^d \setminus I$, we deduce from (4.13), in particular,

$$ -\Delta\psi_E + \nabla\tilde{\Sigma}_E = 0 \quad \text{in } \mathbb{R}^d \setminus I. $$

In view of (2.8), recalling that $\mathbb{R}^d \setminus I$ is almost surely connected, we deduce that the pressure $\tilde{\Sigma}_E$ must coincide with $\Sigma_E$ in $\mathbb{R}^d \setminus I$ up to a global constant. Therefore, $\tilde{\Sigma}_E$ is uniquely determined for instance by the choice $\tilde{\Sigma}_E|_{\mathbb{R}^d \setminus I} = \Sigma_E|_{\mathbb{R}^d \setminus I}$. For this choice, as $\tilde{q}_E$ and $\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I}$ are stationary, uniqueness entails that $\tilde{\Sigma}_E$ is also stationary.

Dividing both sides of (4.14) by $R^{d/\alpha}$, recalling definition (2.20) of $\Lambda$, passing to the limit $R \uparrow \infty$, appealing to the ergodic theorem, in view of the stationarity of $\tilde{q}_E$, $\tilde{\Sigma}_E$, and using the energy bound $\mathbb{E}[[q_E|^2]] \lesssim |E|^2$, we obtain

$$ \|\tilde{q}_E\|_{L^\alpha(\Omega)} + \|\tilde{\Sigma}_E - \mathbb{E}[\tilde{\Sigma}_E]\|_{L^\alpha(\Omega)} \lesssim_{\alpha,r} |E| \left( 1 + \sum_n \mathbb{E}[\mu_r(\rho_n)^{\frac{2\alpha}{2-\alpha}} \mathbb{1}_{0 \in I_n}] \right)^{\frac{2-\alpha}{2\alpha}}. $$

Combined with the moment condition (4.10), with $\mu_r(\rho_n) \le \mu_r(\rho_n; U, \varepsilon)$, this yields the conclusion.

Step 2. Proof of (ii).

Equation (2.13) ensures that the flux $p_\varepsilon = D(u_\varepsilon) \in L^2(\Omega; L^2(U))_{\text{sym}}^{d \times d}$ satisfies $\text{tr}(p_\varepsilon) = 0$ and

$$ 2 \int_U D(g) : p_\varepsilon = \int_{U \setminus I_\varepsilon(U)} g \cdot f \quad \forall g \in C_c^1(U)^d : \text{div}(g) = 0, D(g)|_{I_\varepsilon(U)} = 0. $$

Given $\alpha, r$ as in (4.10), Theorem 4 provides an extension $\tilde{p}_\varepsilon \in L^\alpha(\Omega; L^\alpha(U))_{\text{sym}}^{d \times d}$ with $\text{tr}(\tilde{p}_\varepsilon) = 0$, and an associated pressure field $\tilde{S}_\varepsilon \in L^\alpha(\Omega; L^\alpha(U))$, such that

$$ \tilde{p}_\varepsilon|_{U \setminus I_\varepsilon(U)} = p_\varepsilon|_{U \setminus I_\varepsilon(U)}, \quad \text{and} \quad -\text{div}(2\tilde{p}_\varepsilon - \tilde{S}_\varepsilon \text{Id}) = f \mathbb{1}_{U \setminus I_\varepsilon(U)} \text{ in } U, $$
*Proof.* We split the proof into two steps.

Step 1. Proof of (i)

The corrector equation (2.8) ensures that the flux $q_E = D(\psi_E) + E \in L^2(\Omega; L_{\text{loc}}^2(\mathbb{R}^d)^{d \times d})$ satisfies $\text{tr}(q_E) = 0$ and

$$ \int_{\mathbb{R}^d} D(g) : q_E = 0 \quad \forall g \in C_c^1(\mathbb{R}^d)^d : \text{div}(g) = 0, D(g)|_{I} = 0. $$

Given $\alpha, r$ as in (4.10), Theorem 4 provides an extension $\tilde{q}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d)^{d \times d})$ with $\text{tr}(\tilde{q}_E) = 0$, and an associated pressure field $\tilde{\Sigma}_E \in L^\alpha(\Omega; L_{\text{loc}}^\alpha(\mathbb{R}^d))$, such that

$$ \tilde{q}_E|_{\mathbb{R}^d \setminus I} = q_E|_{\mathbb{R}^d \setminus I}, \quad \text{and} \quad -\text{div}(2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}) = 0 \text{ in } \mathbb{R}^d, \quad (4.13) $$

and such that the following estimate holds, for all $R \ge 1$:

$$ \|\tilde{q}_E\|_{L^\alpha(B_R)} + \left\|\tilde{\Sigma}_E - \int_{B_R} \tilde{\Sigma}_E\right\|_{L^\alpha(B_R)} \lesssim_{\alpha,r} \Lambda(B_R; r, \frac{2\alpha}{2-\alpha}) \|q_E\|_{L^2(\hat{B}_R \setminus I)}. \quad (4.14) $$

In addition, the construction in the proof of Theorem 4 ensures that $\tilde{q}_E$ can be chosen stationary. Since $\tilde{q}_E$ coincides with $q_E = D(\psi_E) + E$ on $\mathbb{R}^d \setminus I$, we deduce from (4.13), in particular,

$$ -\Delta\psi_E + \nabla\tilde{\Sigma}_E = 0 \quad \text{in } \mathbb{R}^d \setminus I. $$

In view of (2.8), recalling that $\mathbb{R}^d \setminus I$ is almost surely connected, we deduce that the pressure $\tilde{\Sigma}_E$ must coincide with $\Sigma_E$ in $\mathbb{R}^d \setminus I$ up to a global constant. Therefore, $\tilde{\Sigma}_E$ is uniquely determined for instance by the choice $\tilde{\Sigma}_E|_{\mathbb{R}^d \setminus I} = \Sigma_E|_{\mathbb{R}^d \setminus I}$. For this choice, as $\tilde{q}_E$ and $\Sigma_E \mathbb{1}_{\mathbb{R}^d \setminus I}$ are stationary, uniqueness entails that $\tilde{\Sigma}_E$ is also stationary.

Dividing both sides of (4.14) by $R^{d/\alpha}$, recalling definition (2.20) of $\Lambda$, passing to the limit $R \uparrow \infty$, appealing to the ergodic theorem, in view of the stationarity of $\tilde{q}_E$, $\tilde{\Sigma}_E$, and using the energy bound $\mathbb{E}[[q_E|^2]] \lesssim |E|^2$, we obtain

$$ \|\tilde{q}_E\|_{L^\alpha(\Omega)} + \|\tilde{\Sigma}_E - \mathbb{E}[\tilde{\Sigma}_E]\|_{L^\alpha(\Omega)} \lesssim_{\alpha,r} |E| \left( 1 + \sum_n \mathbb{E}[\mu_r(\rho_n)^{\frac{2\alpha}{2-\alpha}} \mathbb{1}_{0 \in I_n}] \right)^{\frac{2-\alpha}{2\alpha}}. $$

Combined with the moment condition (4.10), with $\mu_r(\rho_n) \le \mu_r(\rho_n; U, \varepsilon)$, this yields the conclusion.

Step 2. Proof of (ii).

Equation (2.13) ensures that the flux $p_\varepsilon = D(u_\varepsilon) \in L^2(\Omega; L^2(U))_{\text{sym}}^{d \times d}$ satisfies $\text{tr}(p_\varepsilon) = 0$ and

$$ 2 \int_U D(g) : p_\varepsilon = \int_{U \setminus I_\varepsilon(U)} g \cdot f \quad \forall g \in C_c^1(U)^d : \text{div}(g) = 0, D(g)|_{I_\varepsilon(U)} = 0. $$

Given $\alpha, r$ as in (4.10), Theorem 4 provides an extension $\tilde{p}_\varepsilon \in L^\alpha(\Omega; L^\alpha(U))_{\text{sym}}^{d \times d}$ with $\text{tr}(\tilde{p}_\varepsilon) = 0$, and an associated pressure field $\tilde{S}_\varepsilon \in L^\alpha(\Omega; L^\alpha(U))$, such that

$$ \tilde{p}_\varepsilon|_{U \setminus I_\varepsilon(U)} = p_\varepsilon|_{U \setminus I_\varepsilon(U)}, \quad \text{and} \quad -\text{div}(2\tilde{p}_\varepsilon - \tilde{S}_\varepsilon \text{Id}) = f \mathbb{1}_{U \setminus I_\varepsilon(U)} \text{ in } U, $$
and by scaling, the estimate (2.19) takes on the following guise:

$$
\| \tilde{p}_{\varepsilon} \|_{L^{\alpha}(U)} + \left\| \tilde{S}_{\varepsilon} - \int_{U} \tilde{S}_{\varepsilon} \right\|_{L^{\alpha}(U)} \lesssim_{U, \alpha, r} \left( |U| + \varepsilon^d \sum_{n \in \mathcal{N}_{\varepsilon}(U)} \mu_r (\rho_n; U, \varepsilon)^{\frac{2\alpha}{2-\alpha}} \right)^{\frac{2-\alpha}{2\alpha}} \\
\times ( \|f\|_{L^{\frac{2d}{d+2}}(U)}^{\frac{2d}{d+2}} + \|p_{\varepsilon}\|_{L^2(U \setminus I_{\varepsilon}(U))} ).
$$

Combined with the energy bound $\|p_\epsilon\|_{L^2(U \setminus I_\epsilon(U))} \lesssim \|f\|_{L^{2d/(d+2)}(U)}$, and with the moment condition (4.10), this yields the conclusion. ■

**Remark 4.2.** In view of the construction in the proof of Theorem **4**, it is easily checked that the above-constructed extended fluxes $\tilde{q}_E$, $\tilde{p}_\epsilon$ can be viewed as limiting fluxes for corresponding Stokes problems with a suspension of droplets with diverging shear viscosity. More precisely, for all $\kappa > 0$, we consider the following corrector problem:

$$
\begin{cases}
-\operatorname{div}(2(1_{\mathbb{R}^d \setminus I} + \kappa \mathbb{1}_I)(\mathrm{D}(\psi_E^\kappa) + E)) + \nabla \Sigma_E^\kappa = 0 & \text{in } \mathbb{R}^d, \\
\operatorname{div}(\psi_E^\kappa) = 0 & \text{in } \mathbb{R}^d.
\end{cases}
$$

Under the assumptions of Corollary 4.1, in the limit $\kappa \uparrow \infty$, there holds $\mathrm{D}(\psi_E^\kappa) \to \mathrm{D}(\psi_E)$ in $L^2(\Omega)$ and corresponding fluxes converge,

$$
2(1_{\mathbb{R}^d \setminus I} + \kappa \mathbb{1}_I)(\mathrm{D}(\psi_E^\kappa) + E) - \Sigma_E^\kappa \mathrm{Id} \to \tilde{q}_E \quad \text{in } L^\alpha(\Omega),
$$

and a similar result holds for $\tilde{p}_\epsilon$. We skip the detail for brevity.

Next we compute $\mathbb{E}[\tilde{q}_E]$ and $\mathbb{E}[\tilde{\Sigma}_E]$, which happen to provide alternative definitions of the effective constants $\bar{B}, \bar{b}$. Note in particular that these ensemble averages do not depend on the actual choice of the extension $\tilde{q}_E$ in Corollary 4.1 (i).

**Lemma 4.3** (Effective constants). On top of Assumptions (H°) and (H'°), let $(\tilde{q}_E, \tilde{\Sigma}_E)$ be defined as in Corollary 4.1 (i) for some $\alpha > 1$. Then we have almost surely, as $\varepsilon \downarrow 0$,

$$
\tilde{q}_E(\cdot) \to \mathbb{E}[\tilde{q}_E] = \bar{B}E, \quad \tilde{\Sigma}_E(\cdot) \to \mathbb{E}[\tilde{\Sigma}_E] = -\bar{b} : E, \quad \text{weakly in } L_{\text{loc}}^{\alpha}(\mathbb{R}^d). \quad (4.15)
$$

In addition, provided $\alpha \ge \frac{2d}{d+2}$, these convergences are almost surely strong in $H_{\text{loc}}^{-1}(\mathbb{R}^d)$.

*Proof.* We split the proof into two steps.

*Step 1.* Proof of weak convergences (4.15).

As $\tilde{q}_E$ and $\tilde{\Sigma}_E$ are stationary, the ergodic theorem implies almost surely the weak convergences $\tilde{q}_E(\cdot|_E) \to \mathbb{E}[\tilde{q}_E]$ and $\tilde{\Sigma}_E(\cdot|_E) \to \mathbb{E}[\tilde{\Sigma}_E]$ in $L_{\text{loc}}^{\alpha}(\mathbb{R}^d)$, and it remains to compute these two expectations. For that purpose, up to an approximation argument as in the proof of Theorem **1**, we may assume without loss of generality $\alpha > \frac{2d}{d+2}$. We split the proof into two further substeps.

*Substep 1.1.* Proof that $\bar{B}E = \mathbb{E}[\tilde{q}_E]$.

For all $R \ge 1$, we set $\chi_R := R^{-d} \chi(\frac{1}{R})$, for some smooth averaging function $\chi \in$
and by scaling, the estimate (2.19) takes on the following guise:

$$
\| \tilde{p}_{\varepsilon} \|_{L^{\alpha}(U)} + \left\| \tilde{S}_{\varepsilon} - \int_{U} \tilde{S}_{\varepsilon} \right\|_{L^{\alpha}(U)} \lesssim_{U, \alpha, r} \left( |U| + \varepsilon^d \sum_{n \in \mathcal{N}_{\varepsilon}(U)} \mu_r (\rho_n; U, \varepsilon)^{\frac{2\alpha}{2-\alpha}} \right)^{\frac{2-\alpha}{2\alpha}} \\
\times ( \|f\|_{L^{\frac{2d}{d+2}}(U)}^{\frac{2d}{d+2}} + \|p_{\varepsilon}\|_{L^2(U \setminus I_{\varepsilon}(U))} ).
$$

Combined with the energy bound $\|p_\epsilon\|_{L^2(U \setminus I_\epsilon(U))} \lesssim \|f\|_{L^{2d/(d+2)}(U)}$, and with the moment condition (4.10), this yields the conclusion. ■

**Remark 4.2.** In view of the construction in the proof of Theorem **4**, it is easily checked that the above-constructed extended fluxes $\tilde{q}_E$, $\tilde{p}_\epsilon$ can be viewed as limiting fluxes for corresponding Stokes problems with a suspension of droplets with diverging shear viscosity. More precisely, for all $\kappa > 0$, we consider the following corrector problem:

$$
\begin{cases}
-\operatorname{div}(2(1_{\mathbb{R}^d \setminus I} + \kappa \mathbb{1}_I)(\mathrm{D}(\psi_E^\kappa) + E)) + \nabla \Sigma_E^\kappa = 0 & \text{in } \mathbb{R}^d, \\
\operatorname{div}(\psi_E^\kappa) = 0 & \text{in } \mathbb{R}^d.
\end{cases}
$$

Under the assumptions of Corollary 4.1, in the limit $\kappa \uparrow \infty$, there holds $\mathrm{D}(\psi_E^\kappa) \to \mathrm{D}(\psi_E)$ in $L^2(\Omega)$ and corresponding fluxes converge,

$$
2(1_{\mathbb{R}^d \setminus I} + \kappa \mathbb{1}_I)(\mathrm{D}(\psi_E^\kappa) + E) - \Sigma_E^\kappa \mathrm{Id} \to \tilde{q}_E \quad \text{in } L^\alpha(\Omega),
$$

and a similar result holds for $\tilde{p}_\epsilon$. We skip the detail for brevity.

Next we compute $\mathbb{E}[\tilde{q}_E]$ and $\mathbb{E}[\tilde{\Sigma}_E]$, which happen to provide alternative definitions of the effective constants $\bar{B}, \bar{b}$. Note in particular that these ensemble averages do not depend on the actual choice of the extension $\tilde{q}_E$ in Corollary 4.1 (i).

**Lemma 4.3** (Effective constants). On top of Assumptions (H°_δ) and (H'_δ), let $(\tilde{q}_E, \tilde{\Sigma}_E)$ be defined as in Corollary 4.1 (i) for some $\alpha > 1$. Then we have almost surely, as $\varepsilon \downarrow 0$,

$$
\tilde{q}_E(\cdot) \to \mathbb{E}[\tilde{q}_E] = \bar{B}E, \quad \tilde{\Sigma}_E(\cdot) \to \mathbb{E}[\tilde{\Sigma}_E] = -\bar{b} : E, \quad \text{weakly in } L_{\text{loc}}^{\alpha}(\mathbb{R}^d). \quad (4.15)
$$

In addition, provided $\alpha \ge \frac{2d}{d+2}$, these convergences are almost surely strong in $H_{\text{loc}}^{-1}(\mathbb{R}^d)$.

*Proof.* We split the proof into two steps.

*Step 1.* Proof of weak convergences (4.15).

As $\tilde{q}_E$ and $\tilde{\Sigma}_E$ are stationary, the ergodic theorem implies almost surely the weak convergences $\tilde{q}_E(\cdot|_E) \to \mathbb{E}[\tilde{q}_E]$ and $\tilde{\Sigma}_E(\cdot|_E) \to \mathbb{E}[\tilde{\Sigma}_E]$ in $L_{\text{loc}}^{\alpha}(\mathbb{R}^d)$, and it remains to compute these two expectations. For that purpose, up to an approximation argument as in the proof of Theorem 1, we may assume without loss of generality $\alpha > \frac{2d}{d+2}$. We split the proof into two further substeps.

*Substep 1.1.* Proof that $\bar{B}E = \mathbb{E}[\tilde{q}_E]$.

For all $R \ge 1$, we set $\chi_R := R^{-d} \chi(\frac{1}{R})$, for some smooth averaging function $\chi \in$
$C_c^\infty(\mathbb{R}^d; \mathbb{R}^+)$ such that $\chi$ is constant in $B$, vanishes outside $B_2$, and satisfies $\int_{\mathbb{R}^d} \chi = 1$. Given $E' \in \mathbb{M}_0^{\text{sym}}$, as $q_E = \text{D}(\psi_E) + E$ is stationary, definition (2.7) of $\bar{B}$ and the ergodic theorem yield almost surely,

$$E' : \bar{B} E = \mathbb{E}[q_{E'} : q_E] = \lim_{R \uparrow \infty} \int_{\mathbb{R}^d} \chi_R q_{E'} : q_E. \quad (4.16)$$

Since $q_{E'}$ vanishes in $\mathcal{I}$, cf. (2.8), and since $q_E$ coincides with its extension $\tilde{q}_E$ in $\mathbb{R}^d \setminus \mathcal{I}$, we find

$$q_{E'} : q_E = q_{E'} : \tilde{q}_E = E' : \tilde{q}_E + \text{D}(\psi_{E'}) : \tilde{q}_E.$$

Inserting this identity into (4.16), and noting that the ergodic theorem implies the almost sure convergence $\int_{\mathbb{R}^d} \chi_R \tilde{q}_E \to \mathbb{E}[\tilde{q}_E]$, we find

$$E' : \bar{B} E = E' : \mathbb{E}[\tilde{q}_E] + \lim_{R \uparrow \infty} \int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E. \quad (4.17)$$

In order to prove the claim $\bar{B} E = \mathbb{E}[\tilde{q}_E]$, it remains to show that the last limit vanishes:

$$\lim_{R \uparrow \infty} \int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E = 0. \quad (4.18)$$

Integrating by parts, using properties (4.11) of the extensions $(\tilde{q}_E, \tilde{\Sigma}_E)$, and using the constraint $\text{div}(\psi_{E'}) = 0$, we find

$$\begin{aligned}
\int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E &= \int_{\mathbb{R}^d} \text{D}(\chi_R \psi_{E'}) : \tilde{q}_E - \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_{E'}) : \tilde{q}_E \\
&= \frac{1}{2} \int_{\mathbb{R}^d} \tilde{\Sigma}_E \text{div}(\chi_R \psi_{E'}) - \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_{E'}) : \tilde{q}_E \\
&= -\frac{1}{2} \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_{E'}) : (2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}).
\end{aligned}$$

The relation $\text{div}(\psi_{E'}) = 0$ entails $\int_{\mathbb{R}^d} \nabla \chi_R \cdot \psi_{E'} = 0$, which allows us to add any constant to the pressure $\tilde{\Sigma}_E$ in the right-hand side. In view of the properties of the averaging function $\chi_R$, Hölder's inequality leads to

$$\begin{aligned}
\left| \int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E \right| & \lesssim \int_{B_2} \left| \frac{1}{R} \psi_{E'}(R \cdot) \right| \left( \left| \tilde{q}_E(R \cdot) \right| + \left| \tilde{\Sigma}_E(R \cdot) - \int_{B_2} \tilde{\Sigma}_E(R \cdot) \right| \right) \quad (4.19) \\
& \lesssim \left\| \frac{1}{R} \psi_{E'}(R \cdot) \right\|_{L^{\alpha'}(B_2)} \left( \left\| \tilde{q}_E(R \cdot) \right\|_{L^{\alpha}(B_2)} + \left\| \tilde{\Sigma}_E(R \cdot) - \int_{B_2} \tilde{\Sigma}_E(R \cdot) \right\|_{L^{\alpha}(B_2)} \right).
\end{aligned}$$

As the choice $\alpha > \frac{2d}{d+2}$ entails $\alpha' < \frac{2d}{d-2}$, we can use the sublinearity of $\psi_{E'}$ in $L^{\alpha'}$, cf. (2.9), together with the boundedness of $\{(\tilde{q}_E, \tilde{\Sigma}_E(R \cdot) - f_{B_2} \tilde{\Sigma}_E(R \cdot))\}_R$ in $L^\alpha(B_2)$, cf. Corollary 4.1 (i), and claim (4.18) follows.
$C_c^\infty(\mathbb{R}^d; \mathbb{R}^+)$ such that $\chi$ is constant in $B$, vanishes outside $B_2$, and satisfies $\int_{\mathbb{R}^d} \chi = 1$. Given $E' \in \mathbb{M}_0^{\text{sym}}$, as $q_E = \text{D}(\psi_E) + E$ is stationary, definition (2.7) of $\bar{B}$ and the ergodic theorem yield almost surely,

$$E' : \bar{B} E = \mathbb{E}[q_{E'} : q_E] = \lim_{R \uparrow \infty} \int_{\mathbb{R}^d} \chi_R q_{E'} : q_E. \quad (4.16)$$

Since $q_{E'}$ vanishes in $\mathcal{I}$, cf. (2.8), and since $q_E$ coincides with its extension $\tilde{q}_E$ in $\mathbb{R}^d \setminus \mathcal{I}$, we find

$$q_{E'} : q_E = q_{E'} : \tilde{q}_E = E' : \tilde{q}_E + \text{D}(\psi_{E'}) : \tilde{q}_E.$$

Inserting this identity into (4.16), and noting that the ergodic theorem implies the almost sure convergence $\int_{\mathbb{R}^d} \chi_R \tilde{q}_E \to \mathbb{E}[\tilde{q}_E]$, we find

$$E' : \bar{B} E = E' : \mathbb{E}[\tilde{q}_E] + \lim_{R \uparrow \infty} \int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E. \quad (4.17)$$

In order to prove the claim $\bar{B} E = \mathbb{E}[\tilde{q}_E]$, it remains to show that the last limit vanishes:

$$\lim_{R \uparrow \infty} \int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E = 0. \quad (4.18)$$

Integrating by parts, using properties (4.11) of the extensions $(\tilde{q}_E, \tilde{\Sigma}_E)$, and using the constraint $\text{div}(\psi_{E'}) = 0$, we find

$$\begin{align*}
\int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E &= \int_{\mathbb{R}^d} \text{D}(\chi_R \psi_{E'}) : \tilde{q}_E - \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_{E'}) : \tilde{q}_E \\
&= \frac{1}{2} \int_{\mathbb{R}^d} \tilde{\Sigma}_E \text{div}(\chi_R \psi_{E'}) - \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_{E'}) : \tilde{q}_E \\
&= -\frac{1}{2} \int_{\mathbb{R}^d} (\nabla \chi_R \otimes \psi_{E'}) : (2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}).
\end{align*}$$

The relation $\text{div}(\psi_{E'}) = 0$ entails $\int_{\mathbb{R}^d} \nabla \chi_R \cdot \psi_{E'} = 0$, which allows us to add any constant to the pressure $\tilde{\Sigma}_E$ in the right-hand side. In view of the properties of the averaging function $\chi_R$, Hölder's inequality leads to

$$\begin{align}
\left| \int_{\mathbb{R}^d} \chi_R \text{D}(\psi_{E'}) : \tilde{q}_E \right| & \lesssim \int_{B_2} \left| \frac{1}{R} \psi_{E'}(R \cdot) \right| \left( |\tilde{q}_E(R \cdot)| + \left| \tilde{\Sigma}_E(R \cdot) - \int_{B_2} \tilde{\Sigma}_E(R \cdot) \right| \right) \tag{4.19} \\
& \lesssim \left\| \frac{1}{R} \psi_{E'}(R \cdot) \right\|_{L^{\alpha'}(B_2)} \left( \|\tilde{q}_E(R \cdot)\|_{L^{\alpha}(B_2)} + \left\| \tilde{\Sigma}_E(R \cdot) - \int_{B_2} \tilde{\Sigma}_E(R \cdot) \right\|_{L^{\alpha}(B_2)} \right). \nonumber
\end{align}$$

As the choice $\alpha > \frac{2d}{d+2}$ entails $\alpha' < \frac{2d}{d-2}$, we can use the sublinearity of $\psi_{E'}$ in $L^{\alpha'}$, cf. (2.9), together with the boundedness of $\{ (\tilde{q}_E, \tilde{\Sigma}_E(R \cdot) - \int_{B_2} \tilde{\Sigma}_E(R \cdot)) \}_R$ in $L^\alpha(B_2)$, cf. Corollary 4.1 (i), and claim (4.18) follows.
Substep 1.2. Proof that $\tilde{\mathbf{b}}: E = -\mathbb{E}[\tilde{\Sigma}_E]$.

In terms of the cutoff function $w_n$ that we constructed in Lemma 3.3, integrating by parts, and recalling that the corrector equation (2.8) yields $\text{div}(\sigma(\psi_E + Ex, \Sigma_E)) = 0$ in $I_n^+ \setminus I_n$, the definition (2.12) of $\tilde{\mathbf{b}}$ becomes

$$ 
\begin{aligned}
\tilde{\mathbf{b}}: E &= \frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{\partial I_n} (x-x_n) \cdot \sigma(\psi_E + Ex, \Sigma_E) v \right] \\
&= -\frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n^+ \setminus I_n} \text{div}(w_n \sigma(\psi_E + Ex, \Sigma_E)(x-x_n)) \right] \\
&= -\frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n^+ \setminus I_n} D((x-x_n)w_n) : \sigma(\psi_E + Ex, \Sigma_E) \right].
\end{aligned}
$$

Writing $\sigma(\psi_E + Ex, \Sigma_E) = 2q_E - \Sigma_E \text{Id in } I_n^+ \setminus I_n$, and using the extensions $\tilde{q}_E$ and $\tilde{\Sigma}_E$ as in (4.11), we are led to

$$ \tilde{\mathbf{b}}: E = \frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n} D((x-x_n)w_n) : (2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}) \right]. $$

Since $D((x-x_n)w_n) = \text{Id in } I_n$ and since $\text{tr}(\tilde{q}_E) = 0$, we deduce

$$ \tilde{\mathbf{b}}: E = -\mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n} \tilde{\Sigma}_E \right], $$

and the claim $\tilde{\mathbf{b}}: E = -\mathbb{E}[\tilde{\Sigma}_E]$ easily follows by stationarity.

Step 2. Proof of strong convergences in $H_{\text{loc}}^{-1}(\mathbb{R}^d)$.

For $\alpha > \frac{2d}{d+2}$, strong convergences in $H_{\text{loc}}^{-1}(\mathbb{R}^d)$ follow from (4.15) and the compact Rellich embedding. It remains to consider the critical case $\alpha = \frac{2d}{d+2}$, for which we appeal to a two-scale argument inspired by [21, Lemma 1.15]. By stationarity, it suffices to prove $\tilde{q}_E(\cdot|_\epsilon) \to \bar{\mathbf{B}}E$ and $\tilde{\Sigma}_E(\cdot|_\epsilon) \to -\tilde{\mathbf{b}}: E$ strongly in $H^{-1}(B)$ almost surely as $\epsilon \downarrow 0$. As the argument is the same for $\tilde{q}_E$ and for $\tilde{\Sigma}_E$, we may focus on the former.

Let $h \in H^1(B)$ be momentarily fixed. Given $\eta > 0$, we choose a partition $\{Q_i\}_i$ of $B$ into Lipschitz subsets with $|Q_i| \simeq \eta^d$. In these terms, we can decompose

$$ 
\begin{aligned}
\int_B h(\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E) &= \sum_i \left( \int_{Q_i} h \right) \int_{Q_i} (\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E) \\
&\quad + \sum_i \int_{Q_i} \left( h - \int_{Q_i} h \right) (\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E). \quad (4.20)
\end{aligned}
$$

On the one hand, for all $s \in (1, \infty)$, noting that $1_{Q_i}$ belongs to $W^{\frac{1}{s}, s}(B)$, we can bound

$$ \left| \int_{Q_i} (\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E) \right| \le \left\| \frac{1}{|Q_i|} \right\|_{W^{\frac{1}{s}, s}(B)} \left\| \tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E \right\|_{W^{-\frac{1}{s}, s'}(B)}, $$
Substep 1.2. Proof that $\bar{\mathbf{b}}: E = -\mathbb{E}[\tilde{\Sigma}_E]$.

In terms of the cutoff function $w_n$ that we constructed in Lemma 3.3, integrating by parts, and recalling that the corrector equation (2.8) yields $\text{div}(\sigma(\psi_E + Ex, \Sigma_E)) = 0$ in $I_n^+ \setminus I_n$, the definition (2.12) of $\bar{\mathbf{b}}$ becomes

$$
\begin{align*}
\bar{\mathbf{b}}: E &= \frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{\partial I_n} (x-x_n) \cdot \sigma(\psi_E + Ex, \Sigma_E) v \right] \\
&= -\frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n^+ \setminus I_n} \text{div}(w_n \sigma(\psi_E + Ex, \Sigma_E)(x-x_n)) \right] \\
&= -\frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n^+ \setminus I_n} D((x-x_n)w_n) : \sigma(\psi_E + Ex, \Sigma_E) \right].
\end{align*}
$$

Writing $\sigma(\psi_E + Ex, \Sigma_E) = 2q_E - \Sigma_E \text{Id in } I_n^+ \setminus I_n$, and using the extensions $\tilde{q}_E$ and $\tilde{\Sigma}_E$ as in (4.11), we are led to

$$
\bar{\mathbf{b}}: E = \frac{1}{d} \mathbb{E} \left[ \sum_n \frac{1}{|I_n|} \int_{I_n} D((x-x_n)w_n) : (2\tilde{q}_E - \tilde{\Sigma}_E \text{Id}) \right].
$$

Since $D((x-x_n)w_n) = \text{Id in } I_n$ and since $\text{tr}(\tilde{q}_E) = 0$, we deduce

$$
\bar{\mathbf{b}}: E = -\mathbb{E}\left[\sum_n \frac{1}{|I_n|} \int_{I_n} \tilde{\Sigma}_E\right],
$$

and the claim $\bar{\mathbf{b}}: E = -\mathbb{E}[\tilde{\Sigma}_E]$ easily follows by stationarity.

Step 2. Proof of strong convergences in $H_{\text{loc}}^{-1}(\mathbb{R}^d)$.

For $\alpha > \frac{2d}{d+2}$, strong convergences in $H_{\text{loc}}^{-1}(\mathbb{R}^d)$ follow from (4.15) and the compact Rel-
lich embedding. It remains to consider the critical case $\alpha = \frac{2d}{d+2}$, for which we appeal to
a two-scale argument inspired by [21, Lemma 1.15]. By stationarity, it suffices to prove
$\tilde{q}_E(\cdot) \to \bar{\mathbf{B}}E$ and $\tilde{\Sigma}_E(\cdot) \to -\bar{\mathbf{b}}: E$ strongly in $H^{-1}(B)$ almost surely as $\varepsilon \downarrow 0$. As the
argument is the same for $\tilde{q}_E$ and for $\tilde{\Sigma}_E$, we may focus on the former.

Let $h \in H^1(B)$ be momentarily fixed. Given $\eta > 0$, we choose a partition $\{Q_i\}_i$ of B
into Lipschitz subsets with $|Q_i| \simeq \eta^d$. In these terms, we can decompose

$$
\int_B h(\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E) = \sum_i \left( \int_{Q_i} h \right) \int_{Q_i} (\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E) \\
+ \sum_i \int_{Q_i} \left( h - \int_{Q_i} h \right) (\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E). \quad (4.20)
$$

On the one hand, for all $s \in (1, \infty)$, noting that $1_{Q_i}$ belongs to $W^{\frac{1}{s}, s}(B)$, we can bound

$$
\left| \int_{Q_i} (\tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E) \right| \le \left\| \frac{1}{|Q_i|} \right\|_{W^{\frac{1}{s},s}(B)} \left\| \tilde{q}_E(\cdot|_\epsilon) - \bar{\mathbf{B}}E \right\|_{W^{-\frac{1}{s},s'}(B)},
$$
and thus, further using the Sobolev embedding in the form of $\|h\|_{L^1(B)} \lesssim \|h\|_{H^1(B)}$, we deduce for the first right-hand-side term in (4.20),

$$ \left| \sum_i \left( \int_{Q_i} h \right) \int_{Q_i} \left( \tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E \right) \right| \\ \lesssim \|h\|_{H^1(B)} \left( \sup_i \left\| \frac{1}{|Q_i|} \right\|_{W^{\frac{1}{s},s}(B)} \right) \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{W^{-\frac{1}{s},s'}(B)}. \quad (4.21) $$

On the other hand, using Hölder's inequality and the Poincaré–Sobolev embedding, the second right-hand-side term in (4.20) can be estimated as

$$ \begin{aligned} \left| \sum_i \int_{Q_i} \left(h - \int_{Q_i} h\right) \left(\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\right) \right| & \le \sum_i \left\|h - \int_{Q_i} h\right\|_{L^{\frac{2d}{d-2}}(Q_i)} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)} \\ & \lesssim \sum_i \|\nabla h\|_{L^2(Q_i)} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)} \\ & \le \|\nabla h\|_{L^2(B)} \left( \sum_i \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)}^2 \right)^{\frac{1}{2}}. \end{aligned} $$

Combining this with (4.20) and (4.21), and taking the supremum over test functions $h \in H^1(B)$, we conclude for all $s \in (1, \infty)$,

$$ \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{H^{-1}(B)} \lesssim \left( \sup_i \left\| \frac{1}{|Q_i|} \right\|_{W^{\frac{1}{s},s}(B)} \right) \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{W^{-\frac{1}{s},s'}(B)} \\ + \left( \sum_i \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)}^2 \right)^{\frac{1}{2}}. $$

We now pass to the limit $\epsilon \downarrow 0$ in this estimate. Choosing $2^{\frac{d-1}{d-2}} < s < \infty$, the compact Rellich embedding ensures that $L^{\frac{2d}{d+2}}(B)$ is compactly embedded in $W^{-\frac{1}{s},s'}(B)$. Therefore, in view of (4.15) with $\alpha = \frac{2d}{d+2}$, we deduce $\tilde{q}_E(\frac{\cdot}{\epsilon}) \to \bar{B}E$ strongly in $W^{-\frac{1}{s},s'}(B)$ almost surely as $\epsilon \downarrow 0$. Further using the stationarity and the boundedness of $\tilde{q}_E(\frac{\cdot}{\epsilon}) - 2\bar{B}E$ in $L^{\frac{2d}{d+2}}(\Omega)$, cf. Corollary 4.1 (i), we get almost surely,

$$ \limsup_{\epsilon \downarrow 0} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{H^{-1}(B)} \lesssim \left( \sum_i |Q_i|^{\frac{d+2}{d}} \right)^{\frac{1}{2}}. $$

Using that $\sum_i |Q_i| \lesssim 1$ and $|Q_i| \lesssim \eta^d$, this turns into

$$ \limsup_{\epsilon \downarrow 0} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{H^{-1}(B)} \lesssim \eta. $$

Finally, letting the mesh $\eta$ of the partition $\{Q_i\}_i$ tend to 0, the conclusion follows. ■
and thus, further using the Sobolev embedding in the form of $\|h\|_{L^1(B)} \lesssim \|h\|_{H^1(B)}$, we deduce for the first right-hand-side term in (4.20),

$$ \left| \sum_i \left( \int_{Q_i} h \right) \int_{Q_i} \left( \tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E \right) \right| \\ \lesssim \|h\|_{H^1(B)} \left( \sup_i \left\| \frac{1}{|Q_i|} \right\|_{W^{\frac{1}{s},s}(B)} \right) \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{W^{-\frac{1}{s},s'}(B)}. \quad (4.21) $$

On the other hand, using Hölder's inequality and the Poincaré–Sobolev embedding, the second right-hand-side term in (4.20) can be estimated as

$$ \begin{aligned} \left| \sum_i \int_{Q_i} \left(h - \int_{Q_i} h\right) \left(\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\right) \right| & \le \sum_i \left\|h - \int_{Q_i} h\right\|_{L^{\frac{2d}{d-2}}(Q_i)} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)} \\ & \lesssim \sum_i \|\nabla h\|_{L^2(Q_i)} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)} \\ & \le \|\nabla h\|_{L^2(B)} \left( \sum_i \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)}^2 \right)^{\frac{1}{2}}. \end{aligned} $$

Combining this with (4.20) and (4.21), and taking the supremum over test functions $h \in H^1(B)$, we conclude for all $s \in (1, \infty)$,

$$ \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{H^{-1}(B)} \lesssim \left( \sup_i \left\| \frac{1}{|Q_i|} \right\|_{W^{\frac{1}{s},s}(B)} \right) \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{W^{-\frac{1}{s},s'}(B)} \\ + \left( \sum_i \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{L^{\frac{2d}{d+2}}(Q_i)}^2 \right)^{\frac{1}{2}}. $$

We now pass to the limit $\epsilon \downarrow 0$ in this estimate. Choosing $2^{\frac{d-1}{d-2}} < s < \infty$, the compact Rellich embedding ensures that $L^{\frac{2d}{d+2}}(B)$ is compactly embedded in $W^{-\frac{1}{s},s'}(B)$. Therefore, in view of (4.15) with $\alpha = \frac{2d}{d+2}$, we deduce $\tilde{q}_E(\frac{\cdot}{\epsilon}) \to \bar{B}E$ strongly in $W^{-\frac{1}{s},s'}(B)$ almost surely as $\epsilon \downarrow 0$. Further using the stationarity and the boundedness of $\tilde{q}_E(\frac{\cdot}{\epsilon}) - 2\bar{B}E$ in $L^{\frac{2d}{d+2}}(\Omega)$, cf. Corollary 4.1 (i), we get almost surely,

$$ \limsup_{\epsilon \downarrow 0} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{H^{-1}(B)} \lesssim \left( \sum_i |Q_i|^{\frac{d+2}{d}} \right)^{\frac{1}{2}}. $$

Using that $\sum_i |Q_i| \lesssim 1$ and $|Q_i| \lesssim \eta^d$, this turns into

$$ \limsup_{\epsilon \downarrow 0} \|\tilde{q}_E(\frac{\cdot}{\epsilon}) - \bar{B}E\|_{H^{-1}(B)} \lesssim \eta. $$

Finally, letting the mesh $\eta$ of the partition $\{Q_i\}_i$ tend to 0, the conclusion follows. ■
### 4.3. Proof of Theorem 3

The moment condition (2.11) amounts to the following: for some $2 \le r \ne \frac{2d}{d-2}$ and $\frac{2d}{d+2} \le \alpha \le 2 \wedge \frac{2dr}{r(d-2)+2d}$, with $\alpha < \frac{d}{d-1}$ if $r=2$, the interparticle distances satisfy almost surely,

$$ \limsup_{\epsilon \downarrow 0} \epsilon^d \sum_{n \in \mathcal{N}_\epsilon(U)} \mu_r(\rho_{n;U,\epsilon})^{\frac{2\alpha}{2-\alpha}} < \infty. \quad (4.22) $$

We split the proof into two steps. First, we establish the convergence of the velocity field by a direct div-curl argument inspired by the work of Jikov ([24, 25]) on homogenization problems with stiff inclusions (see also [20, Section 3.2]), and then we turn to the convergence of the pressure.

Step 1. Div-curl argument: we prove that almost surely, as $\epsilon \downarrow 0$,

$$ \begin{aligned} & u_{\epsilon} \rightarrow \bar{u} \quad \text{weakly in } H_0^1(U), \\ & \tilde{p}_{\epsilon} \rightarrow \bar{B} \, \mathrm{D}(\bar{u}) \quad \text{weakly in } L^{\alpha}(U), \\ & \tilde{S}_{\epsilon} - \int_{U} \tilde{S}_{\epsilon} \rightarrow \bar{S} \quad \text{weakly in } L^{\alpha}(U), \end{aligned} \quad (4.23) $$

where $(\bar{u}, \bar{S})$ is the solution of the homogenized equation (2.14). By a standard energy argument as e.g. in [7, Step 8.1 of the proof of Proposition 2.1], provided that $f \in L^p(U)$ for some $p > d$, this weak convergence result easily implies the following corresponding corrector result, almost surely:

$$ \begin{aligned} & p_{\epsilon} - \sum_{E \in \mathcal{E}} q_E(\dot{\epsilon}) \nabla_E \bar{u} \rightarrow 0 \quad \text{strongly in } L^2(U), \\ & u_{\epsilon} - \bar{u} - \sum_{E \in \mathcal{E}} \epsilon \psi_E(\dot{\epsilon}) \nabla_E \bar{u} \rightarrow 0 \quad \text{strongly in } H_0^1(U), \end{aligned} \quad (4.24) $$

where we recall the shorthand notation $\nabla_E \bar{u} = E : D(\bar{u})$ and where $\mathcal{E}$ stands for an orthonormal basis of $M_0^{\text{sym}}$. We omit the proof of this standard consequence (4.24) and rather focus on the proof of (4.23).

For $\kappa > 0$ we set for abbreviation $U^\kappa := \{x \in U : \text{dist}(x, \partial U) > \kappa\}$. Since $q_E|_I = 0$ and $p_\epsilon|_{I_\epsilon(U)} = 0$, since $\tilde{q}_E$ and $q_E$ coincide on $\mathbb{R}^d \setminus I$, since $\tilde{p}_\epsilon$ and $p_\epsilon$ coincide on $U \setminus I_\epsilon(U)$, and since definition (2.10) of $I_\epsilon(U)$ entails $I_\epsilon(U) \cap U^\kappa = (\epsilon I) \cap U^\kappa$ whenever $\epsilon < \frac{\kappa}{2}$, we deduce the following identity on $U^\kappa$ for $\epsilon < \frac{\kappa}{2}$:

$$ \tilde{q}_E(\dot{\epsilon}) : p_\epsilon = q_E(\dot{\epsilon}) : \tilde{p}_\epsilon, \quad (4.25) $$

and we aim at passing to the limit in both sides. Since the energy bound entails that $(u_\epsilon)_\epsilon$ is almost surely bounded in $H_0^1(U)$, since Corollary 4.1 (ii) ensures that $(\tilde{p}_\epsilon, \tilde{S}_\epsilon)_\epsilon$ is almost surely bounded in $L^\alpha(U)$, further recalling (2.9) and Lemma 4.3, we find almost surely,
### 4.3. Proof of Theorem 3

The moment condition (2.11) amounts to the following: for some $2 \le r \ne \frac{2d}{d-2}$ and $\frac{2d}{d+2} \le \alpha \le 2 \wedge \frac{2dr}{r(d-2)+2d}$, with $\alpha < \frac{d}{d-1}$ if $r=2$, the interparticle distances satisfy almost surely,

$$ \limsup_{\epsilon \downarrow 0} \epsilon^d \sum_{n \in \mathcal{N}_\epsilon(U)} \mu_r(\rho_{n;U,\epsilon})^{\frac{2\alpha}{2-\alpha}} < \infty. \quad (4.22) $$

We split the proof into two steps. First, we establish the convergence of the velocity field by a direct div-curl argument inspired by the work of Jikov ([24, 25]) on homogenization problems with stiff inclusions (see also [20, Section 3.2]), and then we turn to the convergence of the pressure.

Step 1. Div-curl argument: we prove that almost surely, as $\epsilon \downarrow 0$,

$$ \begin{align} u_{\epsilon} &\to \bar{u} \quad \text{weakly in } H_0^1(U), \nonumber \\ \tilde{p}_{\epsilon} &\to \bar{B} \, \mathrm{D}(\bar{u}) \quad \text{weakly in } L^{\alpha}(U), \tag{4.23} \\ \tilde{S}_{\epsilon} - \int_{U} \tilde{S}_{\epsilon} &\to \bar{S} \quad \text{weakly in } L^{\alpha}(U), \nonumber \end{align} $$

where $(\bar{u}, \bar{S})$ is the solution of the homogenized equation (2.14). By a standard energy argument as e.g. in [7, Step 8.1 of the proof of Proposition 2.1], provided that $f \in L^p(U)$ for some $p > d$, this weak convergence result easily implies the following corresponding corrector result, almost surely:

$$ \begin{align} p_{\epsilon} - \sum_{E \in \mathcal{E}} q_E(\dot{\epsilon}) \nabla_E \bar{u} &\to 0 \quad \text{strongly in } L^2(U), \tag{4.24} \\ u_{\epsilon} - \bar{u} - \sum_{E \in \mathcal{E}} \epsilon \psi_E(\dot{\epsilon}) \nabla_E \bar{u} &\to 0 \quad \text{strongly in } H_0^1(U), \nonumber \end{align} $$

where we recall the shorthand notation $\nabla_E \bar{u} = E : D(\bar{u})$ and where $\mathcal{E}$ stands for an orthonormal basis of $M_0^{\text{sym}}$. We omit the proof of this standard consequence (4.24) and rather focus on the proof of (4.23).

For $\kappa > 0$ we set for abbreviation $U^\kappa := \{x \in U : \text{dist}(x, \partial U) > \kappa\}$. Since $q_E|_I = 0$ and $p_\epsilon|_{I_\epsilon(U)} = 0$, since $\tilde{q}_E$ and $q_E$ coincide on $\mathbb{R}^d \setminus I$, since $\tilde{p}_\epsilon$ and $p_\epsilon$ coincide on $U \setminus I_\epsilon(U)$, and since definition (2.10) of $I_\epsilon(U)$ entails $I_\epsilon(U) \cap U^\kappa = (\epsilon I) \cap U^\kappa$ whenever $\epsilon < \frac{\kappa}{2}$, we deduce the following identity on $U^\kappa$ for $\epsilon < \frac{\kappa}{2}$:

$$ \tilde{q}_E(\dot{\epsilon}) : p_\epsilon = q_E(\dot{\epsilon}) : \tilde{p}_\epsilon, \quad (4.25) $$

and we aim at passing to the limit in both sides. Since the energy bound entails that $(u_\epsilon)_\epsilon$ is almost surely bounded in $H_0^1(U)$, since Corollary 4.1 (ii) ensures that $(\tilde{p}_\epsilon, \tilde{S}_\epsilon)_\epsilon$ is almost surely bounded in $L^\alpha(U)$, further recalling (2.9) and Lemma 4.3, we find almost surely,