Derivation of Darcy's law in randomly punctured domains

We consider the homogenization of a Poisson problem or a Stokes system in a randomly punctured domain with Dirichlet boundary conditions. We assume that the holes are spherical and have random centres and radii. We impose that the average distance between the balls is of size $\eps$ and their average radius is $\varepsilon^{\alpha}$, $\alpha \in (1; 3)$. We prove that, as in the periodic case [G. Allaire, ``Homogenization of the Navier-Stokes equations in domains perforated with tiny holes. II''], the solutions converge to the solution of Darcy's law (or its scalar analogue in the case of Poisson). In the same spirit of [A. Giunti, R. H\"ofer and J. Vel\'azquez, ``Homogenization of the Poisson equation in randomly perforated domains under minimal assumptions on the size of the holes''], we work under minimal conditions on the integrability of the random radii. These ensure that the problem is well-defined but do not rule out the onset of clusters of holes.

B ε α ρz (εz). (0.1) Here, the set of centres Φ is a Poisson point process of intensity λ > 0 and the set 1 ε D := {x ∈ R 3 : εx ∈ D}. The radii R = {ρ z } z∈Φ ⊆ [1; +∞) are independent and identically distributed random variables satisfying for a constant C < +∞ This condition is minimal in order to ensure that, P-almost surely and when ε is small, the set H ε does not fully cover the domain D, hence implying that D ε = ∅ (see Lemma 1.1). However, condition (0.2) does not prevent that, with high probability, the balls in H ε do overlap.
For ε > 0 and D ε as above, we consider the (weak) solution to either In the case of the Stokes system, we further assume that We refer to the next section for a more detailed discussion on what conditions (0.2) and (0.5) entail in terms of the geometric properties of the set H ε .
It is easy to see that in the case of spherical periodic holes having distance ε and radius ε α , α ∈ (1, 3], the density of harmonic capacity of H ε is asymptotically of order ε −3+α ; The same is true in the case of the Stokes capacity. When α = 3 these limits are thus finite. In the case of the Poisson problem, the solutions to (0.3) thus converge to the solution u ∈ H 1 0 (D) to −∆u + µu = f in D, where the constant µ > 0 is the limit of the capacity density [7]. Similarly, the limit problem for (0.4) is given by a Brinkmann system, namely a Stokes system in D with no-slip boundary conditions and with the additional termμu in the system of equations [1]. The termμ > 0 is as well strictly related to the limit of the Stokes capacity density. We also mention that, for holes that are periodic but not spherical, the termμ is a positive-definite matrix. For α ∈ (1; 3) as in the present paper, the solutions to (0.3) or (0.4) need to be rescaled by the factor ε −3+α in order to converge to a non-trivial limit. The effective equations, in this case, are either u = kf in D or Darcy's law u = K(f − ∇p) in D [2]. Here, k, K are related to the rescaled limit of the density of capacity and admit a representation in terms of a corrector problem solved in the exterior domain R 3 \B 1 (0).
When α = 1, namely when the distance between holes and their size have the same order ε, the effective equations for (0.3) and (0.4) are as in the case α ∈ (1, 3); the effective constants k, K obtained in the limit, however, are determined by a corrector problem of different nature. In this case indeed, there is only one microscopic scale ε and the relative distance between the connected components of the holes H ε does not tends to infinity for ε → 0. This yields that the corrector equations are solved in the periodic cell and not in the exterior domain R 3 \B 1 (0) [3].
For holes that are not periodic, the extremal regimes α ∈ {1, 3} have been rigorously studied both in deterministic and random settings. For α = 3 we mention, for instance [6,9,16,17,18,21,22,23] and refer to the introductions in [12] and [14] for a detailed overview of these results. We stress that the homogenization of (0.3) and (0.4) when H ε is as in (0.1) with α = 3 has been studied in the series of papers [12,13,14]. These works prove the convergence to the effective equation under the minimal assumption that H ε has finite averaged capacity density. There is no additional condition on the minimal distance between the balls in the set of H ε .
There are many works devoted also to the regime α = 1. We refer, in particular, to [4] where (0.3) and (0.4) are studied for a very general class of stationary and ergodic punctured domains. For these domains, the formulation of the corrector equation for the the effective quantities k, K is solved in the probability space (Ω, F, P) generating the holes.
There is fewer mathematical literature concerning the homogenization of (0.3) or (0.4) in the regime α ∈ (1; 3). For periodic holes, this has been studied in [2]. These results have been extended for certain regimes to compressible Navier-Stokes systems [15] or to elliptic systems in the context of linear elasticity [19]. We are not aware of analogous results when the holes H ε are not periodic. The present paper considers this problem when H ε is random and, in the same spirit of [12,14], allows that the balls in H ε overlap and cluster.
The main result of this paper is the following: Theorem 0.1. Let σ ε := ε − 3−α 2 and let H ε and D ε be the random sets defined in (0.1).
Here, and in the rest of the paper, E · denotes the expectation under the probability measure for (Φ, R).
As mentioned above, condition (0.2) is minimal in order to ensure that the set D ε is non-empty for P-almost every realization. A lower stochastic integrability assumption for the radii, indeed, yields that, in the limit ε ↓ 0 and P-almost surely H ε , covers the full set D (see Lemma 1.1 in the next section). By the Strong Law of the Large Numbers, condition (0.2) implies that the density of capacity is almost surely of order ε −3+α as in the periodic case. As already remarked in [12] in the case α = 3, with (0.4) we require that the radii satisfy the slightly stronger assumption (0.5). While (0.2) seems to be the optimal condition in order to control the density of harmonic capacity, the lack of subadditivity of the Stokes capacity calls for a better control on the geometry of the set H ε .
The ideas used in the proof of Theorem 0.1 are an adaptation of the techniques used in [2,7] for the periodic case. They are combined with the tools developed in [12,14] to tackle the case of domains having holes that may overlap. As shown in [2], the uniform bounds on the sequences {σ 2 ε u ε } ε>0 , {σ ε ∇u ε } ε>0 are obtained by means of a Poincaré's inequality for functions that vanish on ∂D ε . If v ∈ H 1 (D ε ), since the function vanishes on the holes H ε , the constant in the Poincaré' s inequality is of order σ −1 ε << 1.
, this would instead be of order 1 (dependent on the domain D). Note that, as for α = 3 we have σ ε = 1, there is no gain in using a Poincar'e's inequality in H 1 0 (D ε ) instead of in H 1 0 (D) in this regime. In the case of centres of H ε that are distributed like a Poisson point process, the is a low probability that some regions of D ε have few holes, thus leading to a worse Poincaré's constant. This causes the lack of uniform bounds for the family {σ 2 ε u ε } ε>0 in L 2 (D).
Equipped with uniform bounds for the rescaled solutions of (0.3), one may prove Theorem 0.1, (a) by constructing suitable oscillating test functions {w ε } ε>0 . These allow to pass to the limit in the equation and identify the effective problem. We stress that a crucial ingredient in these arguments is given by the quantitative bounds obtained in [11] in the case α = 3. These bounds may indeed also be extended to the current setting sot that the rate of convergence of the measures −σ −2 ε ∆w ε ∈ H −1 (D) is quantified. This allows to control the convergence of the duality term −∆w ε ; u ε H −1 (D);H 1 0 (D) . There is a fine balance the convergence of −σ 2 ε ∆w ε with the right space where we have uniform bounds for {σ 2 ε u ε } ε>0 . In contrast with the periodic case, the unboundedness of {σ 2 ε u ε } ε>0 in L 2 (D) requires for a careful study of the duality term above. For the precise statements, we refer to (3.3) in Lemma 3.1 and Lemma 3.3. The same ideas sketched here apply also to the case of solutions to (0.4). This time, the oscillating test functions {w ε } ε>0 are replaced by the reduction operator R ε of Lemma 4.1.
Remark 0.2. We comment below on some variations and corollaries of Theorem 0.1: then the convergence of Theorem 0.1 holds also with p = 2. In this case, indeed, we may drop the logarithmic factor in the bounds of Lemma 2.1.
(ii) A careful inspection of the proof of Theorem 0.1 yields that, under assumption (0.5) and for a source f ∈ W 1,∞ , the convergences in both (a) and (b) may be upgraded to for an exponent κ > 0 depending on α, β.
(iii) The quenched version of Theorem 0.1, namely the P-almost sure convergence of the families in L p (D), holds as well provided that we restrict to any vanishing sequence {ε j } j∈N that converges fast enough. For instance, it suffices that j 1 3 +ǫ ε j → 0, ǫ > 0. It is a technical but easy argument to observe that, under this assumption, limits (3.3) of Lemma 3.1 and (4.1)-(4.2) of Lemma 4.1 vanish also P-almost surely. From these, the quenched version of Theorem 0.1 may be shown as done in the annealed case. To control the limits in (3.3), (4.1) and (4.2) without taking the expectation, one may follow the same lines of the current proof and control most of the terms by the Strong Law of Large Numbers. Condition j 1 3 +ǫ ε j → 0 on the speed of the convergence for {ε j } j∈N is needed in order to obtain quenched bounds for the term in (3.31) by means of Borel-Cantelli's Lemma. The paper is structured as follows: In the next section we describe the setting and introduce the notation that we use throughout the proofs. Subsection 1.2 is devoted to discussing the minimality of assumption (0.2) and what condition (0.5) implies on the geometry of the holes H ε . In Section 2, we show the uniform bounds on the family {σ 2 ε u ε } ε>0 , with u ε solving (0.3) or (0.4). In Section 3 we argue Theorem 0.1 in case (a), while in Section 4 we adapt it to case (b). The proof of case (b) is conceptually similar to the one for (a), but it is technically more challenging. It heavily relies on the geometric properties of the holes implied by condition (0.5). Finally, Section 5 contains the proof of the main auxiliary results used throughout the paper.

Setting and notation
Let D ⊆ R 3 be an open set having C 1,1 -boundary. We assume that D is star-shaped with respect to a point x 0 ∈ R 3 . This assumption is purely technical and allows us to give an easier formulation for the set of holes H ε . With no loss of generality we assume that x 0 = 0. The process (Φ; R) is a stationary marked point process on R 3 having identically and independent distributed marks on [1 + ∞). In other words, (Φ; R) may be seen as a Poisson point process on the space R 3 × [1; +∞), having intensityλ(x, ρ) = λf (ρ). The expectation in (0.2) or (0.5) is therefore taken with respect to the measure f (ρ)dρ. We denote by (Ω; F, P) the probability space associated to (Φ, R), so that the random sets in (0.1) and the random fields solving (0.3) or (0.4) may be written as H ε = H ε (ω), D ε = D ε (ω) and u ε (ω; ·), respectively. The set of realizations Ω may be seen as the set of atomic measures n∈N δ (zn,ρn) in R 3 × [1; +∞) or, equivalently, as the set of (unordered) collections We choose as F the smallest σ-algebra such that the random variables N (B) : Ω → N, ω → #{ω ∩ B} are measurable for every set B ⊆ R 4 the Borel σ-algebra B R 4 . Here and throughout the paper, # stands for the cardinality of the set considered. For every p ∈ [1; +∞) we define the space L p (Ω) as the space of (F-measurable) random variables F : Ω → R endowed with the norm E |F (ω)| p 1 p . For p = +∞, we set L ∞ (Ω) as the space of P-essentially bounded random variables. We denote by L p (Ω × D), p ∈ [1; +∞), the space of random fields F : Ω × R 3 → R that are measurable with respect to the product σ-algebra are separable for p ∈ [1, +∞) and reflexive for p ∈ (1, +∞) (see e.g. [5][ Section 13,4]). The same definition, with obvious modifications, holds in the case of the target space R replaced by R 3 .
We often appeal to the Strong Law of Large Numbers (SSLN) for averaged sums of the form where {X z } z∈Φ ε (D) are identically distributed random variables that have sufficiently decaying correlations. Here, we send the radius of the ball B R to infinity. It is well-known that such results hold and we refer to [14][Section 5] for a detailed proof of the result that is tailored to the current setting.

Notation
We use the notation or for C or C where the constant depends only on α, λ, D and, in case (b), also on β in (0.5). Given a parameter p ∈ R, we use the notation p if the implicit constant also depends on the value p. For r > 0, we write B r for the ball of radius r centred in the origin of R 3 . We denote by · ; · the duality bracket between the spaces H −1 (D) and H 1 0 (D). When no ambiguity occurs, we skip the argument ω ∈ Ω in all the random objects considered in the paper. If (Φ; R) is as in the previous subsection, for a set A ⊆ R d , we define

On the assumptions on the radii
In this subsection we discuss the choice of assumptions (0.2) and (0.5) in Theorem 0.1. We postpone to the Appendix the proofs of the statements. The next result states that assumption (0.2) is sufficient to have only microscopic holes whose size vanishes in the limit ε ↓ 0. Moreover, it is also necessary in order to have that holes H ε do not cover the full domain D.
Lemma 1.1. The following conditions are equivalent: (i) The process satisfies (0.2); (ii) For P-almost every realization and for every ε small enough the set D ε = ∅.
In the following result we provide the geometric information on H ε that may be inferred by strengthening condition (0.2) to (0.5). Roughly speaking, the next lemma tells that, under condition (0.5), we have a control on the maximum number of holes of comparable size that intersect. More precisely, we may discretize the range of the size of the radii {ρ z } z∈Φ ε (D) and partition the set of centres Φ ε (D) according to the order of magnitude of the associated radii. The next statement says that there exists an M ∈ N (that is independent from the realization ω ∈ Ω) such that, provided that the step-size of the previous discretization is small enough, each sub-collection contains at most M holes that overlap when dilated by a factor 4. This result allows to treat also the case of the Stokes system in Theorem 0.1, (b) and motivates the need of the stronger assumption (0.5) in that setting.
such that for P-almost every realization and for every ε small enough it holds and we may rewrite

Uniform bounds
In this section we provide uniform bounds for the family {σ 2 ε u ε } ε>0 and {σ ε ∇u ε } ε>0 . We stress that, as in [2], this is done by relying on a Poincaré's inequality for functions that vanish in the holes H ε . The order of magnitude of the typical size (i.e. ε α ) and distance (i.e. ε) of the holes yields that the Poincaré's constant scales as the factor σ ε introduced in Theorem 0.1. This, combined with the energy estimate for (0.3) or (0.4), allows to obtain the bounds on the rescaled solutions. We mention that the next results contain both annealed and quenched uniform bounds. The quenched versions are not needed to prove Theorem 0.1, but may be used to prove the quenched analogue described in Remark 0.2, (iii).

Lemma 2.1. Let u ε be is as in Theorem 0.1. Then for every
Furthermore, for P-almost every realization, the sequences This, in turn, is a consequence of Proof of Lemma 2.2. As first step, we argue that the following Poincaré's inequality holds: Let V be a convex domain. Assume that V ⊆ B r for some r > 0. Let s < r. Then, for every q ∈ [1; 2] and u ∈ H 1 (V \B s ) such that u = 0 on ∂B s it holds The proof of this result is standard and may be easily proven by writing the integrals in spherical coordinates. We stress that the assumptions on V allows to write the domain V \B s as As second step, we construct an appropriate random tesselation for D: We consider the Voronoi tesselation {V z } z∈Φ associated to the point process Φ, namely the sets Note that, by the previous rescaling, we have that, if diam(V z ) := r z , then diam(V ε,z ) = εr z .
It is immediate to see that, for every realization ω ∈ Ω, the sets {V ε,z } z∈A ε are essentially disjoint, convex and cover the set D. Since Φ is stationary, the random variables {r z } z∈Φ are identically distributed. Furthermore, they are distributed as a generalized Gamma distribution having intensity g(r) = C(λ)r 8 and that there exists a constant c = c(λ) > 0 such that for every function F : R + → R (that is integrable with respect to the measure g(r)dr) Since p ∈ [1, 2), we may appeal to Hölder's inequality and conclude that i.e. inequality (2.6). This concludes the proof of (2.2) in the case p ∈ [1; 2).
To tackle the case p = 2 we need a further manipulation: we distinguish between points z ∈ A ε having r z > − log ε or r z − log ε:ˆD (2.8) We apply Poincaré's inequality in H 1 0 (D) on every integral of the second sum above. This implies that so that Chebyschev's inequality and (2.5) yield where we set C ε (2) := ε 3 z∈A ε exp r 2 z . Note that, again by (2.4)-(2.5) and the Law of Large Numbers, this definition of C ε (2) satisfies (2.3). Inserting the previous display into (2.10) implies that We now apply Lemma 2 in the remaining sum and obtain (2.7) with p = 2, where the sum is restricted to the points z ∈ A ε such that r z − log ε. From this, we infer that (2.10) By redefining C ε (2) 2 = min ε 3 z∈A ε exp r 2 z ; 1 , the above inequality immediately implies (2.2) for p = 2. The proof of Lemma 2.2 is complete.
Proof of Lemma 2.1. We prove Lemma 2.1 for u ε solving (0.3). The case (0.4) is analogous. Since f ∈ L q (D) with q ∈ (2; +∞], we may test (0.3) with u ε and use Hölder's inequality to control We thus appeal to (2.2) with p = q q−1 and obtain that Thanks to (2.3) of Lemma 2.2, this yields that the sequence {σ ε ∇u ε } ε>0 is bounded in L 2 (D) for P-almost every realization. Similarly, we infer (2.1) by taking the expectation and applying Hölder's inequality.
We argue the remaining bounds for the terms of u ε in a similar way: We combine Lemma 2.2 with the same calculation above for (2.11) and apply Hölder's inequality. This establishes Lemma 2.1.

3
Proof of Theorem 0.1, (a) Lemma 3.1. There exists an ε 0 = ε 0 (d) such that for every ε < ε 0 and P-almost every realization there exists a family In addition, and for every φ ∈ C ∞ 0 (D) and v ε ∈ H 1 0 (D ε ) satisfying the bounds of Lemma 2.1 and such that σ 2 Here, the constant k is as in Theorem 0.1, (a).
Proof of Theorem 0.1, (a). The proof is similar to the one in [2]. We first show that σ 2 ε u ε ⇀ u in L p (D × Ω), p ∈ [1, 2). By the uniform bounds of Lemma 2.1, we have that, up to a subsequence, there exists a weak limit u * ∈ L p (Ω × R d ), p ∈ [1, 2). We prove that, P-almost surely, the function u * = kf in D. This, in particular, also implies that the full family {σ 2 ε u ε } ε>0 weakly converges to u * . We restrict to the converging subsequence {σ 2 ε j u ε j } j∈N . However, for the sake of a lean notation, we forget about the subsequence {ε j } j∈N and continue using the notation u ε and ε ↓ 0. Let ε 0 and {w ε } ε>0 be as in Lemma 3.1. For every ε < ε 0 , χ ∈ L ∞ (Ω) and φ ∈ C ∞ 0 (D) we test equation (0.3) with χw ε φ and take the expectation: Using Leibniz's rule, integration by parts and the bounds for u ε and w ε in Lemma 2.1 and 3.1 we reduce to We now appeal to (3.3) in Lemma 3.1 applied to the converging subsequence {u ε } ε>0 and conclude that Since χ ∈ L ∞ (Ω) and φ ∈ C ∞ 0 (D) are arbitrary, we infer that for P-almost every realization u * = kf for (Lebesgue-)almost every x ∈ D. We stress that in this last statement we used the separability of L p (D), p ∈ [1, ∞). This establishes that the full family To conclude Theorem 0.1, (a) it remains to upgrade the previous convergence from weak to strong. We fix p ∈ [1, 2). By the assumption on f , the function u * ∈ L q (D), for some q ∈ (2; +∞]. Let {u n } n∈N ⊆ C ∞ 0 (D) be an approximating sequence for u * in L q (D). Since w ε ∈ W 1,∞ (D), the function w ε u n ∈ H 1 0 (D). Hence, by Lemma 2.2 applied to u ε − w ε u n we obtain We claim that Provided this holds, we establish Theorem 0.1, (a), as follows: By the triangle inequality we have that Since u * and u n ∈ C ∞ 0 (D) are deterministic, we take the expectation and use Lemma 3.1 with (3.6) to get This implies the statement of Theorem 0.1, (a), since p < 2 and {u n } n∈N converges to u in L 2 (D).
We thus turn to (3.5): We skip the lower index n ∈ N and write u instead of u n . If we expand the inner square, we write For first term in the right-hand we use (0.3) and the fact that We focus on the remaining two terms in (3.7): Using Leibniz's rule and an integration by parts we have that Thanks to Lemma 2.1, Lemma 3.1 and since u ∈ C ∞ 0 (D), the first and second term vanish in the limit ε ↓ 0. Hence, By Lemma 2.1 and since u ε ⇀ u * , we may apply (3.3) of Lemma 3.1 with φ = u and v ε = u ε to the limit on the right-hand side above. This yields We now turn to the last term in (3.7). Also here, we use Leibniz rule to compute By an argument similar to the one for (3.10), we reduce to We now apply (3.3) of Lemma 3.1 to v ε = w ε u and φ = u. This implies that Inserting (3.8), (3.10) and (3.11) into (3.7) we have that Since u * = kf , it is easy to see the the right-hand side above equals the right-hand side of (3.5). This establishes (3.5) and concludes the proof of Theorem 0.1, case (a).

Proof of Lemma 3.1
Lemma 3.1 may be proven in a way that is similar to [14][Lemma 3.1]. The first crucial ingredient is the following lemma, that allows to find a suitable partition of the holes H ε by dividing this set into a part containing well separated holes and another one containing the clusters. The next result is the analogue of [14][ Lemma 4.2] with the different rescaling of the radii of the balls generating the set H ε .
For every x ∈ R 3 , we recall the definition of R ε,x in (1.1). We have: . Then there exists a partition H ε := H ε g ∪ H ε b , with the following properties: • There exists a subset of centres n ε (D) ⊆ Φ ε (D) such that Finally, we have that Let γ in Lemma 3.2 be fixed. We construct w ε as done in [14, ]: we set where for each z ∈ n ε (D), the function w ε,z vanishes in the hole B ε α ρz (εz) and solves We also define of the measure We stress that all the previous objects depend on the choice of the parameter γ in Lemma 3.2. The next result states that this parameter may be chosen in so that the norm µ ε − 4πλE ρ H −1 (D) is suitably small. This, together with Lemma 3.2, provides the crucial tool to show Lemma 3.1: There exists γ ∈ (0, α − 1) such that if µ ε is as in (3.17) there exists κ > 0 such that for every random field v ∈ H 1 0 (D) Proof of Lemma 3.1. By construction, it is clear that, for P-almost every realization, the functions w ε ∈ W 1,∞ (R 3 ) ∩ H 1 (R 3 ), vanish in H ε and are such that w ε L ∞ (R 3 ) = 1.
We now turn to (3.1). Using the definitions of w ε g and w ε b and Lemma 3.2 we have that Thanks to (0.2) and the Strong law of Large numbers, for P-a.e. realization the right-hand side vanishes in the limit ε ↓ 0.
We now turn to the second term: Since by the maximum principle |w b ε − 1| 1, we may use the definition of D ε b to bound Thanks to (3.14) in Lemma 3.2, the right-hand side vanishes in the limit ε ↓ 0 for P-almost every realization. Combining this with (3.19) and (3.18) yields (3.1) for w ε − 1. Inequality (3.1) for σ −1 ε ∇w ε follows by Lemma 3.2 and the definition (3.15) of w ε as done in [14][Lemma 3.1]. Limit (3.2) may be argued as done above for (3.1), this time appealing to the bound (0.2) and the stationarity of (Φ, R).
Proof of Lemma 3.3. We divide the proof into steps. The strategy of this proof is similar to the one for [11][Theorem 2.1, (b)].
Step 2. For k ∈ N fixed, let {K ε,x,k } x∈N k,ε be the covering of D constructed in the previous step. We define the random variables and construct the random step function Let v be as in the statement of the lemma and m ε (k) as above. The triangle and Cauchy-Schwarz inequalities imply that so that the proof of the lemma reduces to estimating the norms We now claim that there exists a γ > 0, k ∈ N for a positive exponent κ > 0. Combining these two inequalities with (3.24) establishes Lemma 3.3.
In the remaining part of the proof we tackle inequalities (3.25). We follow the same lines of [11][Theorem 1.1, (b)]. and thus only sketch the main steps for the argument.
Step 3. We claim that We first argue that that This follows by Lemma 5.1 applied to the measure σ −2 ε µ ε : In this case, the random set of centres is Z =Φ ε (D), the random radii R = {R ε,z } z∈Φ(D) , the functions g i = σ −2 ε ∇ ν w ε,z , z ∈Φ ε (D) and the partition {K ε,k,x } x∈N ε,k of the previous step. Note that, by construction, this partition satisfies the assumptions of Lemma 5.1. The explicit formulation of the harmonic functions {w ε,z } z∈n ε (D) defined in (3.16) (c.f. also [11][(2.24)]) implies that for every z ∈ n ε (D) Therefore, Lemma 5.1 and the bounds (3.28) yield that which implies (3.27) thanks to (3.22).
It thus remains to pass from (3.27) to (3.26): We do this by taking the expectation and arguing as for [11][Inequality (4.22)]. We rely on the stationarity of (φ, R), the properties of the Poisson point process and the fact that z ∈ n ε implies that ε α ρ z ε 1+γ and R ε,z ε 1+ 1 2 γ .
Step 4. We now turn to the left-hand side in the second inequality of (3.25) and show that The proof of this step is similar to [11][Theorem 2.1, (b)]: Using the explicit formulation of m ε (k) we reduce to Since ∂D is C 1 and compact, for ε small enough (depending on D) we have By stationarity, the second term in (3.30) is controlled by The remaining term on the righ-hand side may be controlled by the right-hand side in (3.29) by means of standard CLT arguments as done in [11][Inequality (4.23)] for the analogous term. We stress that the crucial observation is that the random variables S k,ε,x − λE ρ are centred up to an error term. We mention that in this case the set K ε,k,x has been defined in a different way from [11] and we use properties (3.22) instead of [11][ (4.13)]. This yields (3.29) Step 5. We show that, given (3.26) and (3.29) of the previous two steps, we may pick γ and k ∈ N such that inequalities (3.25) hold: Thanks to the definition of σ ε and since α ∈ (1; 3), we may find γ close enough to α − 1, e.g. γ = 20 21 (α − 1), and a k ∈ N, e.g. k = − 9 20 (α − 1), such that This, thanks to (3.26), implies that the first inequality in (3.25) holds with the choice κ = 1 2 0(α − 1) > 0. The same values of γ and κ yield that also the right hand side of (3.29) is bounded by ε 1− 1 2 (α−1) . This yields also the remaining inequality in (3.25) and thus concludes the proof of Lemma 3.3.

Proof of Lemma 3.2.
The proof of this lemma follows the same construction implemented in the proof of [11][Lemma 4.1] with d = 3, δ = γ and with the radii {ρ z } z∈Φ ε (D) rescaled by ε α instead of ε 3 . Note that the constraint for γ is due to this different rescaling. In the current setting, we replace ε 2 by ε 1+ 1 2 γ in the definition of the set K ε b in [11][(4.7)]. Estimate (3.14) may be argued as [11][Lemma 4.4] by relying on (0.2).

Proof of Theorem 0.1, (b)
The next lemma is the analogue of Lemma 3.1: Lemma 4.1. For every δ > 0, there exists an ε 0 > 0 and a set A δ ∈ F, having P(A δ ) 1 − δ, such that for every ω ∈ A δ and ε ε 0 there exists a linear map Furthermore, if v ε satisfies the bounds of Lemma 2.1 and Proof of Theorem 0.1, (b). The proof of this statement is very similar to the one for case (a) and we only emphasize the few technical differences: Using the bounds of Lemma 2.1, we have that, up to a subsequence, u ε j ⇀ u * in L p (Ω × D), 1 p < 2. We prove that Identity (4.3) also implies that the full {u ε } ε>0 converges to u * .
As for the proof of Theorem 0.1, case (a), we restrict to the converging subsequence {u ε j } j∈N but we skip the index j ∈ N in the notation. We start by noting that, using the divergence-free condition for u ε and that u ε vanishes on ∂D, we have that for every φ ∈ C ∞ (D) and χ ∈ L ∞ (Ω) E χˆD ∇φ · u * = 0. (4.5) Let χ ∈ L ∞ (Ω) and φ ∈ C ∞ 0 (D) with ∇ · φ = 0 in D be fixed. For every δ > 0, we appeal to Lemma 4.1 to infer that there exists an ε δ > 0 and a set A δ ∈ F, having P(A δ ) 1 − δ, such that for every ω ∈ A δ and for every ε ε δ we may consider the function R ε φ ∈ H 1 0 (D ε ) of Lemma 4.1. Testing equation (0.4) with R ε (ρ), and using that the vector field R ε v is divergence-free, we infer that Using Lemma 4.1 and the bounds of Lemma 2.1 this implies that in the limit ε ↓ 0 we have We now send δ ↓ 0 and appeal to the Dominated Convergence Theorem to infer that Since D has C 1,1 -boundary and is simply connected, the spaces L p (D), p ∈ (1, +∞) admit an L p -Helmoltz decomposition L p (D) = L p div (D) ⊕ L p curl (D) [10][Section III.1]. This, the separability of L p (D), p ∈ [1, +∞), and the arbitrariness of χ and φ in (4.6), allows us to infer that for P-almost realization the function u * satisfies u * = Kf + ∇p(ω; ·) for p(ω; ·) ∈ W 1,p (D), p ∈ [1; 2). By a similar argument, we may use (4.5) to infer that for P-almost every realization and for every v ∈ W 1.q (D), q > 2 we havê D (∇p(·; ω) + Kf ) · ∇v = 0.
Since (4.4) admits a unique mean-zero solution, we conclude that p(ω, ·) does not depend on ω. Finally, since D is regular enough and f ∈ L q (D), standard elliptic regularity yields that p ∈ H 1 (D). This concludes the proof of (4.3).
We now upgrade the convergence of the family {u ε } ε>0 to u * from weak to strong: We claim that for every δ > 0 we may find a set A δ ⊆ Ω with P(A δ ) > 1 − δ such that Here, q ∈ [1, 2). The proof of this inequality follows the same lines of the proof for (3.12) in case (a): In this case, we rely on Lemma 4.1 instead of Lemma 3.1 and use that, thanks to the definition (4.4), it holdsˆD From (4.7), the statement of Theorem 0.1, (b) easily follows: Let, indeed, q ∈ [1, 2) be fixed. For every δ > 0, let A δ be as above. We rewrite and, given an exponent p ∈ (q; 2), we use Hölder's inequality and the assumption on A δ to control Since by Lemma 2.1 the family σ 2 ε u ε is uniformly bounded in every L p (Ω × D) for p ∈ [1, 2), we establish Since δ is arbitrary, we conclude the proof of Theorem 0.1, (b).

Proof of Lemma 4.1
This section is devoted to arguing Lemma 4.1 by leveraging on the geometric information on the clusters of holes H ε contained Lemma 1.2. The idea behind these proof is in spirit very similar to the one for Lemma 3.1 in case (a): As in that setting, indeed, we aim at partitioning the holes of H ε into a subset H ε g of disjoint and "small enough" holes and H ε b where the clustering occurs. The main difference with case (a), however, is due to the fact that we need to ensure that the so-called Stokes capacity of the set H ε b , namely the vector vanishes in the limit ε ↓ 0. The divergence-free constraint implies that, in contrast with the harmonic capacity of case (a), the Stokes capacity is not subadditive. This yields that, if H ε b is constructed as in Lemma 3.2, then we cannot simply control its Stokes-capacity by the sum of the capacity of each ball of H ε b . We circumvent this issue by relying on the information on the length of the clusters given by Lemma 1.2. We do this by adopting the exact same strategy used to tackle the same issue in the case of the Brinkmann scaling in [12]. The following result is a simple generalization of [12][Lemma 3.2] and upgrades the partition of Lemma 3.2 in such a way that we may control the Stokes-capacity of the clustering holes in H ε b . For a detailed discussion on the main ideas behind this construction, we refer to [12][Subsection 2.3].

Lemma 4.2.
Let γ > 0 be as chosen in Lemma 3.3. For every δ > 0 there exists ε 0 > 0 and A δ ⊆ Ω with P(A δ ) > 1 − δ such that for every ω ∈ Ω and ε ε 0 we may choose H ε g , H ε b of Lemma 3.2 as follows: • There exist Λ(β) > 0, a sub-collection J ε ⊆ I ε and constants {λ ε l } z l ∈J ε ⊆ [1, Λ] such that • There exists k max = k max (β, d) > 0 such that we may partition with I ε k ⊆ J ε k for all k = 1, · · · , k max and • For all k = −3, · · · , k max and every • For each k = −3, · · · , k max and z i ∈ I ε k and for all z j ∈ k−1 l=−3 J ε l we have Finally, the set D ε b of Lemma 3.2 may be chosen as The same statement is true for P-almost every ω ∈ Ω for every ε > ε 0 (with ε 0 depending, in this case, also on the realization ω).
Proof of Lemma 4.2. The proof of this result follows the exact same lines of of [12][Lemma 3.2]. We thus refer to it for the proof and to [12][Subsection 3.1] for a sketch of the ideas behind the quite technical argument. We stress that the different scaling of the radii does not affect the argument since the necessary requirement is that ε α << ε. This holds for every choice of α ∈ (1, 3). We also emphasize that in the current setting, Lemma 1.2 plays the role of [12][ Lemma 5.1]. This result is crucial as it provides information on the length of the overlapping balls of H ε . For every δ > 0, we thus select the set A δ of Lemma 1.2 containing those realizations where the partition of H ε satisfies (1.2) and (1.4). Once restricted to the set A δ , the construction of the set H ε b is as in [12][Lemma 3.1].
Equipped with the previous result, we may now proceed to prove Lemma 4.1: Proof of Lemma 4.1. The proof of this is similar to the one in [12][Lemma 2.5] for the analogous operator and we sketch below the main steps and the main differences in the argument. For δ > 0, let ε 0 > 0 and A δ ⊆ Ω be the set of Lemma 4.2; From now on, we restrict to the realization ω ∈ A δ . For every ε < ε 0 we appeal to Lemma 3.2 and Lemma 4.2 to partition H ε = H ε b ∪ H ε g . We recall the definitions of the set n ε ⊆ Φ ε (D) in (3.13) in Lemma 3.2 and of the subdomain D ε b ⊆ D in (4.10) of Lemma 4.2.
It is clear that the previous quantities also depend on ε. However, in order to keep a leaner notation, we skip it in the notation. We use the same understanding for the function φ ε b and the sets {I ε,i } kmax i=−3 and {J ε,i } kmax i=−3 of Lemma 4.2. We define φ b by solving a finite number of boundary value problems in the annuli We stress that, thanks to Lemma 4.2, for every k = −3, · · · , k max , each one of the above collections contains only disjoint annuli. Let φ (kmax+1) = φ. Starting from k = k max , at every iteration step k = k max , · · · , −3, we solve for every z ∈ I ε,k the Stokes system We then extend φ (k) to φ (k+1) outside z∈I k B θ,z and to zero in z∈I k B z .
The analogue of inequalities of [12][(4.12)-(4.14)], this time with the factor ε d−2 d replaced by ε α and with d = 3, is ∇φ (k) 2 (4.13) and (4.14) Moreover, (4.15) These inequalities may be proven exactly as in [12]. We stress that condition (4.9) in Lemma 4.2 is crucial in order to ensure that this construction satisfies the right-boundary conditions. In other words, the main role of Lemma 4.2 is to ensure that, if at step k the function φ (k) vanishes on a certain subset of H ε b , then φ (k+1) also vanishes in that set (and actually vanishes on a bigger set).
We set φ ε b = φ (−3) obtained by the previous iteration. The first property in (4.11) is an easy consequence of (4.14) and the first identity in (4.15). We recall, indeed, that thanks to Lemma 4.2 we have that The second property in (4.11) follows immediately from (4.14). The third line in (4.11) is an easy consequence of the first line in (4.11) and the second inequality in (4.13). Finally, the last inequality in (4.11) follows by multiplying the last inequality in (4.15) with the factor σ ε and using that, since φ ∈ C ∞ , we have that Thanks to Lemma 4.2 and the definition of the set n ε in Lemma 3.2, the previous inequality yields the last bound in (4.11).
Step 3. (Construction of φ ε g ) Equipped with φ ε b satisfying (4.11), we now turn to the construction of φ ε g . Also in this case, we follow the same lines of [12][Proof of Lemma 2.5, Step 3] and exploit the fact that the set H ε g is only made by balls that are disjoint and have radii ε α ρ that are sufficiently small. We define the function φ ε g exactly as in [12] (εz), B 2,i := B dε,z (εz), C z := B z \T z , D z := B 2,z \B z .
We now turn to show the remaining part of (4.12): We remark that, since z ∈ n ε (D), Lemma 3.2 and definition (4.16) yield that where γ > 0 is as in Lemma 4.2 and β > 0 is as in (0.5). Equipped with the previous bounds, the analogue of estimates [12][(4.26)-(4.30)] yield that for every z ∈ n ε (D) Since B 2,z = D z ∪ C z ∪ T z and the function φ ε g − φ is supported only on z∈n ε (D) B 2,z , we infer that for every z ∈ n ε (D), it holds Summing over z ∈ n ε we obtain the last two inequalities in (4.12). We thus established (4.12) and completed the proof of Step 1.
Step 4. (Properties of R ε ) We now argue that R ε defined in Step 1. satisfies all the properties enumerated in Lemma 4.2. It is immediate to see from (4.12) and (4.11) that R ε φ vanishes on H ε and is divergence-free in D. Inequalities (4.1) also follow easily from the inequalities in (4.12) and (4.11) and arguments analogous to the ones in Lemma 3.1. We stress that, in this case, we appeal to condition (0.5) and, in the expectation, we need to restrict to the subset A δ ⊆ Ω of the realizations for which R ε may be constructed as in Step 1.
To conclude the proof, it only remains to tackle (4.2). We do this by relying on the same ideas used in Lemma 3.1 in the case of the Poisson equation. We use the same notation introduced in Step 2. We begin by claiming that (4.2) reduces to show that for every i = 1, · · · , 3 where We use the definition of R ε φ to rewrite for every ω We claim that, after multiplying by 1 A δ and taking the expectation, the last two integrals on the righthand side vanish in the limit. In fact, using the triangle and Cauchy-Schwarz's inequalities and combining them with (4.11) and the uniform bounds for {v ε } ε>0 we have that (3.14) = 0.
Hence, we show (4.2) provided that We further reduce this to (4.19) if An argument analogous to the one outlined in [12] to pass from the left-hand side of [12][(4.34)] to the one in [12][(4.39)] yields that We stress that in the current setting we use again the uniform bounds on the sequence σ −1 ε ∇u ε and we rely on estimates (4.18) instead of [12][(4.26)-(4.30)]. To pass from (4.20) to (4.19) it suffices to use the smoothness of φ and, again, the bounds on the family {v ε } ε>0 . We thus established that (4.2) reduces to (4.19).
Note that, up to a relabelling of the indices k = −3, · · · , k max , the previous partition satisfies (1.3) of Lemma 1.2.
For any set χ ⊆ Φ ε (D), we say that A contains a chain of length M ∈ N, M 2, if there exist z 1 , · · · , z M ∈ χ such that B 4ε α ρz i (εz i ) ∩ B 4ε α ρz j (εz j ) = ∅, for all i, j = 1, · · · M . We say that A contains a chain of size 1 if and only if A = ∅.
Equipped with this notation, (1.2) follows provided we argue that for κ suitably chosen, there exists k 0 < k max − 1 such that, P-almost surely and for ε small, the sets {I ε,k ∪ I ε,k+1 } kmax k=k 0 are empty. This is equivalent to prove that they do not contain any chain of size at least 1. Similarly, (1.4) is obtained if we find an M ∈ N such that P-almost surely and for ε small enough, all the sets {I ε,k ∪ I ε,k+1 } k 0 k=−3 contain chains of length at most M − 1.
For M ∈ N and k = −3, · · · , k max , we define the events A k,ε,M := I ε,k ∪ I ε,k+1 contains a chain of length at least M .
Using these two bounds, we may apply Borel-Cantelli to the family of events {A k,ε j ,M } and conclude (5.3) and (5.4).
We now turn to case (ii). Identity (5.4) may be rewritten as P( This implies that for every δ > 0, we may pick ε 0 > 0 such that the set P( ε<ε 0 ( kmax k=−3 A k,ε,M ) c ) > 1 − δ. The statement of (ii) immediately follows if we set A δ := ε<ε 0 ( kmax k=−3 A k,ε,M ) c . The same argument applied to (5.3) implies the same statement for (1.2). Lemma 5.1. Let Z := {z i } i∈I ⊆ D be a collection of points and let R := {r i } i∈I ⊆ R + such that the balls {B r i (z i )} i∈I are disjoint. We define the measure where g i ∈ L 2 (∂B r i (z i )). Then, there exists a constant C < +∞ such that for every Lipschitz and (essentially) disjoint covering {K j } j∈J of D such that