Inverse of Divergence and Homogenization of Compressible Navier–Stokes Equations in Randomly Perforated Domains

Abstract

We analyze the behavior of weak solutions to compressible viscous fluid flows in a bounded domain in ${{\,\mathrm{{\mathbb {R}}}\,}}^3$, randomly perforated by tiny balls with random size. Assuming the radii of the balls scale like $\varepsilon ^\alpha$, $\alpha > 3$, with $\varepsilon$ denoting the average distance between the balls, the problem homogenize to the same limiting equation. Our main contribution is a construction of the Bogovskiĭ operator, uniformly in $\varepsilon$, without any assumptions on the minimal distance between the balls.

1 Introduction

The goal of this paper is to analyze the effective behavior of a compressible viscous fluid in randomly perforated domains. We consider a bounded domain $D \subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ which for $\varepsilon > 0$ is perforated by random balls $B_{\varepsilon ^{\alpha } r_i}(\varepsilon z_i)$ with $\alpha > 3$, and show that weak solutions to the Navier–Stokes equations in these perforated domains converge as $\varepsilon \rightarrow 0$ to a weak solution to the same equation in D. Compared to the previous results, mostly restricted to the periodic arrangement of the holes [9, 15, 28] or at least assuming minimal distance between the holes being $\varepsilon$ [30], we do not require any such assumptions on the minimal distance between the holes.

As in the previous works on this topic, the key step in the proof of the homogenization result is the construction of the Bogovskiĭ operator (inverse of the divergence), bounded independently of $\varepsilon$. This operator is then used in a classical way to construct a test function, hence providing uniform estimates on the density and velocity of the fluid.

The question of how small perforations in the domain influence the original equation has a long history. In the case of a perforation by periodically arranged balls in d dimensions, Cioranescu and Murat [7] (see also [8, 25]) studied effective behavior of the Poisson equation with zero Dirichlet boundary conditions on the balls. Denoting the distance between balls by $\varepsilon$ and assuming the radii scale like $\varepsilon ^{\frac{d}{d-2}}$, they identified an additional Brinkman term (“A strange term coming from nowhere”) in the limiting equation.

Heuristically, the homogeneous boundary conditions on the perforation push the solution u towards 0, with the strength related to the size of holes. Focusing on the case $d=3$, holes of size $\varepsilon ^\alpha$, $\alpha > 3$ are too tiny to make any difference. If the balls are larger than $\varepsilon ^3$ (hence push u stronger to 0), solutions converge to 0 as $\varepsilon \rightarrow 0$, and only a rescaling by some negative power of $\varepsilon$ may lead to a reasonable limiting problem.

For the incompressible stationary Stokes and Navier–Stokes equations with periodic distribution of holes, Allaire [3, 4] gave a full description to all cases $\alpha >1$ and all dimensions $d\ge 2$. More precisely, for $d=3$ if $\alpha \in (1,3)$, which corresponds to the supercritical case of large particles, he obtained Darcy’s law; for the critical case $\alpha =3$, an additional friction term occurs and gives rise to Brinkman’s law. The subcritical case $\alpha >3$ corresponding to small particles leads to the same system of Stokes and Navier–Stokes equations. The case $\alpha =1$ was studied in [2] for the steady incompressible Stokes system.

The results on the effective behavior of compressible fluids in perforated domains are much more recent. Masmoudi [31, 32] considered compressible Navier–Stokes equations in the domain perforated by periodic balls with $\alpha = 1$, and obtained in the limit Darcy’s law. This result was later generalized by Feireisl et al. [16] also to the case of the full Navier–Stokes–Fourier system, which besides density and velocity of the fluid takes also into account the fluid temperature. In the case of slightly smaller balls with $1< \alpha < 3$, assuming simultaneous rescaling of the pressure (Low-Mach number limit) to avoid the need to study the “cell problem” with unknown density, Höfer et al. [23] showed convergence of the rescaled solution to Darcy’s law. Finally, for tiny balls with $\alpha > 3$ (the subcritical case) Feireisl and Lu [15] considered stationary Navier–Stokes equations and showed convergence to the same equations in the domain without holes. This result was later improved to the case of more general adiabatic exponent in the pressure [9] as well as to the time-dependent case [30]. In all these works the perforation is assumed to be periodic, or at least the minimal distance between the holes is comparable with $\varepsilon$.

In this work we also consider only the subcritical case $\alpha > 3$, but with random arrangement of holes instead of the periodic one. Unless one additionally assumes that the holes in the random case may not lie close to each other, the key argument in the previous works, the Bogovkiĭ operator, can not be constructed as before. Instead, we show that while the holes can be close to each other, there exists a fixed number N such that there are at most N balls clustered together. Since N is fixed, the construction of the Bogovkiĭ operator can be done.

Our inspiration how to approach the case of randomly perforated domains comes from a recent work of Giunti, Höfer, and Velázquez. In [22], they considered a Poisson equation in a domain perforated by random balls of critical size, thus obtaining the Brinkman law in the limit. The main challenge is to understand possible clustering of the holes and control the capacity of those. In subsequent works, they also considered the incompressible Stokes problem [20, 21] as well as convergence to the Darcy’s model in the supercritical situation [19]. Compared to these works our situation is simpler, since in the subcritical situation we can prove a deterministic upper bound on the size of clusters. Let us also mention that previously Beliaev and Kozlov [5] also considered homogenization of Stokes equations in a supercritically perforated domain.

A different, though quite related, setting is of the flow of a colloidal suspension, that is, of a fluid mixed with moving obstacles. This question traces back to one part of Einstein’s PhD thesis [13], where he formally derives effective viscosity of such suspension, assuming low-volume fraction of the obstacles. With the recent progress in the field of stochastic homogenization, this question was rigorously approached by several groups, first under the assumption of uniform separation of balls [12, 33] and very recently under less restrictive assumptions [11, 18].

2 Setting and the Main Results

In this section we define the perforated domain, formulate the Navier–Stokes equations governing the fluid motion, and state the main results. We consider $D \subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ being a bounded domain with a $C^2$-boundary. To simplify the probabilistic argument we farther assume the domain D is star-shaped with respect to the origin, that is, for any $x \in D$ the segment $\{ \lambda x: \lambda \in [0,1] \} \subset D$.

We model the perforation of D using the Poisson point process, though the arguments can be easily generalized to a larger class of point processes. For an intensity parameter $\lambda > 0$, the Poisson point process is defined as a random collection of points $\Phi = \{z_j\}$ in ${{\,\mathrm{{\mathbb {R}}}\,}}^3$ characterized by the following two properties:

• for any two measurable and disjoint sets $S_1, S_2 \subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$, the random sets $S_1 \cap \Phi$ and $S_2 \cap \Phi$ are independent;

• for any measurable set $S \subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ and $k \in {\mathbb {N}}$ holds. ${\mathbb {P}}(N(S) = k) = \frac{(\lambda |S|)^ke^{-\lambda |S|}}{k!}$

Here $N(S)=\# (S\cap \Phi )$ counts the number of points $z_j\in S$ and |S| denotes measure of S. In addition to the random locations of the balls, modeled by the above Poisson point process, we also assume the balls have random size. For that, let ${\mathcal {R}} = \{r_i\} \subset [0,\infty )$ be another random process of independent identically distributed random variables with finite mth moment, that is,

$$\begin{aligned} {{\,\mathrm{{\mathbb {E}}}\,}}(r_i^m) < \infty \hbox { for some}\ m>0, \end{aligned}$$

and which are independent of $\Phi$. In other words, to each point $z_j \in \Phi$ (center of a ball) we associate also a radius of the ball $r_j \in [0,\infty )$. The exact range of m we can work with will be specified in Theorem 1 below. The random process $(\Phi ,{\mathcal {R}})$ on ${{\,\mathrm{{\mathbb {R}}}\,}}^3 \times {{\,\mathrm{{\mathbb {R}}}\,}}_+$ is called marked Poisson point process, and can be viewed as a random variable $\omega \in \Omega \mapsto (\Phi (\omega ),{\mathcal {R}}(\omega ))$, defined on an abstract probability space $(\Omega ,{\mathcal {F}},{\mathbb {P}})$.

To define the perforated domain $D_\varepsilon$, for $\alpha > 2$ and $\varepsilon > 0$ we set

$$\begin{aligned} \Phi ^{\varepsilon }(D):= \left\{ z \in \Phi \cap \frac{1}{\varepsilon }D: {{\,\textrm{dist}\,}}(\varepsilon z,\partial D) > \varepsilon \right\} , \quad D_\varepsilon := D \setminus \bigcup _{z_j \in \Phi ^{\varepsilon }(D)} B_{\varepsilon ^{\alpha } r_j}(\varepsilon z_j). \end{aligned}$$

(1)

To simplify the exposition and to avoid the need to analyze behavior near the boundary, we only removed those balls from D which are not too close to the boundary $\partial D$. This is also a common assumption in the periodic situation, see, for example, [15, relation (1.3)]. The domain D being star-shaped implies that $\Phi ^\varepsilon (D)$ are monotonically increasing as $\varepsilon \rightarrow 0$.

Our main result is the following existence result for a uniformly bounded Bogovskiĭ operator:

Theorem 1

Let $\alpha > 2$, $D\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ be a bounded star-shaped domain with respect to the origin with $C^2$-boundary, and $(\Phi ,{\mathcal {R}})=(\{z_j\},\{r_j\})$ be a marked Poisson point process with intensity $\lambda > 0$. We assume the radii $r_j \geqq 0$ fulfil ${{\,\mathrm{{\mathbb {E}}}\,}}(r_j^m)<\infty$ for some $m>3/(\alpha -2)$. Then for all $1<q<3$ which fulfil

$$\begin{aligned} \alpha -\frac{3}{m}>\frac{3}{3-q}, \end{aligned}$$

(2)

there exists a random almost surely positive $\varepsilon _0 = \varepsilon _0(\omega )$ such that for $0<\varepsilon \le \varepsilon _0$ there exists a bounded linear operator

$$\begin{aligned} {\mathcal {B}}_\varepsilon :L^q(D_\varepsilon ) / {\mathbb {R}} \rightarrow W_0^{1,q}(D_\varepsilon ;{{\,\mathrm{{\mathbb {R}}}\,}}^3) \end{aligned}$$

with $D_\varepsilon$ defined in (1), such that for all $f\in L^q(D_\varepsilon )$ with $\int _{D_\varepsilon } f = 0$

$$\begin{aligned} {{\,\textrm{div}\,}}({\mathcal {B}}_\varepsilon (f))=f \text { in } D_\varepsilon ,\quad \Vert {\mathcal {B}}_\varepsilon (f)\Vert _{W_0^{1,q}(D_\varepsilon )}\le C \, \Vert f\Vert _{L^q(D_\varepsilon )}, \end{aligned}$$

where the constant $C>0$ is independent of $\omega$ and $\varepsilon$.

In other words, Theorem 1 provides a solution $\textbf{u}$ to the equation ${{\,\textrm{div}\,}}\textbf{u}= f$ in $D_\varepsilon$ with $\textbf{u}|_{\partial D_\varepsilon }=0$ such that $\Vert \nabla {\textbf {u}}\Vert _{L^q(D_\varepsilon )} \leqq \Vert {\textbf {u}}\Vert _{W^{1,q}_0(D_\varepsilon )} \leqq C \Vert f \Vert _{L^q(D_\varepsilon )}$.

As an application for this result we show homogenization of compressible Navier–Stokes equations in perforated domains $D_\varepsilon$. For $\varepsilon > 0$, the unknown density $\varrho _\varepsilon$ and velocity $\textbf{u}_\varepsilon$ of a viscous compressible fluid are described by

$$\begin{aligned} \begin{aligned} \partial _t \varrho _\varepsilon + {{\,\textrm{div}\,}}(\varrho _\varepsilon \textbf{u}_\varepsilon )&=0{} & {} \text { in } (0,T) \times D_\varepsilon ,\\\partial _t(\varrho _\varepsilon \textbf{u}_\varepsilon ) + {{\,\textrm{div}\,}}(\varrho _\varepsilon \textbf{u}_\varepsilon \otimes \textbf{u}_\varepsilon ) + \nabla p(\varrho _\varepsilon )&= {{\,\textrm{div}\,}}{\mathbb {S}}(\nabla \textbf{u}_\varepsilon ) + \varrho _\varepsilon \textbf{f}+\textbf{g}{} & {} \text { in } (0,T) \times D_\varepsilon ,\\\textbf{u}_\varepsilon&=0{} & {} \text { on } (0,T) \times \partial D_\varepsilon , \end{aligned} \end{aligned}$$

(3)

where ${\mathbb {S}}$ denotes the Newtonian viscous stress tensor of the form

$$\begin{aligned} {\mathbb {S}}(\nabla {\textbf {u}})=\mu \bigg (\nabla {\textbf {u}}+\nabla ^T {\textbf {u}}-\frac{2}{3} {{\,\text {div}\,}}({\textbf {u}}){\mathbb {I}}\bigg )+\eta {{\,\text {div}\,}}({\textbf {u}}){\mathbb {I}},\quad \mu > 0, \quad \eta \geqq 0,\end{aligned}$$

$p(\varrho ) = a\varrho ^{\gamma }$ denotes the pressure with $a > 0$ and the adiabatic exponent $\gamma \geqq 1$, and $\textbf{f}$ and $\textbf{g}$ are external forces, which are for simplicity assumed to satisfy $\Vert {\textbf {f}}\Vert _{L^\infty ((0,T)\times {{\,\mathrm {{\mathbb {R}}}\,}}^3;{{\,\mathrm {{\mathbb {R}}}\,}}^3)} + \Vert {\textbf {g}}\Vert _{L^\infty ((0,T)\times {{\,\mathrm {{\mathbb {R}}}\,}}^3;{{\,\mathrm {{\mathbb {R}}}\,}}^3)} \leqq C$. We also fix the total mass $\int _{D_\varepsilon } \varrho _\varepsilon (0,\cdot )={\mathfrak {m}}>0$ independently of $\varepsilon >0$, and supplement the equations with the initial conditions for $\varrho _\varepsilon$ and $\varrho _\varepsilon \textbf{u}_\varepsilon$.

While the existence of classical solutions to (3) is known only in some special cases, the existence theory for weak solutions is quite developed [14, 27, 34]. In particular, for fixed $\varepsilon > 0$ the domain $D_\varepsilon$ is smooth enough to grant existence of global weak solutions. The key point here is to overcome the lack of uniform estimates in $\varepsilon$ on the smoothness of $D_\varepsilon$, in particular to obtain uniform bounds on the solution, which in this setting are usually obtained using a bounded Bogovkiĭ operator. More precisely, Theorem 1 together with a simple existence result for the cut-off function (see Lemma 4) are the only points in the arguments [9, 15, 30], where the information on the structure of the perforation in $D_\varepsilon$ is being used.

In the following we state one of the implications of Theorem 1, the corresponding precise formulation as well as the definition of the finite energy weak solutions in the case of periodic perforation being [30, Definition 1.3, Theorem 1.6]:

Theorem 2

Assume $\alpha >3$. Let $D\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ be a bounded star-shaped domain with respect to the origin with $C^2$-boundary and let $(\Phi ,{\mathcal {R}})=(\{z_j\},\{r_j\})$ be a marked Poisson point process with intensity $\lambda > 0$, and $r_j \geqq 0$ with ${\mathbb {E}}(r_j^M) < \infty$, $M=\max \{3,m\}$, where $m>3/(\alpha -3)$. Furthermore let

$$\begin{aligned} {\mathfrak {m}}>0,\; \gamma >6. \end{aligned}$$

For $0< \varepsilon < 1$ let $[\varrho _\varepsilon ,\textbf{u}_\varepsilon ]$ be a family of finite energy weak solutions for the no-slip compressible Navier–Stokes equations (3) in $(0,T)\times D_\varepsilon$, with $D_\varepsilon$ as in (1). Assume that the initial conditions

$$\begin{aligned} \varrho _\varepsilon (0,\cdot ) = \varrho _{\varepsilon , 0}\quad \text {and}\quad (\varrho _\varepsilon \textbf{u}_\varepsilon )(0,\,\cdot )=\textbf{q}_{\varepsilon ,0} \end{aligned}$$

satisfy

$$\begin{aligned} \varrho _{\varepsilon , 0}\in L^\gamma (D_\varepsilon ),\quad \textbf{q}_{\varepsilon ,0}=0 \text { whenever } \varrho _{\varepsilon ,0}=0,\quad \bigg \Vert \frac{|\textbf{q}_{\varepsilon ,0}|^2}{\varrho _{\varepsilon ,0}}\bigg \Vert _{L^1(D_\varepsilon )}\le C,\\ {\tilde{\varrho }}_{\varepsilon ,0}\rightharpoonup \varrho _0 \text { weakly in } L^\gamma (D),\quad \tilde{\textbf{q}}_{\varepsilon ,0}\rightharpoonup \textbf{q}_0 \text { weakly in } L^\frac{2\gamma }{\gamma +1}(D), \end{aligned}$$

where $C>0$ is independent of $\varepsilon$. Then for almost every $\omega \in \Omega$ there exists $\varepsilon _0=\varepsilon _0(\omega ) > 0$ such that

$$\begin{aligned} \sup _{\varepsilon \in (0,\varepsilon _0)} ( \Vert \varrho _\varepsilon \Vert _{L^\infty (0,T;L^{\gamma }(D_\varepsilon ))}+ \Vert \varrho _\varepsilon \Vert _{L^{\frac{5\gamma }{3}-1}((0,T) \times D_\varepsilon )} + \Vert \textbf{u}_\varepsilon \Vert _{L^2(0,T;W^{1,2}_0(D_\varepsilon ))}) \le C \end{aligned}$$

and, up to a subsequence, the zero extensions satisfy

$$\begin{aligned} {\tilde{\varrho }}_\varepsilon \overset{*}{\rightharpoonup }\varrho \text { in } L^\infty (0,T;L^\gamma (D)),\quad \tilde{\textbf{u}}_\varepsilon \rightharpoonup \textbf{u}\text { in } L^2(0,T;W^{1,2}_0(D)), \end{aligned}$$

where the limit $[\varrho ,\textbf{u}]$ is a renormalized finite energy weak solution to the problem (3) in the limit domain D provided $\frac{\gamma -6 }{2\gamma - 3}\big (\alpha - \frac{3}{m} \big ) > 3$.

The restriction $M\ge 3$ is made in order to construct suitable cut-off functions, as will be clear from the proof of Lemma 4 later on. Using Theorem 1 and Lemma 4, which are the only two spots in the proof where the structure of the perforation plays any role, the proof of Theorem 2 follows verbatim as in [30]. To manifest how Theorem 1 and Lemma 4 are actually applied, in Sect. 6 we will formulate a similar (but simpler) homogenization statement for the stationary case and sketch its proof.

2.1 Notation

Through the whole paper, we use the following notation:

• $(\Omega , {\mathcal {F}},{\mathbb {P}})$ is the probability space associated to the marked point process $(\Phi ,{\mathcal {R}})$.

• $L_0^p(D):=\{f\in L^p(D): \int _D f=0\}$

• |S| denotes the Lebesgue measure of a measurable set $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$.

• For a function f with domain of definition D or $D_\varepsilon$, we denote by ${\tilde{f}}$ the zero extension to ${{\,\mathrm{{\mathbb {R}}}\,}}^3$, that is, we define

  $$\begin{aligned} {\tilde{f}}=f \text { in } D \text { or } D_\varepsilon ,\quad {\tilde{f}}=0 \text { in } {{\,\mathrm{{\mathbb {R}}}\,}}^3\setminus D \text { or } {{\,\mathrm{{\mathbb {R}}}\,}}\setminus D_\varepsilon . \end{aligned}$$

• Boxes are sets of the form $A_x \times A_y \times A_z$, where $A_x,A_y,A_z \subset {{\,\mathrm{{\mathbb {R}}}\,}}$ are intervals.

• For a factor $\lambda >0$ and a set $M\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$ we define $\lambda M:= \{ \lambda x: x \in M \}$.

• For two sets $M,N\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$, we set $\displaystyle {{\,\textrm{dist}\,}}(M,N)=\inf _{x\in M, y\in N} |x-y|$, where |x| is the usual Euclidean norm, and $\displaystyle {{\,\textrm{dist}\,}}_\infty (M,N)=\inf _{x\in M, y\in N} \Vert x-y\Vert _\infty =\inf _{x\in M, y\in N} \max _{1\le i\le 3} |x_i-y_i|$.

• We write $a\lesssim b$ whenever there is a constant $C>0$ that does not depend on $\varepsilon , a$, and b such that $a\le C\, b$. The constant C might change its value whenever it occurs.

Moreover, if no ambiguity occurs, we denote the function spaces as in the scalar case even if the functions are vector- or matrix-valued, for example, we write $L^p(D)$ instead of $L^p(D;{{\,\mathrm{{\mathbb {R}}}\,}}^3)$.

Organization of the Paper: The rest of this paper is organized as follows. In the next section we formulate the probabilistic statements (Theorem 3 and Proposition 1) as well as the analytical framework (Lemma 2) needed for the construction of the Bogovskiĭ operator (Theorem 1). The proofs of these results are content of Sect. 4 (probabilistic part) and Sect. 5 (analytical part). The last section is devoted to a quick sketch of the homogenization result.

3 Ingredients for the Proof of Theorem 1

The proof of Theorem 1 consists of two parts: stochastic and analytical. The stochastic result, Theorem 3, states that for small enough (depending on $\omega \in \Omega$) $\varepsilon > 0$ the balls with radii $\varepsilon ^{\alpha } r_j$ are disjoint and actually little bit separated. The previous construction of the uniformly bounded Bogovkiĭ operator in the $L^2$-setting requires a boundary layer of size third root of the radius of the balls without hitting other balls—this is where the condition $\alpha \geqq 3$ enters. Since in a generic random arrangement of balls the balls are not that much separated, we relax this assumption by replacing one ball with finitely many balls. More precisely, we show that there exists a deterministic number $N=N(\alpha )$ such that we can group balls into clusters of size at most N so that the clusters stay separated from each other.

Theorem 3

Let $\alpha >2$ and $\lambda >0$ be the intensity of a marked Poisson point process $(\Phi ,{\mathcal {R}}) = (\{z_j\},\{r_j\})$ with $r_j \geqq 0$ and ${{\,\mathrm{{\mathbb {E}}}\,}}(r_j^m)<\infty$, where $m>0$ satisfies

$$\begin{aligned} m>\frac{3}{\alpha -2}. \end{aligned}$$

Let $0<\delta < \alpha -1-\frac{3}{m}$, $\kappa \in (\max (1,\delta ),\alpha -1-\frac{3}{m})$, and $\tau \geqq 1$. Then there exists a random variable $\varepsilon _0$, which is almost surely positive, satisfying the following:

• 1. For every $0<\varepsilon \le \varepsilon _0$ it holds that

  $$\begin{aligned} \max _{z_i\in \Phi ^\varepsilon (D)} \tau \varepsilon ^\alpha r_i\le \varepsilon ^{1+\kappa }, \end{aligned}$$

  and for every $z_i,z_j\in \Phi ^\varepsilon (D),\, z_i\ne z_j$,

  $$\begin{aligned} B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z_i)\cap B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z_j)=\emptyset . \end{aligned}$$

• 2. Let

  

  (4)

  Then for each $0 < \varepsilon \leqq \varepsilon _0$, there are finitely many open boxes $\{I_i^\varepsilon \} \subset D$ satisfying that

  • (a) the boxes $I_i^\varepsilon$ cover the balls, that is, for any $z \in \Phi ^\varepsilon (D)$ we have $B_{\varepsilon ^{1+\kappa }}(\varepsilon z)\subset \bigcup _i I_i^\varepsilon$;

  • (b) any box $I_i^\varepsilon$ contains at most N points from $\varepsilon \Phi ^\varepsilon (D)$;

  • (c) balls are well inside the box: for $\varepsilon z \in I_i^\varepsilon \cap \varepsilon \Phi ^\varepsilon (D)$ holds ${{\,\textrm{dist}\,}}(B_{\varepsilon ^{1+\kappa }}(\varepsilon z),\partial I_i^\varepsilon )\ge \frac{1}{16N}\varepsilon ^{1+\delta }$;

  • (d) any two distinct boxes $I_i^\varepsilon$ and $I_j^\varepsilon$ are well separated: ${{\,\text {dist}\,}}_\infty (I_i^\varepsilon ,I_j^\varepsilon ) \geqq \frac{1}{4N}\varepsilon ^{1+\delta }$;

  • (e) the shortest side of $I^\varepsilon _i$ is at least $\frac{1}{2N}\varepsilon ^{1+\delta }$ while the longest side is at most $\varepsilon ^{1+\delta }$.

The proof of the second part of Theorem 3 uses that for $0 < \varepsilon \leqq \varepsilon _0$ any cube with side length $\varepsilon ^{1+\delta }$ contains at most N points from the Poisson point process. This can hold only if $\delta > 0$, since the number of points in a cube of size $\varepsilon ^{1+0}$ is Poisson distributed, that is, any number of points appears there with small but positive probability.

Proposition 1

Let $d\ge 1$, $\delta >0$ be fixed, and let $\{z_j\}\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$ be points generated by a Poisson point process of intensity $\lambda >0$. In addition, let $D \subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$ be a bounded star-shaped domain. Then there exists a deterministic constant $N(\delta , d)\in {{\,\mathrm{{\mathbb {N}}}\,}}$ and a random variable $\varepsilon _0(\omega ,\lambda ,D)$, which is almost surely positive, such that for all $0<\varepsilon \le \varepsilon _0$ and any $x \in {{\,\mathrm{{\mathbb {R}}}\,}}^d$ the cube $x + [0,\varepsilon ^{1+\delta }]^d$ contains at most N points from $D \cap \varepsilon \Phi$.

To construct the Bogovskiĭ operator ${\mathcal {B}}_\varepsilon$ in $D_\varepsilon$ from Theorem 1 we use local Bogovskiĭ operators for each box $I_i^\varepsilon$ to modify the Bogovskiĭ operator in D. Instead of making explicit construction in each box $I_i^\varepsilon$, we invoke a general result on the existence of Bogovskiĭ operator [1] for a class of domains (so-called John domains) and show that each box $I_i^\varepsilon$ minus the balls is a John domain—for this the outcomes of Theorem 3 will be crucial. In particular, we need that there are at most N balls in one box, the balls are not close to each other and they are tiny, compared to the size of the box.

Definition 1

For a constant $c>0$, a domain $U\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$ is said to be a c-John domain if there exists a point $x_0\in U$ such that for any point $x\in U$ there is a rectifiable path $\Gamma :[0,\ell ]\rightarrow U$ which is parametrized by arc length with

$$\begin{aligned} \Gamma (0)=x,\quad \Gamma (\ell )=x_0,\quad \forall t\in [0,\ell ]: |\Gamma (t)-x|\le c\, {{\,\textrm{dist}\,}}(\Gamma (t),\partial U). \end{aligned}$$

The following lemma states that any $I_i^\varepsilon$ is a c-John domain:

Lemma 2

Under the assumptions of Theorem 3 for fixed

$$\begin{aligned} 0<\delta <\frac{\alpha -2-\frac{3}{m}}{2}, \end{aligned}$$

(5)

let $0 < \varepsilon \leqq \varepsilon _0$. Then, for every box $I_i^\varepsilon$ constructed in Theorem 3, the domain

$$\begin{aligned} U:=I^\varepsilon _i \setminus \bigcup \limits _{z_j\in \varepsilon ^{-1} I_i^\varepsilon \cap \Phi ^\varepsilon (D)} B_{\varepsilon ^\alpha r_j}(\varepsilon z_j) \end{aligned}$$

(6)

is a c-John domain with $c=c(N)$, where N is defined in (4).

In particular, for any $1< q < \infty$ there exists a uniformly bounded Bogovkiĭ operator ${\mathcal {B}}_U: L^q_0(U) \rightarrow W^{1,q}_0(U)$, that is, there exists a constant C, independent of $\varepsilon$, such that for any $f \in L^q_0(U)$

$$\begin{aligned} {{\,\textrm{div}\,}}({\mathcal {B}}_{U}(f))=f,\quad ||{\mathcal {B}}_U(f)||_{W_0^{1,q}(U)} \le C ||f||_{L_0^q(U)}. \end{aligned}$$

4 The Probabilistic Results

The goal of this section is to prove the stochastic part of the result, Theorem 3, second part of which is based on Proposition 1 about the distribution of the random points, modeled by the Poisson point process. Fixing $\delta > 0$, this proposition states that for $\varepsilon$ small enough, for any cube of side length $\varepsilon ^{1+\delta }$ inside a fixed domain D there are at most $N=N(\delta ,d)$ of the rescaled points $\varepsilon z$ in the cube. The heuristic explanation of this is as follows: assuming we only need to consider a disjoint set of cubes and fixing $\varepsilon >0$, the number of cubes in D which we have to consider scales like $\frac{1}{\varepsilon ^{(1+\delta )d}}$. At the same time, the probability of one cube of side length $\varepsilon ^{1+\delta }$ having more than N points scales in the case of the Poisson point process like $(\frac{\varepsilon ^{(1+\delta )d}}{\varepsilon ^d})^N = \varepsilon ^{\delta N d}$. Hence, choosing N large enough so that $\frac{1}{\varepsilon ^{(1+\delta )d}} \varepsilon ^{\delta N d} \ll 1$ should lead to the result.

Proof of Proposition 1

We start with a special case, which will be later used to prove the general case.

Claim: There exists $N_1 \in {\mathbb {N}}$ and an almost surely positive random variable $\varepsilon _0(\omega )$ such that for any dyadic $\varepsilon = 2^{-l}$ smaller than $\varepsilon _0$, any half-closed cube $Q_{\varepsilon ,z} = \varepsilon ^{1+\delta } z + [0,\varepsilon ^{1+\delta })^d$, $z \in {\mathbb {Z}}^d$, contains at most $N_1$ points from $\frac{\varepsilon }{2}\Phi (\omega ) \cap D$.

If rescaled up by a factor 2 the claim says that in a cube with side length $(2\varepsilon )^{1+\delta }$ there are at most $N_1$ points, and we are considering points (more precisely cubes) inside 2D instead of D only. The reason for this choice will be clear later in the proof.

For $l\in {{\,\mathrm{{\mathbb {N}}}\,}}$ and $\varepsilon =2^{-l}$, we define

$$\begin{aligned} B_l:=\{\omega \in \Omega :&\text{ one } \text{ of } \text{ the } \text{ dyadic } \text{ cubes } Q_{2^{-l},z} \text{ contains }\\&\text{ at } \text{ least } N_1 \text{ points } \text{ from } 2^{-l-1} \Phi \cap D\}. \end{aligned}$$

In order to proof the claim, it is enough to show $\sum _{l\ge 0} {\mathbb {P}}(B_l)<\infty$ and apply the Borel–Cantelli lemma. Recall that for any measurable bounded set $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$, we denote by $N(S)=\#(S\cap \Phi )$ the number of random points in S. First, by rescaling, we see that

$$\begin{aligned} B_l= \{\omega \in \Omega :&\text { there exists a dyadic cube } Q_{2^{-l},z}\text { such that }N(2^{l+1} (Q_{2^{-l},z}\cap D))\ge N_1\}. \end{aligned}$$

Since we can cover D with at most $C\,|D|\, (2^l)^d$ cubes $Q_{2^{-l},z}$ and due to the stationarity of the process $\Phi$, we estimate

$$\begin{aligned} \begin{aligned} {\mathbb {P}}(B_l)&\le C(D)\, 2^{ld}\, {\mathbb {P}}\big (N(2^{l+1}Q_{2^{-l},0})\ge N_1\big )=C(D)\, 2^{ld}\, {\mathbb {P}}\big (N(2^{l+1}\cdot 2^{-l(1+\delta )}[0,1)^d)\ge N_1\big )\\&\le C(D,N_1)\, 2^{ld}\cdot 2^{(1-l\delta )N_1d}=C(D,N_1)\, 2^{ld(1-N_1\delta )}, \end{aligned} \end{aligned}$$

(7)

where in the last inequality we used that the points in $\Phi$ are Poisson-distributed, that is,

$$\begin{aligned} {\mathbb {P}}(N(S)=n)=e^{-\lambda |S|}\frac{(\lambda |S|)^n}{n!} \end{aligned}$$

for any $n\in {{\,\mathrm{{\mathbb {N}}}\,}}$, and that for any $x>0$ we have that

$$\begin{aligned} e^{-x}\sum _{k\ge n} \frac{x^k}{k!}=\frac{x^n}{n!} e^{-x}\sum _{k\ge 0} \frac{x^kn!}{(n+k)!}\le \frac{x^n}{n!}e^{-x}\sum _{k\ge 0} \frac{x^k}{k!}=\frac{x^n}{n!}. \end{aligned}$$

(8)

Choosing now $N_1=1+\lceil \frac{1}{\delta }\rceil$ in (7), we have ${\mathbb {P}}(B_l)\le C\, 2^{-ql}$ for some $q>0$, meaning that $\sum _{l\ge 0} {\mathbb {P}}(B_l)\le \sum _{l\ge 0}2^{-ql}<\infty$. The Borel–Cantelli lemma now implies that

$$\begin{aligned} {\mathbb {P}}\big (\limsup _{l\rightarrow \infty }B_l\big )=0, \end{aligned}$$

meaning that almost surely there is an $\varepsilon _0(\omega )>0$ such that for all $0<\varepsilon _l=2^{-l}\le \varepsilon _0$, any cube $Q_{2^{-l},z}$ contains not more than $N_1$ points from $\frac{\varepsilon }{2} \Phi \cap D$, thus proving the claim.

To show the general case, for $\omega \in \Omega$ we consider $\varepsilon _0(\omega )$ coming from the claim. Without loss of generality we assume $\varepsilon _0 = 2^{-l_0}$ for some $l_0 \in {{\,\mathrm{{\mathbb {N}}}\,}}$ (otherwise replace $\varepsilon _0$ with the largest smaller power of 2). To finish the proof we need to show that for any $0 < \varepsilon \leqq \varepsilon _0$, any cube $Q_\varepsilon = x + [0,\varepsilon ^{1+\delta }]^d$, $x \in {{\,\mathrm{{\mathbb {R}}}\,}}^d$, contains at most N points from $D \cap \varepsilon \Phi (\omega )$. Let $0< \varepsilon < \varepsilon _0$ and $Q_\varepsilon = x + [0,\varepsilon ^{1+\delta }]^d$ be chosen arbitrary, and let $N:= 2^d N_1$. Let $l \geqq l_0$ be the unique l such that $2^{-(l+1)} \leqq \varepsilon < 2^{-l}$.

Observe that for $\lambda > 0$ we have $\# (Q_\varepsilon \cap \varepsilon \Phi ) = \# (\lambda Q_\varepsilon \cap \lambda \varepsilon \Phi )$, where $\lambda Q_\varepsilon = \{ \lambda x: x \in Q_\varepsilon \}$, which together with that star-shapedness of D yields, for $\lambda = \frac{2^{-(l+1)}}{\varepsilon } \in (0,1]$, that

$$\begin{aligned} \# (Q_\varepsilon \cap \varepsilon \Phi \cap D) = \# (\lambda Q_\varepsilon \cap \lambda \varepsilon \Phi \cap \lambda D) \leqq \# (\lambda Q_\varepsilon \cap 2^{-(l+1)} \Phi \cap D). \end{aligned}$$

We now cover $\lambda Q_\varepsilon$ with (at most) $2^d$ cubes $Q_{2^{-l},z}$. Observe that even if $\lambda Q_\varepsilon$ is closed and $Q_{2^{-l},z}$ are only half-closed, the covering is possible since $\lambda \varepsilon = 2^{-(l+1)} < 2^{-l}$. In particular, the claim implies that any $Q_{2^{-l},z}$ contains at most $N_1$ points from $\frac{2^{-l}}{2}\Phi \cap D$, thus implying that $\lambda Q_\varepsilon$, being covered by at most $2^d$ cubes $Q_{2^{-l},z}$, contains at most $2^d N_1$ points from $2^{-(l+1)}\Phi \cap D$. This together with the last display implies $\# (Q_\varepsilon \cap \varepsilon \Phi \cap D) \leqq 2^d N_1 = N$, thus concluding the proof of the proposition.

The first part of Theorem 3 is based on the following Strong Law of Large Numbers, a proof of which can be found in [22, Lemma 6.1] (see also [26, Theorem 8.14]):

Lemma 3

Let $d\ge 1$ and $(\Phi ,{\mathcal {R}})=(\{z_j\},\{r_j\})$ be a marked Poisson point process with intensity $\lambda >0$. Assume that the marks $\{r_j\}$ are non-negative independent identically distributed random variables independent of $\Phi$ such that ${{{\,\mathrm{{\mathbb {E}}}\,}}(r_j^m) <\infty }$ for some $m > 0$. Then, for every bounded set $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$ which is star-shaped with respect to the origin, we have, almost surely,

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0} \varepsilon ^d N(\varepsilon ^{-1} S)=\lambda |S|, \qquad \lim _{\varepsilon \rightarrow 0} \varepsilon ^d \sum _{z_j \in \varepsilon ^{-1} S} r_j^m=\lambda {{\,\mathrm{{\mathbb {E}}}\,}}( r^m)|S|. \end{aligned}$$

Remark 4

Assuming the boundary of the set S from the previous lemma is not too large, the same argument also shows

$$\begin{aligned} \lim _{\varepsilon \rightarrow 0} \varepsilon ^d \sum _{z_j \in \Phi ^{\varepsilon }(S)} r_j^m=\lambda {{\,\mathrm{{\mathbb {E}}}\,}}( r^m)|S|. \end{aligned}$$

(9)

In particular, it is enough that S has as D a $C^2$-boundary.

Using this remark as well as Proposition 1, we can prove Theorem 3.

Proof of Theorem 3

Part (1): We start with the first part of the theorem, which actually holds for any dimension $d \ge 1$, $\alpha > 2$, $m > \frac{d}{\alpha - 2}$, and $\kappa \in (1,\alpha - 1 - \frac{d}{m})$.

Using (9) and the choice of $\kappa$, we have, for almost all $\omega$ that

$$\begin{aligned} \limsup _{\varepsilon \rightarrow 0} \varepsilon ^{\frac{d}{m}} \max _{z_i\in \Phi ^\varepsilon (D)} r_i \le \limsup _{\varepsilon \rightarrow 0} \varepsilon ^\frac{d}{m} \biggl ( \sum _{z_i\in \Phi ^\varepsilon (D)} r_i^m \biggr )^{\frac{1}{m}} \le [\lambda {{\,\mathrm{{\mathbb {E}}}\,}}(r^m)|D|]^{\frac{1}{m}}. \end{aligned}$$

This implies for $\varepsilon >0$ small enough that

$$\begin{aligned} \max _{z_i\in \Phi ^\varepsilon (D)} \tau \varepsilon ^\alpha r_i \leqq 2\tau \varepsilon ^{\alpha -\frac{d}{m}}[\lambda {{\,\mathrm {{\mathbb {E}}}\,}}(r^m)|D|]^{\frac{1}{m}} \leqq \varepsilon ^{1+\kappa }, \end{aligned}$$

(10)

with the last inequality coming from $\alpha - \frac{d}{m} > \kappa +1$, and therefore being true for $\varepsilon$ possibly even smaller.

To show two balls do not intersect we consider an event

$$\begin{aligned} A_\tau ^\varepsilon :=\{\omega \in \Omega : \text { there are 2 intersecting balls in } \{ B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z)\}_{z \in \Phi ^\varepsilon (D)} \}, \end{aligned}$$

and it is enough to show that

$$\begin{aligned} {\mathbb {P}}\bigg (\bigcap _{\varepsilon _0>0}\bigcup _{\varepsilon \le \varepsilon _0} A_\tau ^\varepsilon \bigg )=0. \end{aligned}$$

(11)

We reduce this to the case of dyadic $\varepsilon$, by showing that

$$\begin{aligned} {\mathbb {P}}\bigg (\bigcap _{l_0\ge 1}\bigcup _{l\ge l_0} A_{{\bar{\tau }}}^{\varepsilon _l}\bigg )=0, \end{aligned}$$

(12)

where $\varepsilon _l=2^{-l}$ and ${\bar{\tau }}=2^{1+\kappa }\tau$.

Indeed, let $l\in {{\,\mathrm{{\mathbb {N}}}\,}}$ be such that $\varepsilon _{l+1}\le \varepsilon <\varepsilon _l$. Now suppose $z_i,z_j\in \Phi ^\varepsilon (D)$, $z_i \ne z_j$ such that

$$\begin{aligned} B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z_i)\cap B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z_j)\ne \emptyset . \end{aligned}$$

Then

$$\begin{aligned} \varepsilon _{l+1}|z_i-z_j|\le \varepsilon |z_i-z_j|\le 2\tau \varepsilon ^{1+\kappa }\le 2\tau \varepsilon _l^{1+\kappa } = 2 \tau (2\varepsilon _{l+1})^{1+\kappa } = 2\cdot 2^{1+\kappa }\tau \varepsilon _{l+1}^{1+\kappa }, \end{aligned}$$

which means that

$$\begin{aligned} B_{2^{1+\kappa }\tau \varepsilon _{l+1}^{1+\kappa }}(\varepsilon _{l+1} z_i)\cap B_{2^{1+\kappa }\tau \varepsilon _{l+1}^{1+\kappa }}(\varepsilon _{l+1} z_j)\ne \emptyset . \end{aligned}$$

The domain D being star-shaped implies monotonicity of $\Phi ^{\varepsilon }(D)$ in $\varepsilon$, in particular $\Phi ^\varepsilon (D)\subset \Phi ^{\varepsilon _{l+1}}(D)$, which combined with the previous display yields

$$\begin{aligned} A_\tau ^\varepsilon \subset A_{{\bar{\tau }}}^{\varepsilon _{l+1}}, \end{aligned}$$

thus showing that (12) implies (11).

It remains to show (12). Let $\varepsilon > 0$ and $\tau \geqq 1$ be fixed. Observe that if for $z_i,z_j\in \Phi ^\varepsilon (D)$ we have $B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z_i)\cap B_{\tau \varepsilon ^{1+\kappa }}(\varepsilon z_j)\ne \emptyset$, then $\varepsilon |z_i-z_j|\le 2\tau \varepsilon ^{1+\kappa }$ and after simplifying $|z_i-z_j|\le 2\tau \varepsilon ^\kappa$. In other words,

$$\begin{aligned} A_\tau ^\varepsilon \subset \{\omega \in \Omega : \exists x\in \frac{1}{\varepsilon }D:\#(\Phi ^\varepsilon (D)\cap B_{2\tau \varepsilon ^\kappa }(x))\ge 2\}. \end{aligned}$$

(13)

Recall that for $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$, we denote by $N(S)=\#(S\cap \Phi )$ the random variable providing the number of points of the process which lie inside S. Let us also note that the points are distributed according to a Poisson distribution with intensity $\lambda >0$. We now recall a basic estimate from [20, Proof of Lemma 6.1]: For $0<\eta <1$, define the set of cubes with side length $\eta$ centered at the grid $\eta {{\,\mathrm{{\mathbb {Z}}}\,}}^d$ by

$$\begin{aligned} {\mathcal {Q}}_\eta :=\{y+[-\eta /2,\eta /2]^d:y\in \eta {{\,\mathrm{{\mathbb {Z}}}\,}}^d\}. \end{aligned}$$

Since it is not true that any ball of radius $\frac{\eta }{4}$ is contained in one of these cubes, we need to add (finitely many) shifted copies of ${\mathcal {Q}}_\eta$. For that let $S_\eta$ be the vertices of the cube $[0,\eta /2]^d$, that is,

$$\begin{aligned} S_\eta =\{z=(z_1,\dots ,z_d)\in {{\,\mathrm{{\mathbb {R}}}\,}}^d: z_k\in \{0,\eta /2\} \text { for } k=1,\dots ,d\}. \end{aligned}$$

Observe that for any $x\in {{\,\mathrm{{\mathbb {R}}}\,}}^d$, there exist $z\in S_\eta$ and a cube $Q\in {\mathcal {Q}}_\eta$ such that $B_\frac{\eta }{4}(x)\subset z+Q$, which immediately implies

$$\begin{aligned} {\mathbb {P}}(\exists x\in \frac{1}{\varepsilon }D&: N(B_\frac{\eta }{4}(x))\geqq 2)\\ {}&\leqq {\mathbb {P}}(\exists Q\in {\mathcal {Q}}_\eta , z\in S_\eta : (z+Q)\cap \frac{1}{\varepsilon }D\ne \emptyset , N(z+Q)\geqq 2). \end{aligned}$$

Since $S_\eta$ has $2^d$ elements and the number of cubes $Q\in {\mathcal {Q}}_\eta$ that intersect $\varepsilon ^{-1} D$ is bounded by $C(D)(\varepsilon \eta )^{-d}$, we use the distribution of Poisson point process to conclude that

$$\begin{aligned} {\mathbb {P}}(\exists x\in \frac{1}{\varepsilon }D: N(B_\frac{\eta }{4}(x))\geqq 2)&\leqq \sum _{z \in S_\eta } \sum _{ Q } {\mathbb {P}}( N(z+Q) \geqq 2) \\ {}&\leqq 2^d C(D)(\varepsilon \eta )^{-d}e^{-\lambda \eta ^d}\sum _{k=2}^\infty \frac{(\lambda \eta ^d)^k}{k!} \\ {}&\leqq C(D) 2^d(\varepsilon \eta )^{-d}(\lambda \eta ^d)^2,\end{aligned}$$

where the last inequality follows from (8). Letting $\eta _\varepsilon :=8\tau \varepsilon ^\kappa$, this together with (13) and the fact that $\#(\Phi ^\varepsilon (D)\cap S)\le N(S)$ for any $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$, yields

$$\begin{aligned} {\mathbb {P}}(A_\tau ^\varepsilon )\le C(\varepsilon ^{1+\kappa })^{-d} \varepsilon ^{2d\kappa }=C\varepsilon ^{d(\kappa -1)}. \end{aligned}$$

To show (12) we take a sum over l with $\varepsilon =\varepsilon _l = 2^{-l}$, which using $\kappa >1$ can be estimated as

$$\begin{aligned} \sum _{l=0}^\infty {\mathbb {P}}(A_{{\bar{\tau }}}^{\varepsilon _l})\le C\sum _{l=0}^\infty 2^{-l d (\kappa -1)}<\infty , \end{aligned}$$

and (12) follows from direct application of the Borel–Cantelli lemma.

Part (2): We now turn to the second part of the theorem, that is, the construction of boxes $I^\varepsilon _i$. Fixing $\varepsilon$, the first step is to construct a finite collection ${{\,\mathrm{{\mathcal {I}}}\,}}= \{ {{\tilde{I}}}_i \}$ of auxiliary boxes such that

• these boxes cover the points, that is, $\bigcup _{i} {{\tilde{I}}}_i \supset \varepsilon \Phi ^{\varepsilon }(D)$,

• ${{\,\text {dist}\,}}_{\infty }({{\tilde{I}}}_i,{{\tilde{I}}}_j) \geqq \frac{1}{2N} \varepsilon ^{1+\delta }$,

• $s({{\tilde{I}}}_i) \leqq \varepsilon ^{1+\delta }$, where s(I) of a box I denotes the size of its longest side,

• each box ${{\tilde{I}}}_i$ satisfies $|{{\tilde{I}}}_i \cap \varepsilon \Phi ^\varepsilon (D)| \leqq N$.

Here the crucial condition is the second one, that is, that the boxes are well-separated. Let $l = \frac{1}{2N}\varepsilon ^{1+\delta }$. We will grow the boxes ${{\tilde{I}}}$ from the collection ${{\,\mathrm{{\mathcal {I}}}\,}}$ step by step, starting with cubes of side length l. At every moment of this growth process, every box ${{\tilde{I}}} \in {{\,\mathrm{{\mathcal {I}}}\,}}$ will satisfy the following conditions:

• i. ${{\tilde{I}}} = [a_xl,b_xl) \times [a_yl,b_yl) \times [a_zl,b_zl)$ for some $a_x,b_x,a_y,b_y,a_z,b_z \in {\mathbb {Z}}$, that is, each box is a union of many small cubes;

• ii. for each $a \in [a_x,b_x) \cap {\mathbb {Z}}$ holds $[al,(a+1)l) \times [a_yl,b_yl) \times [a_zl,b_zl) \cap \varepsilon \Phi ^\varepsilon (D) \ne \emptyset$, and similarly for y and z, that is, in every slice there is some point from $\varepsilon \Phi ^\varepsilon (D)$;

• iii. $\#({{\tilde{I}}} \cap \varepsilon \Phi ^\varepsilon (D)) \leqq N$.

At the beginning, let ${{\,\mathrm{{\mathcal {I}}}\,}}$ consist of all cubes $[a_xl,(a_x+1)l)\times [a_yl,(a_y+1)l)\times [a_zl,(a_z+1)l)$ which have a point from $\varepsilon \Phi ^\varepsilon (D)$ in it. Since D is bounded, ${{\,\mathrm{{\mathcal {I}}}\,}}$ consists of finitely many boxes (cubes). We then repeat the following procedure:

If there exist two different boxes ${{\tilde{I}}}, {{\tilde{J}}} \in {{\,\mathrm{{\mathcal {I}}}\,}}$ such that ${{\,\textrm{dist}\,}}({{\tilde{I}}},{{\tilde{J}}}) = 0$, we fix them and merge them together. That means, we remove ${{\tilde{I}}} = [a_xl,b_xl) \times [a_yl,b_yl) \times [a_zl,b_zl)$ and ${{\tilde{J}}} = [a_x'l,b_x'l) \times [a_y'l,b_y'l) \times [a_z'l,b_z'l)$ from ${{\,\mathrm{{\mathcal {I}}}\,}}$ and add

$$\begin{aligned} {{\tilde{K}}}&= [A_xl,B_xl) \times [A_yl,B_yl) \times [A_zl,B_zl)\\ :\!\!&= [(a_x \wedge a_x')l,(b_x \vee b_x')l) \times [(a_y \wedge a_y')l,(b_y \vee b_y')l) \times [(a_z \wedge a_z')l,(b_z \vee b_z')l) \end{aligned}$$

to ${{\,\mathrm{{\mathcal {I}}}\,}}$ instead. Here $\wedge$ and $\vee$ stand as usual for minimum and maximum, respectively.

First, observe that (i) trivially follows from the definition of ${{\tilde{K}}}$. Next, to verify that ${{\tilde{K}}}$ satisfies (ii), let us fix $i \in \{ x,y,z \}$, and observe that ${{\,\textrm{dist}\,}}({{\tilde{I}}},{{\tilde{J}}}) = 0$ implies $[a_i,b_i] \cap [a_i',b_i'] \ne \emptyset$. Hence, for any ${a \in [\min (a_i,a_i'),\max (b_i,b_i'))}$ either $a \in [a_i,b_i)$, in which case (ii) for ${{\tilde{I}}}$ implies (ii) for ${{\tilde{K}}}$, or $a \in [a_i',b_i')$, in which case (ii) for ${{\tilde{J}}}$ implies (ii) for ${{\tilde{K}}}$.

It remains to argue that ${{\tilde{K}}}$ satisfies also (iii). Since ${{\tilde{I}}}$ satisfies both (ii) and (iii), in particular to each $a \in [a_i,b_i)\cap {{\,\mathrm{{\mathbb {Z}}}\,}}$ there is assigned at least one point from $\varepsilon \Phi ^\varepsilon (D)$ and there are at most N such points, it follows that $[a_i l,b_i l)$ has length at most Nl. The same argument applies verbatim to ${{\tilde{J}}}$, and so the union of $[a_i l,b_i l)$ and $[a_i'l,b_i'l)$ has length at most 2Nl. Hence, each side of ${{\tilde{K}}}$ has length at most $s(\tilde{K})\le 2N l = 2N \frac{1}{2N} \varepsilon ^{1+\delta } = \varepsilon ^{1+\delta }$. In addition ${{\tilde{K}}}$ satisfies (i), and so there exists a (closed) cube $Q_{{{\tilde{K}}}}$ with side length $\varepsilon ^{1+\delta }$ such that ${{\tilde{K}}} \subset Q_{{{\tilde{K}}}}$. By Proposition 1, the number of points in $Q_{{{\tilde{K}}}}$ is at most N, which implies the same for ${{\tilde{K}}}$, that is,

$$\begin{aligned} \#({{\tilde{K}}} \cap \varepsilon \Phi ^\varepsilon (D)) \leqq \#(Q_{{{\tilde{K}}}} \cap \varepsilon \Phi ^\varepsilon (D)) \leqq N, \end{aligned}$$

which shows (iii) for ${{\tilde{K}}}$; moreover, since ${{\tilde{K}}}$ also fulfils (ii), this shows that ${{\tilde{K}}}$ has length at most $s({{\tilde{K}}})\le Nl$.

Since the collection ${{\,\mathrm{{\mathcal {I}}}\,}}$ was finite at the beginning, and in each iteration we decrease the number of boxes in ${{\,\mathrm{{\mathcal {I}}}\,}}$ by one (we remove ${{\tilde{I}}}$ and ${{\tilde{J}}}$ and add ${{\tilde{K}}}$), this process has to terminate. In particular, at the end ${{\,\mathrm{{\mathcal {I}}}\,}}$ consists of boxes which have positive distance from each other, since otherwise the process would not terminate at this point. Since all boxes in ${{\,\mathrm{{\mathcal {I}}}\,}}$ satisfy (i), this in particular implies that this positive distance has to be at least $l = \frac{1}{2N} \varepsilon ^{1+\delta }$. Moreover, since each box has side length at most $Nl = \frac{1}{2} \varepsilon ^{1+\delta }$, and each point in $\varepsilon \Phi ^{\varepsilon }(D)$ is at least distance $\varepsilon$ to $\partial D$, we see that each box (and actually also its small neighborhood) lies inside D.

Using boxes from ${{\,\mathrm{{\mathcal {I}}}\,}}$ we define boxes $I^\varepsilon _i$: for each auxiliary box ${{\tilde{I}}}_i \in {{\,\mathrm{{\mathcal {I}}}\,}}$ set $I^\varepsilon _i:= \{ x \in {{\,\mathrm {{\mathbb {R}}}\,}}^3: {{\,\text {dist}\,}}_\infty (x,{{\tilde{I}}}_i) \leqq \frac{1}{8N} \varepsilon ^{1+\delta } \}$, and it remains to show that $\{ I^\varepsilon _i \}$ satisfy (2a)–(2e). First, by the assumption $\kappa > \delta$, and so for small enough $\varepsilon$ we have $\varepsilon ^{1+\kappa } \leqq \frac{1}{16N} \varepsilon ^{1+\delta }$. Therefore, by the triangle inequality we have for any $\varepsilon z \in {{\tilde{I}}}_i$ that ${{\,\text {dist}\,}}_{\infty } (B_{\varepsilon ^{1+\kappa }}(\varepsilon z),\partial I^\varepsilon _i) \geqq \frac{1}{8N}\varepsilon ^{1+\delta } - \varepsilon ^{1+\kappa } \geqq \frac{1}{16N}\varepsilon ^{1+\delta }$, thus (2a) and (2c) hold. Since by the construction the auxiliary boxes satisfy ${{\,\text {dist}\,}}_\infty ({{\tilde{I}}}_i,{{\tilde{I}}}_j) \geqq \frac{1}{2N}\varepsilon ^{1+\delta }$, and all the points from $\varepsilon \Phi ^{\varepsilon }(D)$ are inside these boxes, we see that $I^\varepsilon _i {\setminus } {{\tilde{I}}}_i$ contains no point from $\varepsilon \Phi ^\varepsilon (D)$. Therefore (iii) for ${{\tilde{I}}}_i \in {{\,\mathrm{{\mathcal {I}}}\,}}$ implies (2b) for $I^\varepsilon _i$. Finally, (2d) trivially follows from the definition of $I^\varepsilon _i$ and the separation of elements in ${{\,\mathrm{{\mathcal {I}}}\,}}$ in form of ${{\,\text {dist}\,}}_\infty ({{\tilde{I}}}_i,{{\tilde{I}}}_j) \geqq \frac{1}{2N} \varepsilon ^{1+\delta }$, and (2e) uses that $I^\varepsilon _i$ consist in each direction of at least one cube and of at most N of them.

5 Proofs of Lemma 2 and Theorem 1

Before proving that a box from which we remove finitely many small well-separated balls is a John domain, let us recall how John domains are defined.

For a constant $c>0$, a domain $U\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$ is said to be a c-John domain if there exists a point $x_0\in U$ such that for any point $x\in U$ there is a rectifiable path $\Gamma :[0,\ell ]\rightarrow U$ which is parametrized by arc length with

$$\begin{aligned} \Gamma (0)=x,\quad \Gamma (\ell )=x_0,\quad \forall t\in [0,\ell ]: |\Gamma (t)-x|\le c\, {{\,\textrm{dist}\,}}(\Gamma (t),\partial U). \end{aligned}$$

(14)

John domains may have fractal boundaries or internal cusps, whereas external cusps are forbidden. For instance, the interior of Koch’s snowflake as well as any convex domain are John domains. In the case of bounded domains, there are several equivalent definitions of John domains, see [37, Section 2.17]. We state the following characterisation, which is used in [10, Section 3.1]: A bounded domain U is a c-John domain in the sense of Definition 1 if and only if there is a $c_1(c)>0$ and a point $x_0\in U$ such that any point $x\in U$ can be connected to $x_0$ by a rectifiable path $\Gamma :[0,\ell ]\rightarrow U$ which is parametrized by arc length and

$$\begin{aligned} \bigcup \limits _{t\in [0,\ell ]} B\big (\Gamma (t), t/c_1\big )\subset U. \end{aligned}$$

One way how to prove Lemma 2 is inductively by showing, that under some assumption on a ball one can remove it from a John domain while changing the John constant of the domain by at most a fixed factor—for that we would need to modify arcs which run close to (or through) this removed ball while estimating how much does this change the situation. For a similar argument with small balls replaced with points, see [24, Theorem 1.4]. Assuming this, since we have to remove at most N balls and at the beginning the domain is rectangle with proportional sides, in particular a John domain, this would lead to the conclusion.

Instead of this we provide a direct constructing argument.

Proof of Lemma 2

To start, we use Theorem 3, part (1), twice: once with $\kappa = \kappa _1:= 1+\delta$ and second time with $\kappa = \kappa _2:= \alpha -1-\frac{3}{m}-\delta$. Observe that both values of $\kappa$ are within the admissible range $(\max (1,\delta ),\alpha -1-\frac{3}{m})$, and therefore the theorem yields the following:

there exists an almost surely positive $\varepsilon _0(\omega )$, obtained as the smaller of the two $\varepsilon _0$, such that, for $0 < \varepsilon \leqq \varepsilon _0$, it holds that

$$\begin{aligned} \max _{z_j \in \Phi ^{\varepsilon }(D)} \varepsilon ^\alpha r_j \leqq \varepsilon ^{1+\kappa _2} \,\,\text{ and } \,\, |\varepsilon z_j - \varepsilon z_k| \geqq 2\varepsilon ^{1+\kappa _1} \quad \text{ for } \text{ any }\,z_j, z_k \in \Phi ^\varepsilon (D).\qquad \end{aligned}$$

(15)

Assume we have $0 < \varepsilon \leqq \varepsilon _0$ small enough and recall that we want to show that

$$\begin{aligned} U:=I^\varepsilon _i \setminus \bigcup _k B_k \end{aligned}$$

is a $c(N)-$John domain in the sense of Definition 1, where $\{ B_k \}_k = \{ B_{\varepsilon ^\alpha r_k}(\varepsilon z_k): \varepsilon z_k \in I^{\varepsilon }_i \}$. For brevity we set $I:= I^{\varepsilon }_i = p + (-l_1/2,l_1/2)\times (-l_2/2,l_2/2)\times (-l_3/2,l_3/2) \subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$, where p is the center and $l_i$ are the side lengths of I. Since (14) is scale-invariant,we can assume $l_1 \geqq l_2 \geqq l_3$. The set

$$\begin{aligned} L:= \{ x \in I: {{\,\textrm{dist}\,}}_\infty (x,\partial I) = \frac{1}{32N}\varepsilon ^{1+\delta }\} \end{aligned}$$

will serve as a “highway” in the set U, and for the specific point $x_0$ from Definition 1 we choose $x_0:= p + (0,0,l_3/2-\frac{1}{32N}\varepsilon ^{1+\delta })$. We also denote the ring around L by $G:= \{ x \in I: {{\,\textrm{dist}\,}}_{\infty }(x,L) < \frac{1}{32N}\varepsilon ^{1+\delta }\}$.

To show that U is a John domain, for each $X \in U$ we need to construct a path from X to $x_0$ along which $|\cdot - X| \leqq c {{\,\text {dist}\,}}(\cdot ,\partial U)$. The idea is first to go from X to L, and then run along L to $x_0$. Observe that for points $x \in L$ the condition is easy to satisfy: for each $x \in L$ we have ${{\,\textrm{dist}\,}}(x,\partial U) = {{\,\textrm{dist}\,}}(x,\partial G) = \frac{1}{32N}\varepsilon ^{1+\delta }$ and $|x-X| \leqq {{\,\text {diam}\,}}(U) \le \sqrt{3}\, l_1$, and so using (2e) to see $l_1 \leqq C(N)\varepsilon ^{1+\delta }$ we get that $|x-X| \leqq c(N) {{\,\text {dist}\,}}(x,\partial U)$ as required.

It remains to describe the path from X to L. For points $X \in G$ this is straightforward (see Fig. 1a): we just choose the shortest path from X to L and observe that any point x on that path satisfies ${{\,\text {dist}\,}}(x,\partial U) \geqq {{\,\text {dist}\,}}(x,\partial G) \ge 3^{-1/2} |X-x|$. The $\sqrt{3}$ is optimal as can be seen from points in corners.

Fig. 1

a The point $X\in G$, first connected to $x_1\in L$ (red) and then to $x_0$ while not leaving L (blue). b The projections (blue) of the balls $B_1$ and $B_k$ onto the sphere S with midpoint X. The cone C illustrated by the red area hits none of the balls and serves as the “outgoing” sector from X to L

In the rest of the proof we deal with the points from the “interior” $U \setminus G$. For $X \in U \setminus G$ we need to construct a path from X to L, while not going too close near the balls $B_k$. We will use two important properties of these balls: the size of the balls is much smaller than their mutual distance (see (15)), and there are at most N of them. We fix $X \in U{\setminus } G$ and show that we can actually use a straight line to connect X with L. Along this line we should be able to move a growing ball without hitting $\{B_k\}$, what is equivalent to an existence of a cone with an opening c(N) which avoids all the balls. For that let S be a unit sphere centered at X, and let ${\mathcal {P}}$ denote the orthogonal projection on S. Farther, let $P:= {\mathcal {P}}(\bigcup _k B_k)$ denote the projection of balls on S. Observe that if we find a disc on S of fixed radius (depending on N) which does not overlap with P, then we are done since such disc corresponds to a cone at X avoiding the balls $B_k$ (see Fig. 1b).

Hence, we reduced our task to a problem of finding a not too small disc in $S {\setminus } P$, with P being a union of at most N discs with some additional properties. First, it can happen that X lies very close to one of the balls, so that the projection of this particular ball on S covers (almost) half of the sphere S. For this reason w.l.o.g. let $B_1$ denote the ball whose center is closest to X, which we treat separately: let $S' \subset S$ be a half-sphere with the pole being exactly opposite to the center of ${\mathcal {P}}(B_1)$, in particular ${\mathcal {P}}(B_1)$ and $S'$ are disjoint. Since $B_1$ was the closest ball to X, it follows from the second estimate in (15) that X is at least $\varepsilon ^{1+\kappa _1}$ away from centers of the remaining balls $B_k, k \geqq 2$. On the other hand, the first relation in (15) bounds the radii of these balls with $\varepsilon ^{1+\kappa _2}$. Therefore, the projections of these remaining balls are discs of radius at most $C\frac{\varepsilon ^{1+\kappa _2}}{\varepsilon ^{1+\kappa _1}} = C\varepsilon ^{\kappa _2 - \kappa _1}$. Since $\kappa _2 - \kappa _1 = \alpha -2-\frac{3}{m}-2\delta >0$ by the choice of $\delta$ in (5), we see that for $\varepsilon$ small these (at most $N-1$) projections are tiny discs (almost points). We can now find a radius $r=r(N)$ with the following property: there exist N discs $D_1,\ldots ,D_N$ of radius r in $S'$ such that the distance between any two discs is at least r as well. One option is to arrange them along the boundary of $S'$ with necessary spacing between them, thus achieving $r \sim N^{-1}$. Provided now $\varepsilon$ is small enough so that the radii of ${\mathcal {P}}(B_k)$, which are bounded by $C \varepsilon ^{\kappa _2 - \kappa _1}$, are smaller than r, we are done: there are N discs $D_1,\ldots ,D_N$ and at most $N-1$ projections ${\mathcal {P}}(B_k)$ where each projection can touch at most one $D_k$, so that one disc will not overlap with any of the projections ${\mathcal {P}}(B_k)$, thus defining the cone we are searching for.

This solution to the last question is naturally far from optimal (in r): consider a well-studied question of finding an optimal cover of a sphere (more precisely half of it) with N identical discs of smallest radius. If $\varrho$ denotes the smallest such radius, then for any configuration of $N-1$ points in $S'$ there exists a disc in $S'$ of radius $\varrho$ which avoids them, thus also providing a solution to our problem [36].

Since the perforated boxes U from Lemma 2 are uniform John domains, in particular we have a Bogovkiĭ operator on each U, Theorem 1 can be proved along the lines of the proof in [9]. First, using a Bogovskiĭ operator on the whole D we obtain a function $\textbf{u}$ with the correct divergence, but which naturally does not vanish on the holes. To achieve that, we modify $\textbf{u}$ in each box $I^{\varepsilon }_i$. More precisely, near $\partial I^\varepsilon _i$ in a boundary layer of size $\frac{1}{16N}\varepsilon ^{1+\delta }$ we change $\textbf{u}$ to its average value over this layer, and then inside the box (where also the balls are removed) cut off this constant function near each hole over a scale $\varepsilon ^\alpha$. Since by this modification we also change the divergence of the function, we employ Bogovskiĭ’s operator both on each box as well as near each hole to fix the divergence.

Proof of Theorem 1

Let us recall definition of $D_\varepsilon =D{\setminus }\bigcup _{z_j\in \Phi ^\varepsilon (D)} B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)$ (see (1)). To prove the theorem we construct a linear operator of Bogovskiĭ type, bounded independently of $\varepsilon$:

$$\begin{aligned} {\mathcal {B}}_\varepsilon :L_0^q(D_\varepsilon )\rightarrow W_0^{1,q}(D_\varepsilon ), \end{aligned}$$

satisfying

$$\begin{aligned} {{\,\textrm{div}\,}}{\mathcal {B}}_\varepsilon (f)=f \text { in } D_\varepsilon ,\quad ||{\mathcal {B}}_\varepsilon (f)||_{W_0^{1,q}(D_\varepsilon )}\le C\, ||f||_{L_0^q(D_\varepsilon )}. \end{aligned}$$

(16)

To this end, we will first use a Bogovskiĭ operator on the whole domain D and then correct this function first to its mean value over a large scale and then to zero near each hole without changing the divergence. We will give some remarks on this procedure after the proof.

For $1<q<\infty$ and $f\in L_0^q(D_\varepsilon )$ we denote by ${\tilde{f}}\in L_0^q(D)$ its zero extension in the holes. Using classical Bogovskiĭ’s operator in Lipschitz domain D [17, Chapter 3], the norm of which depends on the Lipschitz character of D, we can find a function $\textbf{u}={\mathcal {B}}_D({\tilde{f}})\in W_0^{1,q}(D)$ satisfying

$$\begin{aligned} {{\,\textrm{div}\,}}\textbf{u}={\tilde{f}}\text { in } D,\quad \Vert \textbf{u}\Vert _{W_0^{1,q}(D)}\le C\, \Vert {\tilde{f}}\Vert _{L_0^q(D)}=C\, \Vert f\Vert _{L_0^q(D_\varepsilon )}, \end{aligned}$$

with $C=C(D,q)$.

Since $\alpha -3/m>2$, by applying Theorem 3 we find for every $\varepsilon > 0$ small enough a finite collection of boxes $I^{\varepsilon }_i$ such that for any point $z_j \in \Phi ^\varepsilon (D)$ there is i such that

$$\begin{aligned} B_{\varepsilon ^{\alpha } r_j }(\varepsilon z_j)\subset B_{2 \varepsilon ^{\alpha } r_j}(\varepsilon z_j)\subset B_{\varepsilon ^{1+\kappa }}(\varepsilon z_j) \subset I^{\varepsilon ,\text {in}}_i, \end{aligned}$$

where

$$\begin{aligned} I^{\varepsilon ,\text{ in }}_i&:= \{ x \in I^{\varepsilon }_i : {{\,\text {dist}\,}}_\infty (x,\partial I^\varepsilon _i) \geqq \frac{1}{16N}\varepsilon ^{1+\delta } \}. \end{aligned}$$

For any box $I^\varepsilon _i$ and any ball $B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)$, consider the following corresponding cut-off functions:

$$\begin{aligned}&\chi _{\varepsilon ,i} \in C_c^\infty (I^\varepsilon _i),\quad \chi _{\varepsilon ,i} \restriction _{I^{\varepsilon ,\text {in}}_i}=1,\quad \Vert \nabla \chi _{\varepsilon ,i}\Vert _{L^\infty (D)}\lesssim \varepsilon ^{-(1+\delta )},\\&\zeta _{\varepsilon ,j}\in C_c^\infty \Bigl (B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j)\Bigr ),\quad \zeta _{\varepsilon ,j}\restriction _{B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)}=1,\quad \Vert \nabla \zeta _{\varepsilon ,j}\Vert _{L^\infty (B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j))}\lesssim \frac{1}{r_j}\varepsilon ^{-\alpha }. \end{aligned}$$

(17)

Furthermore, we denote the mean value of $\textbf{u}$ over a measurable set $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ by

$$\begin{aligned} \langle \textbf{u}\rangle _S:=\frac{1}{|S|}\int _S \textbf{u}, \end{aligned}$$

and set that

$$\begin{aligned} A^{\varepsilon }_i&:= I^{\varepsilon }_i \setminus I^{\varepsilon ,\text {in}}_i = \{ x \in I^{\varepsilon }_i : {{\,\textrm{dist}\,}}_\infty (x,\partial I^{\varepsilon }_i) < \frac{1}{16N}\varepsilon ^{1+\delta } \},\\ \textbf{b}_{\varepsilon ,i}(\textbf{u})&:=\chi _{\varepsilon ,i} (\textbf{u}-\langle \textbf{u}\rangle _{A^\varepsilon _i})\in W_0^{1,q} (I^\varepsilon _i),\\ \beta _{\varepsilon ,j}(\textbf{u})&:=\zeta _{\varepsilon ,j}\, \langle \textbf{u}\rangle _{A^\varepsilon _i}\in W_0^{1,q} \Bigl (B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j)\Bigr ), \end{aligned}$$

where, as before, i and j are related through $\varepsilon z_j \in I^\varepsilon _i$.

Since all the lengths in the set $A^\varepsilon _i$ are proportional to $\varepsilon ^{1+\delta }$ (with the proportionality depending on N), Poincaré’s inequality implies

$$\begin{aligned} \Vert \textbf{u}-\langle \textbf{u}\rangle _{A^\varepsilon _i}\Vert _{L^q(A^\varepsilon _i)} \lesssim \varepsilon ^{1+\delta }\, \Vert \nabla \textbf{u}\Vert _{L^q(A^\varepsilon _i)}, \end{aligned}$$

and by (17) we get that

$$\begin{aligned} \begin{aligned} \Vert \nabla \textbf{b}_{\varepsilon ,i}(\textbf{u})\Vert _{L^q(A^\varepsilon _i)}&\le \Vert \chi _{\varepsilon ,i} \nabla (\textbf{u}-\langle \textbf{u}\rangle _{A_i^\varepsilon })\Vert _{L^q(A_i^\varepsilon )}+\Vert \nabla \chi _\varepsilon (\textbf{u}-\langle \textbf{u}\rangle _{A_i^\varepsilon })\Vert _{L^q(A_i^\varepsilon )}\\&\lesssim \Vert \nabla (\textbf{u}-\langle \textbf{u}\rangle _{A^\varepsilon _i})\Vert _{L^q(A^\varepsilon _i)}+\varepsilon ^{-(1+\delta )}\, \Vert \textbf{u}-\langle \textbf{u}\rangle _{A^\varepsilon _i}\Vert _{L^q(A^\varepsilon _i)}\\&\lesssim \Vert \nabla \textbf{u}\Vert _{L^q(A^\varepsilon _i)}. \end{aligned} \end{aligned}$$

(18)

Similarly, by (18) and Jensen’s inequality, we get that

$$\begin{aligned} \begin{aligned} \Vert \nabla \beta _{\varepsilon ,j}(\textbf{u})\Vert _{L^q(B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j))}&=\Vert \nabla \zeta _{\varepsilon ,j}\cdot \langle \textbf{u}\rangle _{A^\varepsilon _i}\Vert _{L^q(B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j))}\\&\lesssim r_j^{\frac{3}{q}-1}\varepsilon ^{\big (\frac{3}{q}-1\big )\alpha }\, |\langle \textbf{u}\rangle _{A^\varepsilon _i}| \lesssim r_j^{\frac{3}{q}-1}\varepsilon ^{\big (\frac{3}{q}-1\big )\alpha }\, |A^\varepsilon _i|^{-\frac{1}{q}} \Vert \textbf{u}\Vert _{L^q(A^\varepsilon _i)}\\&\lesssim r_j^{\frac{3}{q}-1}\varepsilon ^{\big (\frac{3}{q}-1\big )\alpha -\frac{3(1+\delta )}{q}}\, \Vert \textbf{u}\Vert _{L^q(A^\varepsilon _i)}. \end{aligned} \end{aligned}$$

(19)

Since $B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)\subset D$, we have $r_j\le \varepsilon ^{1+\kappa _2-\alpha }=\varepsilon ^{-(\frac{3}{m}+\delta )}$ by (15) and the choice of $\kappa _2=\alpha -1-\frac{3}{m}-\delta$. This yields

$$\begin{aligned} r_j^{\frac{3}{q}-1}\varepsilon ^{\big (\frac{3}{q}-1\big )\alpha -\frac{3(1+\delta )}{q}}\le \varepsilon ^{(\frac{3}{q}-1)(\alpha -\frac{3}{m}-\delta )-\frac{3}{q} (1+\delta )}. \end{aligned}$$

Thus, for all $1<q<3$ which satisfy (2), we can choose $\delta$ such that

$$\begin{aligned} 0<\delta \le \frac{(3-q)(\alpha -\frac{3}{m})-3}{6-q} \end{aligned}$$

(20)

to get uniform bounds on $\Vert \beta _{\varepsilon ,j}(\textbf{u})\Vert _{L^q(B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j))}$.

Since $\beta _{\varepsilon ,i}$ as well as $\textbf{b}_{\varepsilon ,j}$ do not have vanishing divergence, we need to correct them using Bogovkiĭ operators on perforations of $B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j)$ and $I^\varepsilon _i$, respectively. In the first case one can construct the Bogovkiĭ operator

$$\begin{aligned} {{\widetilde{{\mathcal {B}}}}}_{\varepsilon ,j}:L_0^q\big (B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j)\setminus B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)\big ) \rightarrow W_0^{1,q}\big (B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j)\setminus B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)\big ) \end{aligned}$$

(21)

as in [17, Chapter III and Theorem III.3.1], which mimics the original proof from Bogovskiĭ in [6]. Alternatively, one can also construct it by observing that $B_{2 \varepsilon ^\alpha r_j}(\varepsilon z_j){\setminus } B_{\varepsilon ^\alpha r_j}(\varepsilon z_j)$ is a uniform John domain independent of $\varepsilon$ and $z_j$ and use [10, Theorem 3.8 and Theorem 5.2]. In the second situation, the existence of the Bogovskiĭ operator ${\mathcal {B}}_{\varepsilon ,i}$ for the set $I^\varepsilon _i {\setminus } \bigcup _{\varepsilon z_j \in I^{\varepsilon }_i} B_{\varepsilon ^{\alpha }r_j}(\varepsilon z_j)$ is content of Lemma 2, provided we choose $\delta$ from (21) possibly even smaller to satisfy also (5). We are now ready to define the restriction operator from D to $D_\varepsilon$ via

$$\begin{aligned} R_\varepsilon (\textbf{u})\!:=\!\textbf{u}-\!\!\!\!\!\!\!\sum _{z_j\in \Phi ^\varepsilon (D)}\!\!\!\! (\beta _{\varepsilon ,j}(\textbf{u})-{{\widetilde{{\mathcal {B}}}}}_{\varepsilon ,j}({{\,\textrm{div}\,}}\beta _{\varepsilon ,j}(\textbf{u})))\!-\!\sum _{i} (\textbf{b}_{\varepsilon ,i}(\textbf{u})\!-\!{\mathcal {B}}_{\varepsilon ,i}({{\,\textrm{div}\,}}\textbf{b}_{\varepsilon ,i}(\textbf{u}))), \end{aligned}$$

where the last sum runs over all boxes $I^\varepsilon _i$ and the functions were extended by 0 outside their domain of definition. This definition is essentially the same as in [9]; note that we just replaced their operators ${\mathcal {B}}_{E_{\varepsilon ,\! n}}$ by our operators ${\mathcal {B}}_{\varepsilon ,i}$. Repeating the arguments shown in [9, Section 3], by (19), (20) and the fact that the boxes $I^\varepsilon _i$ are disjoint, we see that $R_\varepsilon (\textbf{u}) \in W_0^{1,q}(D_\varepsilon )$ is well defined and satisfies

$$\begin{aligned} {{\,\textrm{div}\,}}R_\varepsilon (\textbf{u})&=f \text { in } D_\varepsilon ,\quad \Vert R_\varepsilon (\textbf{u})\Vert _{W_0^{1,q}(D_\varepsilon )}\nonumber \\&\le C\, \bigg (\varepsilon ^{(\frac{3}{q}-1)(\alpha -\frac{3}{m}-\delta )-\frac{3}{q} (1+\delta )}+1\bigg )\Vert \textbf{u}\Vert _{W_0^{1,q}(D)}, \end{aligned}$$

(22)

where the constant $C>0$ is independent of $\varepsilon >0$. Note that due to the choice of $\delta$, the exponent of $\varepsilon$ on the right hand-site is non-negative, so we may bound $R_\varepsilon$ uniformly with respect to $\varepsilon$. For $f\in L_0^q(D_\varepsilon )$ we define

$$\begin{aligned} {\mathcal {B}}_\varepsilon (f):=(R_\varepsilon \circ {\mathcal {B}}_D)({\tilde{f}}), \end{aligned}$$

and observe that we get the desired operator, namely ${\mathcal {B}}_\varepsilon (f)\in W_0^{1,q}(D_\varepsilon )$,

$$\begin{aligned} {{\,\textrm{div}\,}}{\mathcal {B}}_\varepsilon (f)=f \text { in } D_\varepsilon ,\quad \text {and}\quad \Vert {\mathcal {B}}_\varepsilon (f)\Vert _{W_0^{1,q}(D_\varepsilon )}\le C\, \Vert f\Vert _{L^q(D_\varepsilon )}. \end{aligned}$$

This finishes the proof of Theorem 1.

Remark 5

As holes are well separated, one might think that the construction of the Bogovskiĭ operator is possible in just two steps: first in the whole domain and second with a cut-off argument near each hole. This construction would follow the one from [3] and its $L^q$-generalization in [29, Section 5]. However, following their proof, one recognizes that we would get a worse exponent of $\varepsilon$: the term $\frac{3}{q} (1+\delta )$ would change to $\frac{3}{q} (2+\delta )$. This is due to the fact that in our case, we do not have that the holes have mutual distance of order $\varepsilon$ but rather (more than) $\varepsilon ^{2+\delta }$ due to the random distribution of centers.

Remark 6

We note that the $\varepsilon$-dependence in (23) seems not to be optimal but “close to optimal” in the sense of capacity (see also [29, Remark 2.4]): Recall that for $1<q<\infty$ and $S\subset {{\,\mathrm{{\mathbb {R}}}\,}}^d$, the q-capacity is defined as

$$\begin{aligned} \text {Cap}_q(S)=\inf \{\Vert f\Vert _{W^{1,q}({{\,\mathrm{{\mathbb {R}}}\,}}^d)}^q: f\in W^{1,q}({{\,\mathrm{{\mathbb {R}}}\,}}^d), S\subset \{f\ge 1\}\}. \end{aligned}$$

We will here focus on the case $d=3$. For a ball of radius $r>0$, it is known that for any $1<q<3$ there exists a constant $C=C(q)>0$ such that

$$\begin{aligned} \text {Cap}_q(B_r(0))=C r^{3-q}. \end{aligned}$$

Since the capacity is an outer measure, the fact that (for $\varepsilon >0$ small enough) the balls are well separated and the expected number of holes inside D is of order $\varepsilon ^{-3}$, together with (10) and the choice of $\kappa _2=\alpha -1-\frac{3}{m}-\delta$ we have that

$$\begin{aligned} \begin{aligned} \text {Cap}_q\bigg (\bigcup _{z_i\in \Phi ^\varepsilon (D)} B_{r_i\varepsilon ^\alpha }(\varepsilon z_i)\bigg )&\le \sum _{z\in \Phi ^\varepsilon (D)} \text {Cap}_q(B_{r_i\varepsilon ^\alpha }(\varepsilon z_i))\\ \le C\varepsilon ^{-3} \big (\max _{z_i\in \Phi ^\varepsilon (D)} r_i\varepsilon ^\alpha \big )^{3-q}&\le C\varepsilon ^{(1+\kappa _2)(3-q)-3}=C\varepsilon ^{(3-q)(\alpha -\frac{3}{m}-\delta )-3}. \end{aligned} \end{aligned}$$

(23)

We see that the essential quantity $(3-q)(\alpha -\frac{3}{m}-\delta )-3$ arising here is almost the same as in (23). A possible explanation for the connection between the capacity estimates and the estimate for the Bogovskiĭ operator is as follows. Let $u\in W_0^{1,q}(D_\varepsilon )$ and $\varphi \in C_c^\infty ({{\,\mathrm{{\mathbb {R}}}\,}}^d)$ with $\varphi =1$ in D. Then $\varphi (1-u)$ is an admissible function in the definition of the q-capacity for the union of all holes, that is,

$$\begin{aligned} \varphi (1-u)\in W^{1,q}({{\,\mathrm{{\mathbb {R}}}\,}}^d),\quad \varphi (1-u)=1 \text { on } \bigcup _{z_i\in \Phi ^\varepsilon (D)} B_{r_i\varepsilon ^\alpha }(\varepsilon z_i). \end{aligned}$$

Direct calculations yield

$$\begin{aligned} \Vert \varphi (1-u)\Vert _{W^{1,q}({{\,\mathrm{{\mathbb {R}}}\,}}^d)}\le C(\varphi )\,(1+\Vert u\Vert _{W^{1,q}(D_\varepsilon )}), \end{aligned}$$

as well as

$$\begin{aligned} C(\varphi )\,(1+\Vert u\Vert _{W^{1,q}(D_\varepsilon )}^q)\ge \Vert \varphi (1-u)\Vert _{W^{1,q}({{\,\mathrm{{\mathbb {R}}}\,}}^d)}^q \ge \text {Cap}_q\bigg (\bigcup _{z_i\in \Phi ^\varepsilon (D)} B_{r_i\varepsilon ^\alpha }(\varepsilon z_i)\bigg ). \end{aligned}$$

If $\alpha$ is large and the radii $r_i$ are almost constant, meaning that the holes inside D should be very well separated, one might expect that the inequality (24) is close to an equality, yielding

$$\begin{aligned} \Vert u\Vert _{W^{1,q}(D_\varepsilon )}^q\ge C\varepsilon ^{(3-q)(\alpha -\frac{3}{m}-\delta )-3}. \end{aligned}$$

For the Bogovskiĭ operator obtained in Theorem 1, we have ${\mathcal {B}}_\varepsilon (f)\in W_0^{1,q}(D_\varepsilon )$, so the optimal general estimate on $\Vert u\Vert _{W^{1,q}(D_\varepsilon )}$ may be of size $\varepsilon ^{(\frac{3}{q}-1)(\alpha -\frac{3}{m}-\delta )-\frac{3}{q}}$. The suboptimal factor $\varepsilon ^{-\frac{3}{q} (1+\delta )}$ in (23) is due to the fact our construction does not enable us to have a better estimate on $\nabla \beta _{\varepsilon ,j}(\textbf{u})$ in (20).

6 Application to the Navier–Stokes Equations

In this section, we will show the homogenization result for Navier–Stokes equations in a randomly perforated domain in the subcritical case $\alpha >3$. The proof of such result in the case of periodically arranged holes was developed in a series of works [9, 15, 30], and can be split in two parts. First, using Bogovskiĭ’s operator we construct a good test function for the momentum equation, which leads to uniform in $\varepsilon$ estimates on the density as well as the velocity, subsequently providing the compactness. To identify the limiting “effective” equation, we need to construct a suitable cut-off function in order to compare the limiting equation with the equation in $D_\varepsilon$. Since the rest of the proof does not refer in any way to location or structure of the holes, in particular it applies verbatim in our context, for that remaining part of the proof we only sketch the main steps. To shorten the exposition we sketch the argument only in the stationary case, following [15]. An analogous homogenization result holds also in the time-dependent setting, and for stationary Navier–Stokes–Fourier equations—see the statements and the proofs of [30, Theorem 1.6] and [35, Theorem 2], respectively.

6.1 Test Functions

Before we formulate and show the homogenization result, we prove a modification of [30, Lemma 2.1] in the random setting as the last ingredient in the proof of Theorem 2, which makes a reference to the randomness in the structure of the holes.

Lemma 4

Let $\alpha > 2$, $D\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ be a bounded $C^2$ star-shaped domain with $0 \in D$, and $(\Phi ,{\mathcal {R}})=(\{z_i\},\{r_i\})$ be a marked Poisson point process with intensity $\lambda > 0$ and $r_i \geqq 0$ with ${\mathbb {E}}(r_i^M) < \infty$ for $M=\max \{3,m\}$, where $m>3/(\alpha -2)$. Then for any $1< r < 3$ such that $(3-r)\alpha - 3 > 0$ and for almost every $\omega$ there exist a positive $\varepsilon _0(\omega )$ and a family of functions $\{ g_\varepsilon \}_{\varepsilon > 0} \subset W^{1,r}(D)$ such that, for $0 < \varepsilon \le \varepsilon _0$,

$$\begin{aligned} g_\varepsilon = 0 \quad \text { in } \bigcup _{z_j \in \Phi ^\varepsilon (D)} B_{\varepsilon ^\alpha r_j}(\varepsilon z_j), \qquad g_\varepsilon \rightarrow 1 \quad \text { in } W^{1,r}(D) \text { as } \varepsilon \rightarrow 0, \end{aligned}$$

(24)

and there is a constant $C>0$ such that

$$\begin{aligned} \Vert g_\varepsilon - 1 \Vert _{W^{1,r}(D)} \leqq C\varepsilon ^{\sigma } \qquad \text{ with } \sigma := ((3-r)\alpha -3)/r.\end{aligned}$$

(25)

Proof

By $M>3/(\alpha -2)$ and Theorem 3, all the balls $\{B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j)\}_{z_j\in \Phi ^\varepsilon (D)}$ are disjoint. Thus, there exist functions $g_\varepsilon \in C^\infty (D)$ such that

$$\begin{aligned}&0\le g_\varepsilon \le 1,\quad g_\varepsilon = 0 \text { in } \bigcup _{z_j \in \Phi ^\varepsilon (D)} B_{\varepsilon ^\alpha r_j}(\varepsilon z_j),\quad g_\varepsilon = 1 \text { in } D\setminus \bigcup _{z_j \in \Phi ^\varepsilon (D)} B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j),\\&\Vert \nabla g_\varepsilon \Vert _{L^\infty (B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j))}\le C (\varepsilon ^\alpha r_j)^{-1} \text { for all } z_j\in \Phi ^\varepsilon (D), \end{aligned}$$

where the constant $C>0$ is independent of $\varepsilon$ and $r_j$. Moreover, since $M\ge 3$, (9) yields $\lim \limits _{\varepsilon \rightarrow 0} \varepsilon ^3 \sum _{z_j \in \Phi ^\varepsilon (D)} r_j^3 = C$, thus implying that

$$\begin{aligned} \biggl | \bigcup _{z_j \in \Phi ^\varepsilon (D)} B_{2\varepsilon ^\alpha r_j}(\varepsilon z_j) \biggr | \leqq |B_2| \varepsilon ^{3\alpha } \sum _{z_j \in \Phi ^\varepsilon (D)} r_j^3 \leqq C \varepsilon ^{3(\alpha -1)} \end{aligned}$$

(26)

for $\varepsilon > 0$ small enough. This, together with direct calculation, yields that, for any $1<r<3$,

$$\begin{aligned} \Vert 1-g_\varepsilon \Vert _{L^r(D)}\le C\varepsilon ^\frac{3(\alpha -1)}{r},\quad \Vert \nabla g_\varepsilon \Vert _{L^r(D)}\le C\varepsilon ^{\frac{(3-r)\alpha -3}{r}}, \end{aligned}$$

which finally leads to

$$\begin{aligned} \Vert 1-g_\varepsilon \Vert _{W^{1,r}(D)}=\Vert 1-g_\varepsilon \Vert _{L^r(D)}+\Vert \nabla g_\varepsilon \Vert _{L^r(D)}\le C\varepsilon ^\sigma . \end{aligned}$$

$\square$

6.2 The Navier–Stokes Equations, Weak Solutions, and the Convergence Result

For $\varepsilon > 0$, in the domain $D_\varepsilon$ as in (1) we consider the stationary Navier–Stokes equations for compressible viscous fluids

$$\begin{aligned} {{\,\textrm{div}\,}}(\varrho _\varepsilon \textbf{u}_\varepsilon )&=0{} & {} \text { in } D_\varepsilon , \end{aligned}$$

(27)

$$\begin{aligned} {{\,\textrm{div}\,}}(\varrho _\varepsilon \textbf{u}_\varepsilon \otimes \textbf{u}_\varepsilon ) + \nabla p(\varrho _\varepsilon )&= {{\,\textrm{div}\,}}{\mathbb {S}}(\nabla \textbf{u}_\varepsilon ) + \varrho _\varepsilon \textbf{f}+\textbf{g}{} & {} \text { in } D_\varepsilon , \end{aligned}$$

(28)

$$\begin{aligned} \textbf{u}_\varepsilon&=0{} & {} \text { on } \partial D_\varepsilon , \end{aligned}$$

(29)

where ${\mathbb {S}}$ denotes the Newtonian viscous stress tensor of the form

$$\begin{aligned} {\mathbb {S}}(\nabla {\textbf {u}})=\mu \bigg (\nabla {\textbf {u}}+\nabla ^T {\textbf {u}}-\frac{2}{3} {{\,\text {div}\,}}({\textbf {u}}){\mathbb {I}}\bigg )+\eta {{\,\text {div}\,}}({\textbf {u}}){\mathbb {I}},\quad \mu > 0, \quad \eta \geqq 0, \end{aligned}$$

$p(\varrho ) = a\varrho ^{\gamma }$ denotes the pressure with $a > 0$ and the adiabatic exponent $\gamma \geqq 1$, and $\textbf{f}$ and $\textbf{g}$ are external forces satisfying $\Vert {\textbf {f}}\Vert _{L^\infty ({{\,\mathrm {{\mathbb {R}}}\,}}^3;{{\,\mathrm {{\mathbb {R}}}\,}}^3)} + \Vert {\textbf {g}}\Vert _{L^\infty ({{\,\mathrm {{\mathbb {R}}}\,}}^3;{{\,\mathrm {{\mathbb {R}}}\,}}^3)} \leqq C$. We also fix the total mass

$$\begin{aligned} \int _{D_\varepsilon } \varrho _\varepsilon ={\mathfrak {m}}>0 \end{aligned}$$

independently of $\varepsilon >0$.

Definition 2

We call a couple $[\varrho ,\textbf{u}]$ a renormalized finite energy weak solution to equations (28)–(30) if

$$\begin{aligned} \varrho \ge 0 \text { almost everywhere in }D_\varepsilon ,\quad \int _{D_\varepsilon }\varrho ={\mathfrak {m}},\quad \varrho \in L^{2\gamma }(D_\varepsilon ),\quad \textbf{u}\in W^{1,2}_0(D_\varepsilon );\\ \int _{D_\varepsilon } p(\varrho ){{\,\textrm{div}\,}}\varphi +\varrho \textbf{u}\otimes \textbf{u}:\nabla \varphi -{\mathbb {S}}(\nabla \textbf{u}):\nabla \varphi +(\varrho \textbf{f}+\textbf{g})\cdot \varphi =0 \end{aligned}$$

for all all test functions $\varphi \in C_c^\infty (D_\varepsilon ;{{\,\mathrm{{\mathbb {R}}}\,}}^3),$ the energy inequality

$$\begin{aligned} \int _{D_\varepsilon } {\mathbb {S}}(\nabla \textbf{u}):\nabla \textbf{u}\le \int _{D_\varepsilon } (\varrho \textbf{f}+\textbf{g})\cdot \textbf{u}\end{aligned}$$

(30)

holds, and the zero extension $[{\tilde{\varrho }},\tilde{\textbf{u}}]$ satisfies, in ${\mathcal {D}}^\prime ({{\,\mathrm{{\mathbb {R}}}\,}}^3)$, that

$$\begin{aligned} {{\,\textrm{div}\,}}({\tilde{\varrho }}\tilde{\textbf{u}})=0,\quad {{\,\textrm{div}\,}}(b({\tilde{\varrho }})\tilde{\textbf{u}})+({\tilde{\varrho }}b^\prime ({\tilde{\varrho }})-b({\tilde{\varrho }})){{\,\textrm{div}\,}}\tilde{\textbf{u}}=0 \end{aligned}$$

for any $b\in C([0,\infty ))\cap C^1((0,\infty ))$ such that there are constants

$$\begin{aligned} c>0,\quad \lambda _0<1,\quad -1<\lambda _1\le \gamma -1 \end{aligned}$$

with

$$\begin{aligned} b^\prime (s)\le cs^{-\lambda _0} \text { for } s\in (0,1],\quad b^\prime (s)\le cs^{\lambda _1} \text { for } s\in [1,\infty ). \end{aligned}$$

As announced above we show how to apply our results only in the stationary case. For that let us first formulate the actual result as follows:

Theorem 7

Assume $\alpha >3$. Let $D\subset {{\,\mathrm{{\mathbb {R}}}\,}}^3$ be a bounded star-shaped domain with respect to the origin with $C^2$-boundary and let $(\Phi ,{\mathcal {R}})=(\{z_j\},\{r_j\})$ be a marked Poisson point process with intensity $\lambda > 0$, and $r_j \geqq 0$ with ${\mathbb {E}}(r_j^M) < \infty$, $M=\max \{3,m\}$, $m>3/(\alpha -3)$. Furthermore, let

$$\begin{aligned} {\mathfrak {m}}>0,\; \gamma >3. \end{aligned}$$

Then, for almost every $\omega \in \Omega$, there exists $\varepsilon _0=\varepsilon _0(\omega ) > 0$, such that the following holds: For $0< \varepsilon < 1$ let $D_\varepsilon$ be as in (1) and let $[\varrho _\varepsilon ,\textbf{u}_\varepsilon ]$ be a family of renormalized finite energy weak solutions to (28)–(30). Then there is a constant $C>0$, which is independent of $\varepsilon$, such that

$$\begin{aligned} \sup _{\varepsilon \in (0,\varepsilon _0)} \Vert {\tilde{\varrho }}_\varepsilon \Vert _{L^{2\gamma }(D)}+\Vert \tilde{\textbf{u}}_\varepsilon \Vert _{W_0^{1,2}(D)} \le C, \end{aligned}$$

and, up to a subsequence,

$$\begin{aligned} {\tilde{\varrho }}_\varepsilon \rightharpoonup \varrho \text { in } L^{2\gamma }(D),\quad \tilde{\textbf{u}}_\varepsilon \rightharpoonup \textbf{u}\text { in } W_0^{1,2}(D), \end{aligned}$$

where the limit $[\varrho ,\textbf{u}]$ is a renormalized finite energy weak solution to the problem (28)–(30) in the limit domain D.

Remark 8

The condition $M>3/(\alpha -3)$ on the size of radii of the perforations is not just needed for technical purposes, but it is in a sense an optimal assumption. Indeed, one can show that in the case $m=3/(\alpha -3)$, for almost every realization of points and radii there is a sequence of $\varepsilon \rightarrow 0$ such that for each such $\varepsilon$ the rescaled radius $\varepsilon ^\alpha r_j$ of the largest ball in D is of size $\varepsilon ^3$, that is, $r_j \sim \varepsilon ^{3-\alpha }$. While one large ball of size $\varepsilon ^3$ does not necessarily mean that the system should behave as in the critical case (which is expected to lead to a law of Brinkman type), nevertheless the presence of such large ball might change some of the properties of the system. Moreover, in the case $m < 3/(\alpha -3)$, the size of the largest ball would scale like $\varepsilon ^\nu$ with $\nu <3$, and there might be many balls of size at least $\varepsilon ^3$.

Sketch of the proof of Theorem 7

We want to give uniform bounds for the velocity $\textbf{u}_\varepsilon$ and the density $\varrho _\varepsilon$ arising in the Navier–Stokes equations (28)–(30). First, by the energy inequality (31), Korn’s and Hölder’s inequality, we have that

$$\begin{aligned} \Vert \nabla \textbf{u}_\varepsilon \Vert _{L^2(D_\varepsilon )}^2\lesssim \Vert \textbf{f}\Vert _{L^\infty (D_\varepsilon )}\Vert \varrho _\varepsilon \Vert _{L^\frac{6}{5}(D_\varepsilon )}\Vert \textbf{u}_\varepsilon \Vert _{L^6(D_\varepsilon )}+\Vert \textbf{g}\Vert _{L^\infty (D_\varepsilon )} \Vert \textbf{u}_\varepsilon \Vert _{L^6(D_\varepsilon )}. \end{aligned}$$

Since $\textbf{u}_\varepsilon \in W_0^{1,2}(D_\varepsilon )$, we can use Poincaré’s inequality and Sobolev embedding to obtain $\Vert \textbf{u}_\varepsilon \Vert _{L^6(D_\varepsilon )}\lesssim \Vert \nabla \textbf{u}_\varepsilon \Vert _{L^2(D_\varepsilon )}$, which combined with the previous display yields

$$\begin{aligned} \begin{aligned} \Vert \nabla \textbf{u}_\varepsilon \Vert _{L^2(D_\varepsilon )}+\Vert \textbf{u}_\varepsilon \Vert _{L^6(D_\varepsilon )}&\lesssim \Vert \textbf{f}\Vert _{L^\infty (D_\varepsilon )}\Vert \varrho _\varepsilon \Vert _{L^\frac{6}{5}(D_\varepsilon )}+\Vert \textbf{g}\Vert _{L^\infty (D_\varepsilon )} \\&\lesssim \Vert \varrho _\varepsilon \Vert _{L^\frac{6}{5} (D_\varepsilon )}+1. \end{aligned} \end{aligned}$$

(31)

We define a test function

$$\begin{aligned} \varphi :={\mathcal {B}}_\varepsilon \bigl (\varrho _\varepsilon ^\gamma -\langle \varrho _\varepsilon ^\gamma \rangle _{D_\varepsilon }\bigr ), \end{aligned}$$

where $\langle \varrho _\varepsilon ^\gamma \rangle _{D_\varepsilon }:=|D_\varepsilon |^{-1}\int _{D_\varepsilon } \varrho _\varepsilon ^\gamma$ is the mean value of $\varrho _\varepsilon ^\gamma$ over the domain $D_\varepsilon$ and ${\mathcal {B}}_\varepsilon$ is the Bogovskiĭ operator constructed in Theorem 1. We remark that $\varphi$ is well-defined due to the fact $\varrho _\varepsilon ^\gamma \in L^2(D_\varepsilon )$. By the properties of ${\mathcal {B}}_\varepsilon$, we obtain ${{\,\textrm{div}\,}}\varphi =\varrho _\varepsilon ^\gamma -\langle \varrho _\varepsilon ^\gamma \rangle _{D_\varepsilon }$ in $D_\varepsilon$ and

$$\begin{aligned} \Vert \varphi \Vert _{W_0^{1,2}(D_\varepsilon )}\lesssim \Vert \varrho _\varepsilon ^\gamma \Vert _{L^2(D_\varepsilon )}+\Vert \varrho _\varepsilon ^\gamma \Vert _{L^1(D_\varepsilon )}\lesssim \Vert \varrho _\varepsilon \Vert _{L^{2\gamma }(D_\varepsilon )}^\gamma . \end{aligned}$$

Testing (29) with $\varphi$ yields

$$\begin{aligned} \int _{D_\varepsilon } p(\varrho _\varepsilon )\varrho _\varepsilon ^\gamma =\sum _{j=1}^4 I_j, \end{aligned}$$

where the integrals $I_j$ are defined as

$$\begin{aligned}&I_1:=\int _{D_\varepsilon } p(\varrho _\varepsilon )\langle \varrho _\varepsilon ^\gamma \rangle _{D_\varepsilon },\quad{} & {} I_2:=\int _{D_\varepsilon }\mu \nabla \textbf{u}:\nabla \varphi +\bigg (\frac{\mu }{3}+\eta \bigg ){{\,\textrm{div}\,}}\textbf{u}_\varepsilon {{\,\textrm{div}\,}}\varphi ,\\&I_3:=-\int _{D_\varepsilon } \varrho _\varepsilon \textbf{u}_\varepsilon \otimes \textbf{u}_\varepsilon :\nabla \varphi ,\quad{} & {} I_4:=-\int _{D_\varepsilon } (\varrho _\varepsilon \textbf{f}+\textbf{g})\cdot \varphi . \end{aligned}$$

By interpolation of Lebesgue spaces, we estimate $I_1$ by

$$\begin{aligned} |I_1|\le \frac{1}{|D_\varepsilon |} \Vert \varrho _\varepsilon \Vert _{L^\gamma (D_\varepsilon )}^{2\gamma }\le \frac{1}{|D_\varepsilon |}\bigg (\Vert \varrho _\varepsilon \Vert _{L^1(D_\varepsilon )}^{\theta _1}\Vert \varrho _\varepsilon \Vert _{L^{2\gamma }(D_\varepsilon )}^{1-\theta _1}\bigg )^{2\gamma }\lesssim \Vert \varrho _\varepsilon \Vert _{L^{2\gamma }(D_\varepsilon )}^{2\gamma (1-\theta _1)}, \end{aligned}$$

where $\theta _1\in (0,1)$ is determined by

$$\begin{aligned} \frac{1}{\gamma }=\frac{\theta _1}{1}+\frac{1-\theta _1}{2\gamma }. \end{aligned}$$

The estimates for the remaining integrals are the same as in [15], so we refer to [15, page 386] for details. Finally, we obtain $\Vert \varrho _\varepsilon \Vert _{L^{2\gamma }(D_\varepsilon )}^{2\gamma }\lesssim 1+\Vert \varrho _\varepsilon \Vert _{L^{2\gamma }(D_\varepsilon )}^{2\gamma (1-\beta )}$ for some $\beta >0$, which yields $\Vert \varrho _\varepsilon \Vert _{L^{2\gamma }(D_\varepsilon )}\le C$. In view of (32), we also have $\Vert \textbf{u}_\varepsilon \Vert _{W_0^{1,2}(D_\varepsilon )}\le C$, where the constant $C>0$ does not depend on $\varepsilon$. This completes the proof for the uniform bounds.

In the following we want to identify the limiting equations. First, using the fact that $[\varrho _\varepsilon , \textbf{u}_\varepsilon ]$ is a renormalized weak solution in $D_\varepsilon$ we get that the zero extensions of $\varrho _\varepsilon$ and $\textbf{u}_\varepsilon$ solve

$$\begin{aligned} {{\,\textrm{div}\,}}({\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon )=0,\quad {{\,\textrm{div}\,}}(b({\tilde{\varrho }}_\varepsilon )\tilde{\textbf{u}}_\varepsilon )+({\tilde{\varrho }}_\varepsilon b^\prime ({\tilde{\varrho }}_\varepsilon )-b({\tilde{\varrho }}_\varepsilon )){{\,\textrm{div}\,}}\tilde{\textbf{u}}_\varepsilon =0 \text { in } {\mathcal {D}}^\prime , \end{aligned}$$

where $b\in C([0,\infty ))\cap C^1((0,\infty ))$ is as in Definition 2.

Considering the momentum equation in the whole domain, we get an error $F_\varepsilon$ on the right-hand side of the equation. Since the balls are tiny ($\alpha > 3$), this friction term is in the limit negligible. More precisely, the zero prolongations of the density and velocity satisfy

$$\begin{aligned} \nabla p({\tilde{\varrho }}_\varepsilon )+{{\,\textrm{div}\,}}({\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \otimes \tilde{\textbf{u}}_\varepsilon )-{{\,\textrm{div}\,}}{\mathbb {S}}(\nabla \tilde{\textbf{u}}_\varepsilon )={\tilde{\varrho }}_\varepsilon \textbf{f}+\textbf{g}+F_\varepsilon , \end{aligned}$$

(32)

where $F_\varepsilon$ is a distribution satisfying for all $\varphi \in C_c^\infty (D)$

$$\begin{aligned} |\langle F_\varepsilon ,\varphi \rangle _{{\mathcal {D}}^\prime , {\mathcal {D}}}|\lesssim \varepsilon ^\sigma \Vert \varphi \Vert _{L^r(D)}+\varepsilon ^\frac{3(\alpha -1)\sigma _0}{2(2+\sigma _0)}\Vert \nabla \varphi \Vert _{L^{2+\sigma _0}(D)} \end{aligned}$$

with

$$\begin{aligned} \sigma :=\frac{\alpha -3}{4},\quad r:=\frac{12(\alpha -1)}{\alpha -3},\quad \sigma _0\in (0,\infty ). \end{aligned}$$

To show this, we will use the cut-off function $g_\varepsilon$ from Lemma 4. For any test function $\varphi \in C_c^\infty (D)$, we test the momentum equation in $D_\varepsilon$ with $g_\varepsilon \varphi$ to get:

$$\begin{aligned}&\int _D {\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \otimes \tilde{\textbf{u}}_\varepsilon : \nabla \varphi +p({\tilde{\varrho }}_\varepsilon ) {{\,\textrm{div}\,}}\varphi -{\mathbb {S}}(\nabla \tilde{\textbf{u}}_\varepsilon ) :\nabla \varphi +({\tilde{\varrho }}_\varepsilon \textbf{f}+\textbf{g})\cdot \varphi \, dx\\&\quad =I_\varepsilon +\int _D {\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \otimes \tilde{\textbf{u}}_\varepsilon : \nabla (g_\varepsilon \varphi )+p({\tilde{\varrho }}_\varepsilon ) {{\,\textrm{div}\,}}(g_\varepsilon \varphi )-{\mathbb {S}}(\nabla \tilde{\textbf{u}}_\varepsilon ) :\nabla (g_\varepsilon \varphi ) \\&\qquad \quad +({\tilde{\varrho }}_\varepsilon \textbf{f}+\textbf{g})\cdot (g_\varepsilon \varphi )\, dx\\&\quad =I_\varepsilon , \end{aligned}$$

where we used that $g_\varepsilon \varphi \in C_c^\infty (D)$ is an appropriate test function, and the term $I_\varepsilon$ is given by

$$\begin{aligned} I_\varepsilon :=\sum _{j=1}^4 I_{\varepsilon , j}:=&\int _D {\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \otimes \tilde{\textbf{u}}_\varepsilon : (1-g_\varepsilon )\nabla \varphi -{\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \otimes \tilde{\textbf{u}}_\varepsilon : (\nabla g_\varepsilon \otimes \varphi )\, dx\\&+\int _D p({\tilde{\varrho }}_\varepsilon ) (1-g_\varepsilon ){{\,\textrm{div}\,}}\varphi -p({\tilde{\varrho }}_\varepsilon )\nabla g_\varepsilon \cdot \varphi \, dx\\&+\int _D -{\mathbb {S}}(\nabla \tilde{\textbf{u}}_\varepsilon ):(1-g_\varepsilon )\nabla \varphi +{\mathbb {S}}(\nabla \tilde{\textbf{u}} _\varepsilon ):(\nabla g_\varepsilon \otimes \varphi )\, dx\\&+\int _D ({\tilde{\varrho }}_\varepsilon \textbf{f}+\textbf{g})\cdot (1-g_\varepsilon )\varphi \, dx. \end{aligned}$$

Using the bounds on the cut-off function $g_\varepsilon$, we combine the previous estimates on the density and velocity to prove (33) (for details see [15, Proof of Proposition 2.2]).

By the uniform estimates on $\varrho _\varepsilon$ and $\textbf{u}_\varepsilon$, we get a subsequence (not relabeled) such that

$$\begin{aligned} {\tilde{\varrho }}_\varepsilon \rightharpoonup \varrho \text { in } L^{2\gamma }(D),\quad \tilde{\textbf{u}}_\varepsilon \rightharpoonup \textbf{u}\text { in } W_0^{1,2}(D). \end{aligned}$$

By compact Sobolev embedding, this yields

$$\begin{aligned}&\tilde{\textbf{u}}_\varepsilon \rightarrow \textbf{u}\text { strongly in }L^q(D) \text { for all } 1\le q<6,\\&{\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \rightharpoonup \varrho \textbf{u}\text { weakly in }L^q(D) \text { for any } 1< q< \frac{6\gamma }{\gamma +3},\\&{\tilde{\varrho }}_\varepsilon \tilde{\textbf{u}}_\varepsilon \otimes \tilde{\textbf{u}}_\varepsilon \rightharpoonup \varrho \textbf{u}\otimes \textbf{u}\text { weakly in }L^q(D) \text { for all } 1<q<\frac{6\gamma }{2\gamma +3}. \end{aligned}$$

Letting $\varepsilon \rightarrow 0$ in equations (28) and (29), we get the following equations in ${\mathcal {D}}^\prime (D)$:

$$\begin{aligned} {{\,\textrm{div}\,}}(\varrho \textbf{u})&=0,\\ {{\,\textrm{div}\,}}(\varrho \textbf{u}\otimes \textbf{u})+\overline{p(\varrho )}&={{\,\textrm{div}\,}}{\mathbb {S}}(\nabla \textbf{u})+\varrho \textbf{f}+\textbf{g}. \end{aligned}$$

Here $\overline{p(\varrho )}$ is the weak limit of $p({\tilde{\varrho }}_\varepsilon )$ in $L^2(D)$. Moreover, the couple $[\varrho ,u]$ satisfies the renormalized equations. To finish the proof of Theorem 2, we have to prove that $\overline{p(\varrho )}=p(\varrho )$, arguing as in [15, Section 2.4.2].