Identifying equal powers of η implies $w_{i,0}(0) = w_{i,1}(0) = 0$ and (3.20). Considering now $x \in \partial P$ and using the far field expansions of $\mathbf{\Phi}_i$ and $\mathbf{v}_{i,p}$, we find

$$
\begin{align*}
\mathcal{X}_i(x) &= \eta^{2-d} \left( \mathcal{F}^{\infty} \mathbf{\Phi}_i - \delta_{d=2} \frac{\log \eta}{4\pi} + \sum_{p=0}^{+\infty} \eta^{d-2+p} \nabla^p G(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i \right) + \sum_{p=2}^{+\infty} \sum_{k=1}^{+\infty} \eta^{d-2+p+k} \nabla^k G(x) \mathcal{F}^k \mathbf{v}_{i,p} + \sum_{p=0}^{+\infty} \eta^p \mathbf{w}_{i,p}(x) \\
&= \eta^{2-d} \left( \mathcal{F}^{\infty} \mathbf{\Phi}_i - \delta_{d=2} \frac{\log \eta}{4\pi} \right) + \sum_{p=0}^{+\infty} \eta^p \nabla^p G(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i + \sum_{p=d+1}^{+\infty} \sum_{k=d-1}^{p-2} \eta^p \nabla^{k-d+2} G(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} \\
&\quad + \sum_{p=0}^{+\infty} \eta^p \mathbf{w}_{i,p}.
\end{align*}
$$

Identifying again equal powers of η yields the boundary conditions of (3.19) (the constant being absorbed by the periodicity condition). Finally, let us verify that the systems (3.19) and (3.20) are well-posed by checking that their respective compatibility conditions are satisfied. The compatibility condition for (3.19) reads

$$
\int_{\partial P} \mathbf{n} \cdot \left( \nabla^p G \cdot \mathcal{F}^p \mathbf{\Phi}_i + \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} \right) \, d\sigma = 0,
$$

which is satisfied due to div(G) = 0 in P. The compatibility condition for (3.20) reads

$$
- \int_{\partial T} \left( \sum_{k=0}^{p} \nabla^k \mathbf{w}_{i,p}(0) \cdot \mathbf{y}^k \right) \cdot \mathbf{n} \, d\sigma(y) = 0.
$$

This condition is also satisfied because div($\mathbf{w}_{i,p}$) = 0 implies

$$
\operatorname{div}(\nabla^k \mathbf{w}_{i,p}(0) \cdot \mathbf{y}^k) = \sum_{1 \le i_1, \dots, i_k \le d} \sum_{j=1}^{d} \partial_{y_{i_1} \dots y_{i_k}} w_{i,p,j}(0) \partial_{y_j} (y_{i_1} \dots y_{i_k}) = k(\nabla^{k-1} \operatorname{div}(\mathbf{w}_{i,p}))(0) \cdot \mathbf{y}^{k-1} = 0.
$$

□

The next proposition proves that the asymptotic expansions (3.16) hold in the energy norm.

**Proposition 3.5.** For any $N \in \mathbb{N}$, denote by $(\mathcal{X}_i^N, \alpha_i^N)$ the approximations of $(\mathcal{X}_i, \alpha_i)$ obtained by truncating the expansions (3.16) at rank $N$:

$$
\left\{
\begin{array}{l}
\mathcal{X}_i^N(x) := \eta^{2-d} \mathbf{\Phi}_i(x/\eta) + \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \mathbf{w}_{i,p}(x), \\
\\
\alpha_i^N(x) := \eta^{1-d} \sigma_i(x/\eta) + \sum_{p=2}^N \eta^{p-1} \xi_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \zeta_{i,p}(x).
\end{array}
\right.
\qquad (3.21)
$$

There exists a uniform constant $C > 0$ such that for any $N \ge 0$, $(\mathcal{X}_i^N, \alpha_i^N)$ is an approximation of the solution $(\mathcal{X}_i, \alpha_i)$ to (1.3) at order $O(\eta^{N+1})$ according to the following error bounds:

$$
\begin{align}
& \|\mathcal{X}_i - \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d})} + \|\nabla\mathcal{X}_i - \nabla\mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \nonumber \\
& \qquad + \|\alpha_i - \alpha_i^N\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \le C\eta^{N+1} \quad \text{if } d \ge 3 \text{ or } (d=2 \text{ and } N=0), \tag{3.22}
\end{align}
$$

$$
\begin{align}
& |\log \eta|^{-\frac{1}{2}} \|\mathcal{X}_i - \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} + \|\nabla \mathcal{X}_i - \nabla \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \nonumber \\
& \qquad + \|\alpha_i - \alpha_i^N\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \le C\eta^{N+1} \quad \text{if } d=2 \text{ and } N \ge 1. \tag{3.23}
\end{align}
$$

*Proof.* Let us assume without loss of generality that $N \ge 1$. By utilizing the summability properties (3.17) and (3.18), we find the existence of a constant $C$ uniform in $N$ and $\eta$ such that

$$
\begin{align}
\sum_{p=N+1}^{+\infty} \eta^{p-N-1} \|\nabla(\mathbf{v}_{i,p}(\cdot/\eta))\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} &= \eta^{\frac{d}{2}-1} \sum_{p=N+1}^{+\infty} \eta^{p-N-1} \|\nabla \mathbf{v}_{i,p}\|_{L^2(\eta^{-1}P\setminus T, \mathbb{R}^{d \times d})} \tag{3.24} \\
&\le \eta^{\frac{d}{2}-1} \sum_{p=N+1}^{+\infty} \eta^{p-N-1} \|\nabla \mathbf{v}_{i,p}\|_{L^2(\mathbb{R}^d \setminus T, \mathbb{R}^{d \times d})} \le C\eta^{\frac{d}{2}-1}, \nonumber
\end{align}
$$
$$ \begin{aligned} \|\mathcal{X}_i - \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} &\le \eta^{N+1} \sum_{p=N+1}^{+\infty} \eta^{p-N-1} (\|\mathbf{v}_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d})} + \|\mathbf{w}_{i,p}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d})}) \\ &\le \begin{cases} C|\log \eta|^{\frac{1}{2}}\eta^{N+1} & \text{if } d=2, \\ C\eta^{N+1} & \text{if } d=3, \end{cases} \end{aligned} $$

and

$$ \|\alpha_i - \alpha_i^N\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \le \eta^{N+1} \sum_{p=N+1}^{+\infty} \left( \eta^{p-N-2}\|\xi_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T))/\mathbb{R}} + \|\zeta_{i,p}\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \right) \le C\eta^{N+1}. $$

□

**Remark 3.2.** The proof of Proposition 3.5 uses crucially the fact that the series arising in (3.16) are summable, which has been obtained from the layer potential representation. A different proof that uses only variational estimates obtained directly from the definition of the interior and exterior problems (3.19) and (3.20) is given in Appendix A. However, the proof yields a bounding constant which is a priori not uniform in the truncation order N.

**Remark 3.3.** The $|\log \eta|^{-\frac{1}{2}}$ factor in dimension $d=2$ cannot be removed. Using $\mathbf{v}_{i,p}(y) = O(|y|^{-1})$ as $|y| \to +\infty$, we find that the corresponding term in the series (3.16) has indeed a magnitude of order $\|\mathbf{v}_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} = O(|\log \eta|^{\frac{1}{2}})$.

We conclude this section by refining the results of Proposition 3.4: some simplification arise by exploiting the cubic cell symmetries of $P$ and the radial symmetries of the fundamental solution $(G, \psi)$. We also provide a more explicit formula for the functions $(\mathbf{w}_{i,p})_{p \ge 0}$.

**Proposition 3.6.** If the dimension $d$ is an odd integer, then for any $p \in \{0, 1, \dots, d\}$:

* $\mathbf{w}_{i,p}(-x) = \mathbf{w}_{i,p}(x)$ and $\zeta_{i,p}(-x) = -\zeta_{i,p}(x)$ for all $x \in P$ if $p$ is even;

* $\mathbf{w}_{i,p}(-x) = -\mathbf{w}_{i,p}(x)$ and $\zeta_{i,p}(-x) = \zeta_{i,p}(x)$ for all $x \in P$ if $p$ is odd;

* $\mathbf{v}_{i,p} = 0$ if $p$ is odd.

If $d$ is even, these properties actually hold for all $p \in \mathbb{N}$.

**Proof.** Let us consider an integer $p \in \mathbb{N}$ satisfying $0 \le p \le d$. Since $G$ and $\psi$ are respectively even and odd, it holds $\nabla^p\Gamma(-x) = (-1)^p\Gamma(x)$ and $\nabla^p\psi(-x) = (-1)^{p+1}\nabla^p\psi(x)$. Using the boundary conditions for $\mathbf{w}_{i,p}$ for $0 \le p \le d$, we deduce the existence of constant vectors $C_{i,p}$ such that for all $0 \le i, p \le d$:

$$ \mathbf{w}_{i,p}(-x) = (-1)^p \mathbf{w}_{i,p}(x) + C_{i,p} \quad \text{and} \quad \zeta_{i,p}(-x) = (-1)^{p+1} \zeta_{i,p}(x) \quad \forall x \in P. \quad (3.26) $$

If $p$ is even, it is clear that $C_{i,p} = 0$ and thus $\mathbf{w}_{i,p}(-x) = \mathbf{w}_{i,p}(x)$ when $p$ is even. If $p$ is odd, differentiating (3.26) $k$ times implies $\nabla^k \mathbf{w}_{i,p-k}(0) = (-1)^p \nabla^k \mathbf{w}_{i,p-k}(0)$, and thus $\nabla^k \mathbf{w}_{i,p-k}(0) = 0$ for any $k \ge 1$. Consequently, the boundary condition for $\mathbf{v}_{i,p}$ in (3.20) with $p \in \{2, \dots, d\}$ reduces to

$$ \mathbf{v}_{i,p}(y) = -\mathbf{w}_{i,p}(0) \quad \text{for all } y \in \partial T, $$

which implies $\mathbf{v}_{i,p} = -\mathbf{w}_{i,p}(0)$ if $d=2$ or $\mathbf{v}_{i,p} = -\sum_{j=1}^d (e_j^T \mathbf{w}_{i,p}(0)) \Psi_j$ if $d \ge 3$. Due to the far field condition $\mathbf{v}_{i,p}(y) = O(|y|^{1-d})$, it is then necessary that $\mathbf{w}_{i,p}(0) = 0$. It follows that $\mathbf{v}_{i,p} = 0$ and $C_{i,p} = 0$, proving the result for all $p$ satisfying $0 \le p \le d$.

Let us now assume that $d$ is even and that the properties of the proposition are proved until rank $p-1$ with $p > d$. We rewrite the sum arising in the boundary condition for $\mathbf{w}_{i,p}$ as

$$ \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} = \sum_{k=2}^{p-d+1} \nabla^{p-k-d+2} G \cdot \mathcal{F}^{p-k-d+2} \mathbf{v}_{i,k} = \sum_{k=1}^{\lfloor \frac{p-d+1}{2} \rfloor} \nabla^{p-2k-d+2} G \cdot \mathcal{F}^{p-2k-d+2} \mathbf{v}_{i,2k}. $$

When $d$ is even, the exponent $p-2k-d+2$ in $\nabla^{p-2k-d+2}\Gamma$ has the same parity than the integer $p$. This implies again (3.26). The same reasoning as for the case $0 \le p \le d$ can then be applied, showing that result holds at rank $p$. Thus, the proof is complete by induction. □

An explicit formula is actually available for the functions $\mathbf{w}_{i,p}$ once one has computed the functions $\mathcal{R}_i$ introduced in (2.48).

**Proposition 3.7.** For any $p \ge 0$, it holds

$$ \mathbf{w}_{i,p}(x) = \mathbf{w}_{i,p}(0) + (\nabla^p R(x) - \nabla^p R(0)) \cdot \mathcal{F}^p \Phi_i + \sum_{k=d-1}^{p-2} (\nabla^{k-d+2} R(x) - \nabla^{k-d+2} R(0)) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k}, \quad (3.27) $$
$$
\begin{align*}
& \sum_{p=N+1}^{+\infty} \eta^{p-N-1} (\|\mathbf{w}_{i,p}\|_{L^2(P\setminus(\eta T),\mathbb{R}^d)} + \|\nabla \mathbf{w}_{i,p}\|_{L^2(P\setminus(\eta T),\mathbb{R}^{d\times d})} + \|\zeta_{i,p}\|_{L^2(P\setminus(\eta T))/\mathbb{R}}) \\
& \le \sum_{p=N+1}^{+\infty} \eta^{p-N-1} (\|\mathbf{w}_{i,p}\|_{H^1(P,\mathbb{R}^d)} + \|\zeta_{i,p}\|_{L^2(P)/\mathbb{R}}) \le C. \tag{3.25} \\
& \sum_{p=N+1}^{+\infty} \eta^{p-N-2} \|\xi_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T))/\mathbb{R}} = \sum_{p=N+1}^{+\infty} \eta^{p-N-1} \eta^{\frac{d}{2}-1} \|\xi_{i,p}\|_{L^2(\eta^{-1}P\setminus T)} \\
& \le \eta^{\frac{d}{2}-1} \sum_{p=N+1}^{+\infty} \eta^{p-N-1} \|\xi_{i,p}\|_{L^2(\mathbb{R}^d\setminus T)} \le C\eta^{\frac{d}{2}-1}.
\end{align*}
$$

Combining (3.24) and the Poincaré inequality (A.3), we also obtain

$$
\sum_{p=N+1}^{+\infty} \eta^{p-N-1} \|\mathbf{v}_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d\times d})} \le \begin{cases} C |\log \eta|^{\frac{1}{2}} & \text{if } d=2, \\ C & \text{if } d \ge 3. \end{cases}
$$

By using these bounds, we can directly estimate the norm of $\mathcal{X}_i - \mathcal{X}_i^N$ and $\alpha_i - \alpha_i^N$:

$$
\begin{align*}
\|\nabla \mathcal{X}_i - \nabla \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d\times d})} &\le \eta^{N+1} \sum_{p=N+1}^{+\infty} \eta^{p-N-1} (\|\nabla(\mathbf{v}_{i,p}(\cdot/\eta))\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d\times d})} + \|\nabla \mathbf{w}_{i,p}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d\times d})}) \\
&\le C\eta^{N+1}, \\
\|\mathcal{X}_i - \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d\times d})} &\le \eta^{N+1} \sum_{p=N+1}^{+\infty} \eta^{p-N-1} (\|\mathbf{v}_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d})} + \|\mathbf{w}_{i,p}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d})}) \\
&\le \begin{cases} C |\log \eta|^{\frac{1}{2}} \eta^{N+1} & \text{if } d=2, \\ C\eta^{N+1} & \text{if } d=3, \end{cases}
\end{align*}
$$

and

$$
\|\alpha_i - \alpha_i^N\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \le \eta^{N+1} \sum_{p=N+1}^{+\infty} \left( \eta^{p-N-2} \|\xi_{i,p}(\cdot/\eta)\|_{L^2(P\setminus(\eta T))/\mathbb{R}} + \|\zeta_{i,p}\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \right) \le C\eta^{N+1}.
$$

□

Remark 3.2. The proof of Proposition 3.5 uses crucially the fact that the series arising in (3.16) are summable, which has been obtained from the layer potential representation. A different proof that uses only variational estimates obtained directly from the definition of the interior and exterior problems (3.19) and (3.20) is given in Appendix A. However, the proof yields a bounding constant which is a priori not uniform in the truncation order N.

Remark 3.3. The $|\log \eta|^{-\frac{1}{2}}$ factor in dimension $d=2$ cannot be removed. Using $\mathbf{v}_{i,p}(y) = O(|y|^{-1})$ as $|y| \to +\infty$, we find that the corresponding term in the series (3.16) has indeed a magnitude of order $\|\mathbf{v}_{i,p}(\cdot, \eta)\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} = O(|\log \eta|^{\frac{1}{2}})$.

We conclude this section by refining the results of Proposition 3.4: some simplification arise by exploiting the cubic cell symmetries of P and the radial symmetries of the fundamental solution $(G, \psi)$. We also provide a more explicit formula for the functions $(\mathbf{w}_{i,p})_{p \ge 0}$.

**Proposition 3.6.** If the dimension $d$ is an odd integer, then for any $p \in \{0, 1, \dots, d\}$:

* $\mathbf{w}_{i,p}(-x) = \mathbf{w}_{i,p}(x)$ and $\zeta_{i,p}(-x) = -\zeta_{i,p}(x)$ for all $x \in P$ if $p$ is even;

* $\mathbf{w}_{i,p}(-x) = -\mathbf{w}_{i,p}(x)$ and $\zeta_{i,p}(-x) = \zeta_{i,p}(x)$ for all $x \in P$ if $p$ is odd;

* $\mathbf{v}_{i,p} = 0$ if $p$ is odd.

If $d$ is even, these properties actually hold for all $p \in \mathbb{N}$.

*Proof.* Let us consider an integer $p \in \mathbb{N}$ satisfying $0 \le p \le d$. Since $G$ and $\psi$ are respectively even and odd, it holds $\nabla^p \Gamma(-x) = (-1)^p \Gamma(x)$ and $\nabla^p \psi(-x) = (-1)^{p+1} \nabla^p \psi(x)$. Using the boundary conditions for $\mathbf{w}_{i,p}$ for $0 \le p \le d$, we deduce the existence of constant vectors $C_{i,p}$ such that for all $0 \le i, p \le d$:

$$
\mathbf{w}_{i,p}(-x) = (-1)^p \mathbf{w}_{i,p}(x) + \mathbf{C}_{i,p} \quad \forall x \in P. \quad (3.26)
$$

If $p$ is even, it is clear that $C_{i,p} = 0$ and thus $\mathbf{w}_{i,p}(-x) = \mathbf{w}_{i,p}(x)$ when $p$ is even. If $p$ is odd, differentiating (3.26) $k$ times implies $\nabla^k \mathbf{w}_{i,p-k}(0) = (-1)^p \nabla^k \mathbf{w}_{i,p-k}(0)$, and thus $\nabla^k \mathbf{w}_{i,p-k}(0) = 0$ for any $k \ge 1$. Consequently, the boundary condition for $\mathbf{v}_{i,p}$ in (3.20) with $p \in \{2, \dots, d\}$ reduces to

$$
\mathbf{v}_{i,p}(y) = -\mathbf{w}_{i,p}(0) \text{ for all } y \in \partial T,
$$
where $\mathbf{w}_{i,p}(0)$ can be computed inductively from the explicit formula

$$ \mathbf{w}_{i,p}(0) = -\sum_{k=1}^{p} (-1)^k \mathcal{F}^k \Phi \cdot \nabla^k \mathbf{w}_{i,p-k}(0), \qquad (3.28) $$

and where $\mathcal{F}^k \Phi$ denotes the matrix valued tensor $\mathcal{F}^k \Phi = (e_i \cdot \mathcal{F}^k \Phi_j)_{1\le i,j \le d}$.

*Proof.* It is straightforward to observe that formula (3.27) solves the interior Stokes problem (3.19) because $\nabla^p G_{\text{per}} = \nabla^p R + \nabla^p G$ is a smooth periodic function in $P\setminus\{0\}$. In order to obtain the value of $\mathbf{w}_{i,p}(0)$, we notice that, if $d \ge 3$, the condition $\mathbf{v}_{i,p}(y) = O(|y|^{1-d})$ as $|y| \to +\infty$ is equivalent to $\mathcal{F}^0 \mathbf{v}_{i,p} = 0$, i.e.

$$ \mathcal{F}^0 \mathbf{v}_{i,p} \cdot \mathbf{e}_j = \mathbf{e}_j \cdot \int_{\partial T} \sigma(\mathbf{v}_{i,p}, \xi_{i,p}) \mathbf{n} \, d\sigma = \int_{\partial T} \mathbf{\Psi}_j \cdot \sigma(\mathbf{v}_{i,p}, \xi_{i,p}) \mathbf{n} \, d\sigma = 0 \text{ for all } 1 \le j \le d. $$

Integrating by parts on $\mathbb{R}^d$, this equation reduces to $\int_{\partial T} \mathbf{v}_{i,p} \cdot \sigma(\mathbf{\Psi}_j, \tau_j) \mathbf{n} \, d\sigma = 0$ for all $1 \le j \le d$. Using the boundary condition (3.20) for $\mathbf{v}_{i,p}$, this implies

$$ 0 = \int_{\partial T} \left( -\mathbf{w}_{i,p}(0) - \sum_{k=1}^{p} \frac{1}{k!} \nabla^k \mathbf{w}_{i,p-k}(0) \cdot \mathbf{y}^k \right) \cdot \sigma(\mathbf{\Psi}_j, \tau_j) \mathbf{n} \, d\sigma = \mathbf{w}_{i,p}(0) \cdot M \mathbf{e}_j - \sum_{k=1}^{p} (-1)^k \nabla^k \mathbf{w}_{i,p-k}(0) \cdot \mathcal{F}^k \mathbf{\Psi}_j. $$

Remembering that $\mathcal{F}^k \Phi_j = -\sum_{l=1}^d M_{jl}^{-1} \mathcal{F}^k \Psi_l$ (equation (2.28)), we thus obtain

$$ \mathbf{e}_j \cdot \mathbf{w}_{i,p}(0) = \sum_{l=1}^{d} M_{jl}^{-1} \sum_{k=1}^{p} (-1)^k \mathbf{w}_{i,p-k}(0) \cdot \mathcal{F}^k \mathbf{\Psi}_l = - \sum_{j=1}^{d} \sum_{k=1}^{p} (-1)^k \mathbf{w}_{i,p-k}(0) \cdot \mathcal{F}^k \Phi_j, $$

which yields (3.28).

If $d=2$, the condition $\mathbf{v}_{i,p}(y) = O(|y|^{-1})$ is satisfied if and only if $\mathcal{F}^\infty \mathbf{v}_{i,p} = 0$. Using the identity (2.34), this is equivalent to

$$ \mathcal{F}^\infty \mathbf{v}_{i,p} \cdot \mathbf{e}_j = \int_{\partial T} \mathbf{v}_{i,p} \cdot \sigma(\mathbf{\Phi}_j, \sigma_j) \mathbf{n} \, d\sigma = 0, \quad \text{for all } 1 \le j \le d. $$

Using once more the boundary condition (3.20), we deduce that (3.28) also holds in dimension $d=2$. $\square$

## 4. FULL ASYMPTOTIC EXPANSIONS OF THE PERMEABILITY MATRIX

In this section, we compute the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$ defined in (1.5). We start by noticing that it is rather straightforward to state the existence of a convergent series representation for the permeability matrix from the layer potential representation of $(\mathcal{X}_i, \alpha_i)$. Then, the result of Proposition 3.4 will enable us to identify all the terms of the series.

**Proposition 4.1.** The permeability matrix $\mathcal{X}^*$ admits the following absolutely convergent series representations:

$$ \mathcal{X}^* = \eta^{2-d} M^{-1} - cI + \sum_{p=2}^{+\infty} \eta^p M_p \text{ if } d \ge 3, $$

and

$$ \mathcal{X}^* = -\frac{\log \eta}{4\pi} I - cI + \mathcal{F}^\infty \Phi + \sum_{p=2}^{+\infty} \eta^p M_p \text{ if } d=2, $$

for some family of real $d \times d$ matrices $(M_p)_{p \ge 2}$.

*Proof.* Remembering the definition (3.1) of the potential $S_{\text{per},\eta T}^{\kappa\eta}$ and integrating $\mathcal{X}_i(x) = S_{\text{per},\eta T}^{\kappa\eta}[[\sigma(\mathcal{X}_i, \alpha_i)\mathbf{n}]](x)$ over the cell $P$ yields

$$ \mathbf{e}_j \cdot \int_P \mathcal{X}_i(x) \, dx = \mathcal{X}_{ij}^* + \int_{\eta T} \mathbf{e}_j \cdot \mathcal{X}_i(x) \, dx = \eta^{2-d} \left( -\frac{\log \eta}{4\pi} \delta_{ij} \delta_{d=2} - c\delta_{ij} + \kappa_{ij}^\eta \right), \quad (4.1) $$

where we used $\int_P G_{\text{per}} dx = 0$. For $x \in \eta T$, it can be verified that the asymptotic expansion (3.16) is still valid and yields

$$ \int_{\eta T} \mathcal{X}_i(x) \, dx = \eta^d \left( \sum_{p=2}^{+\infty} \int_T \mathbf{v}_{i,p}(y) \, dy + \sum_{p=0}^{+\infty} \int_T \mathbf{w}_{i,p}(\eta y) \, dy \right). \quad (4.2) $$

Due to the uniform convergence of the series of analytic function $\sum ||\mathbf{w}_{i,p}||_{L^\infty(P, \mathbb{R}^d)} < +\infty$ (equation (3.17)), the series $\sum_{p=0}^{+\infty} \mathbf{w}_{i,p}(\eta y)$ is analytic in $\eta$. Therefore, (4.2) implies that $\int_{\eta T} \mathcal{X}_i(x) \, dx = \eta^d \mathbf{F}_i(\eta)$ for some analytic function $\mathbf{F}_i$. Then, substituting $\kappa^\eta$ by its analytic series representation (3.13) yields the result. $\square$
which implies $v_{i,p} = -w_{i,p}(0)$ if $d=2$ or $v_{i,p} = -\sum_{j=1}^{d}(e_j^T w_{i,p}(0))\Psi_j$ if $d \ge 3$. Due to the far field condition $v_{i,p}(y) = O(|y|^{1-d})$, it is then necessary that $w_{i,p}(0) = 0$. It follows that $v_{i,p} = 0$ and $C_{i,p} = 0$, proving the result for all $p$ satisfying $0 \le p \le d$.

Let us now assume that $d$ is even and that the properties of the proposition are proved until rank $p-1$ with $p > d$. We rewrite the sum arising in the boundary condition for $w_{i,p}$ as

$$ \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G \cdot \mathcal{F}^{k-d+2} v_{i,p-k} = \sum_{k=2}^{p-d+1} \nabla^{p-k-d+2} G \cdot \mathcal{F}^{p-k-d+2} v_{i,k} = \sum_{k=1}^{\lfloor \frac{p-d+1}{2} \rfloor} \nabla^{p-2k-d+2} G \cdot \mathcal{F}^{p-2k-d+2} v_{i,2k}. $$

When $d$ is even, the exponent $p-2k-d+2$ in $\nabla^{p-2k-d+2}\Gamma$ has the same parity than the integer $p$. This implies again (3.26). The same reasoning as for the case $0 \le p \le d$ can then be applied, showing that result holds at rank $p$. Thus, the proof is complete by induction. $\square$

An explicit formula is actually available for the functions $w_{i,p}$ once one has computed the functions $R_i$ introduced in (2.48).

**Proposition 3.7.** For any $p \ge 0$, it holds

$$ w_{i,p}(x) = w_{i,p}(0) + (\nabla^p R(x) - \nabla^p R(0)) \cdot \mathcal{F}^p \Phi_i + \sum_{k=d-1}^{p-2} (\nabla^{k-d+2} R(x) - \nabla^{k-d+2} R(0)) \cdot \mathcal{F}^{k-d+2} v_{i,p-k}, \quad (3.27) $$

where $w_{i,p}(0)$ can be computed inductively from the explicit formula

$$ w_{i,p}(0) = - \sum_{k=1}^{p} (-1)^k \mathcal{F}^k \Phi \cdot \nabla^k w_{i,p-k}(0), \quad (3.28) $$

and where $\mathcal{F}^k \Phi$ denotes the matrix valued tensor $\mathcal{F}^k \Phi = (e_i \cdot \mathcal{F}^k \Phi_j)_{1 \le i,j \le d}$.

*Proof.* It is straightforward to observe that formula (3.27) solves the interior Stokes problem (3.19) because $\nabla^p G_{\text{per}} = \nabla^p R + \nabla^p G$ is a smooth periodic function in $P\setminus\{0\}$. In order to obtain the value of $w_{i,p}(0)$, we notice that, if $d \ge 3$, the condition $v_{i,p}(y) = O(|y|^{1-d})$ as $|y| \to +\infty$ is equivalent to $\mathcal{F}^0 v_{i,p} = 0$, i.e.

$$ \mathcal{F}^0 v_{i,p} \cdot e_j = e_j \cdot \int_{\partial T} \sigma(v_{i,p}, \xi_{i,p}) n \, d\sigma = \int_{\partial T} \Psi_j \cdot \sigma(v_{i,p}, \xi_{i,p}) n \, d\sigma = 0 \text{ for all } 1 \le j \le d. $$

Integrating by parts on $\mathbb{R}^d$, this equation reduces to $\int_{\partial T} v_{i,p} \cdot \sigma(\Psi_j, \tau_j) n \, d\sigma = 0$ for all $1 \le j \le d$. Using the boundary condition (3.20) for $v_{i,p}$, this implies

$$ 0 = \int_{\partial T} \left( -w_{i,p}(0) - \sum_{k=1}^{p} \frac{1}{k!} \nabla^k w_{i,p-k}(0) \cdot y^k \right) \cdot \sigma(\Psi_j, \tau_j) n \, d\sigma = w_{i,p}(0) \cdot M e_j - \sum_{k=1}^{p} (-1)^k \nabla^k w_{i,p-k}(0) \cdot \mathcal{F}^k \Psi_j. $$

Remembering that $\mathcal{F}^k \Phi_j = -\sum_{l=1}^d M_{jl}^{-1} \mathcal{F}^k \Psi_l$ (equation (2.28)), we thus obtain

$$ e_j \cdot w_{i,p}(0) = \sum_{l=1}^{d} M_{jl}^{-1} \sum_{k=1}^{p} (-1)^k w_{i,p-k}(0) \cdot \mathcal{F}^k \Psi_l = - \sum_{j=1}^{d} \sum_{k=1}^{p} (-1)^k w_{i,p-k}(0) \cdot \mathcal{F}^k \Phi_j, $$

which yields (3.28).

If $d=2$, the condition $v_{i,p}(y) = O(|y|^{-1})$ is satisfied if and only if $\mathcal{F}^\infty v_{i,p} = 0$. Using the identity (2.34), this is equivalent to

$$ \mathcal{F}^\infty v_{i,p} \cdot e_j = \int_{\partial T} v_{i,p} \cdot \sigma(\Phi_j, \sigma_j) n \, d\sigma = 0, \quad \text{for all } 1 \le j \le d. $$

Using once more the boundary condition (3.20), we deduce that (3.28) also holds in dimension $d=2$. $\square$

## 4. FULL ASYMPTOTIC EXPANSIONS OF THE PERMEABILITY MATRIX

In this section, we compute the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$ defined in (1.5). We start by noticing that it is rather straightforward to state the existence of a convergent series representation for the permeability matrix from the layer potential representation of $(\mathcal{X}_i, \alpha_i)$. Then, the result of Proposition 3.4 will enable us to identify all the terms of the series.

**Proposition 4.1.** The permeability matrix $\mathcal{X}^*$ admits the following absolutely convergent series representations:

$$ \mathcal{X}^* = \eta^{2-d} M^{-1} - cI + \sum_{p=2}^{+\infty} \eta^p M_p \text{ if } d \ge 3, $$
In the sequel, we aim at identifying explicitly the constant matrices $(M_p)_{p \ge 2}$. The Cauchy-Schwartz inequality implies

$$ \left| \int_{P \setminus (\eta T)} \boldsymbol{\chi}_i \cdot \boldsymbol{e}_j \, dx - \int_{P \setminus (\eta T)} \boldsymbol{\chi}_i^N \cdot \boldsymbol{e}_j \, dx \right| \le (1 - \eta^d |T|)^{\frac{1}{2}} \| \boldsymbol{\chi}_i - \boldsymbol{\chi}_i^N \|_{L^2(P \setminus (\eta T), \mathbb{R}^d)} \quad \text{for all } 1 \le i, j \le d. \quad (4.3) $$

This inequality together with the bounds (3.22) and (3.23) imply that the asymptotic expansion of the permeability matrix $\boldsymbol{\chi}^*$ can be obtained by simply integrating (3.16) on $P \setminus (\eta T)$:

$$ \begin{aligned} \boldsymbol{\chi}_{ij}^* = \int_{P \setminus (\eta T)} \boldsymbol{\chi}_i(x) \cdot \boldsymbol{e}_j \, dx = \eta^{2-d} \int_{P \setminus (\eta T)} \boldsymbol{\Phi}_i(x/\eta) \cdot \boldsymbol{e}_j \, dx & \\ & + \sum_{p=2}^{+\infty} \eta^p \int_{P \setminus (\eta T)} \boldsymbol{v}_{i,p}(x/\eta) \cdot \boldsymbol{e}_j \, dx + \sum_{p=0}^{+\infty} \eta^p \int_{P \setminus (\eta T)} \boldsymbol{w}_{i,p}(x) \cdot \boldsymbol{e}_j \, dx. \quad (4.4) \end{aligned} $$

The terms involving the interior solutions $(\boldsymbol{w}_{i,p})_{p \ge 0}$ can be easily expanded in powers of $\eta$ using their Taylor series expansion at zero:

$$ \int_{P \setminus (\eta T)} \boldsymbol{w}_{i,p}(x) \, dx = \int_{P} \boldsymbol{w}_{i,p}(x) \, dx - \eta^d \int_{T} \boldsymbol{w}_{i,p}(\eta y) \, dy = \int_{P} \boldsymbol{w}_{i,p}(x) \, dx - \sum_{k=0}^{+\infty} \eta^{d+k} \frac{1}{k!} \nabla^k \boldsymbol{w}_{i,p}(0) \cdot \int_{T} y^k \, dy. \quad (4.5) $$

In order to obtain the asymptotic expansion of the terms involving the exterior solutions $\boldsymbol{\Phi}_i$ and $(\boldsymbol{v}_{i,p})_{p \ge 2}$, we rely on the following lemma.

**Lemma 4.1.** Let $v \in H_{\text{loc}}^1(\mathbb{R}^d \setminus T)$ be a function with an asymptotic expansion at infinity of the form:

$$ \boldsymbol{v}(y) = \mathcal{F}^\infty \boldsymbol{v} + \sum_{k=0}^{+\infty} \nabla^k G(y) \cdot \mathcal{F}^k \boldsymbol{v} \quad \text{as } |y| \to +\infty. \quad (4.6) $$

Then, the following asymptotic expansion holds as $\eta \to 0$:

$$ \begin{aligned} \int_{P \setminus (\eta T)} \boldsymbol{v}(x/\eta) \, dx \sim \mathcal{F}^\infty \boldsymbol{v} - \frac{\delta_{d=2}}{4\pi} \left( \log \eta + \frac{1}{2} \right) \mathcal{F}^0 \boldsymbol{v} & \\ & + \eta^{d-2} \left( \frac{1}{2} \mathfrak{G}^0 \cdot \mathcal{F}^0 \boldsymbol{v} + \eta^2 \text{ f. p. } \int_{\eta^{-1}P \setminus T} \boldsymbol{v}(y) \, dy - \sum_{k=2}^{+\infty} \frac{\eta^{2k}}{2k-2} \mathfrak{G}^{2k} \cdot \mathcal{F}^{2k} \boldsymbol{v} \right), \end{aligned} \quad (4.7) $$

where $\mathfrak{G}^{2k}$ is the tensor defined in (2.59) and f. p. denotes the Hadamard finite part integral.

*Proof.* Let us denote by $I(\eta)$ the quantity

$$ I(\eta) := \eta^{-d} \int_{P \setminus (\eta T)} \boldsymbol{v}(x/\eta) \, dx = \int_{\eta^{-1}P \setminus T} \boldsymbol{v}(y) \, dy. \quad (4.8) $$

Applying a classical shape differentiation formula (e.g. [24, Theorem 5.2.2]), the derivative of $I(\eta)$ is given by

$$ \frac{dI}{d\eta}(\eta) = \frac{d}{d\eta} \int_{\eta^{-1}P \setminus T} \boldsymbol{v}(y) \, dy = \int_{\eta^{-1}P} \boldsymbol{v}(y)(\boldsymbol{\theta}(y) \cdot \boldsymbol{n}(y)) \, d\sigma(y), \quad (4.9) $$

where $\boldsymbol{\theta}(y)$ is the “velocity” of the moving boundary at $y \in \eta^{-1} \partial P$. Since a point of the moving boundary $\eta^{-1} \partial P$ can be parameterized as $y = \eta^{-1} y_0$ for a fixed $y_0 \in \partial P$, the velocity is given by $\boldsymbol{\theta}(y) = -\eta^{-2} \boldsymbol{y}_0 = -\eta^{-1} \boldsymbol{y}$. Substituting into (4.9) and performing successive changes of variables, we obtain

$$ \frac{dI}{d\eta}(\eta) = -\eta^{-1} \int_{\eta^{-1} \partial P} \boldsymbol{v}(y)(\boldsymbol{y} \cdot \boldsymbol{n}) \, d\sigma(y) = -\eta^{-1} \eta^{1-d} \eta^{-1} \int_{\partial P} \boldsymbol{v}(x/\eta)(\boldsymbol{x} \cdot \boldsymbol{n}) \, d\sigma(x). $$

Using the asymptotic expansion (4.6), $\boldsymbol{x} \cdot \boldsymbol{n} = \frac{1}{2}$ on $\partial P$ and the fact that $\nabla^k \Gamma$ is homogeneous of degree $2-d-k$ (or “associated homogeneous” if $d=2$ and $k=0$), we find

$$ \begin{aligned} \frac{d}{d\eta} I(\eta) & \sim -\eta^{-1-d} \frac{|\partial P|}{2} \mathcal{F}^\infty \boldsymbol{v} + \frac{1}{4\pi} \frac{|\partial P|}{2} \eta^{-3} \log \eta \mathcal{F}^0 \boldsymbol{v} \delta_{d=2} - \eta^{-1-d} \sum_{k=0}^{+\infty} \eta^{d-2+k} \left( \int_{\partial P} \nabla^k G(x) \boldsymbol{x} \cdot \boldsymbol{n} \, d\sigma \right) \cdot \mathcal{F}^k \boldsymbol{v} \\ & \sim -d\eta^{-1-d} \mathcal{F}^\infty \boldsymbol{v} + \frac{\delta_{d=2}}{2\pi} \eta^{-3} \log \eta \mathcal{F}^0 \boldsymbol{v} - \sum_{k=0}^{+\infty} \eta^{2k-3} \mathfrak{G}^{2k} \cdot \mathcal{F}^{2k} \boldsymbol{v}. \end{aligned} $$

Integrating this expansion with respect to $\eta$, we obtain

$$ \begin{aligned} I(\eta) \sim \eta^{-d} \mathcal{F}^\infty \boldsymbol{v} & - \frac{\delta_{d=2}}{4\pi} \left( \eta^{-2} \log \eta + \frac{1}{2} \eta^{-2} \right) \mathcal{F}^0 \boldsymbol{v} + \frac{1}{2} \eta^{-2} \mathfrak{G}^0 \cdot \mathcal{F}^0 \boldsymbol{v} + \text{f. p. } \int_{\eta^{-1}P \setminus T} \boldsymbol{v}(y) \, dy \\ & - \log \eta \mathfrak{G}^2 \cdot \mathcal{F}^2 \boldsymbol{v} - \sum_{k=2}^{+\infty} \frac{\eta^{2k-2}}{2k-2} \mathfrak{G}^{2k} \cdot \mathcal{F}^{2k} \boldsymbol{v}. \end{aligned} $$
and

$$
\chi^* = - \frac{\log \eta}{4\pi} I - cI + \mathcal{F}^\infty \Phi + \sum_{p=2}^{+\infty} \eta^p M_p \text{ if } d=2,
$$

for some family of real $d \times d$ matrices $(M_p)_{p \ge 2}$.

Proof. Remembering the definition (3.1) of the potential $S_{\text{per},\eta T}^{\kappa\eta}$ and integrating $\mathbf{X}_i(x) = S_{\text{per},\eta T}^{\kappa\eta}[[\sigma(\mathbf{X}_i, \alpha_i)\mathbf{n}]](x)$ over the cell $P$ yields

$$
e_j \cdot \int_P \mathbf{X}_i(x) \, dx = \mathcal{X}_{ij}^* + \int_{\eta T} e_j \cdot \mathbf{X}_i(x) \, dx = \eta^{2-d} \left( -\frac{\log \eta}{4\pi} \delta_{ij} \delta_{d=2} - c \delta_{ij} + \kappa_{ij}^\eta \right), \quad (4.1)
$$

where we used $\int_P G_{\text{per}} dx = 0$. For $x \in \eta T$, it can be verified that the asymptotic expansion (3.16) is still valid and yields

$$
\int_{\eta T} \mathbf{X}_i(x) \, dx = \eta^d \left( \sum_{p=2}^{+\infty} \int_T \mathbf{v}_{i,p}(y) \, dy + \sum_{p=0}^{+\infty} \int_T \mathbf{w}_{i,p}(\eta y) \, dy \right). \quad (4.2)
$$

Due to the uniform convergence of the series of analytic function $\sum ||w_{i,p}||_{L^\infty(P, \mathbb{R}^d)} < +\infty$ (equation (3.17)), the series $\sum_{p=0}^{+\infty} w_{i,p}(\eta y)$ is analytic in $\eta$. Therefore, (4.2) implies that $\int_{\eta T} \mathbf{X}_i(x) \, dx = \eta^d \mathbf{F}_i(\eta)$ for some analytic function $\mathbf{F}_i$. Then, substituting $\kappa^\eta$ by its analytic series representation (3.13) yields the result. $\square$

In the sequel, we aim at identifying explicitly the constant matrices $(M_p)_{p \ge 2}$. The Cauchy-Schwartz inequality
implies

$$
\left| \int_{P \setminus (\eta T)} \mathbf{X}_i \cdot \mathbf{e}_j \, dx - \int_{P \setminus (\eta T)} \mathbf{X}_i^N \cdot \mathbf{e}_j \, dx \right| \le (1 - \eta^d |T|)^{\frac{1}{2}} \|\mathbf{X}_i - \mathbf{X}_i^N\|_{L^2(P \setminus (\eta T), \mathbb{R}^d)} \quad \text{for all } 1 \le i, j \le d. \quad (4.3)
$$

This inequality together with the bounds (3.22) and (3.23) imply that the asymptotic expansion of the permeability matrix $\mathcal{X}^*$ can be obtained by simply integrating (3.16) on $P\setminus(\eta T)$:

$$
\begin{equation}
\begin{split}
\chi_{ij}^* = \int_{P \setminus (\eta T)} \mathbf{X}_i(x) \cdot \mathbf{e}_j \, dx ={}& \eta^{2-d} \int_{P \setminus (\eta T)} \mathbf{\Phi}_i(x/\eta) \cdot \mathbf{e}_j \, dx \\
& + \sum_{p=2}^{+\infty} \eta^p \int_{P \setminus (\eta T)} \mathbf{v}_{i,p}(x/\eta) \cdot \mathbf{e}_j \, dx + \sum_{p=0}^{+\infty} \eta^p \int_{P \setminus (\eta T)} \mathbf{w}_{i,p}(x) \cdot \mathbf{e}_j \, dx.
\end{split}
\tag{4.4}
\end{equation}
$$

The terms involving the interior solutions $(w_{i,p})_{p \ge 0}$ can be easily expanded in powers of $\eta$ using their Taylor
series expansion at zero:

$$
\int_{P \setminus (\eta T)} \mathbf{w}_{i,p}(x) \, dx = \int_{P} \mathbf{w}_{i,p}(x) \, dx - \eta^d \int_{T} \mathbf{w}_{i,p}(\eta y) \, dy = \int_{P} \mathbf{w}_{i,p}(x) \, dx - \sum_{k=0}^{+\infty} \eta^{d+k} \frac{1}{k!} \nabla^k \mathbf{w}_{i,p}(0) \cdot \int_{T} y^k \, dy. \quad (4.5)
$$

In order to obtain the asymptotic expansion of the terms involving the exterior solutions $\mathbf{\Phi}_i$ and $(\mathbf{v}_{i,p})_{p \ge 2}$, we
rely on the following lemma.

**Lemma 4.1.** Let $v \in H_{\text{loc}}^1(\mathbb{R}^d \setminus T)$ be a function with an asymptotic expansion at infinity of the form:

$$
\mathbf{v}(y) = \mathcal{F}^{\infty} \mathbf{v} + \sum_{k=0}^{+\infty} \nabla^k G(y) \cdot \mathcal{F}^k \mathbf{v} \quad \text{as } |y| \to +\infty. \qquad (4.6)
$$

Then, the following asymptotic expansion holds as $\eta \to 0$:

$$
\int_{P \setminus (\eta T)} \mathbf{v}(x/\eta) \, dx \sim \mathcal{F}^\infty \mathbf{v} - \frac{\delta_{d=2}}{4\pi} \left( \log \eta + \frac{1}{2} \right) \mathcal{F}^0 \mathbf{v} \\
+ \eta^{d-2} \left( \frac{1}{2} \mathfrak{G}^0 \cdot \mathcal{F}^0 \mathbf{v} + \eta^2 \text{ f. p. } \int_{\eta^{-1}P \setminus T} \mathbf{v}(y) \, dy - \sum_{k=2}^{+\infty} \frac{\eta^{2k}}{2k-2} \mathfrak{G}^{2k} \cdot \mathcal{F}^{2k} \mathbf{v} \right), \quad (4.7)
$$

where $\mathfrak{G}^{2k}$ is the tensor defined in (2.59) and f. p. denotes the Hadamard finite part integral.

*Proof.* Let us denote by $I(\eta)$ the quantity

$$
I(\eta) := \eta^{-d} \int_{P \setminus (\eta T)} \mathbf{v}(x/\eta) \, dx = \int_{\eta^{-1}P \setminus T} \mathbf{v}(y) \, dy. \quad (4.8)
$$

Applying a classical shape differentiation formula (e.g. [25, Theorem 5.2.2]), the derivative of $I(\eta)$ is given by

$$
\frac{dI}{d\eta}(\eta) = \frac{d}{d\eta} \int_{\eta^{-1}P \setminus T} \mathbf{v}(y) \, dy = \int_{\eta^{-1}\partial P} \mathbf{v}(y)(\boldsymbol{\theta}(y) \cdot \mathbf{n}(y)) \, d\sigma(y), \quad (4.9)
$$
The result follows by using (4.8) and $\mathfrak{G}^2 = 0$ (Lemma 2.2). $\square$

Remark 4.1. As visible from the proof, it came very close to having a logarithmic term in the asymptotic expansion (4.7): its absence is due to the property $\mathfrak{G}^2 = 0$ of Proposition 2.11.

Following the classical definition for one dimensional integrals [14], we define the Hadamard finite part integral $f.p. \eta \to 0^+$ $\int_{\eta^{-1}P \setminus T} v(y) dy$ arising in (4.7) as the constant $c_{0,0} \in \mathbb{R}^d$ such that

$$ \int_{\eta^{-1}P \setminus T} v \, dy = \sum_{p=0}^{P} \sum_{k=-M}^{0} c_{p,k} \eta^k \log \eta^p + o(1) \text{ as } \eta \to 0^+. $$

The next lemma shows that an explicit expression for the Hadamard finite part integral can be found in some particular situations.

**Lemma 4.2.** (i) If $v = S_T[\phi]$ with $\phi \in H^{-\frac{1}{2}}(\partial T, \mathbb{R}^d)$. Then,

$$ \text{f.p.} \int_{\eta^{-1}P \setminus T} S_T[\phi](y) dy = -\frac{1}{2} \mathfrak{R}^2 \cdot \int_{\partial T} \phi(y)y^2 d\sigma(y) - \int_T S_T[\phi](y) dy. \quad (4.10) $$

$$ (ii) \text{ f.p.} \int_{\eta^{-1}P \setminus T} dx = -|T|. $$

*Proof.* (i) By the definition (2.16) of the single layer potential,

$$ \begin{aligned} \int_{\eta^{-1}P \setminus T} S_T[\phi](y) dy &= \int_{\eta^{-1}P} S_T[\phi](y) dy - \int_T S_T[\phi](y) dy \\ &= \int_{\partial T} \left( \int_{\eta^{-1}P} G(y - y') dy \right) \phi(y') d\sigma(y') - \int_T S_T[\phi](y) dy. \end{aligned} \quad (4.11) $$

Euler's identity $\mathrm{div}(xG) = 2G(x) + \delta_{d=2} \frac{1}{4\pi} I$ and the divergence theorem enable us to rewrite the integral on $\eta^{-1}P$ as:

$$ \begin{aligned} \int_{\eta^{-1}P} G(y - y') dy &= \int_{\eta^{-1}P} \left( \frac{1}{2} \mathrm{div}_y((y - y')G(y - y')) - \frac{1}{8\pi} \delta_{d=2} I \right) dy \\ &= \frac{1}{2} \int_{\eta^{-1}\partial P} (y - y') \cdot n(y) G(y - y') d\sigma(y) - \frac{1}{8\pi} \delta_{d=2} \eta^{-d} I. \end{aligned} \quad (4.12) $$

Then, using a change of variable and the far field expansion of $G$, the first integral of the above equality reads

$$ \begin{aligned} \frac{1}{2} \int_{\eta^{-1}\partial P} (y - y') \cdot n(y) G(y - y') d\sigma(y) &= \frac{1}{2} \eta^{1-d} \int_{\partial P} (x/\eta - y') \cdot n(x) G(x/\eta - y') d\sigma(x) \\ &= \frac{1}{2} \eta^{1-d} \int_{\partial P} (x/\eta - y') \cdot n(x) \times \sum_{k=0}^{+\infty} \eta^{d-2+k} \frac{(-1)^k}{k!} \nabla^k G(x) \cdot y'^k d\sigma(x) \\ &= \frac{1}{2} \sum_{k=0}^{+\infty} \frac{(-1)^k}{k!} \eta^{-2+k} \int_{\partial P} (x \cdot n) \nabla^k G(x) \cdot y'^k d\sigma(x) - \frac{1}{2} \sum_{k=0}^{+\infty} \frac{(-1)^k}{k!} \eta^{-1+k} \int_{\partial P} (y' \cdot n(x)) \nabla^k G(x) \cdot y'^k d\sigma(x). \end{aligned} \quad (4.13) $$

Using (4.12) and (4.13), the finite part integral corresponds to the zero-th order term in the expansion of $\int_{\eta^{-1}P \setminus T} S_T[\phi](y) dy$ as $\eta \to 0$, which is

$$ \text{f.p.} \int_{\eta^{-1}P \setminus T} S_T[\phi](y) dy = \frac{1}{8} \int_{\partial T} \int_{\partial P} \nabla^2 G(x) d\sigma(x) \cdot y'^2 \phi(y') d\sigma(y') \\ + \frac{1}{2} \int_{\partial T} \int_{\partial P} (y' \cdot n(x)) \nabla G(x) d\sigma(x) \cdot y' \phi(y') d\sigma(y') - \int_T S_T[\phi](y) dy. $$

The first term is null due to $\mathfrak{G}^2 = 0$ while the second term is $-\frac{1}{2} \mathfrak{R}^2 \cdot \int_{\partial T} \phi(y) y^2 d\sigma(y)$ (Proposition 2.11). This yields the result.

$$ (ii) \text{ This is immediate since } \int_{\eta^{-1}P \setminus T} dx = \eta^{-d} - |T|. \square $$

By applying Lemma 4.1 to (4.4) and using (4.5), we obtain the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$.
where $\theta(y)$ is the “velocity” of the moving boundary at $y \in \eta^{-1}\partial P$. Since a point of the moving boundary $\eta^{-1}\partial P$ can be parameterized as $y = \eta^{-1}y_0$ for a fixed $y_0 \in \partial P$, the velocity is given by $\theta(y) = -\eta^{-2}\mathbf{y}_0 = -\eta^{-1}\mathbf{y}$. Substituting into (4.9) and performing successive changes of variables, we obtain

$$ \frac{dI}{d\eta}(\eta) = -\eta^{-1} \int_{\eta^{-1}\partial P} \mathbf{v}(y)(\mathbf{y} \cdot \mathbf{n}) \, d\sigma(y) = -\eta^{-1}\eta^{1-d}\eta^{-1} \int_{\partial P} \mathbf{v}(x/\eta)(\mathbf{x} \cdot \mathbf{n}) \, d\sigma(x). $$

Using the asymptotic expansion (4.6), $\mathbf{x} \cdot \mathbf{n} = \frac{1}{2}$ on $\partial P$ and the fact that $\nabla^k \Gamma$ is homogeneous of degree $2-d-k$ (or “associated homogeneous” if $d=2$ and $k=0$), we find

$$ \begin{aligned} \frac{d}{d\eta} I(\eta) &\sim -\eta^{-1-d} \frac{|\partial P|}{2} \mathcal{F}^{\infty} \mathbf{v} + \frac{1}{4\pi} \frac{|\partial P|}{2} \eta^{-3} \log \eta \mathcal{F}^0 \mathbf{v} \delta_{d=2} - \eta^{-1-d} \sum_{k=0}^{+\infty} \eta^{d-2+k} \left( \int_{\partial P} \nabla^k G(x) \mathbf{x} \cdot \mathbf{n} \, d\sigma \right) \cdot \mathcal{F}^k \mathbf{v} \\ &\sim -d\eta^{-1-d} \mathcal{F}^{\infty} \mathbf{v} + \frac{\delta_{d=2}}{2\pi} \eta^{-3} \log \eta \mathcal{F}^0 \mathbf{v} - \sum_{k=0}^{+\infty} \eta^{2k-3} \mathfrak{G}^{2k} \cdot \mathcal{F}^{2k} \mathbf{v}. \end{aligned} $$

Integrating this expansion with respect to $\eta$, we obtain

$$ \begin{aligned} I(\eta) &\sim \eta^{-d} \mathcal{F}^{\infty} \mathbf{v} - \frac{\delta_{d=2}}{4\pi} \left( \eta^{-2} \log \eta + \frac{1}{2} \eta^{-2} \right) \mathcal{F}^0 \mathbf{v} + \frac{1}{2} \eta^{-2} \mathfrak{G}^0 \cdot \mathcal{F}^0 \mathbf{v} + \text{f. p.} \int_{\eta^{-1}P \setminus T} \mathbf{v}(y) \, dy \\ &\quad - \log \eta \mathfrak{G}^2 \cdot \mathcal{F}^2 \mathbf{v} - \sum_{k=2}^{+\infty} \frac{\eta^{2k-2}}{2k-2} \mathfrak{G}^{2k} \cdot \mathcal{F}^{2k} \mathbf{v}. \end{aligned} $$

The result follows by using (4.8) and $\mathfrak{G}^2 = 0$ (Lemma 2.2). $\square$

**Remark 4.1.** As visible from the proof, it came very close to having a logarithmic term in the asymptotic expansion (4.7): its absence is due to the property $\mathfrak{G}^2 = 0$ of Proposition 2.11.

Following the classical definition for one dimensional integrals [15], we define the Hadamard finite part integral $\text{f. p.}_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} \mathbf{v}(y) \, dy$ arising in (4.7) as the constant $c_{0,0} \in \mathbb{R}^d$ such that

$$ \int_{\eta^{-1}P \setminus T} \mathbf{v} \, dy = \sum_{p=0}^{P} \sum_{k=-M}^{0} c_{p,k} \eta^k \log \eta^p + o(1) \text{ as } \eta \to 0^+. $$

The next lemma shows that an explicit expression for the Hadamard finite part integral can be found in some particular situations.

**Lemma 4.2.** (i) If $\mathbf{v} = S_T[\phi]$ with $\phi \in H^{-\frac{1}{2}}(\partial T, \mathbb{R}^d)$. Then,

$$ \text{f. p.} \int_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} S_T[\phi](y) \, dy = -\frac{1}{2} \mathfrak{R}^2 \cdot \int_{\partial T} \phi(y) y^2 \, d\sigma(y) - \int_T S_T[\phi](y) \, dy. \quad (4.10) $$

$$ (ii) \text{ f. p.} \int_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} dx = -|T|. $$

*Proof.* (i) By the definition (2.16) of the single layer potential,

$$ \begin{aligned} \int_{\eta^{-1}P \setminus T} S_T[\phi](y) \, dy &= \int_{\eta^{-1}P} S_T[\phi](y) \, dy - \int_T S_T[\phi](y) \, dy \\ &= \int_{\partial T} \left( \int_{\eta^{-1}P} G(y - y') \, dy \right) \phi(y') \, d\sigma(y') - \int_T S_T[\phi](y) \, dy. \end{aligned} \quad (4.11) $$

Euler's identity $\text{div}(\mathbf{x}G) = 2G(x) + \delta_{d=2} \frac{1}{4\pi} I$ and the divergence theorem enable us to rewrite the integral on $\eta^{-1}P$ as:

$$ \begin{aligned} \int_{\eta^{-1}P} G(y - y') \, dy &= \int_{\eta^{-1}P} \left( \frac{1}{2} \text{div}_y((\mathbf{y} - \mathbf{y}')G(y - y')) - \frac{1}{8\pi} \delta_{d=2} I \right) \, dy \\ &= \frac{1}{2} \int_{\eta^{-1}\partial P} (\mathbf{y} - \mathbf{y}') \cdot \mathbf{n}(y) G(y - y') \, d\sigma(y) - \frac{1}{8\pi} \delta_{d=2} \eta^{-d} I. \end{aligned} \quad (4.12) $$
**Proposition 4.2.** If $d \ge 3$, the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$ defined in (1.5) reads, as $\eta \to 0^+$:

$$ \begin{aligned} \mathcal{X}_{ij}^* \sim & \eta^{2-d} M_{ij}^{-1} - c \delta_{ij} - \sum_{p=1}^{+\infty} \eta^{2p} e_j \cdot \nabla^{2p} R(0) \cdot \mathcal{F}^{2p} \mathbf{\Phi}_i \\ & + \sum_{p=2}^{+\infty} \eta^p e_j \cdot \mathbf{w}_{i,p}(0) - \sum_{p=2}^{+\infty} \eta^{d+p} \sum_{k=0}^{p} \left( \frac{1}{k!} \int_T y^k dy \cdot \nabla^k \mathbf{w}_{i,p-k}(0) \right) \cdot e_j \tag{4.14} \\ & - \sum_{p=2}^{+\infty} \eta^{p+d} \int_T e_j \cdot \mathbf{v}_{i,p}(y) dy - \sum_{p=2}^{+\infty} \eta^{p+d} \sum_{k=2}^{p} e_j \cdot \nabla^k R(0) \cdot \mathcal{F}^k \mathbf{v}_{i,p-k+2}. \end{aligned} $$

for any $1 \le i, j \le d$, where $c$ is the universal constant arising in (2.47).

If $d = 2$, the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$ reads instead

$$ \begin{aligned} \mathcal{X}_{ij}^* \sim & -\frac{\log \eta}{4\pi} \delta_{ij} - c\delta_{ij} + \mathcal{F}^\infty \mathbf{\Phi}_i \cdot e_j - \sum_{p=1}^{+\infty} \eta^{2p} e_j \cdot \nabla^{2p} R(0) \cdot \mathcal{F}^{2p} \mathbf{\Phi}_i \\ & + \sum_{p=1}^{+\infty} \eta^{2p} e_j \cdot \mathbf{w}_{i,2p}(0) - \sum_{p=1}^{+\infty} \eta^{2p+2} \sum_{k=0}^{2p} \left( \frac{1}{k!} \int_T y^k dy \cdot \nabla^k \mathbf{w}_{i,2p-k}(0) \right) \cdot e_j \tag{4.15} \\ & - \sum_{p=1}^{+\infty} \eta^{2p+2} \int_T e_j \cdot \mathbf{v}_{i,2p}(y) dy - \sum_{p=1}^{+\infty} \eta^{2p+2} \sum_{k=1}^{p} e_j \cdot \nabla^{2k} R(0) \cdot \mathcal{F}^{2k} \mathbf{v}_{i,2p-2k+2}. \end{aligned} $$

*Proof.* Applying Lemma 4.1 to (4.4) using $\mathcal{F}^0 \mathbf{\Phi}_i = e_i$ and $\mathfrak{G}^{2k}/(2k-2) = \mathfrak{R}^{2k}$ (Proposition 2.11), we obtain,

$$ \begin{aligned} \mathcal{X}_{ij}^* \sim & \eta^{2-d} \mathcal{F}^\infty \mathbf{\Phi}_i \cdot e_j - \frac{\delta_{d=2}}{4\pi} \left( \log \eta \delta_{ij} + \frac{1}{2} \right) \delta_{ij} + \frac{1}{2} \mathfrak{G}_{ij}^0 + \eta^2 \operatorname{f.p.} \int_{\eta^{-1}P\setminus T} \mathbf{\Phi}_i(y) \cdot e_j \, dy \\ & - \sum_{p=2}^{+\infty} \eta^{2p} e_j \cdot \mathfrak{R}^{2p} \cdot \mathcal{F}^{2p} \mathbf{\Phi}_i \\ & + \sum_{p=2}^{+\infty} \eta^{p+d} \left( \operatorname{f.p.} \int_{\eta^{-1}P\setminus T} e_j \cdot \mathbf{v}_{i,p}(y) \, dy - \sum_{k=2}^{+\infty} \eta^{2k-2} e_j \cdot \mathfrak{R}^{2k} \cdot \mathcal{F}^{2k} \mathbf{v}_{i,p} \right) \\ & + \sum_{p=0}^{+\infty} \eta^p \int_P e_j \cdot \mathbf{w}_{i,p} \, dx - \sum_{p=0}^{+\infty} \eta^{d+p} \sum_{k=0}^{p} \frac{1}{k!} \int_T y^k dy \cdot \nabla^k \mathbf{w}_{i,p-k}(0) \cdot e_j. \end{aligned} \tag{4.16} $$

Using (3.27) and Proposition 2.11, we find that the average of $\mathbf{w}_{i,p}$ on the cell $P$ can be expressed as

$$ \int_P \mathbf{w}_{i,p} \, dx = \mathbf{w}_{i,p}(0) + (\mathfrak{R}^p - \nabla^p R(0)) \cdot \mathcal{F}^p \mathbf{\Phi}_i + \sum_{k=d-1}^{p-2} (\mathfrak{R}^{k-d+2} - \nabla^{k-d+2} R(0)) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k}. \tag{4.17} $$

Moreover, the finite part integrals of $\mathbf{\Phi}_i = S_T[\sigma(\mathbf{\Phi}_i, \sigma_i)\mathbf{n}|_+] + \mathcal{F}^\infty \mathbf{\Phi}_i$ and $\mathbf{v}_{i,p} = S_T[[\sigma(\mathbf{v}_{i,p}, \xi_{i,p})\mathbf{n}]]$ can be computed by using Lemma 4.2:

$$ \operatorname{f.p.} \int_{\eta^{-1}P\setminus T} \mathbf{\Phi}_i(y) \, dy = -|T| \mathcal{F}^\infty \mathbf{\Phi}_i - \mathfrak{R}^2 \cdot \mathcal{F}^2 \mathbf{\Phi}_i - \int_T (\mathbf{\Phi}_i - \mathcal{F}^\infty \mathbf{\Phi}_i) \, dy = -\mathfrak{R}^2 \cdot \mathcal{F}^2 \mathbf{\Phi}_i, \tag{4.18} $$

and

$$ \operatorname{f.p.} \int_{\eta^{-1}P\setminus T} \mathbf{v}_{i,p}(y) \, dy = -\mathfrak{R}^2 \cdot \mathcal{F}^2 \mathbf{v}_{i,p} - \int_T \mathbf{v}_{i,p}(y) \, dy. \tag{4.19} $$

Then, we observe that, since $\mathcal{F}^0 \mathbf{v}_{i,p} = 0$ for any $p \ge 2$ and $\mathfrak{R}^1 = 0$:

$$ \begin{aligned} & \sum_{p=0}^{+\infty} \eta^p \sum_{k=d-1}^{p-2} \mathfrak{R}^{k-d+2} \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} = \sum_{p=d}^{+\infty} \eta^p \sum_{k=d-2}^{p-2} \mathfrak{R}^{k-d+2} \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} = \sum_{p=0}^{+\infty} \eta^{p+d} \sum_{k=0}^{p} \mathfrak{R}^k \cdot \mathcal{F}^k \mathbf{v}_{i,p-k+2} \\ & = \sum_{p=0}^{+\infty} \sum_{k=0}^{+\infty} \eta^{p+d+k} \mathfrak{R}^k \cdot \mathcal{F}^k \mathbf{v}_{i,p+2} = \sum_{p=2}^{+\infty} \eta^{p+d} \sum_{k=1}^{+\infty} \eta^{2k-2} \mathfrak{R}^{2k} \cdot \mathcal{F}^{2k} \mathbf{v}_{i,p} \end{aligned} \tag{4.20} $$
Then, using a change of variable and the far field expansion of G, the first integral of the above equality reads

$$
\begin{align}
\frac{1}{2} \int_{\eta^{-1} \partial P} (\mathbf{y} - \mathbf{y}') \cdot \mathbf{n}(y) G(y - y') \, d\sigma(y) &= \frac{1}{2} \eta^{1-d} \int_{\partial P} (\mathbf{x}/\eta - \mathbf{y}') \cdot \mathbf{n}(x) G(x/\eta - y') \, d\sigma(x) \nonumber \\
&= \frac{1}{2} \eta^{1-d} \int_{\partial P} (\mathbf{x}/\eta - \mathbf{y}') \cdot \mathbf{n}(x) \times \sum_{k=0}^{+\infty} \eta^{d-2+k} \frac{(-1)^k}{k!} \nabla^k G(x) \cdot y'^k \, d\sigma(x) \nonumber \\
&= \frac{1}{2} \sum_{k=0}^{+\infty} \frac{(-1)^k}{k!} \eta^{-2+k} \int_{\partial P} (\mathbf{x} \cdot \mathbf{n}) \nabla^k G(x) \cdot y'^k \, d\sigma(x) - \frac{1}{2} \sum_{k=0}^{+\infty} \frac{(-1)^k}{k!} \eta^{-1+k} \int_{\partial P} (\mathbf{y}' \cdot \mathbf{n}(x)) \nabla^k G(x) \cdot y'^k \, d\sigma(x). \tag{4.13}
\end{align}
$$

Using (4.12) and (4.13), the finite part integral corresponds to the zero-th order term in the expansion of
$\int_{\eta^{-1}P\setminus T} \mathcal{S}_T[\phi](y) dy$ as $\eta \to 0$, which is

$$
\begin{align*}
\underset{\eta \to 0^+}{\text{f. p.}} \int_{\eta^{-1}P \setminus T} \mathcal{S}_T[\phi](y) \, dy &= \frac{1}{8} \int_{\partial T} \int_{\partial P} \nabla^2 G(x) \, d\sigma(x) \cdot y'^2 \phi(y') \, d\sigma(y') \\
&\quad + \frac{1}{2} \int_{\partial T} \int_{\partial P} (\mathbf{y}' \cdot \mathbf{n}(x)) \nabla G(x) \, d\sigma(x) \cdot y' \phi(y') \, d\sigma(y') - \int_T \mathcal{S}_T[\phi](y) \, dy.
\end{align*}
$$

The first term is null due to $\mathfrak{G}^2 = 0$ while the second term is $-\frac{1}{2}\mathfrak{R}^2 \cdot \int_{\partial T} \phi(y)y^2 d\sigma(y)$ (Proposition 2.11). This
yields the result.

(ii) This is immediate since $\int_{\eta^{-1}P\setminus T} dx = \eta^{-d} - |T|$.

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

By applying Lemma 4.1 to (4.4) and using (4.5), we obtain the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$.

**Proposition 4.2.** If $d \ge 3$, the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$ defined in (1.5) reads, as $\eta \to 0^+$:

$$
\begin{align}
\mathcal{X}_{ij}^* \sim \eta^{2-d} M_{ij}^{-1} - c \delta_{ij} - \sum_{p=1}^{+\infty} \eta^{2p} e_j \cdot \nabla^{2p} R(0) \cdot \mathcal{F}^{2p} \Phi_i & \nonumber \\
& + \sum_{p=2}^{+\infty} \eta^p e_j \cdot \mathbf{w}_{i,p}(0) - \sum_{p=2}^{+\infty} \eta^{d+p} \sum_{k=0}^{p} \left( \frac{1}{k!} \int_T y^k \, dy \cdot \nabla^k \mathbf{w}_{i,p-k}(0) \right) \cdot e_j \tag{4.14} \\
& - \sum_{p=2}^{+\infty} \eta^{p+d} \int_T e_j \cdot \mathbf{v}_{i,p}(y) \, dy - \sum_{p=2}^{+\infty} \eta^{p+d} \sum_{k=2}^{p} e_j \cdot \nabla^k R(0) \cdot \mathcal{F}^k \mathbf{v}_{i,p-k+2}, \nonumber
\end{align}
$$

for any $1 \le i, j \le d$, where $c$ is the universal constant arising in (2.47).

If $d=2$, the full asymptotic expansion of the permeability matrix $\mathcal{X}^*$ reads instead

$$
\begin{align}
\mathcal{X}_{ij}^* \sim & -\frac{\log \eta}{4\pi} \delta_{ij} - c\delta_{ij} + \mathcal{F}^\infty \Phi_i \cdot e_j - \sum_{p=1}^{+\infty} \eta^{2p} e_j \cdot \nabla^{2p} R(0) \cdot \mathcal{F}^{2p} \Phi_i \nonumber \\
& + \sum_{p=1}^{+\infty} \eta^{2p} e_j \cdot \mathbf{w}_{i,2p}(0) - \sum_{p=1}^{+\infty} \eta^{2p+2} \sum_{k=0}^{2p} \left( \frac{1}{k!} \int_T y^k \, dy \cdot \nabla^k \mathbf{w}_{i,2p-k}(0) \right) \cdot e_j \tag{4.15} \\
& - \sum_{p=1}^{+\infty} \eta^{2p+2} \int_T e_j \cdot \mathbf{v}_{i,2p}(y) \, dy - \sum_{p=1}^{+\infty} \eta^{2p+2} \sum_{k=1}^{p} e_j \cdot \nabla^{2k} R(0) \cdot \mathcal{F}^{2k} \mathbf{v}_{i,2p-2k+2}. \nonumber
\end{align}
$$

*Proof.* Applying Lemma 4.1 to (4.4) using $\mathcal{F}^0 \Phi_i = e_i$ and $\mathfrak{G}^{2k}/(2k-2) = \mathfrak{R}^{2k}$ (Proposition 2.11), we obtain,

$$
\begin{align}
\mathcal{X}_{ij}^* \sim & \eta^{2-d} \mathcal{F}^\infty \Phi_i \cdot e_j - \frac{\delta_{d=2}}{4\pi} \left( \log \eta \delta_{ij} + \frac{1}{2} \right) \delta_{ij} + \frac{1}{2} \mathfrak{G}_{ij}^0 + \eta^2 \int_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} \Phi_i(y) \cdot e_j \, dy \nonumber \\
& - \sum_{p=2}^{+\infty} \eta^{2p} e_j \cdot \mathfrak{R}^{2p} \cdot \mathcal{F}^{2p} \Phi_i \nonumber \\
& + \sum_{p=2}^{+\infty} \eta^{p+d} \left( \text{f. p.} \int_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} e_j \cdot \mathbf{v}_{i,p}(y) \, dy - \sum_{k=2}^{+\infty} \eta^{2k-2} e_j \cdot \mathfrak{R}^{2k} \cdot \mathcal{F}^{2k} \mathbf{v}_{i,p} \right) \tag{4.16} \\
& + \sum_{p=0}^{+\infty} \eta^p \int_P e_j \cdot \mathbf{w}_{i,p} \, dx - \sum_{p=0}^{+\infty} \eta^{d+p} \sum_{k=0}^{p} \frac{1}{k!} \int_T y^k \, dy \cdot \nabla^k \mathbf{w}_{i,p-k}(0) \cdot e_j. \nonumber
\end{align}
$$
Consequently, substituting (4.17) into (4.16) and using the previous identities together with $\mathfrak{R}^0 = -\frac{1}{2}\mathfrak{G}^0 - cI + \frac{1}{8\pi}\delta_{d=2}I$, we find that all the terms involving the tensors $\mathfrak{R}^k$ disappear

$$ \begin{aligned} \mathcal{X}_{ij}^* \sim & \eta^{2-d} \mathcal{F}^\infty \mathbf{\Phi}_i \cdot \mathbf{e}_j - \frac{\delta_{d=2}}{4\pi} \log \eta \delta_{ij} - c \delta_{ij} - \sum_{p=2}^{+\infty} \eta^{p+d} \int_T e_j \cdot \mathbf{v}_{i,p}(y) \, dy + \sum_{p=0}^{+\infty} \eta^p e_j \cdot (\mathbf{w}_{i,p}(0) - \nabla^p R(0) \cdot \mathcal{F}^p \mathbf{\Phi}_i) \\ & - \sum_{p=0}^{+\infty} \eta^{p+d} \sum_{k=0}^{p} e_j \cdot \nabla^k R(0) \cdot \mathcal{F}^k \mathbf{v}_{i,p-k+2} - \sum_{p=0}^{+\infty} \eta^{d+p} \sum_{k=0}^{p} \frac{1}{k!} \int_T y^k \, dy \cdot \nabla^k \mathbf{w}_{i,p-k}(0) \cdot e_j. \end{aligned} \quad (4.21) $$

The final result follows in dimension $d \ge 3$ by using $\mathbf{w}_{i,0}(0) = \mathbf{w}_{i,1}(0) = 0$, $\nabla \mathbf{w}_{i,0}(0) = 0$ and the parity properties of $R$. In dimension $d = 2$, all terms featuring odd powers of $\eta$ vanish due to $\mathbf{v}_{i,2p+1} = 0$ for all $p \ge 2$ (Proposition 3.6). $\square$

**Remark 4.2.** All the coefficients of the asymptotic expansions (4.14) and (4.15) can be obtained numerically from the computation of

* the far field coefficients $(\mathcal{F}^k \mathbf{\Phi}_i)_{0 \le k \le \infty}$ and $(\mathcal{F}^k \mathbf{v}_{i,p})_{k \ge 0}$ for $1 \le i \le d$ and $p \ge 2$. These can be obtained by solving the exterior problem with an integral equation method [29, 31];

* the integrals in the interior domain $\int_T e_j \cdot \mathbf{v}_{i,p} \, dy$, which can be obtained either from the layer potential representation of $\mathbf{v}_{i,p}$ (after having solved the exterior problem) or by solving the interior problem for instance with the finite element method;

* the derivatives $\nabla^k R(0)$ for any $k \ge 0$, which also occur in the evaluation of $\mathbf{w}_{i,p}(0)$ with formulas (3.27) and (3.28). These can be found from the asymptotic expansion at zero of the periodic Green matrix $G_{\text{per}}$:

$$ G_{\text{per}}(x) = G(x) + cI + \sum_{k=2}^{+\infty} \frac{1}{(2k)!} \nabla^{2k} R(0) \cdot x^{2k} \quad \text{as } x \to 0. $$

The coefficients of the expansion can be estimated numerically thanks to Ewald's method described in [15, 6, 23].

To conclude our analysis, we verify that formulas (4.14) and (4.15) agree at least up to the order $O(\eta^2)$ with the expansions (1.6) and (1.7) found by Hasimoto when obstacle $T$ is spherical.

**Proposition 4.3.** For a spherical obstacle $T = B(0, 1)$, the leading order asymptotic expansion of the permeability matrix $\mathcal{X}^*$ reads

$$ \mathcal{X}_{ij}^* = \left( \eta^{d-2} \frac{d-1}{d(d-2)|\partial B(0,1)|} - c + \frac{d-1}{d^2} \eta^2 \right) \delta_{ij} + O(\eta^5) \quad \text{if } d \ge 3. \quad (4.22) $$

and

$$ \mathcal{X}_{ij}^* = \left( -\frac{1}{4\pi} \log \eta \delta_{ij} - c + \frac{1}{8\pi} + \frac{1}{4} \eta^2 \right) \delta_{ij} + O(\eta^4) \quad \text{if } d=2. \quad (4.23) $$

**Proof.** Case $d \ge 3$. Computing the Laplacian of (2.4) yields

$$ \Delta \mathbf{G}_i(x) = -\nabla^2 \Gamma(x) e_i = -\frac{1}{|\partial B(0,1)|} \left( \frac{e_i}{|x|^d} - d \frac{x_i x}{|x|^{d+2}} \right). \quad (4.24) $$

From this identity, it is easily obtained that for a spherical obstacle, the function $\mathbf{\Psi}_j$ solving (2.26) is given by

$$ \mathbf{\Psi}_j(y) = -\frac{d(d-2)}{d-1} |\partial B(0,1)| \left( \mathbf{G}_j(y) + \frac{1}{2d} \Delta \mathbf{G}_j(y) \right), \quad y \in \mathbb{R}^d \setminus \overline{B(0,1)}. \quad (4.25) $$

This implies in particular that $M = \frac{d(d-2)}{d-1}|\partial B(0,1)|I$ and that the function $\mathbf{\Phi}_i$ of (2.28) is given by

$$ \mathbf{\Phi}_i = \frac{d-1}{d(d-2)} \frac{1}{|\partial B(0,1)|} e_i + \mathbf{G}_i + \frac{1}{2d} \Delta \mathbf{G}_i. \quad (4.26) $$

From this expression, it can be verified that

$$ \sigma(\mathbf{\Phi}_i, \sigma_i) \mathbf{n} = \frac{1}{|\partial B(0,1)|} e_i, \quad (4.27) $$

from where it can be deduced that $\mathcal{F}^0 \mathbf{\Phi}_i = e_i$, $\mathcal{F}^1 \mathbf{\Phi}_i = 0$, $\mathcal{F}^2 \mathbf{\Phi}_i = \frac{1}{2d} I e_i$ and $\mathcal{F}^3 \mathbf{\Phi}_i = 0$. From (3.27) and (3.28), we infer $\mathbf{w}_{i,2}(0) = -\frac{1}{2d} \Delta \mathbf{R}_i(0)$. Substituting into (4.14), we obtain

$$ \mathcal{X}_{ij}^* = \eta^{2-d} \frac{d-1}{d(d-2)|\partial B(0,1)|} \delta_{ij} - c \delta_{ij} - \frac{\eta^2}{d} \Delta R_{ij}(0) \eta^2 + O(\eta^4). $$

The result follows by using $\Delta R(0) = -\frac{d-1}{d}I$ (Proposition 2.10).
Using (3.27) and Proposition 2.11, we find that the average of $w_{i,p}$ on the cell $P$ can be expressed as

$$ \int_P w_{i,p} \, dx = w_{i,p}(0) + (\mathfrak{R}^p - \nabla^p R(0)) \cdot \mathcal{F}^p \Phi_i + \sum_{k=d-1}^{p-2} (\mathfrak{R}^{k-d+2} - \nabla^{k-d+2} R(0)) \cdot \mathcal{F}^{k-d+2} v_{i,p-k}. \quad (4.17) $$

Moreover, the finite part integrals of $\Phi_i = S_T[\sigma(\Phi_i, \sigma_i)n_+]+\mathcal{F}^\infty\Phi_i$ and $v_{i,p} = S_T[\llbracket\sigma(v_{i,p}, \xi_{i,p})n\rrbracket]$ can be computed by using Lemma 4.2:

$$ \text{f. p.} \quad \int_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} \Phi_i(y) \, dy = -|T| \mathcal{F}^\infty \Phi_i - \mathfrak{R}^2 \cdot \mathcal{F}^2 \Phi_i - \int_T (\Phi_i - \mathcal{F}^\infty \Phi_i) \, dy = -\mathfrak{R}^2 \cdot \mathcal{F}^2 \Phi_i, \quad (4.18) $$

and

$$ \text{f. p.} \quad \int_{\eta \to 0^+} \int_{\eta^{-1}P \setminus T} v_{i,p}(y) \, dy = -\mathfrak{R}^2 \cdot \mathcal{F}^2 v_{i,p} - \int_T v_{i,p}(y) \, dy. \quad (4.19) $$

Then, we observe that, since $\mathcal{F}^0 v_{i,p} = 0$ for any $p \ge 2$ and $\mathfrak{R}^1 = 0$:

$$ \begin{aligned} \sum_{p=0}^{+\infty} \eta^p \sum_{k=d-1}^{p-2} \mathfrak{R}^{k-d+2} \cdot \mathcal{F}^{k-d+2} v_{i,p-k} &= \sum_{p=d}^{+\infty} \eta^p \sum_{k=d-2}^{p-2} \mathfrak{R}^{k-d+2} \cdot \mathcal{F}^{k-d+2} v_{i,p-k} = \sum_{p=0}^{+\infty} \eta^{p+d} \sum_{k=0}^{p} \mathfrak{R}^k \cdot \mathcal{F}^k v_{i,p-k+2} \\ &= \sum_{p=0}^{+\infty} \sum_{k=0}^{+\infty} \eta^{p+d+k} \mathfrak{R}^k \cdot \mathcal{F}^k v_{i,p+2} = \sum_{p=2}^{+\infty} \eta^{p+d} \sum_{k=1}^{+\infty} \eta^{2k-2} \mathfrak{R}^{2k} \cdot \mathcal{F}^{2k} v_{i,p} \end{aligned} \quad (4.20) $$

Consequently, substituting (4.17) into (4.16) and using the previous identities together with $\mathfrak{R}^0 = -\frac{1}{2}\mathfrak{G}^0 - cI + \frac{1}{8\pi}\delta_{d=2}I$, we find that all the terms involving the tensors $\mathfrak{R}^k$ disappear

$$ \begin{aligned} \mathcal{X}_{ij}^* \sim & \eta^{2-d} \mathcal{F}^\infty \Phi_i \cdot e_j - \frac{\delta_{d=2}}{4\pi} \log \eta \delta_{ij} - c \delta_{ij} - \sum_{p=2}^{+\infty} \eta^{p+d} \int_T e_j \cdot v_{i,p}(y) \, dy + \sum_{p=0}^{+\infty} \eta^p e_j \cdot (w_{i,p}(0) - \nabla^p R(0) \cdot \mathcal{F}^p \Phi_i) \\ & - \sum_{p=0}^{+\infty} \eta^{p+d} \sum_{k=0}^{p} e_j \cdot \nabla^k R(0) \cdot \mathcal{F}^k v_{i,p-k+2} - \sum_{p=0}^{+\infty} \eta^{d+p} \sum_{k=0}^{p} \frac{1}{k!} \int_T y^k dy \cdot \nabla^k w_{i,p-k}(0) \cdot e_j. \end{aligned} \quad (4.21) $$

The final result follows in dimension $d \ge 3$ by using $w_{i,0}(0) = w_{i,1}(0) = 0$, $\nabla w_{i,0}(0) = 0$ and the parity properties of $R$. In dimension $d = 2$, all terms featuring odd powers of $\eta$ vanish due to $v_{i,2p+1} = 0$ for all $p \ge 2$ (Proposition 3.6). $\square$

**Remark 4.2.** All the coefficients of the asymptotic expansions (4.14) and (4.15) can be obtained numerically from the computation of

* the far field coefficients $(\mathcal{F}^k \Phi_i)_{0 \le k \le \infty}$ and $(\mathcal{F}^k v_{i,p})_{k \ge 0}$ for $1 \le i \le d$ and $p \ge 2$. These can be obtained by solving the exterior problem with an integral equation method [32, 35];

* the integrals in the interior domain $\int_T e_j \cdot v_{i,p} \, dy$, which can be obtained either from the layer potential representation of $v_{i,p}$ (after having solved the exterior problem) or by solving the interior problem for instance with the finite element method;

* the derivatives $\nabla^k R(0)$ for any $k \ge 0$, which also occur in the evaluation of $w_{i,p}(0)$ with formulas (3.27) and (3.28). These can be found from the asymptotic expansion at zero of the periodic Green matrix $G_{\text{per}}$:

$$ G_{\text{per}}(x) = G(x) + cI + \sum_{k=2}^{+\infty} \frac{1}{(2k)!} \nabla^{2k} R(0) \cdot x^{2k} \text{ as } x \to 0. $$

The coefficients of the expansion can be estimated numerically thanks to Ewald's method described in [16, 6, 24].

To conclude our analysis, we verify that formulas (4.14) and (4.15) agree at least up to the order $O(\eta^2)$ with the expansions (1.6) and (1.7) found by Hasimoto when obstacle $T$ is spherical.

**Proposition 4.3.** For a spherical obstacle $T = B(0,1)$, the leading order asymptotic expansion of the permeability matrix $\mathcal{X}^*$ reads

$$ \mathcal{X}_{ij}^* = \left( \eta^{d-2} \frac{d-1}{d(d-2)|\partial B(0,1)|} - c + \frac{d-1}{d^2} \eta^2 \right) \delta_{ij} + O(\eta^5) \quad \text{if } d \ge 3. \quad (4.22) $$

and

$$ \mathcal{X}_{ij}^* = \left( -\frac{1}{4\pi} \log \eta \delta_{ij} - c + \frac{1}{8\pi} + \frac{1}{4} \eta^2 \right) \delta_{ij} + O(\eta^4) \quad \text{if } d=2. \quad (4.23) $$
Case $d=2$. Formula (4.24) still holds, which allows to infer that the function $\Phi_i$ solving (2.29) is given by

$$ \Phi_i = G_i + \frac{1}{4}\Delta G_i + \frac{1}{8\pi}e_i. $$

Formula (4.27) remains valid, which implies $\mathcal{F}^\infty\Phi_i = \frac{1}{8\pi}e_i$, $\mathcal{F}^0\Phi_i = e_i$, $\mathcal{F}^2\Phi_i = \frac{1}{4}I$. Consequently, we have again $w_{i,2}(0) = -\frac{1}{4}\Delta R_i(0)$ and substituting into (4.15) yields

$$ \chi_{ij}^* = -\frac{\log \eta}{4\pi} \delta_{ij} - c\delta_{ij} + \frac{1}{8\pi}\delta_{ij} - \frac{\eta^2}{2}\Delta R_{ij}(0) + O(\eta^4). $$

The formula (4.23) is obtained after using $\Delta R(0) = -\frac{1}{2}I$. $\square$

**Remark 4.3.** If $T = B(0,1)$, the far field coefficients $\mathcal{F}^k\Phi_i$ involve the spherical moments $\int_{\partial B(0,1)} y^{2k} dy$ (eq. (4.27)). These can be computed explicitly for all $k \in \mathbb{N}$, see Lemma C.1 in the appendix.

**Remark 4.4.** The reader can verify by dividing these asymptotic expansions for $d=3$ by $6\pi$ and for $d=2$ by $4\pi$ that they agree with the expressions (1.6) and (1.7) derived by Hasimoto.

**Remark 4.5.** The computation of higher order terms of the asymptotic expansion is more involved even for spherical obstacles, because it would require to compute the functions $v_{i,p}$ for $p \ge 2$. It would be possible to solve the defining exterior problem (3.20) by using Lamb-type representations relying on spherical harmonics [27].

## APPENDIX A. HIGH-ORDER ERROR ESTIMATES WITHOUT POTENTIAL THEORY

In this appendix, we give an alternative proof of the bounds of Proposition 3.5 which does not use potential theory, but directly derives the estimates (3.22) and (3.23) from the boundary value problems (3.19) and (3.20). The proof gives a less strong result, however, because it yields the following estimates with a constant $C_N$ depending a priori on $N$:

$$ \begin{aligned} & \|\mathcal{X}_i - \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} + \|\nabla \mathcal{X}_i - \nabla \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \\ & \qquad + \|\alpha_i - \alpha_i^N\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \le C_N \eta^{N+1} \quad \text{if } d \ge 3 \text{ or } (d=2 \text{ and } N=0), \quad \text{(A.1)} \end{aligned} $$

$$ \begin{aligned} & |\log \eta|^{-\frac{1}{2}} \|\mathcal{X}_i - \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} + \|\nabla \mathcal{X}_i - \nabla \mathcal{X}_i^N\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \\ & \qquad + \|\alpha_i - \alpha_i^N\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \le C_N \eta^{N+1} \quad \text{if } d=2 \text{ and } N \ge 1. \quad \text{(A.2)} \end{aligned} $$

Before stating the proof, we need to establish some energy bounds for the solution to periodic Stokes problem with a Dirichlet boundary condition on the small hole $\eta T$. Let us first recall the following lemma from [28, 1].

**Lemma A.1.** There exists a constant $C > 0$ independent of $\eta$ such that the following Poincaré inequalities hold for any $v \in H^1(P\setminus(\eta T))$ vanishing on $\partial(\eta T)$:

$$ \|v\|_{L^2(P\setminus(\eta T))} \le \begin{cases} C|\log \eta|^{\frac{1}{2}} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d=2, \\ C\eta^{1-d/2} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d \ge 3, \end{cases} \quad \text{(A.3)} $$

$$ \left| \int_{P\setminus(\eta T)} v \, dx \right| \le \begin{cases} C |\log \eta|^{\frac{1}{2}} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d=2, \\ C \eta^{1-d/2} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d \ge 3, \end{cases} \quad \text{(A.4)} $$

and

$$ \left\| v - \int_{P\setminus(\eta T)} v \, dx \right\|_{L^2(P\setminus(\eta T))} \le C \|\nabla v\|_{L^2(P\setminus(\eta T))}. \quad \text{(A.5)} $$

This allows us to prove the following result, which is the analogue for the Stokes system of what has already been established in [18, Proposition 5.6] for the periodic perforated Poisson equation. The proof needs some adaptations to account for the presence of the pressure.

**Proposition A.1.** Let $f \in C_{\text{per}}^\infty(P, \mathbb{R}^d)$, $g \in H^{\frac{1}{2}}(\partial(\eta T), \mathbb{R}^d)$, $h \in C^\infty(P\setminus T, \mathbb{R}^d)$ and $q \in C^\infty(P\setminus T)$ be functions satisfying the following assumptions:

$$ (i) \int_{\partial P} h \cdot n \, d\sigma - \int_{\partial(\eta T)} g \cdot n \, d\sigma = 0, $$

$$ (ii) -\Delta h + \nabla q \in C_{\text{per}}^{\infty}(P\setminus T, \mathbb{R}^d). $$
Proof. Case $d \ge 3$. Computing the Laplacian of (2.4) yields

$$ \Delta \mathbf{G}_i(x) = -\nabla^2 \Gamma(x) e_i = - \frac{1}{|\partial B(0,1)|} \left( \frac{e_i}{|x|^d} - d \frac{x_i x}{|x|^{d+2}} \right). \quad (4.24) $$

From this identity, it is easily obtained that for a spherical obstacle, the function $\Psi_j$ solving (2.26) is given by

$$ \Psi_j(y) = - \frac{d(d-2)}{d-1} |\partial B(0,1)| \left( \mathbf{G}_j(y) + \frac{1}{2d} \Delta \mathbf{G}_j(y) \right), \quad y \in \mathbb{R}^d \setminus \overline{B(0,1)}. \quad (4.25) $$

This implies in particular that $M = \frac{d(d-2)}{d-1} |\partial B(0,1)| I$ and that the function $\Phi_i$ of (2.28) is given by

$$ \Phi_i = \frac{d-1}{d(d-2)} \frac{1}{|\partial B(0,1)|} e_i + \mathbf{G}_i + \frac{1}{2d} \Delta \mathbf{G}_i. \quad (4.26) $$

From this expression, it can be verified that

$$ \sigma(\mathbf{\Phi}_i, \sigma_i) \mathbf{n} = \frac{1}{|\partial B(0,1)|} e_i, \quad (4.27) $$

from where it can be deduced that $\mathcal{F}^0 \mathbf{\Phi}_i = e_i$, $\mathcal{F}^1 \mathbf{\Phi}_i = 0$, $\mathcal{F}^2 \mathbf{\Phi}_i = \frac{1}{2d} I e_i$ and $\mathcal{F}^3 \mathbf{\Phi}_i = 0$. From (3.27) and (3.28), we infer $w_{i,2}(0) = -\frac{1}{2d} \Delta \mathbf{R}_i(0)$. Substituting into (4.14), we obtain

$$ \mathcal{X}_{ij}^* = \eta^{2-d} \frac{d-1}{d(d-2)|\partial B(0,1)|} \delta_{ij} - c \delta_{ij} - \frac{\eta^2}{d} \Delta R_{ij}(0) \eta^2 + O(\eta^4). $$

The result follows by using $\Delta R(0) = -\frac{d-1}{d} I$ (Proposition 2.10).

Case $d=2$. Formula (4.24) still holds, which allows to infer that the function $\Phi_i$ solving (2.29) is given by

$$ \mathbf{\Phi}_i = \mathbf{G}_i + \frac{1}{4} \Delta \mathbf{G}_i + \frac{1}{8\pi} e_i. $$

Formula (4.27) remains valid, which implies $\mathcal{F}^\infty \mathbf{\Phi}_i = \frac{1}{8\pi} e_i$, $\mathcal{F}^0 \mathbf{\Phi}_i = e_i$, $\mathcal{F}^2 \mathbf{\Phi}_i = \frac{1}{4} I$. Consequently, we have again $w_{i,2}(0) = -\frac{1}{4} \Delta \mathbf{R}_i(0)$ and substituting into (4.15) yields

$$ \mathcal{X}_{ij}^* = -\frac{\log \eta}{4\pi} \delta_{ij} - c \delta_{ij} + \frac{1}{8\pi} \delta_{ij} - \frac{\eta^2}{2} \Delta R_{ij}(0) + O(\eta^4). $$

The formula (4.23) is obtained after using $\Delta R(0) = -\frac{1}{2} I$. $\square$

Remark 4.3. If $T = B(0,1)$, the far field coefficients $\mathcal{F}^k \mathbf{\Phi}_i$ involve the spherical moments $\int_{\partial B(0,1)} y^{2k} dy$ (eq. (4.27)). These can be computed explicitly for all $k \in \mathbb{N}$, see Lemma C.1 in the appendix.

Remark 4.4. The reader can verify by dividing these asymptotic expansions for $d=3$ by $6\pi$ and for $d=2$ by $4\pi$ that they agree with the expressions (1.6) and (1.7) derived by Hasimoto.

Remark 4.5. The computation of higher order terms of the asymptotic expansion is more involved even for spherical obstacles, because it would require to compute the functions $v_{i,p}$ for $p \ge 2$. It would be possible to solve the defining exterior problem (3.20) by using Lamb-type representations relying on spherical harmonics [29].

## APPENDIX A. HIGH-ORDER ERROR ESTIMATES WITHOUT POTENTIAL THEORY

In this appendix, we give an alternative proof of the bounds of Proposition 3.5 which does not use potential theory, but directly derives the estimates (3.22) and (3.23) from the boundary value problems (3.19) and (3.20). The proof gives a less strong result, however, because it yields the following estimates with a constant $C_N$ depending a priori on $N$:

$$ \begin{aligned} & \|\mathbf{X}_i - \mathbf{X}_i^N\|_{L^2(P \setminus (\eta T), \mathbb{R}^d)} + \|\nabla \mathbf{X}_i - \nabla \mathbf{X}_i^N\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})} \\ & \qquad + \|\alpha_i - \alpha_i^N\|_{L^2(P \setminus (\eta T))/\mathbb{R}} \le C_N \eta^{N+1} \text{ if } d \ge 3 \text{ or } (d=2 \text{ and } N=0), \quad (A.1) \end{aligned} $$

$$ \begin{aligned} & |\log \eta|^{-\frac{1}{2}} \|\mathbf{X}_i - \mathbf{X}_i^N\|_{L^2(P \setminus (\eta T), \mathbb{R}^d)} + \|\nabla \mathbf{X}_i - \nabla \mathbf{X}_i^N\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})} \\ & \qquad + \|\alpha_i - \alpha_i^N\|_{L^2(P \setminus (\eta T))/\mathbb{R}} \le C_N \eta^{N+1} \text{ if } d=2 \text{ and } N \ge 1. \quad (A.2) \end{aligned} $$

Before stating the proof, we need to establish some energy bounds for the solution to periodic Stokes problem with a Dirichlet boundary condition on the small hole $\eta T$. Let us first recall the following lemma from [30, 1].
There exists a unique solution $(\mathbf{u}, p) \in H^1(P\setminus(\eta T), \mathbb{R}^d) \times L^2(P\setminus(\eta T))/\mathbb{R}$ to the boundary value problem

$$ egin{cases} 
-\Delta \mathbf{u} + \nabla p = \mathbf{f} \text{ in } P\setminus(\eta T), \\
\operatorname{div}(\mathbf{u}) = 0 \text{ in } P\setminus(\eta T), \\
\mathbf{u} = \mathbf{g} \text{ on } \partial(\eta T), \\
\mathbf{u} - \mathbf{h} \text{ is P-periodic}, \\
\sigma(\mathbf{u}, p)\mathbf{n}\mathbf{n}^T - \sigma(\mathbf{h}, q)\mathbf{n}\mathbf{n}^T \text{ is P-periodic.}
\end{cases} 
\quad (A.6) $$

Furthermore, for any $k > \frac{d}{2} - 1$, there exists a constant $C_k$ independent of $\eta$, $f$ and $h$ such that if $d \ge 3$,

$$ \begin{aligned} & \|\nabla \mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + \|p\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \\ & \le C_k \left[ \begin{aligned} & |\alpha(\mathbf{f}, \mathbf{h})|\eta^{1-d/2} + \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)}\eta^{d/2-1} \\ & + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \end{aligned} \right], \quad (A.7) \end{aligned} $$

and

$$ \begin{aligned} & \|\mathbf{u}\|_{L^2(P\setminus(\eta T)), \mathbb{R}^d} \le C_k \left[ \begin{aligned} & |\alpha(\mathbf{f}, \mathbf{h})|\eta^{2-d} + \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} \\ & + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \end{aligned} \right], \quad (A.8) \end{aligned} $$

where $\alpha(\mathbf{f}, \mathbf{h})$ is the “default of compatibility” defined by

$$ \alpha(\mathbf{f}, \mathbf{h}) := \int_P \mathbf{f} \, d\mathbf{x} + \int_{\partial P} \sigma(\mathbf{h}, q) \mathbf{n} \, d\sigma. \quad (A.9) $$

If $d=2$, the following bounds hold:

$$ \begin{aligned} & \|\nabla \mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + \|p\|_{L^2(P\setminus(\eta T))/\mathbb{R}} \\ & \le C_k \left[ \begin{aligned} & |\alpha(\mathbf{f}, \mathbf{h})| |\log \eta|^{\frac{1}{2}} + \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} \\ & + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \end{aligned} \right], \quad (A.10) \end{aligned} $$

and

$$ \begin{aligned} & \|\mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} \le C_k \left[ \begin{aligned} & |\alpha(\mathbf{f}, \mathbf{h})| |\log \eta| + \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} |\log \eta|^{\frac{1}{2}} \\ & + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \end{aligned} \right]. \quad (A.11) \end{aligned} $$

**Proof.** Existence and uniqueness of a solution $(\mathbf{u}, p)$ to (A.6) can be obtained by adapting the proof of Proposition 2.7, where no compatibility condition is required because of the Dirichlet boundary condition on $\partial(\eta T)$. We then decompose $\mathbf{u} = \mathbf{u}_1 + \mathbf{u}_2$ and $p = p_1 + p_2$ as the sums of the solutions to two distinct problems: $(\mathbf{u}_1, p_1) \in H^1(P, \mathbb{R}^d) \times L^2(P)/\mathbb{R}$ is defined as the unique solution to

$$ egin{cases} 
-\Delta \mathbf{u}_1 + \nabla p_1 = \mathbf{f} - \alpha(\mathbf{f}, \mathbf{h}) \text{ in } P, \\
\operatorname{div}(\mathbf{u}_1) = 0 \text{ in } P, \\
\mathbf{u}_1 - \mathbf{h} \text{ is P-periodic}, \\
\sigma(\mathbf{u}_1, p_1)\mathbf{n}\mathbf{n}^T - \sigma(\mathbf{h}, q)\mathbf{n}\mathbf{n}^T \text{ is P-periodic}, \\
\mathbf{u}_1(0) = 0.
\end{cases} 
\quad (A.12) $$

The problem (A.12) is of the type of (2.42), whose compatibility condition is satisfied by the definition (A.9) of $\alpha(\mathbf{f}, \mathbf{h})$. Furthermore, the condition $\mathbf{u}_1(0) = 0$ fixes the additive constant determining $\mathbf{u}_1$. By Propositions 2.7 and 2.8, $(\mathbf{u}_1, p_1)$ satisfies the following bound:

$$ \begin{aligned} & \|\mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^d)} + \|\nabla \mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^{d \times d})} + \|p_1\|_{L^\infty(P)/\mathbb{R}} \\ & \le C_k (\|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)}). \end{aligned} \quad (A.13) $$

Then $(\mathbf{u}_2, p_2)$ is defined as the solution to the difference problem

$$ egin{cases} 
-\Delta \mathbf{u}_2 + \nabla p_2 = \alpha(\mathbf{f}, \mathbf{h}) \text{ in } P, \\
\operatorname{div}(\mathbf{u}_2) = 0 \text{ in } P, \\
\mathbf{u}_2 = \mathbf{g} - \mathbf{u}_1 \text{ on } \partial(\eta T), \\
\mathbf{u}_2 \text{ is P-periodic}, \\
\sigma(\mathbf{u}_2, p_2)\mathbf{n}\mathbf{n}^T \text{ is P-periodic.}
\end{cases} 
\quad (A.14) $$

In order to estimate the energy norm of the solution to (A.14), we construct lifts of the boundary conditions $\hat{\mathbf{u}}_1$ and $\hat{\mathbf{g}}$. Let $K$ be an open set containing the obstacle $T: T \in K$. We consider an extension $\tilde{\mathbf{g}}_\eta$ in $K$ of the
**Lemma A.1.** There exists a constant $C > 0$ independent of $\eta$ such that the following Poincaré inequalities hold for any $v \in H^1(P\setminus(\eta T))$ vanishing on $\partial(\eta T)$:

$$ \|v\|_{L^2(P\setminus(\eta T))} \le \begin{cases} C|\log \eta|^{\frac{1}{2}} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d=2, \\ C\eta^{1-d/2} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d \ge 3, \end{cases} \quad (A.3) $$

$$ \left| \int_{P\setminus(\eta T)} v \, dx \right| \le \begin{cases} C |\log \eta|^{\frac{1}{2}} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d=2, \\ C \eta^{1-d/2} \|\nabla v\|_{L^2(P\setminus(\eta T))} & \text{if } d \ge 3, \end{cases} \quad (A.4) $$

and

$$ \left\| v - \int_{P\setminus(\eta T)} v \, dx \right\|_{L^2(P\setminus(\eta T))} \le C \|\nabla v\|_{L^2(P\setminus(\eta T))}. \quad (A.5) $$

This allows us to prove the following result, which is the analogue for the Stokes system of what has already been established in [19, Proposition 5.6] for the periodic perforated Poisson equation. The proof needs some adaptations to account for the presence of the pressure.

**Proposition A.1.** Let $\mathbf{f} \in C_{\text{per}}^{\infty}(P, \mathbb{R}^d)$, $\mathbf{g} \in H^{\frac{1}{2}}(\partial(\eta T), \mathbb{R}^d)$, $\mathbf{h} \in C^{\infty}(P\setminus T, \mathbb{R}^d)$ and $q \in C^{\infty}(P\setminus T)$ be functions satisfying the following assumptions:

$$ (i) \quad \int_{\partial P} \mathbf{h} \cdot \mathbf{n} \, d\sigma - \int_{\partial(\eta T)} \mathbf{g} \cdot \mathbf{n} \, d\sigma = 0, $$

$$ (ii) \quad -\Delta \mathbf{h} + \nabla q \in C_{\text{per}}^{\infty}(P\setminus T, \mathbb{R}^d). $$

There exists a unique solution $(\mathbf{u}, p) \in H^1(P\setminus(\eta T), \mathbb{R}^d) \times L^2(P\setminus(\eta T))/\mathbb{R}$ to the boundary value problem

$$ \left\{ \begin{array}{l} -\Delta \mathbf{u} + \nabla p = \mathbf{f} \text{ in } P\setminus(\eta T), \\ \operatorname{div}(\mathbf{u}) = 0 \text{ in } P\setminus(\eta T), \\ \mathbf{u} = \mathbf{g} \text{ on } \partial(\eta T), \\ \mathbf{u} - \mathbf{h} \text{ is P-periodic,} \\ \sigma(\mathbf{u}, p)\mathbf{n}\mathbf{n}^T - \sigma(\mathbf{h}, q)\mathbf{n}\mathbf{n}^T \text{ is P-periodic.} \end{array} \right. \quad (A.6) $$

Furthermore, for any $k > \frac{d}{2} - 1$, there exists a constant $C_k$ independent of $\eta$, $f$ and $h$ such that if $d \ge 3$,

$$ \begin{aligned} \|\nabla \mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + \|p\|_{L^2(P\setminus(\eta T))/\mathbb{R}} & \le C_k \left[ |\alpha(\mathbf{f}, \mathbf{h})| \eta^{1-d/2} + \|\mathbf{g}(\eta)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} \eta^{d/2-1} \right. \\ & \qquad \left. + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \right], \quad (A.7) \end{aligned} $$

and

$$ \begin{aligned} \|\mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} & \le C_k \left[ |\alpha(\mathbf{f}, \mathbf{h})| \eta^{2-d} + \|\mathbf{g}(\eta)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} \right. \\ & \qquad \left. + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \right], \quad (A.8) \end{aligned} $$

where $\alpha(\mathbf{f}, \mathbf{h})$ is the “default of compatibility” defined by

$$ \alpha(\mathbf{f}, \mathbf{h}) := \int_P \mathbf{f} \, dx + \int_{\partial P} \sigma(\mathbf{h}, q) \mathbf{n} \, d\sigma. \quad (A.9) $$

If $d=2$, the following bounds hold:

$$ \begin{aligned} \|\nabla \mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + \|p\|_{L^2(P\setminus(\eta T))/\mathbb{R}} & \le C_k \left[ |\alpha(\mathbf{f}, \mathbf{h})| |\log \eta|^{\frac{1}{2}} + \|\mathbf{g}(\eta)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} \right. \\ & \qquad \left. + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \right], \quad (A.10) \end{aligned} $$

and

$$ \begin{aligned} \|\mathbf{u}\|_{L^2(P\setminus(\eta T), \mathbb{R}^d)} & \le C_k \left[ |\alpha(\mathbf{f}, \mathbf{h})| |\log \eta| + \|\mathbf{g}(\eta)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} |\log \eta|^{\frac{1}{2}} \right. \\ & \qquad \left. + \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P\setminus T)} \right]. \quad (A.11) \end{aligned} $$

**Proof.** Existence and uniqueness of a solution $(\mathbf{u}, p)$ to (A.6) can be obtained by adapting the proof of Proposition 2.7, where no compatibility condition is required because of the Dirichlet boundary condition on $\partial(\eta T)$.
function $g(\eta \cdot)$ compactly supported in $K$, satisfying $\tilde{g}_\eta = g(\eta \cdot)$ on $\partial T$ and $||\tilde{g}_\eta||_{H^1(K, \mathbb{R}^d)} \le C||g(\eta \cdot)||_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)}$ for a constant $C$ independent of $g$ and $\eta$. We define $\hat{g} := \tilde{g}_\eta(\cdot/\eta)$, which yields an extension of $g$ compactly supported in $P$ (extending by 0 in $P\setminus(\eta K)$) satisfying $\hat{g} = g$ on $\partial(\eta T)$ and

$$ \begin{aligned} ||\hat{g}||_{H^1(P, \mathbb{R}^d)} &= ||\tilde{g}_\eta(\cdot/\eta)||_{H^1(\eta K)} = \eta^{d/2-1}||\nabla\tilde{g}_\eta||_{L^2(K, \mathbb{R}^{d \times d})} + \eta^{d/2}||\tilde{g}_\eta||_{L^2(K, \mathbb{R}^d)} \\ &\le C\eta^{d/2-1}||\tilde{g}_\eta||_{H^1(K, \mathbb{R}^d)} \le C\eta^{d/2-1}||g(\eta \cdot)||_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)}. \end{aligned} $$

Similarly, by setting $\hat{u}_1(x) = u_1(x)\theta(x/\eta)$ for a cut-off function $\theta \in C_c^\infty(K)$ satisfying $\theta = 1$ on $\partial T$, we construct an extension $\hat{u}_1 \in H^1(P)$ compactly supported in $P$ such that $\hat{u}_1 = u_1$ on $\partial(\eta T)$ and

$$ ||\hat{u}_1||_{H^1(P, \mathbb{R}^d)} = \eta^{d/2}||\theta u_1(\eta \cdot)||_{L^2(K, \mathbb{R}^d)} + \eta^{d/2-1}||\nabla(\theta u_1(\eta \cdot))||_{L^2(K, \mathbb{R}^{d \times d})} \le C\eta\eta^{d/2-1}||\nabla u_1||_{L^\infty(P, \mathbb{R}^{d \times d})}, $$

where we have used $u_1(0) = 0$ and $||u_1(\eta \cdot)||_{L^\infty(K, \mathbb{R}^d)} \le C\eta||\nabla u_1||_{L^\infty(P, \mathbb{R}^{d \times d})}$ in the last estimate.

We then consider the function $v_2 := u_2 - (\hat{g} - \hat{u}_1) + B_\eta(\text{div}(\hat{g} - \hat{u}_1))$ where $B_\eta$ is the Bogovskii's operator provided by Lemma A.2 below. Note that $\text{div}(\hat{g} - \hat{u}_1)$ is in the domain of $B_\eta$ thanks to the assumption (i). Since by construction $v_2 = 0$ on $\partial(\eta T)$, $v_2$ is P-periodic and $\text{div}(v_2) = 0$, an application of the divergence theorem yields, for $d \ge 3$,

$$ \begin{aligned} \int_{P\setminus(\eta T)} |\nabla v_2|^2 dx &= \int_{P\setminus(\eta T)} (\nabla u_2 : \nabla v_2 - p_2 \text{div}(v_2)) dx - \int_{P\setminus(\eta T)} \nabla v_2 \cdot (\nabla \hat{g} + \nabla \hat{u}_1 - B_\eta(\text{div}(\hat{g} - \hat{u}_1))) dx \\ &= \int_{P\setminus(\eta T)} \alpha(f, h) \cdot v_2 dx - \int_{P\setminus(\eta T)} \nabla v_2 \cdot (\nabla \hat{g} + \nabla \hat{u}_1 - B_\eta(\text{div}(\hat{g} - \hat{u}_1))) dx \\ &\le C(|\alpha(f, h)|\eta^{1-d/2} + ||\nabla \hat{g}||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + ||\nabla \hat{u}_1||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})})||\nabla v_2||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \\ &\le C(|\alpha(f, h)|\eta^{1-d/2} + \eta^{d/2-1}||g(\eta \cdot)||_{H^{\frac{1}{2}}(\partial T)} + \eta^{d/2}||\nabla u_1||_{L^\infty(P, \mathbb{R}^{d \times d})})||\nabla v_2||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})}. \end{aligned} \quad (A.15) $$

Using (A.13), we obtain thus if $d \ge 3$,

$$ \begin{aligned} ||\nabla u_2||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} &\le ||\nabla v_2||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + C\eta^{d/2-1}||g(\eta \cdot)||_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} + \eta^{d/2}||\nabla u_1||_{L^\infty(P, \mathbb{R}^{d \times d})} \\ &\le C_k(|\alpha(f, h)|\eta^{1-d/2} + \eta^{d/2-1}||g(\eta \cdot)||_{H^{\frac{1}{2}}(\partial T)} \\ &\qquad + \eta^{d/2}||f - \alpha(f, h)||_{H^k(P, \mathbb{R}^d)} + \eta^{d/2}||h||_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \eta^{d/2}||q||_{H^{k+1}(P\setminus T)}). \end{aligned} \quad (A.16) $$

We then deduce from the Poincaré inequality (A.3)

$$ \begin{aligned} ||u_2||_{L^2(P\setminus(\eta T), \mathbb{R}^d)} \le C_k(|\alpha(f, h)\eta^{2-d} + ||g(\eta \cdot)||_{H^{\frac{1}{2}}(\partial T)} \\ + \eta||f - \alpha(f, h)||_{H^k(P, \mathbb{R}^d)} + \eta||h||_{H^{k+2}(P\setminus T, \mathbb{R}^d)} + \eta||q||_{H^{k+1}(P\setminus T)}). \end{aligned} \quad (A.17) $$

It remains to estimate $||p_2||_{L^2(P\setminus(\eta T))/\mathbb{R}}$. Let the constant of $p_2$ be chosen such that $\int_{P\setminus(\eta T)} p_2 dx = 0$. Then, multiplying (A.14) and integrating by parts yields

$$ \int_{P\setminus(\eta T)} |p_2|^2 dx = \int_{P\setminus(\eta T)} p_2 \text{div}(B_\eta(p_2)) dx = \int_{P\setminus(\eta T)} \nabla u_2 : \nabla (B_\eta p_2) dx - \int_{P\setminus(\eta T)} \alpha(f, h) \cdot B_\eta(p_2) dx. $$

Using the Cauchy-Schwartz inequality and Lemma A.2, we thus find

$$ ||p_2||_{L^2(P\setminus(\eta T))/\mathbb{R}} \le ||\nabla u_2||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + ||\alpha(f, h)||_{L^2(P\setminus(\eta T))}. \quad (A.18) $$

The final result for $d \ge 3$ follows by summing the estimate (A.13) with (A.16) to (A.18). In the case $d = 2$, the estimates (A.10) and (A.11) are similarly obtained by replacing the quantities $\eta^{1-d/2}$ with $|\log \eta|^{1/2}$ coming from the Poincaré inequality (A.3) in (A.15) to (A.17). $\square$

In the proof of the previous proposition, we have used the following lemma stated in [2, Lemma 2.2.4].

**Lemma A.2.** For any $\eta > 0$, there exists a linear “Bogovskii’s” operator $B_\eta : L^2(P\setminus(\eta T)) \to H^1(P\setminus(\eta T), \mathbb{R}^d)$ satisfying for any $\phi \in L^2(P\setminus(\eta Y))$ such that $\int_{P\setminus(\eta T)} \phi dy = 0$:

(i) $\text{div}(B_\eta \phi) = \phi$,

(ii) $B_\eta \phi = 0$ on $\partial(\eta T)$,

(iii) $B_\eta \phi$ is P-periodic,

(iv) $||\nabla(B_\eta \phi)||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \le C||\phi||_{L^2(P\setminus(\eta T))}$ for a constant $C$ independent of $\eta$ and $\phi$.

We can now give an alternate proof of Proposition 3.5 that uses only the a priori energy estimates derived in the previous proposition.
We then decompose $\mathbf{u} = \mathbf{u}_1 + \mathbf{u}_2$ and $p = p_1 + p_2$ as the sums of the solutions to two distinct problems: $(\mathbf{u}_1, p_1) \in H^1(P, \mathbb{R}^d) \times L^2(P)/\mathbb{R}$ is defined as the unique solution to

$$ egin{cases} 
-\Delta \mathbf{u}_1 + \nabla p_1 = \mathbf{f} - \alpha(\mathbf{f}, \mathbf{h}) \text{ in } P, \\
\operatorname{div}(\mathbf{u}_1) = 0 \text{ in } P, \\
\mathbf{u}_1 - \mathbf{h} \text{ is P-periodic}, \\
\sigma(\mathbf{u}_1, p_1) \mathbf{n}\mathbf{n}^T - \sigma(\mathbf{h}, q) \mathbf{n}\mathbf{n}^T \text{ is P-periodic}, \\
\mathbf{u}_1(0) = 0.
\end{cases} 
\quad (A.12) $$

The problem (A.12) is of the type of (2.42), whose compatibility condition is satisfied by the definition (A.9) of $\alpha(\mathbf{f}, \mathbf{h})$. Furthermore, the condition $\mathbf{u}_1(0) = 0$ fixes the additive constant determining $\mathbf{u}_1$. By Propositions 2.7 and 2.8, $(\mathbf{u}_2, p_2)$ satisfies the following bound:

$$ \|\mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^d)} + \|\nabla \mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^{d \times d})} + \|p_1\|_{L^\infty(P)/\mathbb{R}} \\ \le C_k (\|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \|\mathbf{h}\|_{H^{k+2}(P \setminus T, \mathbb{R}^d)} + \|q\|_{H^{k+1}(P \setminus T)}). \quad (A.13) $$

Then $(\mathbf{u}_2, p_2)$ is defined as the solution to the difference problem

$$ egin{cases} 
-\Delta \mathbf{u}_2 + \nabla p_2 = \alpha(\mathbf{f}, \mathbf{h}) \text{ in } P, \\
\operatorname{div}(\mathbf{u}_2) = 0 \text{ in } P, \\
\mathbf{u}_2 = \mathbf{g} - \mathbf{u}_1 \text{ on } \partial(\eta T), \\
\mathbf{u}_2 \text{ is P-periodic}, \\
\sigma(\mathbf{u}_2, p_2) \mathbf{n}\mathbf{n}^T \text{ is P-periodic.}
\end{cases} 
\quad (A.14) $$

In order to estimate the energy norm of the solution to (A.14), we construct lifts of the boundary conditions $\hat{\mathbf{u}}_1$ and $\hat{\mathbf{g}}$. Let $K$ be an open set containing the obstacle $T: T \subset K$. We consider an extension $\tilde{\mathbf{g}}_\eta$ in $K$ of the function $\mathbf{g}(\eta\cdot)$ compactly supported in $K$, satisfying $\tilde{\mathbf{g}}_\eta = \mathbf{g}(\eta\cdot)$ on $\partial T$ and $\|\tilde{\mathbf{g}}_\eta\|_{H^1(K, \mathbb{R}^d)} \le C \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)}$ for a constant $C$ independent of $g$ and $\eta$. We define $\hat{\mathbf{g}} := \tilde{\mathbf{g}}_\eta(\cdot/\eta)$, which yields an extension of $\mathbf{g}$ compactly supported in $P$ (extending by 0 in $P\setminus(\eta K)$) satisfying $\hat{\mathbf{g}} = \mathbf{g}$ on $\partial(\eta T)$ and

$$ \begin{aligned} \|\hat{\mathbf{g}}\|_{H^1(P, \mathbb{R}^d)} &= \|\tilde{\mathbf{g}}_\eta(\cdot/\eta)\|_{H^1(\eta K)} = \eta^{d/2-1} \|\nabla\tilde{\mathbf{g}}_\eta\|_{L^2(K, \mathbb{R}^{d \times d})} + \eta^{d/2} \|\tilde{\mathbf{g}}_\eta\|_{L^2(K, \mathbb{R}^d)} \\ &\le C \eta^{d/2-1} \|\tilde{\mathbf{g}}_\eta\|_{H^1(K, \mathbb{R}^d)} \le C \eta^{d/2-1} \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)}. \end{aligned} $$

Similarly, by setting $\hat{\mathbf{u}}_1(x) = \mathbf{u}_1(x)\theta(x/\eta)$ for a cut-off function $\theta \in C_c^\infty(K)$ satisfying $\theta = 1$ on $\partial T$, we construct an extension $\hat{\mathbf{u}}_1 \in H^1(P)$ compactly supported in $P$ such that $\hat{\mathbf{u}}_1 = \mathbf{u}_1$ on $\partial(\eta T)$ and

$$ \|\hat{\mathbf{u}}_1\|_{H^1(P, \mathbb{R}^d)} = \eta^{d/2} \|\theta \mathbf{u}_1(\eta\cdot)\|_{L^2(K, \mathbb{R}^d)} + \eta^{d/2-1} \|\nabla(\theta \mathbf{u}_1(\eta\cdot))\|_{L^2(K, \mathbb{R}^{d \times d})} \le C \eta^{d/2-1} \|\nabla \mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^{d \times d})}, $$

where we have used $\mathbf{u}_1(0) = 0$ and $\|\mathbf{u}_1(\eta\cdot)\|_{L^\infty(K, \mathbb{R}^d)} \le C\eta \|\nabla \mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^{d \times d})}$ in the last estimate.

We then consider the function $\mathbf{v}_2 := \mathbf{u}_2 - (\hat{\mathbf{g}} - \hat{\mathbf{u}}_1) + B_\eta(\operatorname{div}(\hat{\mathbf{g}} - \hat{\mathbf{u}}_1))$ where $B_\eta$ is the Bogovskii's operator provided by Lemma A.2 below. Note that $\operatorname{div}(\hat{\mathbf{g}} - \hat{\mathbf{u}}_1)$ is in the domain of $B_\eta$ thanks to the assumption (i). Since by construction $\mathbf{v}_2 = 0$ on $\partial(\eta T)$, $\mathbf{v}_2$ is P-periodic and $\operatorname{div}(\mathbf{v}_2) = 0$, an application of the divergence theorem yields, for $d \ge 3$,

$$ \begin{aligned} \int_{P \setminus (\eta T)} |\nabla \mathbf{v}_2|^2 dx &= \int_{P \setminus (\eta T)} (\nabla \mathbf{u}_2 : \nabla \mathbf{v}_2 - p_2 \operatorname{div}(\mathbf{v}_2)) dx - \int_{P \setminus (\eta T)} \nabla \mathbf{v}_2 \cdot (\nabla \hat{\mathbf{g}} + \nabla \hat{\mathbf{u}}_1 - B_\eta(\operatorname{div}(\hat{\mathbf{g}} - \hat{\mathbf{u}}_1))) dx \\ &= \int_{P \setminus (\eta T)} \alpha(\mathbf{f}, \mathbf{h}) \cdot \mathbf{v}_2 dx - \int_{P \setminus (\eta T)} \nabla \mathbf{v}_2 \cdot (\nabla \hat{\mathbf{g}} + \nabla \hat{\mathbf{u}}_1 - B_\eta(\operatorname{div}(\hat{\mathbf{g}} - \hat{\mathbf{u}}_1))) dx \\ &\le C (|\alpha(\mathbf{f}, \mathbf{h})|\eta^{1-d/2} + \|\nabla \hat{\mathbf{g}}\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})} + \|\nabla \hat{\mathbf{u}}_1\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})}) \|\nabla \mathbf{v}_2\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})} \\ &\le C (|\alpha(\mathbf{f}, \mathbf{h})|\eta^{1-d/2} + \eta^{d/2-1} \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T)} + \eta^{d/2} \|\nabla \mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^{d \times d})}) \|\nabla \mathbf{v}_2\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})}). \end{aligned} \quad (A.15) $$

Using (A.13), we obtain thus if $d \ge 3$,

$$ \begin{aligned} \|\nabla \mathbf{u}_2\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})} &\le \|\nabla \mathbf{v}_2\|_{L^2(P \setminus (\eta T), \mathbb{R}^{d \times d})} + C \eta^{d/2-1} \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T, \mathbb{R}^d)} + \eta^{d/2} \|\nabla \mathbf{u}_1\|_{L^\infty(P, \mathbb{R}^{d \times d})} \\ &\le C_k (|\alpha(\mathbf{f}, \mathbf{h})|\eta^{1-d/2} + \eta^{d/2-1} \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T)} \\ &\quad + \eta^{d/2} \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \eta^{d/2} \|\mathbf{h}\|_{H^{k+2}(P \setminus T, \mathbb{R}^d)} + \eta^{d/2} \|q\|_{H^{k+1}(P \setminus T)}). \end{aligned} \quad (A.16) $$

We then deduce from the Poincaré inequality (A.3)

$$ \|\mathbf{u}_2\|_{L^2(P \setminus (\eta T), \mathbb{R}^d)} \le C_k (|\alpha(\mathbf{f}, \mathbf{h})\eta^{2-d} + \|\mathbf{g}(\eta\cdot)\|_{H^{\frac{1}{2}}(\partial T)} \\ + \eta \|\mathbf{f} - \alpha(\mathbf{f}, \mathbf{h})\|_{H^k(P, \mathbb{R}^d)} + \eta \|\mathbf{h}\|_{H^{k+2}(P \setminus T, \mathbb{R}^d)} + \eta \|q\|_{H^{k+1}(P \setminus T)}). \quad (A.17) $$
Proof of Proposition 3.5. Let us denote by $(\mathbf{r}_{i,\eta}, s_{i,\eta})$ the residual functions $\mathbf{r}_{i,\eta} := \mathbf{X}_i - \mathbf{X}_i^N$ and $s_{i,\eta} := \alpha_i - \alpha_i^N$. $(\mathbf{r}_{i,\eta}, s_{i,\eta})$ satisfies

$$ \left\{ \begin{array}{l} -\Delta \mathbf{r}_{i,\eta} + s_{i,\eta} = 0 \text{ in } P\setminus(\eta T), \\ \mathrm{div}(\mathbf{r}_{i,\eta}) = 0 \text{ in } P\setminus(\eta T). \end{array} \right. \qquad (A.19) $$

For $y \in \partial T$, we have for any $\eta > 0$:

$$ \begin{aligned} \mathbf{r}_{i,\eta}(\eta y) &= -\mathbf{X}_i^N(\eta y) = -\sum_{p=1}^N \eta^p \mathbf{v}_{i,p}(y) - \sum_{p=0}^N \eta^p \mathbf{w}_{i,p}(\eta y) \\ &= -\sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(y) - \sum_{p=0}^N \eta^p \sum_{k=0}^p \frac{1}{k!} \nabla^k \mathbf{w}_{i,p-k}(0) \cdot y^k + O(\eta^{N+1}) \\ &= -\mathbf{w}_{i,0}(0) - \eta(\mathbf{w}_{i,1}(0) + \nabla \mathbf{w}_{i,0}(0) \cdot y) + O(\eta^{N+1}) = O(\eta^{N+1}), \end{aligned} $$

where we have used the boundary conditions of $\mathbf{v}_{i,p}$ in (3.20). For $x \in \partial P$, we find

$$ \begin{aligned} \mathbf{r}_{i,\eta}(x) &= \mathbf{X}_i(x) - \eta^{2-d} \mathbf{\Phi}_i(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{w}_{i,p}(x) \\ &= -\eta^{2-d} \mathbf{\Phi}_i(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p G(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i \\ &\quad + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} + \mathbf{H}_i^N(x), \end{aligned} $$

$$ \begin{aligned} s_{i,\eta}(x) &= \alpha_i(x) - \eta^{1-d} \sigma_i(x/\eta) - \sum_{p=2}^N \eta^{p-1} \xi_{i,p}(x/\eta) - \sum_{p=2}^N \eta^p \zeta_{i,p}(x) \\ &= -\eta^{1-d} \sigma_i(x/\eta) - \sum_{p=2}^N \eta^{p-1} \xi_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p \psi(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i \\ &\quad + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} \psi(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} + \beta_i^N(x), \end{aligned} $$

where $\mathbf{H}_i^N$ and $\sigma(\mathbf{H}_i^N, \beta_i^N)nn^T$ are $P$-periodic functions due to the boundary conditions of (3.19). Denote by $\mathbf{h}_{i,\eta}$ and $q_{i,\eta}$ the functions defined on $P\setminus\bar{T}$ (a neighborhood of $\partial P$) by

$$ \mathbf{h}_{i,\eta}(x) := -\eta^{2-d} \mathbf{\Phi}_i(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p G(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k}. $$

$$ q_{i,\eta}(x) := -\eta^{1-d} \sigma_i(x/\eta) - \sum_{p=2}^N \eta^{p-1} \xi_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p \psi(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} \psi(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k}. $$

By using the far field expansions of $\mathbf{\Phi}_i$, $\mathbf{v}_{i,p}$, $\xi_{i,p}$ and $\sigma_i$, it can be verified that

$$ \|\mathbf{h}_{i,\eta}\|_{H^{k+2}(P\setminus\bar{T}, \mathbb{R}^d)} + \|q_{i,\eta}\|_{H^{k+1}(P\setminus\bar{T})} = O(\eta^{N+1}), $$

for any $k > d/2 - 1$. Moreover, $-\Delta \mathbf{h}_{i,\eta} + \nabla q_{i,\eta} = 0$ and $\mathrm{div}(\mathbf{h}_{i,\eta}) = 0$. Hence, we aim at using the result of Proposition A.1 to estimate the functions $(\mathbf{r}_{i,\eta}, s_{i,\eta})$ which solve (A.19) and satisfy $\|\mathbf{r}_{i,\eta}(\eta\cdot)\|_{H^{1/2}(\partial T, \mathbb{R}^d)} = O(\eta^{N+1})$ as well as the boundary conditions

$$ \left\{ \begin{array}{l} \mathbf{r}_{i,\eta} - \mathbf{h}_{i,\eta} \text{ is } P\text{-periodic,} \\ \sigma(\mathbf{r}_{i,\eta}, s_{i,\eta})nn^T - \sigma(\mathbf{h}_{i,\eta}, q_{i,\eta})nn^T \text{ is } P\text{-periodic.} \end{array} \right. $$
It remains to estimate $||p_2||_{L^2(P\setminus(\eta T))/\mathbb{R}}$. Let the constant of $p_2$ be chosen such that $\int_{P\setminus(\eta T)} p_2 dx = 0$. Then, multiplying (A.14) and integrating by parts yields

$$ \int_{P\setminus(\eta T)} |p_2|^2 dx = \int_{P\setminus(\eta T)} p_2 \operatorname{div}(B_\eta(p_2)) dx = \int_{P\setminus(\eta T)} \nabla u_2 : \nabla(B_\eta p_2) dx - \int_{P\setminus(\eta T)} \alpha(f, h) \cdot B_\eta(p_2) dx. $$

Using the Cauchy-Schwartz inequality and Lemma A.2, we thus find

$$ ||p_2||_{L^2(P\setminus(\eta T))/\mathbb{R}} \le ||\nabla u_2||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} + ||\alpha(f, h)||_{L^2(P\setminus(\eta T))}. \quad (A.18) $$

The final result for $d \ge 3$ follows by summing the estimate (A.13) with (A.16) to (A.18). In the case $d = 2$, the estimates (A.10) and (A.11) are similarly obtained by replacing the quantities $\eta^{1-d/2}$ with $|\log \eta|^{1/2}$ coming from the Poincaré inequality (A.3) in (A.15) to (A.17). $\square$

In the proof of the previous proposition, we have used the following lemma stated in [2, Lemma 2.2.4].

**Lemma A.2.** For any $\eta > 0$, there exists a linear “Bogovskii’s” operator $B_\eta : L^2(P\setminus(\eta T)) \to H^1(P\setminus(\eta T), \mathbb{R}^d)$ satisfying for any $\phi \in L^2(P\setminus(\eta Y))$ such that $\int_{P\setminus(\eta T)} \phi dy = 0$:

(i) $\operatorname{div}(B_\eta \phi) = \phi$,

(ii) $B_\eta \phi = 0$ on $\partial(\eta T)$,

(iii) $B_\eta \phi$ is P-periodic,

(iv) $||\nabla(B_\eta \phi)||_{L^2(P\setminus(\eta T), \mathbb{R}^{d \times d})} \le C||\phi||_{L^2(P\setminus(\eta T))}$ for a constant $C$ independent of $\eta$ and $\phi$.

We can now give an alternate proof of Proposition 3.5 that uses only the a priori energy estimates derived in the previous proposition.

**Proof of Proposition 3.5.** Let us denote by $(\mathbf{r}_{i,\eta}, s_{i,\eta})$ the residual functions $\mathbf{r}_{i,\eta} := \mathbf{X}_i - \mathbf{X}_i^N$ and $s_{i,\eta} := \alpha_i - \alpha_i^N$. $(\mathbf{r}_{i,\eta}, s_{i,\eta})$ satisfies

$$ \left\{ \begin{array}{l} -\Delta \mathbf{r}_{i,\eta} + s_{i,\eta} = 0 \quad \text{in } P\setminus(\eta T), \\ \operatorname{div}(\mathbf{r}_{i,\eta}) = 0 \quad \text{in } P\setminus(\eta T). \end{array} \right. \quad (A.19) $$

For $y \in \partial T$, we have for any $\eta > 0$:

$$ \begin{aligned} \mathbf{r}_{i,\eta}(\eta y) &= -\mathbf{X}_i^N(\eta y) = -\sum_{p=1}^N \eta^p \mathbf{v}_{i,p}(y) - \sum_{p=0}^N \eta^p \mathbf{w}_{i,p}(\eta y) \\ &= -\sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(y) - \sum_{p=0}^N \eta^p \sum_{k=0}^p \frac{1}{k!} \nabla^k \mathbf{w}_{i,p-k}(0) \cdot y^k + O(\eta^{N+1}) \\ &= -\mathbf{w}_{i,0}(0) - \eta(\mathbf{w}_{i,1}(0) + \nabla \mathbf{w}_{i,0}(0) \cdot y) + O(\eta^{N+1}) = O(\eta^{N+1}), \end{aligned} $$

where we have used the boundary conditions of $\mathbf{v}_{i,p}$ in (3.20). For $x \in \partial P$, we find

$$ \begin{aligned} \mathbf{r}_{i,\eta}(x) &= \mathbf{X}_i(x) - \eta^{2-d} \mathbf{\Phi}_i(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{w}_{i,p}(x) \\ &= -\eta^{2-d} \mathbf{\Phi}_i(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p G(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i \\ &\quad + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} + \mathbf{H}_i^N(x), \\ s_{i,\eta}(x) &= \alpha_i(x) - \eta^{1-d} \sigma_i(x/\eta) - \sum_{p=2}^N \eta^{p-1} \xi_{i,p}(x/\eta) - \sum_{p=2}^N \eta^p \zeta_{i,p}(x) \\ &= -\eta^{1-d} \sigma_i(x/\eta) - \sum_{p=2}^N \eta^{p-1} \xi_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p \psi(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i \\ &\quad + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} \psi(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k} + \beta_i^N(x), \end{aligned} $$

where $\mathbf{H}_i^N$ and $\sigma(\mathbf{H}_i^N, \beta_i^N)nn^T$ are P-periodic functions due to the boundary conditions of (3.19). Denote by $\mathbf{h}_{i,\eta}$ and $q_{i,\eta}$ the functions defined on $P\setminus\bar{T}$ (a neighborhood of $\partial P$) by

$$ \mathbf{h}_{i,\eta}(x) := -\eta^{2-d} \mathbf{\Phi}_i(x/\eta) - \sum_{p=2}^N \eta^p \mathbf{v}_{i,p}(x/\eta) + \sum_{p=0}^N \eta^p \nabla^p G(x) \cdot \mathcal{F}^p \mathbf{\Phi}_i + \sum_{p=d}^N \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} G(x) \cdot \mathcal{F}^{k-d+2} \mathbf{v}_{i,p-k}. $$
In order to apply Proposition A.1, we need to verify its assumption (i). By using $\int_{\partial P} \nabla^p G_i(x) \cdot n \, d\sigma(x) = 0$ for all $p \ge 0$ and $1 \le i \le d$, we find

$$ 
\begin{aligned}
& \int_{\partial P} h_{i,\eta} \cdot n \, d\sigma - \int_{\partial(\eta T)} r_{i,\eta} \cdot n \, d\sigma \\
&= -\eta^{2-d} \left( \int_{\partial P} \Phi_i(x/\eta) \cdot n \, d\sigma(x) - \int_{\partial(\eta T)} \Phi_i(x/\eta) \cdot n \, d\sigma(x) \right) \\
&\quad - \sum_{p=2}^{N} \eta^p \left( \int_{\partial P} v_{i,p}(x/\eta) \cdot n \, d\sigma(x) - \int_{\partial(\eta T)} \eta^p v_{i,p}(x/\eta) \, d\sigma(x) \right) + \sum_{p=0}^{N} \eta^p \int_{\partial(\eta T)} w_{i,p}(x) \cdot n \, d\sigma(x) \\
&= -\eta^{2-d} \int_{P \setminus (\eta T)} \operatorname{div}(\Phi_i(x/\eta)) \, dx - \sum_{p=2}^{N} \eta^p \int_{P \setminus (\eta T)} \operatorname{div}(v_{i,p}(x/\eta)) \, dx + \sum_{p=0}^{N} \eta^p \int_{\eta T} \operatorname{div}(w_{i,p}) \, dx = 0,
\end{aligned}
$$

hence the assumption is satisfied. It remains to evaluate the default of compatibility which reads

$$ 
\begin{aligned}
\alpha(0, h_{i,\eta}) &= \int_{\partial P} \sigma(h_{i,\eta}, q_{i,\eta}) n \, d\sigma \\
&= -\eta^{2-d} \int_{\partial P} \sigma(\Phi_i(x/\eta), \eta^{-1}\sigma_i(x/\eta)) n \, d\sigma - \sum_{p=2}^{N} \eta^p \int_{\partial P} \sigma(v_{i,p}(x/\eta), \eta^{-1}\xi_{i,p}(x/\eta)) n \, d\sigma + \mathcal{F}^0 \Phi_i, \quad (A.20)
\end{aligned}
$$

where we used Lemma 2.1. Using now (2.35) and the fact that $(\Phi_i, \sigma_i)$ and $(v_{i,p}, \xi_{i,p})$ are solutions to Stokes problems, the application of the divergence theorem yields

$$ 
\begin{aligned}
\alpha(0, h_{i,\eta}) &= -\eta^{2-d} \int_{\partial(\eta T)} \sigma(\Phi_i(x/\eta), \eta^{-1}\sigma_i(x/\eta)) n \, d\sigma - \sum_{p=2}^{N} \eta^p \int_{\partial(\eta T)} \sigma(v_{i,p}(x/\eta), \eta^{-1}\xi_{i,p}(x/\eta)) n \, d\sigma + e_i \\
&= -\int_{\partial T} \sigma(\Phi_i, \sigma_i) n \, d\sigma - \sum_{p=2}^{N} \eta^{p+d-2} \int_{\partial T} \sigma(v_{i,p}, \xi_{i,p}) n \, d\sigma + e_i.
\end{aligned}
$$

This term is precisely zero because $-\int_{\partial T} \sigma(\Phi_i, \sigma_i) n \, d\sigma = -\mathcal{F}^0 \Phi_i = e_i$ and $\int_{\partial T} \sigma(v_{i,p}, \xi_{i,p}) n \, d\sigma = \mathcal{F}^0 v_{i,p} = 0$ due to (3.20). The final estimates (3.22) and (3.23) follow by applying Proposition A.1, however with a resulting constant $C_N$ which is a priori not uniform in $N$.

## APPENDIX B. SOME BOX INTEGRALS IDENTITIES

The proof of the point (iii) of Lemma 2.2 relies on the following result.

**Proposition B.1.** Let $P = (-1/2, 1/2)^d$ the $d$-dimensional unit hypercube with $d \ge 1$. The following identity holds for any $p, k \ge 0$ and $1 \le i \le d$:

$$ \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{p}{2k+1} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{2k+d-p}{(2k+1)2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x), \quad (B.1) $$

and in particular for $p = 2k + d$:

$$ \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{2k+d+2}} d\sigma(x) = \frac{2k+1}{2k+d} \int_{\partial P} \frac{x_i^{2k}}{|x|^{2k+d}} d\sigma(x). \quad (B.2) $$

*Proof.* Let us denote by $F_i = \{(x_1, \dots, x_d) | x_i = \frac{1}{2}$ and $-\frac{1}{2} \le x_j \le \frac{1}{2}$ for $1 \le j \ne i \le d\}$ the face of the cube with normal $e_i$. By decomposing the integrals over the face $F_i$ and the other faces $F_j$, $j \ne i$, we find

$$ \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = 2 \times \frac{1}{2^{2k}} \int_{F_i} \frac{1}{|x|^p} d\sigma(x) + 2(d-1) \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x) $$

where $1 \le j \le d$ is any index different from $i$. Since $\int_{\partial P} \frac{1}{|x|^d} d\sigma(x) = 2d \int_{F_i} \frac{1}{|x|^p} d\sigma(x)$, this identity can be rewritten as

$$ \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{1}{2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x) + 2(d-1) \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x), \quad (B.3) $$

which also implies

$$ \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{1}{2(d-1)} \left( \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) - \frac{1}{2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x) \right) \quad (B.4) $$

for any $p, k \ge 0$. We then use the (tangential) divergence theorem on the face $F_j$ to obtain

$$ \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{1}{2k+1} \int_{F_j} \frac{\partial_i(x_i^{2k+1})}{|x|^p} d\sigma(x) = \frac{p}{2k+1} \int_{F_j} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{1}{2k+1} \int_{\partial F_j} \frac{x_i^{2k+1} n_i}{|x|^p} d\sigma(x), $$
$$ q_{i,\eta}(x) := -\eta^{1-d}\sigma_i(x/\eta) - \sum_{p=2}^{N} \eta^{p-1}\xi_{i,p}(x/\eta) + \sum_{p=0}^{N} \eta^p \nabla^p \psi(x) \cdot \mathcal{F}^p \Phi_i + \sum_{p=d}^{N} \eta^p \sum_{k=d-1}^{p-2} \nabla^{k-d+2} \psi(x) \cdot \mathcal{F}^{k-d+2} v_{i,p-k} $$

By using the far field expansions of $\Phi_i$, $v_{i,p}$, $\xi_{i,p}$ and $\sigma_i$, it can be verified that

$$ \|h_{i,\eta}\|_{H^{k+2}(P\setminus\bar{T},\mathbb{R}^d)} + \|q_{i,\eta}\|_{H^{k+1}(P\setminus\bar{T})} = O(\eta^{N+1}), $$

for any $k > d/2 - 1$. Moreover, $-\Delta h_{i,\eta} + \nabla q_{i,\eta} = 0$ and $\text{div}(h_{i,\eta}) = 0$. Hence, we aim at using the result of Proposition A.1 to estimate the functions $(r_{i,\eta}, s_{i,\eta})$ which solve (A.19) and satisfy $\|r_{i,\eta}(\eta)\|_{H^{1/2}(\partial T, \mathbb{R}^d)} = O(\eta^{N+1})$ as well as the boundary conditions

$$ \begin{cases} r_{i,\eta} - h_{i,\eta} \text{ is P-periodic,} \\ \sigma(r_{i,\eta}, s_{i,\eta})nn^T - \sigma(h_{i,\eta}, q_{i,\eta})nn^T \text{ is P-periodic.} \end{cases} $$

In order to apply Proposition A.1, we need to verify its assumption (i). By using $\int_{\partial P} \nabla^p G_i(x) \cdot n \, d\sigma(x) = 0$ for all $p \ge 0$ and $1 \le i \le d$, we find

$$ \begin{aligned} & \int_{\partial P} h_{i,\eta} \cdot n \, d\sigma - \int_{\partial(\eta T)} r_{i,\eta} \cdot n \, d\sigma \\ &= -\eta^{2-d} \left( \int_{\partial P} \Phi_i(x/\eta) \cdot n \, d\sigma(x) - \int_{\partial(\eta T)} \Phi_i(x/\eta) \cdot n \, d\sigma(x) \right) \\ & \quad - \sum_{p=2}^{N} \eta^p \left( \int_{\partial P} v_{i,p}(x/\eta) \cdot n \, d\sigma(x) - \int_{\partial(\eta T)} \eta^p v_{i,p}(x/\eta) \, d\sigma(x) \right) + \sum_{p=0}^{N} \eta^p \int_{\partial(\eta T)} w_{i,p}(x) \cdot n \, d\sigma(x) \\ &= -\eta^{2-d} \int_{P\setminus(\eta T)} \text{div}(\Phi_i(x/\eta)) \, dx - \sum_{p=2}^{N} \eta^p \int_{P\setminus(\eta T)} \text{div}(v_{i,p}(x/\eta)) \, dx + \sum_{p=0}^{N} \int_{\eta T} \text{div}(w_{i,p}) \, dx = 0, \end{aligned} $$

hence the assumption is satisfied. It remains to evaluate the default of compatibility which reads

$$ \begin{aligned} \alpha(0, h_{i,\eta}) &= \int_{\partial P} \sigma(h_{i,\eta}, q_{i,\eta}) n \, d\sigma \\ &= -\eta^{2-d} \int_{\partial P} \sigma(\Phi_i(x/\eta), \eta^{-1}\sigma_i(x/\eta)) n \, d\sigma - \sum_{p=2}^{N} \eta^p \int_{\partial P} \sigma(v_{i,p}(x/\eta), \eta^{-1}\xi_{i,p}(x/\eta)) n \, d\sigma + \mathcal{F}^0 \Phi_i, \quad (A.20) \end{aligned} $$

where we used Lemma 2.1. Using now (2.35) and the fact that $(\Phi_i, \sigma_i)$ and $(v_{i,p}, \xi_{i,p})$ are solutions to Stokes problems, the application of the divergence theorem yields

$$ \begin{aligned} \alpha(0, h_{i,\eta}) &= -\eta^{2-d} \int_{\partial(\eta T)} \sigma(\Phi_i(x/\eta), \eta^{-1}\sigma_i(x/\eta)) n \, d\sigma - \sum_{p=2}^{N} \eta^p \int_{\partial(\eta T)} \sigma(v_{i,p}(x/\eta), \eta^{-1}\xi_{i,p}(x/\eta)) n \, d\sigma + e_i \\ &= -\int_{\partial T} \sigma(\Phi_i, \sigma_i) n \, d\sigma - \sum_{p=2}^{N} \eta^{p+d-2} \int_{\partial T} \sigma(v_{i,p}, \xi_{i,p}) n \, d\sigma + e_i. \end{aligned} $$

This term is precisely zero because $-\int_{\partial T} \sigma(\Phi_i, \sigma_i) n \, d\sigma = -\mathcal{F}^0 \Phi_i = e_i$ and $\int_{\partial T} \sigma(v_{i,p}, \xi_{i,p}) n \, d\sigma = \mathcal{F}^0 v_{i,p} = 0$ due to (3.20). The final estimates (3.22) and (3.23) follow by applying Proposition A.1, however with a resulting constant $C_N$ which is a priori not uniform in $N$.

## APPENDIX B. SOME BOX INTEGRALS IDENTITIES

The proof of the point (iii) of Lemma 2.2 relies on the following result.

**Proposition B.1.** Let $P = (-1/2, 1/2)^d$ the $d$-dimensional unit hypercube with $d \ge 1$. The following identity holds for any $p, k \ge 0$ and $1 \le i \le d$:

$$ \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{p}{2k+1} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{2k+d-p}{(2k+1)2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x), \quad (B.1) $$

and in particular for $p = 2k + d$:

$$ \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{2k+d+2}} d\sigma(x) = \frac{2k+1}{2k+d} \int_{\partial P} \frac{x_i^{2k}}{|x|^{2k+d}} d\sigma(x). \quad (B.2) $$

*Proof.* Let us denote by $F_i = \{(x_1, \dots, x_d) | x_i = \frac{1}{2}$ and $-\frac{1}{2} \le x_j \le \frac{1}{2}$ for $1 \le j \ne i \le d\}$ the face of the cube with normal $e_i$. By decomposing the integrals over the face $F_i$ and the other faces $F_j$, $j \ne i$, we find

$$ \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = 2 \times \frac{1}{2^{2k}} \int_{F_i} \frac{1}{|x|^p} d\sigma(x) + 2(d-1) \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x) $$
where $n_i$ is the normal tangent to $F_j$ pointing outward $\partial F_j$. On the boundary of the face $F_j$ we have $x_i n_i = 1/2$ and $x_i = \pm \frac{1}{2}$ so by using the divergence theorem again,

$$ \int_{\partial F_j} \frac{x_i^{2k+1} n_i}{|x|^p} d\sigma = \frac{1}{2^{2k}} \int_{\partial F_j} \frac{x_i n_i}{|x|^p} d\sigma(x) = \frac{1}{2^{2k}} \int_{F_j} \partial_i \left( \frac{x_i}{|x|^p} \right) d\sigma(x) = \frac{1}{2^{2k}} \int_{F_j} \left( \frac{1}{|x|^p} - p \frac{x_i^2}{|x|^{p+2}} \right) d\sigma(x). $$

Inserting this equality in (B.4) and using (B.3), we thus arrive at

$$ \begin{aligned} \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma &= \frac{p}{2k+1} \int_{F_j} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{1}{(2k+1)2^{2k}} \int_{F_j} \left( \frac{1}{|x|^p} - p \frac{x_i^2}{|x|^{p+2}} \right) d\sigma(x) \\ &= \frac{p}{(2k+1)2(d-1)} \left( \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) - \frac{1}{2^{2k+2}d} \int_{\partial P} \frac{1}{|x|^{p+2}} d\sigma(x) \right) \\ &\quad + \frac{1}{(2k+1)2^{2k}} \left( \frac{1}{2d} \int_{\partial P} \frac{1}{|x|^p} d\sigma - \frac{p}{2(d-1)} \left( \int_{\partial P} \frac{x_i^2}{|x|^{p+2}} d\sigma(x) - \frac{1}{2^2d} \int_{\partial P} \frac{1}{|x|^{p+2}} d\sigma(x) \right) \right) \\ &= \frac{p}{(2k+1)2(d-1)} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{1}{(2k+1)2^{2k}} \left( \frac{1}{2d} - \frac{p}{2d(d-1)} \right) \int_{\partial P} \frac{1}{|x|^p} d\sigma(x). \end{aligned} $$

Inserting this equality back in (B.3), we obtain

$$ \begin{aligned} \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma &= \frac{p}{2k+1} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \left( \frac{1}{2^{2k}d} + \frac{d-1-p}{(2k+1)2^{2k}d} \right) \int_{\partial P} \frac{1}{|x|^p} d\sigma(x) \\ &= \frac{p}{2k+1} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{2k+d-p}{(2k+1)2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x), \end{aligned} $$

which completes the proof of (B.1). □

## APPENDIX C. MOMENTS OF THE UNIT HYPERSPHERE

The asymptotic formula (4.14) and (4.15) involve the moments integrals $\int_T y^k dy$. In the case $T = B(0, 1)$, these integrals can be computed explicitly.

**Lemma C.1.** Let $l \in \mathbb{N}$ and $1 \le i_1, \dots, i_l \le d$. For $1 \le j \le d$, let $l_j$ be the number of indices equal to $j$: $l_j = \text{Card}\{p \mid i_p = j\}$. The components of the tensor $\int_{B(0,1)} y^l dy$ with $y^l = (y_{i_1} \dots y_{i_l})_{1 \le i_1, \dots, i_l \le d}$ are given by

* $\int_{B(0,1)} y_{i_1} \dots y_{i_l} dy = 0$ whenever one of the numbers $l_j$ is odd for $1 \le j \le d$;

* if $l = 2k$ and $l_j = 2k_j$ for all $1 \le j \le d$, then

$$ \int_{B(0,1)} y_{i_1} \dots y_{i_l} dy = \int_{B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} dy = \frac{2}{(2k+d)} \frac{\prod_{i=1}^d \text{Gam}\left(k_i + \frac{1}{2}\right)}{\text{Gam}\left(k + \frac{d}{2}\right)}, $$

where $\text{Gam}(x) = \int_0^{+\infty} t^{x-1}e^{-t} dt$ denotes the Gamma function.

*Proof.* The first bullet is clear due to the symmetry. For the second bullet, we first use an integration in polar coordinates to obtain

$$ \int_{B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} dy = \frac{1}{2k+d} \int_{\partial B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} d\sigma(y). $$

In order to obtain an explicit expression for the integral arising in the right-hand side, we compute with two different ways the integral

$$ I := \int_{\mathbb{R}^d} e^{-|y|^2} y_1^{2k_1} y_2^{2k_2} \dots y_d^{2k_d} dy. $$

On the one hand, an integration using spherical coordinates yields

$$ I = \int_{0}^{+\infty} e^{-r^2} r^{2k+d-1} dr \int_{\partial B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} d\sigma(y) = \frac{1}{2} \text{Gam}\left(k + \frac{d}{2}\right) \int_{\partial B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} d\sigma(y). $$

On the other hand, Fubini's theorem yields

$$ I = \prod_{i=1}^{d} \int_{-\infty}^{+\infty} e^{-t^2} t^{2k_i} dt = \prod_{i=1}^{d} \text{Gam}\left(k_i + \frac{1}{2}\right). $$

**Remark C.1.** For $k = 0$, we retrieve the well-known formula $|B(0,1)| = \frac{2\text{Gam}\left(\frac{1}{2}\right)^d}{d\text{Gam}\left(\frac{d}{2}\right)} = \frac{\pi^{\frac{d}{2}}}{\text{Gam}\left(1+\frac{d}{2}\right)}$ for the volume of the unit sphere. □
where $1 \le j \le d$ is any index different from $i$. Since $\int_{\partial P} \frac{1}{|x|^d} d\sigma(x) = 2d \int_{F_i} \frac{1}{|x|^p} d\sigma(x)$, this identity can be rewritten as

$$ \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{1}{2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x) + 2(d-1) \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x), \quad (B.3) $$

which also implies

$$ \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{1}{2(d-1)} \left( \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma(x) - \frac{1}{2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x) \right) \quad (B.4) $$

for any $p, k \ge 0$. We then use the (tangential) divergence theorem on the face $F_j$ to obtain

$$ \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma(x) = \frac{1}{2k+1} \int_{F_j} \frac{\partial_i(x_i^{2k+1})}{|x|^p} d\sigma(x) = \frac{p}{2k+1} \int_{F_j} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{1}{2k+1} \int_{\partial F_j} \frac{x_i^{2k+1} n_i}{|x|^p} d\sigma(x), $$

where $n_i$ is the normal tangent to $F_j$ pointing outward $\partial F_j$. On the boundary of the face $F_j$ we have $x_i n_i = 1/2$ and $x_i = \pm \frac{1}{2}$ so by using the divergence theorem again,

$$ \int_{\partial F_j} \frac{x_i^{2k+1} n_i}{|x|^p} d\sigma = \frac{1}{2^{2k}} \int_{\partial F_j} \frac{x_i n_i}{|x|^p} d\sigma(x) = \frac{1}{2^{2k}} \int_{F_j} \partial_i \left( \frac{x_i}{|x|^p} \right) d\sigma(x) = \frac{1}{2^{2k}} \int_{F_j} \left( \frac{1}{|x|^p} - p \frac{x_i^2}{|x|^{p+2}} \right) d\sigma(x). $$

Inserting this equality in (B.4) and using (B.3), we thus arrive at

$$ \begin{aligned} \int_{F_j} \frac{x_i^{2k}}{|x|^p} d\sigma &= \frac{p}{2k+1} \int_{F_j} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{1}{(2k+1)2^{2k}} \int_{F_j} \left( \frac{1}{|x|^p} - p \frac{x_i^2}{|x|^{p+2}} \right) d\sigma(x) \\ &= \frac{p}{(2k+1)2(d-1)} \left( \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) - \frac{1}{2^{2k+2}d} \int_{\partial P} \frac{1}{|x|^{p+2}} d\sigma(x) \right) \\ &\quad + \frac{1}{(2k+1)2^{2k}} \left( \frac{1}{2d} \int_{\partial P} \frac{1}{|x|^p} d\sigma - \frac{p}{2(d-1)} \left( \int_{\partial P} \frac{x_i^2}{|x|^{p+2}} d\sigma(x) - \frac{1}{2^2d} \int_{\partial P} \frac{1}{|x|^{p+2}} d\sigma(x) \right) \right) \\ &= \frac{p}{(2k+1)2(d-1)} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{1}{(2k+1)2^{2k}} \left( \frac{1}{2d} - \frac{p}{2d(d-1)} \right) \int_{\partial P} \frac{1}{|x|^p} d\sigma(x). \end{aligned} $$

Inserting this equality back in (B.3), we obtain

$$ \begin{aligned} \int_{\partial P} \frac{x_i^{2k}}{|x|^p} d\sigma &= \frac{p}{2k+1} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \left( \frac{1}{2^{2k}d} + \frac{d-1-p}{(2k+1)2^{2k}d} \right) \int_{\partial P} \frac{1}{|x|^p} d\sigma(x) \\ &= \frac{p}{2k+1} \int_{\partial P} \frac{x_i^{2k+2}}{|x|^{p+2}} d\sigma(x) + \frac{2k+d-p}{(2k+1)2^{2k}d} \int_{\partial P} \frac{1}{|x|^p} d\sigma(x), \end{aligned} $$

which completes the proof of (B.1). □

## APPENDIX C. MOMENTS OF THE UNIT HYPERSPHERE

The asymptotic formula (4.14) and (4.15) involve the moments integrals $\int_T y^k dy$. In the case $T = B(0, 1)$, these integrals can be computed explicitly.

**Lemma C.1.** Let $l \in \mathbb{N}$ and $1 \le i_1, \dots, i_l \le d$. For $1 \le j \le d$, let $l_j$ be the number of indices equal to $j$: $l_j = \text{Card}\{p \mid i_p = j\}$. The components of the tensor $\int_{B(0,1)} y^l dy$ with $y^l = (y_{i_1} \dots y_{i_l})_{1 \le i_1, \dots, i_l \le d}$ are given by

* $\int_{B(0,1)} y_{i_1} \dots y_{i_l} dy = 0$ whenever one of the numbers $l_j$ is odd for $1 \le j \le d$;

* if $l = 2k$ and $l_j = 2k_j$ for all $1 \le j \le d$, then

$$ \int_{B(0,1)} y_{i_1} \dots y_{i_l} dy = \int_{B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} dy = \frac{2}{(2k+d)} \frac{\prod_{i=1}^d \text{Gam}(k_i + \frac{1}{2})}{\text{Gam}(k + \frac{d}{2})}, $$

where $\text{Gam}(x) = \int_0^{+\infty} t^{x-1}e^{-t} dt$ denotes the Gamma function.

*Proof.* The first bullet is clear due to the symmetry. For the second bullet, we first use an integration in polar coordinates to obtain

$$ \int_{B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} dy = \frac{1}{2k+d} \int_{\partial B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} d\sigma(y). $$

In order to obtain an explicit expression for the integral arising in the right-hand side, we compute with two different ways the integral

$$ I := \int_{\mathbb{R}^d} e^{-|y|^2} y_1^{2k_1} y_2^{2k_2} \dots y_d^{2k_d} dy. $$
REFERENCES

[1] ALLAIRE, G. Continuity of the Darcy's law in the low-volume fraction limit. *Annali della Scuola Normale Superiore di Pisa-Clase di Scienze* **18**, 4 (1991), 475-499.

[2] ALLAIRE, G. Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes I. Abstract framework, a volume distribution of holes. *Archive for Rational Mechanics and Analysis* **113**, 3 (Sept. 1991), 209-259.

[3] AMMARI, H., FITZPATRICK, B., KANG, H., RUIZ, M., YU, S., AND ZHANG, H. *Mathematical and Computational Methods in Photonics and Phononics*, vol. 235. American Mathematical Soc., 2018.

[4] AMMARI, H., AND KANG, H. Boundary layer techniques for solving the Helmholtz equation in the presence of small inhomogeneities. *Journal of mathematical analysis and applications* **296**, 1 (2004), 190-208.

[5] AMMARI, H., KANG, H., AND TOUIBI, K. Boundary layer techniques for deriving the effective properties of composite materials. *Asymptotic Analysis* **41**, 2 (2005), 119-140.

[6] BORN, M., AND BRADBURN, M. The thermodynamics of crystal lattices: II. Calculation of certain lattice sums occurring in thermodynamics. In *Mathematical Proceedings of the Cambridge Philosophical Society* (1943), vol. 39, Cambridge University Press, pp. 104-113.

[7] BRINKMAN, H. C. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. *Flow, Turbulence and Combustion* **1**, 1 (Dec. 1949), 27-34.

[8] BURGERS, J. M. Hydrodynamics. — On the influence of the concentration of a suspension upon the sedimentation velocity (in particular for a suspension of spherical particles). In *Selected Papers of J. M. Burgers*, F. T. M. Nieuwstadt and J. A. Steketee, Eds. Springer Netherlands, Dordrecht, 1941, pp. 452-477.

[9] CHILDRESS, S. Viscous Flow Past a Random Array of Spheres. *The Journal of Chemical Physics* **56**, 6 (Mar. 1972), 2527-2539.

[10] COUSSOT, P. *Mudflow Rheology and Dynamics*. Routledge, 2017.

[11] DE CRISTOFORIS, M. L. Singular perturbation problems in potential theory and applications. In *Complex Analysis and Potential Theory*. World Sci. Publ, 2007, pp. 131-139.

[12] DEBYE, P., AND BUECHE, A. M. Intrinsic Viscosity, Diffusion, and Sedimentation Rate of Polymers in Solution. *The Journal of Chemical Physics* **16**, 6 (June 1948), 573-579.

[13] DHONT, J. K. G. *An Introduction to Dynamics of Colloids*. Elsevier, May 1996.

[14] ESTRADA, R., AND KANWAL, R. P. *A Distributional Approach to Asymptotics*. Birkhäuser Boston, Boston, MA, 2002.

[15] EWALD, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. *Annalen der Physik* **369**, 3 (1921), 253-287.

[16] FABES, E. B., KENIG, C. E., AND VERCHOTA, G. C. The Dirichlet problem for the Stokes system on Lipschitz domains. *Duke Mathematical Journal* **57**, 3 (Dec. 1988), 769-793.

[17] FEPPON, F. High Order Homogenization of the Stokes System in a Periodic Porous Medium. *SIAM Journal on Mathematical Analysis* **53**, 3 (Jan. 2021), 2890-2924.

[18] FEPPON, F., AND AMMARI, H. High Order Topological Asymptotics: Reconciling Layer Potentials and Compound Asymptotic Expansions. *Multiscale Modeling & Simulation* **20**, 3 (Sept. 2022), 957-1001.

[19] FREUND, J. B. Numerical Simulation of Flowing Blood Cells. *Annual Review of Fluid Mechanics* **46**, 1 (Jan. 2014), 67-95.

[20] GALDI, G. *An Introduction to the Mathematical Theory of the Navier-Stokes Equations: Steady-State Problems*. Springer Monographs in Mathematics. Springer New York, New York, NY, 2011.

[21] GIRAULT, V., AND RAVIART, P.-A. *Finite Element Methods for Navier-Stokes Equations*, vol. 5 of Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg, 1986.

[22] GREENGARD, L., AND KROPINSKI, M. C. Integral equation methods for Stokes flow in doubly-periodic domains. *Journal of Engineering Mathematics* **48**, 2 (Feb. 2004), 157-170.

[23] HASIMOTO, H. On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres. *Journal of Fluid Mechanics* **5**, 02 (1959), 317.

[24] HENROT, A., AND PIERRE, M. Shape Variation and Optimization. https://ems.press/books/etm/149, Feb. 2018.

[25] HINCH, E. J. An averaged-equation approach to particle interactions in a fluid suspension. *Journal of Fluid Mechanics* **83**, 4 (Dec. 1977), 695-720.

[26] JING, W. A Unified Homogenization Approach for the Dirichlet Problem in Perforated Domains. *SIAM Journal on Mathematical Analysis* **52**, 2 (Jan. 2020), 1192-1220.

[27] KIM, S., AND KARRILA, S. J. *Microhydrodynamics: Principles and Selected Applications*. Butterworth-Heinemann, Sept. 2013.

[28] LADYZHENSKAYA, O. A. The mathematical theory of viscous incompressible flow. *Gordon & Breach* (1969).

[29] McLEAN, W. C. H. *Strongly Elliptic Systems and Boundary Integral Equations*. Cambridge University Press, Cambridge ; New York, 2000.

[30] MUSOLINO, P., AND MITYUSHEV, V. Asymptotic behavior of the longitudinal permeability of a periodic array of thin cylinders. *Electronic Journal of Differential Equations* **2015** (2015), 1-20.

[31] NÉDÉLEC, J.-C. *Acoustic and Electromagnetic Equations*, vol. 144 of Applied Mathematical Sciences. Springer, New York, NY, 2001.

[32] POZRIKIDIS, C. *Boundary Integral and Singularity Methods for Linearized Viscous Flow*, 1 ed. Cambridge University Press, Feb. 1992.

[33] SANCHEZ-PALENCIA, E. Fluid flow in porous media. In *Non-Homogeneous Media and Vibration Theory*, Lecture Notes in Physics. Springer, Berlin, Heidelberg, 1980, pp. 129-157.

[34] SANGANI, A. S., AND ACRIVOS, A. Slow flow through a periodic array of spheres. *International Journal of Multiphase Flow* **8**, 4 (Aug. 1982), 343-360.

[35] TENNETI, S., GARG, R., AND SUBRAMANIAM, S. Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres. *International Journal of Multiphase Flow* **37**, 9 (2011), 1072-1092.

[36] VOGELIUS, M. S., AND VOLKOV, D. Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter. *ESAIM: Mathematical Modelling and Numerical Analysis* **34**, 4 (2000), 723-748.

[37] WRIGHT, M. *Boundary Value Problems for the Stokes System in Arbitrary Lipschitz Domains*. University of Missouri-Columbia, 2008.
On the one hand, an integration using spherical coordinates yields

$$I = \int_{0}^{+\infty} e^{-r^2} r^{2k+d-1} dr \int_{\partial B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} d\sigma(y) = \frac{1}{2} \mathrm{Gam} \left( k + \frac{d}{2} \right) \int_{\partial B(0,1)} y_1^{2k_1} \dots y_d^{2k_d} d\sigma(y).$$

On the other hand, Fubini's theorem yields

$$I = \prod_{i=1}^{d} \int_{-\infty}^{+\infty} e^{-t^2} t^{2k_i} dt = \prod_{i=1}^{d} \mathrm{Gam} \left( k_i + \frac{1}{2} \right).$$

Remark C.1. For $k = 0$, we retrieve the well-known formula $|B(0, 1)| = \frac{2\mathrm{Gam}(\frac{1}{2})^d}{d\mathrm{Gam}(\frac{d}{2})} = \frac{\pi^{\frac{d}{2}}}{\mathrm{Gam}(1+\frac{d}{2})}$ for the volume of the unit sphere.

## REFERENCES

[1] ALLAIRE, G. Continuity of the Darcy's law in the low-volume fraction limit. *Annali della Scuola Normale Superiore di Pisa* **Classe di Scienze 18**, 4 (1991), 475–499.

[2] ALLAIRE, G. Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes I. Abstract framework, a volume distribution of holes. *Archive for Rational Mechanics and Analysis* 113, 3 (Sept. 1991), 209–259.

[3] AMMARI, H., FITZPATRICK, B., KANG, H., RUIZ, M., YU, S., AND ZHANG, H. *Mathematical and Computational Methods in Photonics and Phononics*, vol. 235. American Mathematical Soc., 2018.

[4] AMMARI, H., AND KANG, H. Boundary layer techniques for solving the Helmholtz equation in the presence of small inhomogeneities. *Journal of mathematical analysis and applications* 296, 1 (2004), 190–208.

[5] AMMARI, H., KANG, H., AND TOUIBI, K. Boundary layer techniques for deriving the effective properties of composite materials. *Asymptotic Analysis* 41, 2 (2005), 119–140.

[6] BORN, M., AND BRADBURN, M. The thermodynamics of crystal lattices: II. Calculation of certain lattice sums occurring in thermodynamics. In *Mathematical Proceedings of the Cambridge Philosophical Society* (1943), vol. 39, Cambridge University Press, pp. 104–113.

[7] BRINKMAN, H. C. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. *Flow, Turbulence and Combustion* 1, 1 (Dec. 1949), 27–34.

[8] BURGERS, J. M. Hydrodynamics. — On the influence of the concentration of a suspension upon the sedimentation velocity (in particular for a suspension of spherical particles). In *Selected Papers of J. M. Burgers*, F. T. M. Nieuwstadt and J. A. Steketee, Eds. Springer Netherlands, Dordrecht, 1941, pp. 452–477.

[9] CHILDRESS, S. Viscous Flow Past a Random Array of Spheres. *The Journal of Chemical Physics* 56, 6 (Mar. 1972), 2527–2539.

[10] COUSSOT, P. *Mudflow Rheology and Dynamics*. Routledge, 2017.

[11] DE CRISTOFORIS, M. L. Singular perturbation problems in potential theory and applications. In *Complex Analysis and Potential Theory*. World Sci. Publ, 2007, pp. 131–139.

[12] DEBYE, P., AND BUECHE, A. M. Intrinsic Viscosity, Diffusion, and Sedimentation Rate of Polymers in Solution. *The Journal of Chemical Physics* 16, 6 (June 1948), 573–579.

[13] DHONT, J. K. G. *An Introduction to Dynamics of Colloids*. Elsevier, May 1996.

[14] DUERINCKX, M., AND GLORIA, A. *On Einstein's Effective Viscosity Formula*, vol. 7 of *Memoirs of the European Mathematical Society*. EMS Press, Berlin, 2023.

[15] ESTRADA, R., AND KANWAL, R. P. *A Distributional Approach to Asymptotics*. Birkhäuser Boston, Boston, MA, 2002.

[16] EWALD, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. *Annalen der Physik* 369, 3 (1921), 253–287.

[17] FABES, E. B., KENIG, C. E., AND VERCHOTA, G. C. The Dirichlet problem for the Stokes system on Lipschitz domains. *Duke Mathematical Journal* 57, 3 (Dec. 1988), 769–793.

[18] FEPPON, F. High Order Homogenization of the Stokes System in a Periodic Porous Medium. *SIAM Journal on Mathematical Analysis* 53, 3 (Jan. 2021), 2890–2924.

[19] FEPPON, F., AND AMMARI, H. High Order Topological Asymptotics: Reconciling Layer Potentials and Compound Asymptotic Expansions. *Multiscale Modeling & Simulation* 20, 3 (Sept. 2022), 957–1001.

[20] FREUND, J. B. Numerical Simulation of Flowing Blood Cells. *Annual Review of Fluid Mechanics* 46, 1 (Jan. 2014), 67–95.

[21] GALDI, G. *An Introduction to the Mathematical Theory of the Navier-Stokes Equations: Steady-State Problems*. Springer Monographs in Mathematics. Springer New York, New York, NY, 2011.

[22] GIRAULT, V., AND RAVIART, P.-A. *Finite Element Methods for Navier-Stokes Equations*, vol. 5 of Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg, 1986.

[23] GREENGARD, L., AND KROPINSKI, M. C. Integral equation methods for Stokes flow in doubly-periodic domains. *Journal of Engineering Mathematics* 48, 2 (Feb. 2004), 157–170.

[24] HASIMOTO, H. On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres. *Journal of Fluid Mechanics* 5, 02 (1959), 317.

[25] HENROT, A., AND PIERRE, M. Shape Variation and Optimization. https://ems.press/books/etm/149, Feb. 2018.

[26] HINCH, E. J. An averaged-equation approach to particle interactions in a fluid suspension. *Journal of Fluid Mechanics* 83, 4 (Dec. 1977), 695–720.

[27] JING, W. A Unified Homogenization Approach for the Dirichlet Problem in Perforated Domains. *SIAM Journal on Mathematical Analysis* 52, 2 (Jan. 2020), 1192–1220.

[28] JING, W., LU, Y., AND PRANGE, C. Unified quantitative analysis of the Stokes equations in dilute perforated domains via layer potentials, Sept. 2024.

[29] KIM, S., AND KARRILA, S. J. *Microhydrodynamics: Principles and Selected Applications*. Butterworth-Heinemann, Sept. 2013.

[30] LADYZHENSKAYA, O. A. The mathematical theory of viscous incompressible flow. *Gordon & Breach* (1969).

[31] LU, Y. Homogenization of Stokes equations in perforated domains: A unified approach. *Journal of Mathematical Fluid Mechanics* 22, 3 (2020), 13.
[32] McLEAN, W. C. H. *Strongly Elliptic Systems and Boundary Integral Equations*. Cambridge University Press, Cambridge ; New York, 2000.

[33] MITYUSHEV, V., AND ADLER, P. Longitudinal Permeability of Spatially Periodic Rectangular Arrays of Circular Cylinders I. A Single Cylinder in the Unit Cell. *ZAMM* **82**, 5 (May 2002), 335.

[34] MUSOLINO, P., AND MITYUSHEV, V. Asymptotic behavior of the longitudinal permeability of a periodic array of thin cylinders. *Electronic Journal of Differential Equations* **2015** (2015), 1–20.

[35] NÉDÉLEC, J.-C. *Acoustic and Electromagnetic Equations*, vol. 144 of *Applied Mathematical Sciences*. Springer, New York, NY, 2001.

[36] POZRIKIDIS, C. *Boundary Integral and Singularity Methods for Linearized Viscous Flow*, 1 ed. Cambridge University Press, Feb. 1992.

[37] SANCHEZ-PALENCIA, E. Fluid flow in porous media. In *Non-Homogeneous Media and Vibration Theory*, Lecture Notes in Physics. Springer, Berlin, Heidelberg, 1980, pp. 129–157.

[38] SANGANI, A. S., AND ACRIVOS, A. Slow flow through a periodic array of spheres. *International Journal of Multiphase Flow* **8**, 4 (Aug. 1982), 343–360.

[39] TENNETI, S., GARG, R., AND SUBRAMANIAM, S. Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres. *International Journal of Multiphase Flow* **37**, 9 (2011), 1072–1092.

[40] VOGELIUS, M. S., AND VOLKOV, D. Asymptotic formulas for perturbations in the electromagnetic fields due to the presence of inhomogeneities of small diameter. *ESAIM: Mathematical Modelling and Numerical Analysis* **34**, 4 (2000), 723–748.

[41] WRIGHT, M. *Boundary Value Problems for the Stokes System in Arbitrary Lipschitz Domains*. University of Missouri-Columbia, 2008.