# A reduced Lagrange multiplier finite element method for fluid-particle interaction

Miguel A Fernández, Céline Grandmont, Fabien Lespagnol, Paolo Zunino

► To cite this version:

Miguel A Fernández, Céline Grandmont, Fabien Lespagnol, Paolo Zunino. A reduced Lagrange multiplier finite element method for fluid-particle interaction. 2026. hal-05521377

HAL Id: hal-05521377

https://inria.hal.science/hal-05521377v1

Preprint submitted on 20 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
A reduced Lagrange multiplier finite element method for
fluid-particle interaction

Miguel A. Fernández¹, Céline Grandmont¹, Fabien Lespagnol¹,²,³, Paolo Zunino²

¹ Sorbonne Université, CNRS, LJLL, Inria, Paris, France

² MOX, Department of Mathematics, Politecnico di Milano, Milano, Italy

³ IMAG, Université de Montpellier, CNRS, Montpellier, France

Abstract

We propose a new reduced-order fictitious domain method for the discretization of the motion of small rigid particles immersed in a viscous incompressible flow. The reduced model is based on the projection of Dirichlet boundary constraints onto a finite-dimensional approximation space of Fourier type, thereby yielding a problem with defective interface conditions imposed through a Lagrange multiplier formulation. We analyze the existence of solutions to the reduced Stokes problem and, for arbitrarily small particles and prescribed rigid velocities at the boundary, we prove its convergence toward the full Stokes problem. The convergence rate depends on both the size of the inclusion and the number of modes in the finite-dimensional space. The numerical discretization of the reduced problem is carried out using the finite element method on a computational mesh that does not conform to the particles, within a fictitious domain framework. We propose a stabilized and robust formulation with respect to the particle size, and prove its stability and convergence in velocity, pressure, and interface fluid stress. A key step in the analysis is the derivation of an inf-sup condition involving both the divergence constraint and the interfacial Dirichlet constraints, with particular attention paid to the dependence on the particle size. The properties of the discretization method are further illustrated by numerical experiments.

Introduction

Systems of small particles suspended in a low Reynolds flow are ubiquitous in numerous applications: from aerosols in the respiratory system [1, 2, 3] and micro-circulation [4, 5], to concrete flow [6] and reinforced plastics [7], among others. As regards the simulation of these systems, several numerical approaches have been reported in the literature, either based on fitted or unfitted mesh discretizations.

In the fitted mesh framework, typically based on an ALE formalism of the fluid with moving
meshes, we refer to [8, 9, 10, 11, 12, 13]. Nevertheless, these methods are known to be cumbersome
for large displacements and rotations of the rigid particles. To overcome this difficulty, unfitted mesh
methods have been developed. We can mention the immersed boundary method [14, 15, 16, 17], the
immersed finite element method [18, 19], the fictitious domain method [20, 21, 22], the extended finite
method [23, 24] and penalization-based methods [25, 26]. In the case of small particles and a Stokes
flow, similarly to the slender body theory developed in 3D, methods based on a truncated multipole
expansion solution of the Stokes equations have been reported in [27, 28].

By building on [29, 30, 31], in this paper we introduce a fictitious domain reduced modeling approach for the simulation of rigid elliptical particles immersed in a two-dimensional Stokes flow. The fundamental idea consists in enforcing the interface coupling via projection onto a finite dimensional Lagrange multiplier space. The objective of this reduction is twofold: to reduce the computational cost induced by the Lagrange multiplier treatment of the interface coupling, as an indirect consequence, to mitigate the need for local refinement of the fluid mesh in the vicinity of the interface. Indeed, coupling methods on non-conforming meshes exhibit stability that is strongly correlated with the ratio between the number of degrees of freedom of the Lagrange multipliers and that of the fluid around the interface. In particular, we introduce a stabilization term for the finite element approximation that removes the locking phenomena on the inf-sup constant that arises when the particle size is small and the number of modes is large.

Reduced interface conditions defined through the introduction of a projection operator have al-
ready been used in the context of so-called mixed-dimensional 3D-1D coupling conditions [30]. In par-
ticular, one objective of this paper is to rigorously justify, within a simplified framework, the reduced
fluid-structure interface conditions introduced in [30], where an elastic beam governed by the Timo-shenko equations is coupled with a surrounding incompressible fluid described by the Navier-Stokes equations. In [31], this strategy was analyzed for two-dimensional elliptic equations of Poisson type. Extending this approach to the Stokes case aims to provide a better understanding of the effect of reduced coupling conditions in a fluid-structure interaction framework.

To this end, we provide a complete mathematical and numerical analysis of the reduced model in a pure Stokes flow setting. Well-posedness is established via an inf-sup condition, with particular attention to its dependence on the particle size $\epsilon$ and the number of modes $N$. Compared to [31], the derivation requires a careful treatment of the incompressibility constraint together with the Lagrange multiplier formulation of the kinematic constraint on the interface. Using separation-of-variables techniques, we prove convergence to the reference (full-order) problem, with the convergence rate explicitly expressed in terms of $\epsilon$ and $N$. In particular, we demonstrate that the reduced-order fictitious-domain formulation provides an accurate approximation of the full model solution (velocity, pressure, and fluid stress) as $\epsilon \to 0$ and $N \to \infty$. A key difficulty, compared to [31], is that both the incompressibility and kinematic constraints must be enforced simultaneously, leading to two Lagrange multipliers (pressure and fluid force) with different discretizations (finite elements and Fourier modes, respectively) which may interact. The main theoretical contribution consists in deriving an inf-sup condition and analyzing its behavior with respect to the particle size, the mesh size, and the number of modes. These theoretical findings are complemented by a comprehensive numerical study. In particular, we provide numerical evidence of the accuracy of the reduced-order model by comparison with fully two-dimensional solutions computed either using an ALE method or a fictitious-domain method with Lagrange multipliers.

The rest of the paper is organized as follows. Section 1 presents the full and reduced fluid-particle interaction models. Section 2 is devoted to the mathematical and numerical analysis of the reduced fluid model in a simplified setting. The numerical approximation of the fluid-particle interaction models is discussed in Section 3. Numerical evidence demonstrating the capabilities of the proposed reduced modeling approach is provided in Section 4. Finally, Section 5 summarizes the main results and outlines some directions for future work.

# 1 Problem setting

In this section, we first introduce the full two-dimensional fluid-particle interaction model, which couples a homogeneous incompressible Newtonian fluid (modeled by the Stokes equations) with an immersed small rigid particle. The associated weak formulation is obtained by enforcing the interface coupling via Lagrange multipliers. Finally, we introduce a reduced coupled problem in which the Lagrange multiplier belongs to a finite-dimensional Fourier-type space.

## 1.1 Fluid-particle interaction

We first introduce the equations governing the kinematics and dynamics of the elliptic particle in $\mathbb{R}^2$. Let $\epsilon > 0$ and $a \ge 1$, we assume that the reference configuration of the particle $\hat{\omega}_\epsilon$, taken for convenience as its initial position $\omega_\epsilon(0)$, is given by

$$ \hat{\omega}_{\epsilon} = \epsilon \mathcal{E}(1, a), \quad \mathcal{E}(1, a) = \left\{ (x, y) \in \mathbb{R}^2 : \frac{x^2}{a^2} + y^2 \le 1 \right\}. \qquad (1.1) $$

The current position $\omega_\epsilon(t)$ of the particle is then parametrized as the image $\omega_\epsilon(t) = \phi(\hat{\omega}_\epsilon, t)$ of the motion map $\phi_\epsilon: \hat{\omega}_\epsilon \times \mathbb{R}_+ \to \mathbb{R}^2$ given by

$$ \phi_{\epsilon}(\hat{x}, t) = r_{\epsilon}(t) + \Lambda_{\epsilon}(t)\hat{x} \quad \forall (\hat{x}, t) \in \hat{\omega}_{\epsilon} \times \mathbb{R}_{+}, \qquad (1.2) $$

where $r_\epsilon: \mathbb{R}_+ \to \mathbb{R}^2$ and $\Lambda_\epsilon: \mathbb{R}_+ \to S_2$ represent the displacement and rotation of the particle, respectively. Here, the set $S_2$ denotes for the Special Orthogonal group defined by

$$ S_2 = \{ \Lambda \in \mathbb{R}^{2 \times 2} : \Lambda^T \Lambda = I, \det \Lambda = 1 \}, \qquad (1.3) $$

where $I$ stands for the identity matrix in $\mathbb{R}^{2 \times 2}$. The rotation matrix $\Lambda_\epsilon$ can be expressed in terms of a scalar function $\theta_\epsilon: \mathbb{R}_+ \to \mathbb{R}$ by the relation

$$ \Lambda_{\epsilon}(\theta_{\epsilon}) = \begin{bmatrix} \cos(\theta_{\epsilon}) & -\sin(\theta_{\epsilon}) \\ \sin(\theta_{\epsilon}) & \cos(\theta_{\epsilon}) \end{bmatrix}, \qquad (1.4) $$

where $\theta_\epsilon(t)$ represents the rotation angle with respect to $e_x$ at time $t$, see Figure 1. We deduce that the Eulerian particle velocity at time $t \in \mathbb{R}_+$ and point $x = \phi_\epsilon(\hat{x}, t) \in \omega_\epsilon(t)$ of the particle is given by

$$ u_b(x, t) = \partial_t \phi_\epsilon(\hat{x}, t) = \dot{r}_\epsilon + \dot{\theta}_\epsilon(x - r_\epsilon)^\perp, $$
where for $\boldsymbol{x} = (x_1, x_2)$, we have used the notation $\boldsymbol{x}^{\perp} = (-x_2, x_1)$, and $\dot{g}$ denotes the time derivative of any function of time $g$. Finally, the dynamics of the particle are described by the following equations of conservation of linear and angular momentum:

$$ \begin{cases} \rho_b A_b^{\varepsilon} \ddot{\boldsymbol{r}}_{\varepsilon} = \boldsymbol{f}, \\ \rho_b I_b^{\varepsilon} \ddot{\theta}_{\varepsilon} = m, \end{cases} \qquad (1.5) $$

where $\rho_b$ is the density of the particle, $A_b^{\varepsilon}$ its area, $I_b^{\varepsilon}$ its moment of inertia, and $\boldsymbol{f}$ and $m$ are the force and torque exerted on the particle, respectively.

<figure><img src="image_2.png" /><figcaption>Figure 1: Current configurations of the elliptical particle $\omega_{\varepsilon}(t)$ and of the fluid domain $\Omega_{\varepsilon}(t)$.</figcaption></figure>

Let $\Omega$ be a fixed smooth bounded domain of $\mathbb{R}^2$ with $(0,0) \in \Omega$. We assume that the particle is immersed in a fluid flow occupying at time $t \in \mathbb{R}_+$ the time-dependent domain

$$ \Omega_{\varepsilon}(t) = \Omega \setminus \overline{\omega_{\varepsilon}(t)} \subset \mathbb{R}^2. $$

The trajectory of the fluid domain is hence given by the following non-cylindrical space-time region:

$$ \mathcal{T}_{\varepsilon} = \bigcup_{t \in \mathbb{R}_+} \Omega_{\varepsilon}(t) \times \{t\}. $$

We assume that the fluid is homogeneous, incompressible and Newtonian, and that its flow is governed by the Stokes equations (low-Reynolds regime), in terms of the fluid velocity $\boldsymbol{u}_{\varepsilon} : \mathcal{T}_{\varepsilon} \to \mathbb{R}^2$ and pressure $p_{\varepsilon} : \mathcal{T}_{\varepsilon} \to \mathbb{R}$ fields as:

$$ \left\{ \begin{array}{ll} \text{div } \boldsymbol{\sigma}(\boldsymbol{u}_{\varepsilon}, p_{\varepsilon}) = \mathbf{0} & \text{in } \Omega_{\varepsilon}(t), \\ \boldsymbol{u}_{\varepsilon} = \mathbf{0} & \text{on } \partial\Omega, \end{array} \right. \qquad (1.6) $$

where the stress tensor $\boldsymbol{\sigma}(\boldsymbol{u}_{\varepsilon}, p_{\varepsilon})$ is given by the relations:

$$ \boldsymbol{\sigma}(\boldsymbol{u}_{\varepsilon}, p_{\varepsilon}) = 2\mu_{f}\boldsymbol{\varepsilon}(\boldsymbol{u}_{\varepsilon}) - p_{\varepsilon}\mathbf{1}, \quad \boldsymbol{\varepsilon}(\boldsymbol{u}_{\varepsilon}) = \frac{1}{2}(\boldsymbol{\nabla}\boldsymbol{u}_{\varepsilon} + (\boldsymbol{\nabla}\boldsymbol{u}_{\varepsilon})^{\top}), \qquad (1.7) $$

with $\mu_{f} > 0$ denotes the dynamic viscosity of the fluid. For the sake of simplicity, homogeneous Dirichlet boundary conditions are enforced in (1.6) on the external boundary $\partial\Omega$.

The proposed coupled model reads therefore as follows: Find $(\boldsymbol{r}_{\varepsilon}, \theta_{\varepsilon})$, $(\boldsymbol{u}_{\varepsilon}, p_{\varepsilon})$ solutions of (1.5) and (1.6) respectively, such that the following coupling conditions are satisfied:

$$ \left\{ \begin{array}{ll} \phi_{\varepsilon}(t) = \boldsymbol{r}_{\varepsilon}(t) + \boldsymbol{\Lambda}_{\varepsilon}(t)\boldsymbol{I}, & \omega_{\varepsilon}(t) = \phi_{\varepsilon}(\hat{\omega}_{\varepsilon}, t), \quad \Omega_{\varepsilon}(t) = \Omega \setminus \overline{\omega_{\varepsilon}(t)}, \\ \boldsymbol{u}_{\varepsilon} = \dot{\boldsymbol{r}}_{\varepsilon} + \dot{\theta}_{\varepsilon}(\boldsymbol{I} - \boldsymbol{r}_{\varepsilon})^{\perp} & \text{on } \partial\omega_{\varepsilon}(t), \\ \boldsymbol{f} = -\int_{\partial\omega_{\varepsilon}(t)} \boldsymbol{\sigma}(\boldsymbol{u}_{\varepsilon}, p_{\varepsilon})\boldsymbol{n}_{\varepsilon}, & \qquad (1.8) \\ m = -\int_{\partial\omega_{\varepsilon}(t)} \boldsymbol{\sigma}(\boldsymbol{u}_{\varepsilon}, p_{\varepsilon})\boldsymbol{n}_{\varepsilon} \cdot (\boldsymbol{I} - \boldsymbol{r})^{\perp}, & \qquad (1.8) \end{array} \right. $$

where $\boldsymbol{I}: \mathbb{R}^2 \to \mathbb{R}^2$ denotes the identity map (i.e., $\boldsymbol{I}(\boldsymbol{x}) = \boldsymbol{x}$ for any $\boldsymbol{x} \in \mathbb{R}^2$) and $\boldsymbol{n}_{\varepsilon}$ the outward unit normal to $\partial\Omega_{\varepsilon}(t)$. The relation $(1.8)_1$ enforces the so-called geometric compatibility condition between the particle and fluid domains. The kinematic coupling is enforced in $(1.8)_2$, while the two remain relations $(1.8)_{3,4}$ respectively interfacial balance of force and torque. Finally, the coupled system (1.5), (1.6) and (1.8) must be complemented with initial conditions on the rigid body, $\boldsymbol{r}_{\varepsilon}(0) = \boldsymbol{r}_0$, $\dot{\boldsymbol{r}}_{\varepsilon}(0) = \dot{\boldsymbol{r}}_0$, $\theta_{\varepsilon}(0) = \theta_0$ and $\dot{\theta}_{\varepsilon}(0) = \dot{\theta}_0$.
## 1.2 Weak formulations

In this section, we first introduce a mixed formulation of the coupled problem (1.5), (1.6) and (1.8) and show that it can be expressed, in an equivalent way, as a fictitious domain model formulation, in which the fluid equations are extended to the entire domain $\Omega$, as the solid motions are rigid ones. Both the kinematic and dynamic coupling conditions between the particle and the fluid on the interface $\partial\omega_{\varepsilon}(t)$ are enforced through Lagrange multipliers. In this formulation, referred to as full coupled problem, a new unknown is thus introduced, which is the Lagrange multiplier $\lambda_{\varepsilon}$, defined on the interface $\partial\omega_{\varepsilon}(t)$, associated to the kinematic constraint (1.8)$_2$. Finally, a reduced formulation is introduced where functional space of the Lagrange multiplier is approximated by a finite dimensional space, leading to an approximated problem, referred in the following as reduced coupled problem.

In what follows, we denote by $L_0^2(O)$ the closed subspace of $L^2(O)$ of functions with zero average on $O$. The scalar product on $L^2(O)$ is denoted by $(\cdot, \cdot)_O$ and the associated norm by $\|\cdot\|_O$. The norm of the Sobolev space $H^s(O)$, $s > 0$ is denoted by $\|\cdot\|_{s,O}$. We also use the notation $H^{-\frac{1}{2}}(\partial O)$ to denote the dual space of $H^{\frac{1}{2}}(\partial O)$, with associated duality pairing $\langle\cdot, \cdot\rangle_{-\frac{1}{2}, \frac{1}{2}, \partial O}$. In the case of vector-valued functions, the associated Hilbert spaces will be indicated in boldface.

### 1.2.1 Full coupled problem

We consider the following mixed weak formulation of the full coupled problem (1.5), (1.6) and (1.8): Find $(u_{\varepsilon}, p_{\varepsilon}, \lambda_{\varepsilon}, r_{\varepsilon}, \theta_{\varepsilon})$ such that

$$ (u_{\varepsilon}(t), p_{\varepsilon}(t), \lambda_{\varepsilon}(t), r_{\varepsilon}(t), \theta_{\varepsilon}(t)) \in H_{\partial\Omega}^{1}(\Omega_{\varepsilon}(t)) \times L_{0}^{2}(\Omega_{\varepsilon}(t)) \times H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R} \quad \text{a.e. } t \in \mathbb{R}^{+} $$

and

$$ \begin{cases} \phi_{\varepsilon}(t) = r_{\varepsilon}(t) + \Lambda_{\varepsilon}(t)I, \quad \omega_{\varepsilon}(t) = \phi_{\varepsilon}(\hat{\omega}_{\varepsilon}, t), \quad \Omega_{\varepsilon}(t) = \Omega \setminus \overline{\omega_{\varepsilon}(t)}, \\ a_{\Omega_{\varepsilon}(t)}^{f}((u_{\varepsilon}, p_{\varepsilon}), (v, q)) + \rho_{b}A_{b}^{\varepsilon}\ddot{r}_{\varepsilon} \cdot \delta r + \rho_{b}I_{b}^{\varepsilon}\ddot{\theta}_{\varepsilon}\delta\theta \\ - \langle \lambda_{\varepsilon}, v - \delta r - \delta\theta(I - r_{\varepsilon})^{\perp} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} + \langle \mu, u_{\varepsilon} - \dot{r}_{\varepsilon} - \dot{\theta}_{\varepsilon}(I - r_{\varepsilon})^{\perp} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = 0 \end{cases} \quad (1.9) $$

for any $(v, q, \mu, \delta r, \delta\theta) \in H_{0, \partial\Omega}^{1}(\Omega_{\varepsilon}(t)) \times L_{0}^{2}(\Omega_{\varepsilon}(t)) \times H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R}$. Here, the fluid bilinear form $a_O^f$, for any domain $O \subset \Omega$, is given by the relations

$$ a_O^f((u, p), (v, q)) = 2(\varepsilon(u), \varepsilon(v))_O - b_1(p, v) + b_1(q, u), \quad b_1(p, v) = (p, \text{div } v)_O. $$

Note that, for $u_{\varepsilon}$ and $p_{\varepsilon}$ smooth enough and by using integration by parts in $\Omega_{\varepsilon}(t)$, it follows from (1.9) that the Lagrange multiplier is nothing but the fluid traction on the interface, viz.,

$$ \lambda_{\varepsilon} = \sigma(u_{\varepsilon}, p_{\varepsilon})n_{\varepsilon} \quad \text{on } \partial\omega_{\varepsilon}(t). $$

Due to the rigid-motion condition (1.8)$_2$ at the fluid-structure interface, and since any rigid-motion velocity (with an associated constant pressure) is solution of the Stokes equations in $\omega_{\varepsilon}(t)$, the weak formulation (1.9) can equivalently be reformulated by defining $u_{\varepsilon}(t)$ and $p_{\varepsilon}(t)$ in the whole domain $\Omega$, namely, by taking

$$ (u_{\varepsilon}(t), p_{\varepsilon}(t)) \in H_0^1(\Omega) \times Q_{\varepsilon}(t) \quad \text{a.e. } t \in \mathbb{R}^+, \quad (1.10) $$

with the constrained pressure space

$$ Q_{\varepsilon}(t) = \{q \in L^2(\Omega) : q|_{\Omega_{\varepsilon}(t)} \in L_0^2(\Omega_{\varepsilon}(t)), q|_{\omega_{\varepsilon}(t)} \in L_0^2(\omega_{\varepsilon}(t))\} \quad (1.11) $$

and by enforcing

$$ \left\{ \begin{aligned} & a_{\Omega}^{f}((u_{\varepsilon}, p_{\varepsilon}), (v, q)) + \rho_{b}A_{b}^{\varepsilon}\ddot{r}_{\varepsilon} \cdot \delta r + \rho_{b}I_{b}^{\varepsilon}\ddot{\theta}_{\varepsilon}\delta\theta \\ & - \langle \lambda_{\varepsilon}, v - \delta r - \delta\theta(x - r_{\varepsilon})^{\perp} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} + \langle \mu, u_{\varepsilon} - \dot{r}_{\varepsilon} - \dot{\theta}(x - r_{\varepsilon})^{\perp} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = 0 \end{aligned} \right. \quad (1.12) $$

for any $(v, q, \mu, \delta r, \delta\theta) \in H_0^1(\Omega) \times Q_{\varepsilon}(t) \times H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^2 \times \mathbb{R}$. Note that, in (1.12), the fluid bilinear form is taken in the whole domain $\Omega$. Furthermore, in order to remove the indetermination of the pressure in the two sub-domains (each defined up to an additive constant, one in $\Omega_{\varepsilon}(t)$ and the other in $\omega_{\varepsilon}(t)$), the trial and test functions for the pressure are taken in the closed time-dependent subspace $Q_{\varepsilon}(t)$, given by (1.11). As a result, the pressure $p_{\varepsilon}$ vanishes in $\omega_{\varepsilon}(t)$, due to the rigid-motion nature of the velocity.

From a practical point of view, it is however more convenient to remove the pressure indetermi- nation by looking for the pressure $p_{\varepsilon}(t)$ in $L_0^2(\Omega)$, i.e., by fixing one of the constants in the whole
domain $\Omega$. The remaining free additive constant is fixed by considering the following constrained space for the Lagrange multiplier:

$$ H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) = \{\lambda \in H^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) : \langle\lambda, n_\varepsilon\rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} = 0\}. \quad (1.13) $$

This enables us to introduce the following standard functional spaces for the velocity and the pressure

$$ V = H_0^1(\Omega), \quad Q = L_0^2(\Omega), $$

respectively, and reformulate the coupled problem (1.5), (1.6), (1.8) as follows: Find $(u_\varepsilon, p_\varepsilon, \lambda_\varepsilon, r_\varepsilon, \theta_\varepsilon)$ such that

$$ (u_\varepsilon(t), p_\varepsilon(t), \lambda_\varepsilon(t), r_\varepsilon(t), \theta_\varepsilon(t)) \in V \times Q \times H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) \times \mathbb{R}^2 \times \mathbb{R} \quad \text{a.e. } t \in \mathbb{R}^+ \quad (1.14) $$

and

$$ \begin{aligned} \phi_\varepsilon(t) = \mathbf{r}_\varepsilon(t) + \mathbf{\Lambda}_\varepsilon(t)\mathbf{I}, \quad \omega_\varepsilon(t) = \phi_\varepsilon(\hat{\omega}_\varepsilon, t), \quad \Omega_\varepsilon(t) = \Omega \setminus \overline{\omega_\varepsilon(t)}, \\ \left\{ \begin{array}{l} a_\Omega^f((u_\varepsilon, p_\varepsilon), (v, q)) + \rho_b A_b^\varepsilon \ddot{\mathbf{r}}_\varepsilon \cdot \delta\mathbf{r} + \rho_b I_b^\varepsilon \ddot{\theta}_\varepsilon \delta\theta \\ - \langle \lambda_\varepsilon, v - \delta\mathbf{r} - \delta\theta(\mathbf{I} - \mathbf{r}_\varepsilon)^\perp \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} + \langle \mu, u_\varepsilon - \dot{\mathbf{r}}_\varepsilon - \dot{\theta}_\varepsilon(\mathbf{I} - \mathbf{r}_\varepsilon)^\perp \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} = 0 \end{array} \right. \end{aligned} \quad (1.15) $$

for any $(v, q, \mu, \delta r, \delta \theta) \in V \times Q \times H_0^{-\frac{1}{2}}(\partial \omega_\varepsilon(t)) \times \mathbb{R}^2 \times \mathbb{R}$.

The equivalence between the formulations (1.10)–(1.12) and (1.14)–(1.15) is rigorously established in Section A.1 of the Appendix. For the mathematical and numerical analysis presented in Section 2 below and also for the computer implementation, the weak formulation (1.14)–(1.15) appears to be more appropriate, in particular because $V \times Q$ is the standard functional space for the Stokes problem.

### 1.2.2 Reduced coupled problem

In this section, we introduce a reduced model of the fluid-particle interaction system, alternative to the full weak formulation (1.14)–(1.15). In the spirit of [29, 30, 31], fundamental idea consists in approximating the space of the Lagrange multiplier $H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t))$ by a Fourier-type finite dimensional space $F_N(t)$, for $N \in \mathbb{N}$.

The proposed reduced coupled problem hence reads as follows: Find $(u_\varepsilon^{(N)}, p_\varepsilon^{(N)}, \lambda_\varepsilon^{(N)}, r_\varepsilon^{(N)}, \theta_\varepsilon^{(N)})$ such that

$$ (u_\varepsilon^{(N)}(t), p_\varepsilon^{(N)}(t), \lambda_\varepsilon^{(N)}(t), r_\varepsilon^{(N)}(t), \theta_\varepsilon^{(N)}(t)) \in V \times Q \times F_N(t) \times \mathbb{R}^2 \times \mathbb{R} \quad \text{a.e. } t \in \mathbb{R}^+ $$

and

$$ \phi_\varepsilon^{(N)}(t) = \mathbf{r}_\varepsilon^{(N)}(t) + \mathbf{\Lambda}_\varepsilon^{(N)}(t)\mathbf{I}, \quad \omega_\varepsilon^{(N)}(t) = \phi_\varepsilon^{(N)}(\hat{\omega}_\varepsilon, t), \quad \Omega_\varepsilon^{(N)}(t) = \Omega \setminus \overline{\omega_\varepsilon^{(N)}(t)}, \quad (1.16) $$

$$ \left\{ \begin{array}{l} a_\Omega^f((u_\varepsilon^{(N)}, p_\varepsilon^{(N)}), (v, q)) + \rho_b A_b^\varepsilon \ddot{\mathbf{r}}_\varepsilon^{(N)} \cdot \delta\mathbf{r} + \rho_b I_b^\varepsilon \ddot{\theta}_\varepsilon^{(N)} \delta\theta \\ - (\lambda_\varepsilon^{(N)}, v - \delta\mathbf{r} - \delta\theta(\mathbf{I} - \mathbf{r}_\varepsilon^{(N)})^\perp)_{\varepsilon,t} + (\mu_N, u_\varepsilon^{(N)} - \dot{\mathbf{r}}_\varepsilon^{(N)} - \dot{\theta}_\varepsilon^{(N)}(\mathbf{I} - \mathbf{r}_\varepsilon^{(N)})^\perp)_{\varepsilon,t} = 0 \end{array} \right. \quad (1.17) $$

for any $(v, q, \mu_N, \delta r, \delta \theta) \in V \times Q \times F_N(t) \times \mathbb{R}^2 \times \mathbb{R}$. Here, the notation $(\cdot, \cdot)_{\varepsilon,t}$ stands for a rescaling of the inner product in $L^2(\partial\omega_\varepsilon^{(N)}(t))$, defined by

$$ (\mu, v)_{\varepsilon,t} = \int_{0}^{2\pi} (\mu \circ \phi_{\varepsilon}^{(N)}(t))(\varepsilon) \cdot (v \circ \phi_{\varepsilon}^{(N)}(t))(\varepsilon) = \int_{\partial\omega_{\varepsilon}^{(N)}(t)} (\mu \cdot v) \hat{h} \circ (\phi_{\varepsilon}^{(N)}(t))^{-1}, \quad (1.18) $$

where

$$ \hat{h}(\varepsilon)(\nu) = \left( \varepsilon \sqrt{a^2 \cos^2 \nu + \sin^2 \nu} \right)^{-1}. $$

The rescaling (1.18) is motivated by the fact that one of our aim is to analyze the behavior of this approximation strategy with respect to the size $\varepsilon$ of the rigid particle. With this rescaling, the reduced counterpart of $\lambda_\varepsilon(t)$ is $\hat{h}(\varepsilon) \circ (\phi_\varepsilon^{(N)}(t))^{-1} \lambda_\varepsilon^{(N)}(t)$ (see also Remark 2.1) and thus in the case of a circular structure, namely $a = 1$, the reduced counterpart is $\varepsilon^{-1} \lambda_\varepsilon^{(N)}(t)$.

Finally, the reduced space $F_N(t)$ is constructed by leveraging the elliptical shape of the fluid-structure interface. To this purpose, we first introduce a local coordinate $\nu$ on the boundary of the elliptical particle, in the reference configuration, such that for any point $\hat{\mathbf{x}} \in \partial\hat{\omega}_\varepsilon$ we have

$$ \hat{\mathbf{x}} = \varepsilon(e_x \cos \nu + e_y a \sin \nu), \quad \nu \in (0, 2\pi). $$
For any field vector field $\hat{v} : \partial\hat{\omega}_\varepsilon \to \mathbb{R}^2$, we can hence introduce the re-scaled vector function $\hat{v}(\varepsilon) : (0, 2\pi) \to \mathbb{R}^2$ by

$$ \hat{v}(\varepsilon)(\nu) = \hat{v}(\varepsilon(\cos(\nu)e_x + a\sin(\nu)e_y)) \quad \forall \nu \in (0, 2\pi). $$

Finally, for $N \in \mathbb{N}$, the reduced space $F_N(t)$ is defined as:

$$ F_N(t) = \left\{ \hat{v} \circ (\phi_{\varepsilon}^{(N)}(t))^{-1} \in H_0^{-\frac{1}{2}}(\partial\omega_{\varepsilon}^{(N)}(t)) : \right. \\ \left. \hat{v}(\varepsilon)(\nu) = a_0 + \sum_{k=1}^{N} (a_k \cos(k\nu) + b_k \sin(k\nu)), a_k, b_k \in \mathbb{R}^2 \right\}. \quad (1.19) $$

**Remark 1.1.** The reduction procedure performed to derive the coupled problem (1.16)-(1.17) can be interpreted as a discretization of the space $H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon)$ where the basis functions are taken as the N first Fourier modes.

## 2 Mathematical and numerical analysis of the reduced Stokes problem

In this section, we focus on the fluid counterpart of the full and reduced fluid-particle interaction models introduced in the previous section. We hence simply limit the study to the steady Stokes problem in a fixed domain, with prescribed rigid-body velocity data on the particle boundary. We assume that the immersed particle is a disk, i.e., $a = 1$ in (1.1). Furthermore, in order to facilitate the presentation, the particle is assumed to be centered at the origin and the viscosity coefficient is taken as $\mu_f = \frac{1}{2}$. We analyze the well-posedness of the reduced Stokes model (Section 2.1) and the induced modeling error with respect to both the obstacle size $\varepsilon$ and the number of modes $N$ (Section 2.2). We also introduce and analyze a stabilized spatial approximation of this reduced model (Sections 2.3–2.5).

In what follows, the rigid particle $\omega_\varepsilon \subset \Omega$ is assumed to be fixed, so that the corresponding fluid domain $\Omega_\varepsilon = \Omega \setminus \bar{\omega}_\varepsilon$ is also fixed. The considered full Stokes problem reads therefore as follows: Find $(u_\varepsilon, p_\varepsilon, \lambda_\varepsilon) \in V \times Q \times H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon)$ such that

$$ \left\{ \begin{array}{ll} (\varepsilon(u_\varepsilon), \varepsilon(v))_\Omega - (p_\varepsilon, \operatorname{div} v)_\Omega - \langle \lambda_\varepsilon, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon} = 0 & \forall v \in V, \\ (q, \operatorname{div} u_\varepsilon)_\Omega = 0 & \forall q \in Q, \\ \langle \mu, u_\varepsilon \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon} = \langle \mu, v_b + \varepsilon w_b e_\theta \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon} & \forall \mu \in H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon), \end{array} \right. \quad (2.1) $$

where $e_\theta = -\sin\theta e_x + \cos\theta e_y$ stands for the angular basis vector of the polar coordinate system and $(v_b, w_b) \in \mathbb{R}^2 \times \mathbb{R}$ are given translational and angular velocities.

Similarly to (1.19), for $N \in \mathbb{N}$, we introduce the following reduced Fourier-type space $F_N$ over $\partial\omega_\varepsilon$:

$$ F_N = \left\{ v \in H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon) : v(\varepsilon)(\nu) = a_0 + \sum_{k=1}^{N} (a_k \cos(k\nu) + b_k \sin(k\nu)), a_k, b_k \in \mathbb{R}^2 \right\}, $$

where $v(\varepsilon) : (0, 2\pi) \to \mathbb{R}^2$ is defined as $v(\varepsilon)(\nu) = v(\varepsilon(e_x \cos\nu + e_y \sin\nu))$ for any $\nu \in (0, 2\pi)$. Finally, the re-scaled inner product in $L^2(\partial\omega_\varepsilon)$ is defined as:

$$ (\mu, v)_\varepsilon = \int_0^{2\pi} \mu(\varepsilon) \cdot v(\varepsilon) = \varepsilon^{-1}(\mu, v)_{\partial\omega_\varepsilon}. \quad (2.2) $$

The reduced Stokes problem reads therefore as follows: Find $(u_\varepsilon^{(N)}, p_\varepsilon^{(N)}, \lambda_\varepsilon^{(N)}) \in V \times Q \times F_N$ such that

$$ \left\{ \begin{array}{ll} (\varepsilon(u_\varepsilon^{(N)}), \varepsilon(v))_{\Omega} - (p_\varepsilon^{(N)}, \operatorname{div} v)_{\Omega} - (\lambda_\varepsilon^{(N)}, v)_\varepsilon = 0 & \forall v \in V, \\ (q, \operatorname{div} u_\varepsilon^{(N)})_\Omega = 0 & \forall q \in Q, \\ (\mu_N, v)_\varepsilon = (\mu_N, v_b + \varepsilon w_b e_\theta)_\varepsilon & \forall \mu_N \in F_N. \end{array} \right. \quad (2.3) $$

Finally, we introduce some notation and fundamental results on Fourier series that will be useful for the subsequent mathematical analysis of the reduced coupled problem (2.3). We shall make use of the following Fourier norm $\| \cdot \|_{\frac{1}{2}, \varepsilon}$ on $H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)$ defined as

$$ \|v\|_{\frac{1}{2}, \varepsilon}^{\mathcal{F}} = \left( (|\xi_{0,r}^{\varepsilon,c}|^2 + |\xi_{0,\theta}^{\varepsilon,c}|^2) + \sum_{k=1}^{\infty} (1+k)(|\xi_{k,r}^{\varepsilon,c}|^2 + |\xi_{k,\theta}^{\varepsilon,c}|^2 + |\xi_{k,r}^{\varepsilon,s}|^2 + |\xi_{k,\theta}^{\varepsilon,s}|^2) \right)^{\frac{1}{2}} \quad (2.4) $$
for any $v \in H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$, and where

$$ \xi_{k,r}^{\varepsilon,j} = \mathcal{F}_k^j(v(\varepsilon) \cdot e_r), \quad \xi_{k,\theta}^{\varepsilon,j} = \mathcal{F}_k^j(v(\varepsilon) \cdot e_\theta) \quad \forall (k,j) \in \{0, \dots, N\} \times \{\text{c,s}\}, $$

with the operators $\mathcal{F}_k^j : L^2(\partial\omega_\varepsilon) \to \mathbb{R}$ defined by

$$ \mathcal{F}_0^c v = \frac{1}{2\pi} (v, 1)_{(0,2\pi)}, \quad \mathcal{F}_k^c v = \frac{1}{\pi} (v, \cos(k \cdot))_{(0,2\pi)}, \quad \mathcal{F}_k^s v = \frac{1}{\pi} (v, \sin(k \cdot))_{(0,2\pi)} \quad (2.5) $$

for $k \ge 1$ and $v \in L^2(\partial\omega_\varepsilon)$. The norm $\|\cdot\|_{\frac{1}{2},\varepsilon}$ is well defined on $H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ and is equivalent to the following rescaled norm

$$ \|v\|_{\frac{1}{2},\varepsilon} = \|v(\varepsilon)\|_{\frac{1}{2},(0,2\pi)} \quad \forall v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon). $$

Moreover, due to the chosen rescaling in $\varepsilon$, there exists a constant $c$ independent of $\varepsilon$ such that

$$ c^{-1} \|v\|_{\frac{1}{2},\varepsilon} \le \|v\|_{\frac{1}{2},\varepsilon} \le c \|v\|_{\frac{1}{2},\varepsilon} \quad \forall v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon). \quad (2.6) $$

Let $\Pi_N^\varepsilon : L^2(\partial\omega_\varepsilon) \to F_N$ be the projection operator on $F_N$ associated to the scalar product $(\cdot, \cdot)_\varepsilon$, that is,

$$ (\mu_N, \Pi_N^\varepsilon v)_\varepsilon = (\mu_N, v)_\varepsilon \quad \forall \mu_N \in F_N. $$

Owing to the definition (2.4), we have

$$ \|\Pi_N^\varepsilon v\|_{\frac{1}{2},\varepsilon} \le \|v\|_{\frac{1}{2},\varepsilon} \quad \forall v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon), \quad (2.7) $$

so that, using (2.6), we get

$$ \|\Pi_N^\varepsilon v\|_{\frac{1}{2},\varepsilon} \le \|v\|_{\frac{1}{2},\varepsilon} \quad \forall v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon). $$

The projection operator $\Pi_N^\varepsilon$ is hence continuous in $H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ with respect to the norm $\|\cdot\|_{\frac{1}{2},\varepsilon}$. Moreover the associated operator norm which is uniformly bounded in $\varepsilon$. Furthermore, we notice that $H^{-\frac{1}{2}}(\partial\omega_\varepsilon)$ can be equipped with the duality product $\langle \cdot, \cdot \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}$ and the associated norm $\|\cdot\|_{-\frac{1}{2},\varepsilon}$ respectively defined as

$$ \langle \mu, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \varepsilon} = \varepsilon^{-1} \langle \mu, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon}, \quad \|\mu\|_{-\frac{1}{2}, \varepsilon} = \sup_{v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)} \frac{\langle \mu, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \varepsilon}}{\|v\|_{\frac{1}{2}, \varepsilon}}. \quad (2.8) $$

for $v \in H^{\frac{1}{2}}(\partial\omega_\varepsilon)$ and $\mu \in H^{-\frac{1}{2}}(\partial\omega_\varepsilon)$.

**Remark 2.1.** Note that the duality product $\langle \cdot, \cdot \rangle_{-\frac{1}{2}, \frac{1}{2}, \varepsilon}$ satisfies in particular

$$ \langle \mu, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \varepsilon} = \varepsilon^{-1} (\mu, v)_{\partial\omega_\varepsilon} = (\mu, v)_\varepsilon \quad \forall (\mu, v) \in L^2(\partial\omega_\varepsilon) \times H^{\frac{1}{2}}(\partial\omega_\varepsilon). \quad (2.9) $$

## 2.1 Well-posedness

This section is devoted to the well-posedness analysis of the reduced Stokes problem (2.3). The main result is stated in the following theorem.

**Theorem 2.2.** Problem (2.3) admits a unique solution $(u_\varepsilon^{(N)}, p_\varepsilon^{(N)}, \lambda_\varepsilon^{(N)})$ in $V \times Q \times F_N$. Moreover, the following a priori estimate holds:

$$ \|u_\varepsilon^{(N)}\|_{1,\Omega} + \|p_\varepsilon^{(N)}\|_{\Omega} + \|\lambda_\varepsilon^{(N)}\|_{-\frac{1}{2},\varepsilon} \lesssim |v_b| + \varepsilon|w_b|. \quad (2.10) $$

This result states that the velocity the pressure depend continuously on the data with continuity constant uniformly bounded with respect to the number of Fourier modes $N$ and the size $\varepsilon$ of the rigid particle. It is also the case of the Lagrange multiplier in the rescaled norm. Theorem 2.2 is a consequence of Lemma 2.5 below, which itself derives from the asymptotic properties of the trace operator given by Lemma 2.3. The proof of these preliminary lemmas are given in the Appendix (Section A.3).

The first lemma gives an estimate on the re-scaled trace norm over $\partial\omega_\varepsilon$ of $H^1$ functions defined in $\Omega$.

**Lemma 2.3.** For any $v \in H^1(\Omega)$, there holds

$$ \|v\|_{\frac{1}{2},\varepsilon} \lesssim |\log \varepsilon|^{\frac{1}{2}} \|v\|_{1,\Omega}. $$
Proof. See Section A.3.

The second lemma introduces a divergence free continuous lifting into $\Omega$ of functions belonging to

$$H_0^{\frac{1}{2}}(\partial\omega_\varepsilon) = \{\eta \in H^{\frac{1}{2}}(\partial\omega_\varepsilon) : \int_{\partial\omega_\varepsilon} \eta \cdot n_\varepsilon = 0\},$$

endowed with the rescaled norm $\|\cdot\|_{\frac{1}{2},\varepsilon}$. With this choice the resulting lifting operator has a continuity constant independent of $\varepsilon$.

**Lemma 2.4.** Let $\eta \in H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)$ be given, there exists a unique $(v_\eta^\eta, q_\eta^\eta) \in V \times Q$ such that

$$\left\{
\begin{array}{ll}
-\Delta v_\eta^\eta + \nabla q_\eta^\eta = 0 & \text{in } \Omega_\varepsilon \cup \omega_\varepsilon, \\
\operatorname{div} v_\eta^\eta = 0 & \text{in } \Omega, \\
v_\eta^\eta = \eta & \text{on } \partial\omega_\varepsilon.
\end{array}
\right.
$$

Moreover, the following a priori estimate holds:

$$\|v_{\varepsilon}^{\eta}\|_{1,\Omega} \lesssim \|\eta\|_{\frac{1}{2},\varepsilon}. \qquad (2.11)$$

Proof. See Section A.3.

The third lemma is a variant of the surjectivity theorem for the divergence operator.

**Lemma 2.5.** Let $q \in Q$ be given, there exists $w_\varepsilon^q \in V$ such that

$$\left\{
\begin{array}{ll}
\operatorname{div} w_\varepsilon^q = q & \text{in } \Omega, \\
w_\varepsilon^q = \frac{1}{|\partial\omega_\varepsilon|} \left( \int_{\omega_\varepsilon} q \right) n_\varepsilon & \text{on } \partial\omega_\varepsilon,
\end{array}
\right.
$$

and

$$\|w_{\varepsilon}^{q}\|_{1,\Omega} \lesssim \|q\|_{\Omega}.$$

Proof. See Section A.3.

To prove Theorem 2.2, we use the following variant of the BNB theorem (see, e.g., [32] or [33, Theorem 2.34]) adapted to the case of twofold saddle point problems.

**Theorem 2.6.** Let $Q_1$ and $Q_2$ be two reflexive Banach spaces, $a: X \times X \to \mathbb{R}$, $b_1: Q_1 \times X \to \mathbb{R}$, $b_2: Q_2 \times X \to \mathbb{R}$ three continuous bilinear forms, $c_1: Q_1 \to \mathbb{R}$, $c_2: Q_2 \to \mathbb{R}$, $d: X \to \mathbb{R}$ three continuous linear forms, we consider the twofold saddle point problem: find $(u, \lambda_1, \lambda_2) \in X \times Q_1 \times Q_2$ such that

$$\left\{
\begin{array}{ll}
a(u, v) - b_1(\lambda_1, v) - b_2(\lambda_2, v) = d(v) & \forall v \in X, \\
b_1(\mu_1, u) = c_1(\mu_1) & \forall \mu_1 \in Q_1, \\
b_2(\mu_2, u) = c_2(\mu_2) & \forall \mu_2 \in Q_2.
\end{array}
\right.
\qquad (2.12)$$

Let

$$Z_{b_i} = \{v \in X : b_i(\mu_i, v) = 0 \quad \forall \mu_i \in Q_i\} \subset X, \quad i = 1, 2. \qquad (2.13)$$

We suppose that there exists $\beta_1, \beta_2 > 0$ such that

$$\sup_{v \in Z_{b_2}} \frac{b_1(\lambda_1, v)}{\|v\|_X} \ge \beta_1 \|\lambda_1\|_{Q_1} \quad \forall \lambda_1 \in Q_1,$$

$$\sup_{v \in Z_{b_1}} \frac{b_2(\lambda_2, v)}{\|v\|_X} \ge \beta_2 \|\lambda_2\|_{Q_2} \quad \forall \lambda_2 \in Q_2,$$

and that there exists $\alpha > 0$ such that

$$a(v, v) \ge \alpha \|v\|_X \quad \forall v \in X. \qquad (2.14)$$

Then Problem (2.12) is well posed and the following estimates on $u, \lambda_1$ and $\lambda_2$ holds:

$$\|u\|_{X} \le \alpha^{-1} \|d\| + \beta_{1}^{-1} (1 + \alpha^{-1} \|a\|) \|c_{1}\| + \beta_{2}^{-1} (1 + \alpha^{-1} \|a\|) \|c_{2}\|, \qquad (2.15)$$

and

$$\|\lambda_1\|_{Q_1} \le \beta_1^{-1} (\|d\| + \|a\| \|u\|_X), \quad \|\lambda_2\|_{Q_2} \le \beta_2^{-1} (\|d\| + \|a\| \|u\|_X).$$
Proof. See Section A.3.

We now have all the ingredients to prove Theorem 2.2.

Proof of Theorem 2.2. We first note that problem (2.3) can be recast into the abstract formalism of Theorem 2.6, namely: Find $(\mathbf{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)}, \lambda_{\varepsilon}^{(N)}) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{F}_N$ such that

$$
\begin{align}
& a(\mathbf{u}_{\varepsilon}^{(N)}, \mathbf{v}) - b_1(p_{\varepsilon}^{(N)}, \mathbf{v}) - b_2(\boldsymbol{\lambda}_{\varepsilon}^{(N)}, \mathbf{v}) = 0 \quad \forall \mathbf{v} \in \mathbf{V}, \nonumber \\
& b_1(q, \mathbf{u}_{\varepsilon}^{(N)}) = 0 \quad \forall q \in \mathcal{Q}, \tag{2.16} \\
& b_2(\boldsymbol{\mu}_N, \mathbf{u}_{\varepsilon}^{(N)}) = c(\boldsymbol{\mu}_N) \quad \forall \boldsymbol{\mu}_N \in \mathbf{F}_N, \nonumber
\end{align}
$$

where the bilinear forms $a: \mathbf{V} \times \mathbf{V} \to \mathbb{R}$, $b_1: \mathbf{Q} \times \mathbf{V} \to \mathbb{R}$, $b_2: \mathbf{F}_N \times \mathbf{V} \to \mathbb{R}$ and the linear form $c: \mathbf{F}_N \to \mathbb{R}$ are given by

$$
\begin{align}
a(\mathbf{u}, \mathbf{v}) &= (\varepsilon(\mathbf{u}), \varepsilon(\mathbf{v}))_{\Omega}, \quad b_1(q, \mathbf{v}) = (q, \operatorname{div} \mathbf{v})_{\Omega}, \quad b_2(\boldsymbol{\mu}_N, \mathbf{v}) = (\boldsymbol{\mu}_N, \mathbf{v})_{\varepsilon} \tag{2.17} \\
c(\boldsymbol{\mu}_N) &= (\boldsymbol{\mu}_N, \mathbf{v}_b + \varepsilon w_b e_{\theta})_{\varepsilon}. \notag
\end{align}
$$

Using Korn's inequality in $\Omega$ and Cauchy-Schwarz's inequality, the bilinear form $a$ is continuous and coercive in $\mathbf{V}$, with norm $\|a\|$ and coercivity constant $\alpha > 0$ independent of both $N$ and $\varepsilon$. Similarly, using Cauchy-Schwarz's inequality and Lemma 2.3, the following properties are satisfied by $b_1$ and $b_2$:

$$
\left\{
\begin{array}{l}
|b_1(q, \mathbf{v})| \le \|q\|_{\Omega} \|\operatorname{div} \mathbf{v}\|_{1, \Omega}, \\
|b_2(\boldsymbol{\mu}_N, \mathbf{v})| \le \|\boldsymbol{\mu}_N\|_{-\frac{1}{2}, \varepsilon} \|\mathbf{v}\|_{\frac{1}{2}, \varepsilon} \le |\log(\varepsilon)|^{\frac{1}{2}} \|\boldsymbol{\mu}_N\|_{-\frac{1}{2}, \varepsilon} \|\mathbf{v}\|_{1, \Omega}
\end{array}
\right.
\quad (2.18)
$$

for any $(q, \boldsymbol{\mu}_N, \mathbf{v}) \in \mathbf{Q} \times \mathbf{F}_N \times \mathbf{V}$. This implies that the bilinear forms $b_1$ and $b_2$ are continuous and

$$
\|b_1\| \lesssim 1, \quad \|b_2\| \lesssim |\log(\varepsilon)|^{\frac{1}{2}}. \tag{2.19}
$$

Furthermore, we have

$$
|c(\boldsymbol{\mu}_N)| \lesssim \|\boldsymbol{\mu}_N\|_{-\frac{1}{2}, \varepsilon} \|\mathbf{u}_b\|_{\frac{1}{2}, \varepsilon} \lesssim \|\boldsymbol{\mu}_N\|_{-\frac{1}{2}, \varepsilon} (|\mathbf{v}_b| + \varepsilon|w_b|) \quad \forall \boldsymbol{\mu}_N \in \mathbf{F}_N, \quad (2.20)
$$

which gives the continuity of the linear form $c$. We also deduce from (2.20) that the norm $\|c\|$ of the linear form $c$ is bounded by $|\mathbf{v}_b| + \varepsilon|w_b|$ up to multiplicative constant uniformly independent of $N$ and $\varepsilon$. Let us now prove the inf-sup conditions for the bilinear forms $b_1$, $b_2$ by establishing the existence of two constants $\beta_1$, $\beta_2 > 0$ uniformly independent of $N$ and $\varepsilon$ satisfying

$$
\left\{
\begin{array}{ll}
\displaystyle \sup_{\mathbf{v} \in \mathbf{V}, \Pi_N^\varepsilon(\mathbf{v})=0} & \displaystyle \frac{b_1(q, \mathbf{v})}{\|\mathbf{v}\|_{1,\Omega}} \ge \beta_1 \|q\|_{\Omega} \quad \forall q \in \mathcal{Q}, \\
\\
\displaystyle \sup_{\mathbf{v} \in \mathbf{H}_0^1(\Omega), \operatorname{div} \mathbf{v}=0} & \displaystyle \frac{b_2(\boldsymbol{\mu}_N, \mathbf{v})}{\|\mathbf{v}\|_{1,\Omega}} \ge \beta_2 \|\boldsymbol{\mu}_N\|_{-\frac{1}{2},\varepsilon} \quad \forall \boldsymbol{\mu}_N \in \mathbf{F}_N.
\end{array}
\right.
\quad (2.21)
$$

Let $q \in \mathbf{Q}$ be given. According to Lemma 2.5, there exists $\mathbf{w}_{\varepsilon}^q \in \mathbf{H}_0^1(\Omega)$ and a strictly positive constant independent of $N$ and $\varepsilon$, denoted $\beta_1$, such that

$$
\operatorname{div} \mathbf{w}_{\varepsilon}^{q} = q \quad \text{in} \quad \Omega, \quad \mathbf{w}_{\varepsilon}^{q} = \frac{1}{|\partial \omega_{\varepsilon}|} \left( \int_{\omega_{\varepsilon}} q \right) n_{\varepsilon} \quad \text{on} \quad \partial \omega_{\varepsilon}, \quad \| \mathbf{w}_{\varepsilon}^{q} \|_{1, \Omega} \le \beta_1 \| q \|_{\Omega}.
$$

We deduce that

$$
\|q\|_{\Omega} = \frac{(q, \operatorname{div} \mathbf{w}_{\varepsilon}^{q})_{\Omega}}{\|q\|_{\Omega}} \le \beta_1^{-1} \frac{(q, \operatorname{div} \mathbf{w}_{\varepsilon}^{q})_{\Omega}}{\|\mathbf{w}_{\varepsilon}^{q}\|_{1, \Omega}}.
$$

Due to the definitions of $F_N$ and of $\mathbf{w}_{\varepsilon}^q$ we have $(\boldsymbol{\mu}_N, \mathbf{w}_{\varepsilon}^q)_{\varepsilon} = \frac{1}{|\partial \omega_{\varepsilon}|} \left( \int_{\omega_{\varepsilon}} q \right) (\boldsymbol{\mu}_N, n_{\varepsilon})_{\varepsilon} = 0$, for any $\boldsymbol{\mu}_N \in \mathbf{F}_N$ and thus $\Pi_N^\varepsilon(\mathbf{w}_{\varepsilon}^q) = 0$. Consequently, we deduce

$$
\|q\|_{\Omega} \le \beta_1^{-1} \sup_{\mathbf{v} \in \mathbf{H}_0^1(\Omega), \Pi_N^\varepsilon(\mathbf{v})=0} \frac{b_1(q, \mathbf{v})}{\|\mathbf{v}\|_{1,\Omega}}.
$$

Let $\boldsymbol{\mu}_N \in \mathbf{F}_N$ be given. For any $\eta \in \mathbf{H}_0^{\frac{1}{2}}(\partial \omega_{\varepsilon})$, $\Pi_N^{\varepsilon}(\eta) \in \mathbf{H}_0^{\frac{1}{2}}(\partial \omega_{\varepsilon})$, so according to Lemma 2.4, there exists $\mathbf{v}_{\varepsilon}^{\eta} \in \mathbf{H}_0^1(\Omega)$ and a strictly positive constant independent of $N$ and $\varepsilon$, denoted $\beta_2$, such that

$$
\mathbf{v}_{\varepsilon}^{\eta} = \Pi_{N}^{\varepsilon}(\eta) \quad \text{on} \quad \partial \omega_{\varepsilon}, \quad \operatorname{div} \mathbf{v}_{\varepsilon}^{\eta} = 0 \quad \text{in} \quad \Omega, \quad \|\mathbf{v}_{\varepsilon}^{\eta}\|_{1,\Omega} \lesssim \|\Pi_{N}^{\varepsilon}(\eta)\|_{\frac{1}{2},\varepsilon} \le \beta_2 \|\eta\|_{\frac{1}{2},\varepsilon}.
$$
We deduce that, by definition of $\Pi_N^\epsilon$ and $v_\epsilon^\eta$,

$$ \frac{(\mu_N, \eta)_{-\frac{1}{2}, \epsilon}}{\|\eta\|_{\frac{1}{2}, \epsilon}} = \frac{(\mu_N, v_\epsilon^\eta)_\epsilon}{\|\eta\|_{\frac{1}{2}, \epsilon}} \le \beta_2^{-1} \frac{(\mu_N, v_\epsilon^\eta)_\epsilon}{\|v_\epsilon^\eta\|_{1, \Omega}} \le \beta_2^{-1} \sup_{v \in H_0^1(\Omega), \operatorname{div} v = 0} \frac{(\mu_N, v)_\epsilon}{\|v\|_{1, \Omega}}. $$

We then take the supremum over $\eta \in H^{\frac{1}{2}}(\partial\omega_\epsilon)$ of the left hand side, which yields

$$ \|\mu_N\|_{-\frac{1}{2}, \epsilon} \le \beta_2^{-1} \sup_{v \in H_0^1(\Omega), \operatorname{div} v = 0} \frac{(\mu_N, v)_\epsilon}{\|v\|_{1, \Omega}}. $$

We are now in position to apply Theorem 2.6 and conclude that Problem (2.3) is well-posed and that the solution $(u_\epsilon^{(N)}, p_\epsilon^{(N)}, \lambda_\epsilon^{(N)})$ satisfies the following energy bounds with constants independent of $N$ and $\epsilon$

$$ \begin{aligned} \|\mathbf{u}_\epsilon^{(N)}\|_{1, \Omega} &\lesssim \beta_2^{-1} (1 + \alpha^{-1} \|a\|)(|\mathbf{v}_b| + \epsilon|w_b|) \lesssim |\mathbf{v}_b| + \epsilon|w_b| \\ \|p_\epsilon^{(N)}\|_{\Omega} &\lesssim \beta_1^{-1} \|a\| \|\mathbf{u}_\epsilon^{(N)}\|_{1, \Omega} \lesssim (|\mathbf{v}_b| + \epsilon|w_b|), \\ \|\lambda_\epsilon^{(N)}\|_{-\frac{1}{2}, \epsilon} &\lesssim \beta_2^{-1} \|a\| \|\mathbf{u}_\epsilon^{(N)}\|_{1, \Omega} \lesssim |\mathbf{v}_b| + \epsilon|w_b|, \end{aligned} $$

which can be written equivalently as (2.10). $\square$

## 2.2 Model error analysis: convergence in $\epsilon$, $N$

In this section, we analyse the asymptotic behavior of the solution of Problem (1.17) as $\epsilon \to 0$. In the following we denote by $\epsilon_m$ the upper bound of parameter such that $\omega_\epsilon \subset \Omega$. The convergence result for Problem (2.1) is stated in next theorem.

**Theorem 2.7.** Let $(u_\epsilon, p_\epsilon, \lambda_\epsilon) \in V \times Q \times H_0^{-\frac{1}{2}}(\partial\omega_\epsilon)$ and $(u_\epsilon^{(N)}, p_\epsilon^{(N)}, \lambda_\epsilon^{(N)}) \in V \times Q \times F_N$ be respectively the solutions of (2.1) and (2.3). There exists $0 < \rho < \epsilon_m$ such that for $\epsilon \in (0, \rho)$,

$$ \|\mathbf{u}_\epsilon^{(N)} - \mathbf{u}_\epsilon\|_{1, \Omega} + \|\boldsymbol{\lambda}_\epsilon^{(N)} - \boldsymbol{\varepsilon}\boldsymbol{\lambda}_\epsilon\|_{-\frac{1}{2}, \partial\omega_\epsilon} + \|p_\epsilon^{(N)} - p_\epsilon\|_{\Omega} \lesssim (1 + N) \left(\frac{\varepsilon}{\rho}\right)^{N+1} (|\mathbf{v}_b| + \varepsilon|w_b|). \quad (2.22) $$

This theorem shows that the modeling error behaves as $N\epsilon^{N+1}$ for the velocity and pressure and as $N\epsilon^N$ for the Lagrange multiplier. Moreover, this estimate degenerates when $\epsilon$ is close to $\rho$ and thus it prevents to consider particles too close of the outer boundary of the domain.

To prove Theorem 2.7, we introduce two preliminary lemmas. The first lemma is a standard result which gives an expression of $\lambda_\epsilon^{(N)}$ in terms of the stress jump across $\partial\omega_\epsilon$.

**Lemma 2.8.** Let $(u_\epsilon^{(N)}, p_\epsilon^{(N)}, \lambda_\epsilon^{(N)}) \in V \times Q \times F_N$ be the solution of Problem (2.16), we have in $H^{-\frac{1}{2}}(\partial\omega_\epsilon)$

$$ \boldsymbol{\lambda}_\epsilon^{(N)} = \varepsilon \left( \gamma_{n+}^\epsilon (\boldsymbol{\sigma}(\mathbf{u}_\epsilon^{(N)}, p_\epsilon^{(N)})) + \gamma_{n-}^\epsilon (\boldsymbol{\sigma}(\mathbf{u}_\epsilon^{(N)}, p_\epsilon^{(N)})) \right), \quad (2.23) $$

where $\gamma_{n+}^\epsilon$ and $\gamma_{n-}^\epsilon$ are the normal traces on $\partial\omega_\epsilon$ of functions in $\Omega_\epsilon$ and $\omega_\epsilon$, respectively.

The proof of this lemma is given in Section A.4 of the Appendix. The second lemma describes the behaviour of the Fourier modes of $u_\epsilon^{(N)}$ on $\partial\omega_\epsilon$.

**Lemma 2.9.** Let $u_\epsilon^{(N)} \in V$ be the solution of Problem (2.16) and let $u_b = v_b + \varepsilon w_b e_\theta$ be a given rigid velocity at the boundary, with $v_b \in \mathbb{R}^3$ and $w_b \in \mathbb{R}$. There exist $a_{k,r}^\epsilon, a_{k,\theta}^\epsilon, b_{k,r}^\epsilon, b_{k,\theta}^\epsilon \in \mathbb{R}$ and $\rho \in (0, \epsilon_m)$ such that for $\epsilon \in (0, \rho)$, the error at the boundary $\partial\omega_\epsilon$ writes

$$ \begin{aligned} (u_\epsilon^{(N)}(\varepsilon)(\theta) - u_b(\varepsilon)(\theta)) \cdot e_r &= \left(\frac{\varepsilon}{\rho}\right)^{N+1} \sum_{k=N}^{\infty} (a_{k,r}^\varepsilon \cos(k\theta) + b_{k,r}^\varepsilon \sin(k\theta)), \\ (u_\epsilon^{(N)}(\varepsilon)(\theta) - u_b(\varepsilon)(\theta)) \cdot e_\theta &= \left(\frac{\varepsilon}{\rho}\right)^{N+1} \sum_{k=N}^{\infty} (a_{k,\theta}^\varepsilon \cos(k\theta) + b_{k,\theta}^\varepsilon \sin(k\theta)), \end{aligned} \quad (2.24) $$

Moreover, for $\epsilon \in (0, \rho)$, the following estimate holds:

$$ \left( \sum_{k=N}^{\infty} (1+k)(|a_{k,r}^{\epsilon}|^2 + |a_{k,\theta}^{\epsilon}|^2 + |b_{k,r}^{\epsilon}|^2 + |b_{k,\theta}^{\epsilon}|^2) \right)^{\frac{1}{2}} \lesssim (1+N)(|v_b| + \epsilon|w_b|). $$

For the sake of clarity, the proof of Lemma 2.9 which relies on development of the Stokes solution is postponed in Section A.4 of the Appendix.
Proof of Theorem 2.7. We consider a fixed $\varepsilon \in (0, \varepsilon_m)$. Let us note that $(\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon)|_{\partial\omega_\varepsilon} = \mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_b$ which belongs to $H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)$. Moreover $\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon$ satisfies the Stokes equations in $\Omega_\varepsilon$ and $\omega_\varepsilon$. Therefore, Lemma 2.4 can be applied and it gives

$$ \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon\|_{1,\Omega} \lesssim \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_b\|_{\frac{1}{2},\varepsilon}. $$

The equivalence of the norms $\| \cdot \|_{\frac{1}{2},\varepsilon}$ and $\| \cdot \|_{\frac{\mathcal{F}}{2},\varepsilon}$, uniformly in $\varepsilon$, then implies

$$ \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon\|_{1,\Omega} \lesssim \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_b\|_{\frac{\mathcal{F}}{2},\varepsilon}. $$

Thus, using Lemma 2.9 and the definition of the norm $\| \cdot \|_{\frac{\mathcal{F}}{2},\varepsilon}$, we have that there exists $\rho > 0$ such that for $\varepsilon \in (0, \rho)$,

$$ \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon\|_{1,\Omega} \lesssim (1+N) \left(\frac{\varepsilon}{\rho}\right)^{N+1} (|v_b| + \varepsilon|w_b|). \quad (2.25) $$

To obtain an estimate for $\|\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon\|_{-\frac{1}{2},\varepsilon}$, by subtracting the variational formulations (2.16) and (2.1), and from (2.8), we deduce that $(\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon, p_\varepsilon^{(N)} - p_\varepsilon, \lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon)$ satisfies

$$ 2(\varepsilon(\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon), \varepsilon(\mathbf{v}))_{\Omega} - (p_\varepsilon^{(N)} - p_\varepsilon, \operatorname{div} \mathbf{v})_{\Omega} - \langle \lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon, \mathbf{v} \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon} = 0 \quad \forall \mathbf{v} \in H_0^1(\Omega). \quad (2.26) $$

By definition (2.8) of the $\| \cdot \|_{-\frac{1}{2},\varepsilon}$ norm and since $\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon \in H_0^{-\frac{1}{2}}(\partial\omega_\varepsilon)$

$$ \|\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon\|_{-\frac{1}{2},\varepsilon} = \sup_{\eta \in H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)} \frac{\langle\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon, \eta\rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}}{\|\eta\|_{\frac{1}{2},\varepsilon}}. $$

For any $\eta \in H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)$, we consider $\mathbf{v}_\eta \in H_0^1(\Omega)$ defined in Lemma 2.4. We then obtain

$$ \frac{\langle \lambda_\varepsilon^{(N)} - \varepsilon \lambda_\varepsilon, \eta \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}}{\|\eta\|_{\frac{1}{2},\frac{1}{2},\varepsilon}} \lesssim \frac{\langle \lambda_\varepsilon^{(N)} - \varepsilon \lambda_\varepsilon, \mathbf{v}_\eta \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}}{\|\mathbf{v}_\eta\|_{1,\Omega}} \lesssim \sup_{\mathbf{v} \in \mathcal{V}, \operatorname{div} \mathbf{v} = 0} \frac{\langle \lambda_\varepsilon^{(N)} - \varepsilon \lambda_\varepsilon, \mathbf{v} \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}}{\|\mathbf{v}\|_{1,\Omega}}. $$

From (2.26) we get

$$ \frac{\langle \lambda_\varepsilon^{(N)} - \varepsilon \lambda_\varepsilon, \eta \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}}{\|\eta\|_{\frac{1}{2},\varepsilon}} \lesssim \sup_{\mathbf{v} \in \mathcal{V}, \operatorname{div} \mathbf{v} = 0} \frac{2(\varepsilon(\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon), \varepsilon(\mathbf{v}))_{\Omega}}{\|\mathbf{v}\|_{1,\Omega}} $$

Taking the supremum over $\eta \in H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)$ of the left hand side leads to

$$ \|\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon\|_{-\frac{1}{2},\varepsilon} \lesssim \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon\|_{1,\Omega}, $$

Finally, taking $\mathbf{v}_\varepsilon^{p_\varepsilon - p_\varepsilon^{(N)}}$ defined in Lemma 2.5 and associated to $p_\varepsilon - p_\varepsilon^{(N)}$, enables to bound the pressure error as

$$ \|p_\varepsilon^{(N)} - p_\varepsilon\|_{\Omega} \lesssim \|\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_\varepsilon\|_{1,\Omega}. $$

We conclude thanks to (2.25) that (2.22) holds. $\square$

**Remark 2.10.** Note that, in the fictitious domain method, as expressed in Lemma 2.8 for the solution of Problem (2.3), the force and moment transmitted to the structure include a spurious internal fluid stress. However, for the solution of (2.1) since $(\mathbf{u}_\varepsilon, p_\varepsilon)$ is solution of the Stokes equation in $\omega_\varepsilon$ with rigid body motion on $\partial\omega_\varepsilon$, we have

$$ \gamma_{n-}^{\varepsilon}(\sigma(\mathbf{u}_{\varepsilon}, p_{\varepsilon})) = \gamma_{n-}^{\varepsilon} - p_{\varepsilon} $$

with $p_\varepsilon$ constant. Since it is a constant it does not contribute to the forcing terms and thus there is no spurious force contributions. In particular for any $\delta\mathbf{r} \in \mathbb{R}^2$,

$$ \begin{aligned} |\langle \varepsilon^{-1} \lambda_\varepsilon^{(N)} - \sigma(\mathbf{u}_\varepsilon, p_\varepsilon) \mathbf{n}^+, \delta\mathbf{r} \rangle_{-\frac{1}{2},\frac{1}{2},\partial\omega_\varepsilon}| &= |\langle \lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon, \delta\mathbf{r} \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}| \\ &\lesssim \|\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon\|_{-\frac{1}{2},\varepsilon} \|\delta\mathbf{r}\|_{\frac{1}{2},\varepsilon} \lesssim (1+N) \left(\frac{\varepsilon}{\rho}\right)^{N+1} (|v_b| + \varepsilon|w_b|) \|\delta\mathbf{r}|, \end{aligned} $$

and for any $\delta\theta \in \mathbb{R}$,

$$ \begin{aligned} |\langle \varepsilon^{-1} \lambda_\varepsilon^{(N)} - \sigma(\mathbf{u}_\varepsilon, p_\varepsilon)^+ \mathbf{n}^+, \delta\theta e_\theta \rangle_{-\frac{1}{2},\frac{1}{2},\partial\omega_\varepsilon}| &= |\langle \lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon, \delta\theta e_\theta \rangle_{-\frac{1}{2},\frac{1}{2},\varepsilon}| \\ &\lesssim \|\lambda_\varepsilon^{(N)} - \varepsilon\lambda_\varepsilon\|_{-\frac{1}{2},\varepsilon} \|\delta\theta e_\theta\|_{\frac{1}{2},\varepsilon} \lesssim (1+N) \left(\frac{\varepsilon}{\rho}\right)^{N+1} (|v_b| + \varepsilon|w_b|) \|\delta\theta|, \end{aligned} $$

such that the resultant force and moment applied on $\partial\omega_\varepsilon$ are well approximated by $\varepsilon^{-1}\lambda_\varepsilon^{(N)}$ as $\varepsilon \to 0$.
## 2.3 Finite element approximation

This section is devoted to the numerical approximation of (2.3).

In this section, we propose and analyze the stability of a low-order stabilized finite element method for discretization of Problem 2.3. We recall here the definition of the standard spaces of continuous piecewise affine functions, for a triangularization $\mathcal{T}_h(\Omega)$ of the domain $\Omega$,

$$ 
\begin{aligned}
\mathbf{X}_h^1(\Omega) &= \{v \in \mathbf{H}^1(\Omega) \mid v|_K \in [\mathbb{P}_1(K)]^3, \quad \forall K \in \mathcal{T}_h(\Omega)\}, \\
\mathbf{X}_h^1(\Omega) &= \{v \in \mathbf{H}^1(\Omega) \mid v|_K \in \mathbb{P}_1(K), \quad \forall K \in \mathcal{T}_h(\Omega)\}.
\end{aligned}
$$

We then introduce the discrete spaces $V_h$ and $Q_h$ for the approximation of the fluid velocity $\mathbf{u}_\varepsilon^{(N)}$ and the pressure $p_\varepsilon^{(N)}$ as follows:

$$ V_h = \mathbf{X}_h^1(\Omega) \cap V, \quad Q_h = \mathbf{X}_h^1(\Omega) \cap Q. $$

The discrete Lagrange multiplier space is still the finite dimensional Fourier space $F_N$. In order to overcome the instability of the saddle point problem due to the choice of the Lagrange multiplier space $F_N$, we introduce the following stabilization term (see also [31]):

$$ s_h^\lambda(\lambda_N, \mu_N) = \frac{\gamma_\lambda h}{\varepsilon} (\mu_N, \lambda_N)_\varepsilon \quad \forall (\lambda_N, \mu_N) \in F_N \times F_N, $$

where $\gamma_\lambda > 0$ is a user defined dimensionless parameter that does not depend on $h$, $N$ and $\varepsilon$. Furthermore, the following Brezzi-Pitkaranta stabilization is considered to overcome the lack of inf-sup stability between the velocity and pressure spaces $V_h$ and $Q_h$:

$$ s_h^p(p_h, q_h) = \gamma_p h^2 (\nabla p_h, \nabla q_h)_\Omega \quad \forall (p_h, q_h) \in Q_h \times Q_h, $$

where $\gamma_p > 0$ is a user-defined dimensionless parameter independent of $h$, $N$ and $\varepsilon$. The proposed stabilized finite element approximation of Problem (2.3) reads as follows: Find $(\mathbf{u}_{\varepsilon,h}^{(N)}, p_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon,h}^{(N)}) \in V_h \times Q_h \times F_N$ such that

$$ 
\begin{aligned}
& a(\mathbf{u}_{\varepsilon,h}^{(N)}, \mathbf{v}_h) - b_1(p_{\varepsilon,h}^{(N)}, \mathbf{v}_h) - b_2(\lambda_{\varepsilon,h}^{(N)}, \mathbf{v}_h) = 0 \quad \forall \mathbf{v}_h \in V_h, \\
& b_1(q_h, \mathbf{u}_{\varepsilon,h}^{(N)}) + s_h^p(p_{\varepsilon,h}^{(N)}, q_h) = 0 \quad \forall q_h \in Q_h, \\
& b_2(\mu_N, \mathbf{u}_{\varepsilon,h}^{(N)}) + s_h^\lambda(\lambda_{\varepsilon,h}^{(N)}, \mu_N) = (\mu_N, \mathbf{v}_b + w_b \mathbf{e}_\theta)_\varepsilon \quad \forall \mu_N \in F_N,
\end{aligned}
\quad (2.27) $$

where the bilinear forms $a$, $b_1$ and $b_2$ are defined in (2.17).

To analyze the stability of (2.27), we introduce the following (augmented) bilinear form

$$ B_h : (V_h \times Q_h \times F_N) \times (V_h \times Q_h \times F_N) \to \mathbb{R} $$

defined as

$$ 
\begin{aligned}
B_h((\mathbf{u}_h, p_h, \lambda_N), (\mathbf{v}_h, q_h, \mu_N)) ={}& a(\mathbf{u}_h, \mathbf{v}_h) - b_1(\lambda_N, \mathbf{v}_h) - b_2(p_h, \mathbf{v}_h) \\
& + b_1(\mu_N, \mathbf{u}_h) + b_2(q_h, \mathbf{u}_h) + s_h^\lambda(\lambda_N, \mu_N) + s_h^p(p_h, q_h)
\end{aligned}
\quad (2.28) $$

for any $(\mathbf{u}_h, p_h, \lambda_N), (\mathbf{v}_h, q_h, \mu_N) \in V_h \times Q_h \times F_N$.

Problem (2.27) can then be equivalently reformulated as: Find $(\mathbf{u}_{\varepsilon,h}^{(N)}, p_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon,h}^{(N)}) \in V_h \times Q_h \times F_N$ such that

$$ B_h((\mathbf{u}_{\varepsilon,h}^{(N)}, p_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon,h}^{(N)}), (\mathbf{v}_h, q_h, \mu_N)) = (\mu_N, \mathbf{v}_b + \varepsilon w_b \mathbf{e}_\theta)_\varepsilon \quad (2.29) $$

for any $(\mathbf{v}_h, q_h, \mu_N) \in V_h \times Q_h \times F_N$.

## 2.4 Numerical stability

Before stating the main stability result, we first need to introduce the Scott-Zhang operator $\Pi_h^{\text{SZ}}: V \to V_h$, whose definition can be found in [34]. This operator is well defined for functions in $H^1$. It satisfies the following approximation properties.

**Proposition 2.11.** The Scott-Zhang operator is continuous and consistent, that is

$$ \|\Pi_h^{\text{SZ}} v\|_{1,\Omega} \lesssim \|v\|_{1,\Omega} \quad \forall v \in \mathbf{H}_0^1(\Omega). $$

Let $s \ge 1$, for $K \in \mathcal{T}_h$, let us denote by $S_K$ the domain of the elements neighboring the considered element $K$. Let $v \in \mathbf{H}^s(S_K)$, then

$$ \|\Pi_h^{\text{SZ}} v - v\|_K + h \|\nabla(\Pi_h^{\text{SZ}} v - v)\|_K \lesssim h^{\min\{2,s\}} |v|_{s,S_K}. \quad (2.30) $$
Moreover, we have

$$ \|\Pi_h^{SZ} v - v\|_{\Omega}^2 + h\|v - \Pi_h^{SZ} v\|_{1,\Omega} \lesssim h^{\min\{2,s\}}|v|_{s,\Omega} \quad \forall v \in H_0^s(\Omega). \quad (2.31) $$

The minimum regularity needed for the Scott-Zhang operator introduced in Proposition 2.11 is $H^1$. However, the reduced order Dirichlet boundary conditions on $\partial\omega_\varepsilon$ do not prevent possible discontinuities of the pressure or the velocity gradient respectively. More precisely, we cannot reach the traditional $H^1(\Omega)$ and $H^2(\Omega)$ regularity for the pressure and velocity. In particular, it is a standard result that the maximum regularity for $u_\varepsilon$ and $p_\varepsilon$ on the whole domain $\Omega$ is $H^{\frac{3}{2}-\delta}$ and $H^{\frac{1}{2}-\delta}$, respectively, with $\delta > 0$. Consequently, we also need quasi-interpolation operators designed for the approximation of non smooth functions. Let $\tilde{\Pi}_h : H^{\frac{1}{2}-\delta}(\Omega) \to X_h^1(\Omega)$ be a given quasi-interpolation operator, satisfying for $q \in H^{\frac{1}{2}-\delta}(\Omega)$,

$$ \|q - \tilde{\Pi}_h q\|_{\Omega} \lesssim h^{\frac{1}{2}-\delta} \|q\|_{\frac{1}{2}-\delta, \Omega}, \quad h\|\nabla(\tilde{\Pi}_h q)\|_{\Omega} \lesssim h^{\frac{1}{2}-\delta} \|q\|_{\frac{1}{2}-\delta, \Omega}. \quad (2.32) $$

The first estimate is a consequence of the approximation properties of quasi-interpolators (see, e.g., [35, Theorem 22.6]) while the second is a consequence of an inverse inequality (see, e.g., [36, Lemma 4.3]). The extension of $\tilde{\Pi}_h$ to a quasi-interpolator $\tilde{\Pi}_h : Q \to Q_h$ satisfying (2.32) is then straightforward since $Q_h \subset Q \subset H^1(\Omega)$.

We finally state two preliminary lemmas useful to prove the inf-sup condition of the discrete problem (2.27). These lemmas exhibit test functions that will enable to deal with the constraints thanks to the stabilization terms.

**Lemma 2.12.** Let $\lambda_N \in F_N$ be given. There exists $v_{\varepsilon,h}^{\lambda_N} \in V_h$ such that

$$ \frac{3}{4} \|\lambda_N\|_{-\frac{1}{2},\varepsilon}^2 \le -(\lambda_N, v_{\varepsilon,h}^{\lambda_N})_\varepsilon + C_\varepsilon \frac{h}{\varepsilon} \|\lambda_N\|_\varepsilon^2, \quad (2.33) $$

$$ \|v_{\varepsilon,h}^{\lambda_N}\|_{1,\Omega} \lesssim \|\lambda_N\|_{-\frac{1}{2},\varepsilon}, \quad (2.34) $$

where C is a strictly positive constant that does not depend on $h$, $N$, $\varepsilon$. Moreover, it holds

$$ |(q_h, \operatorname{div} v_{\varepsilon,h}^{\lambda_N})_{\Omega}| \lesssim h \|\nabla q_h\|_{\Omega} \|\lambda_N\|_{-\frac{1}{2},\varepsilon} \quad \forall q_h \in Q_h. \quad (2.35) $$

**Remark 2.13.** Note that inequality (2.33) bounds $b_2(\lambda_N, v_{\varepsilon,h}^{\lambda_N})$ from below but with a residual term that is the stabilization term $s_h^\lambda(\lambda_N, \lambda_N)$. Moreover (2.35) gives an estimate of $b_1(q_h, v_{\varepsilon,h}^{\lambda_N})$ for any $q_h$ and thus “measures” the interaction between the chosen $v_{\varepsilon,h}^{\lambda_N}$ associated with $\lambda_N$ and the divergence constraint.

*Proof of Lemma 2.12.* Let $\lambda_N \in F_N$. Using the definition (2.8) of the norm $\|\cdot\|_{-\frac{1}{2},\varepsilon}$, and since $F_N$ is a closed subspace of $L^2(\partial\omega_\varepsilon)$, there exists $\mu_N \in F_N$ such that

$$ \|\mu_N\|_{\frac{1}{2},\varepsilon} = 1, \quad \|\lambda_N\|_{-\frac{1}{2},\varepsilon} = (\lambda_N, \mu_N)_\varepsilon. $$

Hence, by setting $\underline{\mu}_N = -\mu_N \|\lambda_N\|_{-\frac{1}{2},\varepsilon}$, we have

$$ \|\underline{\mu}_N\|_{\frac{1}{2},\varepsilon} = \|\lambda_N\|_{-\frac{1}{2},\varepsilon}, \quad \text{and} \quad -(\lambda_N, \underline{\mu}_N)_\varepsilon = \|\lambda_N\|_{-\frac{1}{2},\varepsilon}^2. \quad (2.36) $$

On the other hand, owing to Lemma 2.4, since $\underline{\mu}_N \in F_N \subset H_0^{\frac{1}{2}}(\partial\omega_\varepsilon)$ there exists $v_\varepsilon^{\mu_N} \in V$ such that

$$ \operatorname{div} v_\varepsilon^{\mu_N} = 0, \quad v_\varepsilon^{\mu_N}|_{\partial\omega_\varepsilon} = \underline{\mu}_N, \quad \|v_\varepsilon^{\mu_N}\|_{1,\Omega} \lesssim \|\underline{\mu}_N\|_{\frac{1}{2},\varepsilon} = \|\lambda_N\|_{-\frac{1}{2},\varepsilon}, \quad (2.37) $$

where the last identity follows from (2.36). We can now set $v_{\varepsilon,h}^{\lambda_N} = \Pi_h^{SZ} v_\varepsilon^{\mu_N} \in V_h$, so that thanks to (2.36) and (2.37), we obtain

$$ -(\lambda_N, v_{\varepsilon,h}^{\lambda_N})_\varepsilon = -(\lambda_N, v_{\varepsilon,h}^{\lambda_N} - v_\varepsilon^{\mu_N})_\varepsilon - (\lambda_N, \underline{\mu}_N)_\varepsilon = -(\lambda_N, v_{\varepsilon,h}^{\lambda_N} - v_\varepsilon^{\mu_N})_\varepsilon + \|\lambda_N\|_{-\frac{1}{2},\varepsilon}^2. \quad (2.38) $$

Using (2.30) and a scaling argument (see [31, Proof of Lemma 4.2]), it follows that

$$ \|\Pi_h^{SZ} v_\varepsilon^{\mu_N} - v_\varepsilon^{\mu_N}\|_{K \cap \partial\omega_\varepsilon}^2 \lesssim h \|\nabla v_\varepsilon^{\mu_N}\|_{S_K}^2. $$
for each $K \in \mathcal{T}_h$ with $K \cap \partial\omega_\varepsilon \neq \emptyset$. Summing over all $K \in \mathcal{T}_h$ intersecting $\partial\omega_\varepsilon$ and using (2.37), we get

$$ \|\Pi_h^{\text{SZ}} v_{\varepsilon}^{\mu_N} - v_{\varepsilon}^{\mu_N}\|_{\partial\omega_\varepsilon}^2 \lesssim h \|v_{\varepsilon}^{\mu_N}\|_{1,\Omega}^2 \lesssim h \|\lambda_N\|_{-\frac{1}{2},\varepsilon}^2. \quad (2.39) $$

Then, by using Cauchy-Schwarz inequality and the scaling between the norms $\|\cdot\|_\varepsilon$ and $\|\cdot\|_{\partial\omega_\varepsilon}$ resulting from (2.9), we get

$$ |(\lambda_N, v_{\varepsilon}^{\mu_N} - v_{\varepsilon,h}^{\lambda_N})_\varepsilon| \lesssim C \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\lambda_N\|_\varepsilon \|\lambda_N\|_{-\frac{1}{2},\varepsilon} \le \frac{1}{4} \|\lambda_N\|_{-\frac{1}{2},\varepsilon}^2 + C \frac{h}{\varepsilon} \|\lambda_N\|_\varepsilon^2, $$

with $C > 0$ independent of $N$, $h$ and $\varepsilon$. By combining this estimate with (2.38), and by using the continuity of the Scott-Zhang interpolator and (2.37), we finally get the following relations

$$ \begin{aligned} \frac{3}{4} \|\lambda_N\|_{-\frac{1}{2},\varepsilon}^2 &\lesssim -(\lambda_N, v_{\varepsilon,h}^{\lambda_N})_\varepsilon + C \frac{h}{\varepsilon} \|\lambda_N\|_\varepsilon^2, \\ \|v_{\varepsilon,h}^{\lambda_N}\|_{1,\Omega} &\lesssim \|v_{\varepsilon}^{\mu_N}\|_{1,\Omega} \lesssim \|\lambda_N\|_{-\frac{1}{2},\varepsilon}. \end{aligned} \quad (2.40) $$

Finally, since div $v_{\varepsilon}^{\mu_N} = 0$ in $\Omega$, we obtain for any $q_h \in Q_h$,

$$ \begin{aligned} (q_h, \text{div } v_{\varepsilon,h}^{\lambda_N})_{\Omega} &= (q_h, \text{div}(v_{\varepsilon,h}^{\lambda_N} - v_{\varepsilon}^{\mu_N}))_{\Omega}, \\ &= (q_h, \text{div}(\Pi_h^{\text{SZ}} v_{\varepsilon}^{\mu_N} - v_{\varepsilon}^{\mu_N}))_{\Omega} = -(\nabla q_h, \Pi_h^{\text{SZ}} v_{\varepsilon}^{\mu_N} - v_{\varepsilon}^{\mu_N})_{\Omega} \\ &\lesssim h \|\nabla q_h\|_{\Omega} \|v_{\varepsilon}^{\mu_N}\|_{\Omega} \lesssim h \|\nabla q_h\|_{\Omega} \|\lambda_N\|_{-\frac{1}{2},\varepsilon}. \end{aligned} $$

This concludes the proof.

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

**Lemma 2.14.** Let $q_h \in Q_h$ be given. There exists $w_{\varepsilon,h}^{q_h} \in V_h$ such that

$$ \frac{3}{4} \|q_h\|_{\Omega}^2 \lesssim -(q_h, \text{div } w_{\varepsilon,h}^{q_h})_{\Omega} + Ch^2 \|\nabla q_h\|_{\Omega}^2, \quad (2.41) $$

$$ \|w_{\varepsilon,h}^{q_h}\|_{1,\Omega} \lesssim \|q_h\|_{\Omega}, \quad (2.42) $$

with $C > 0$ independent of $h$, $N$ and $\varepsilon$. Moreover, it holds

$$ |(\mu_N, w_{\varepsilon,h}^{q_h})_\varepsilon| \lesssim \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\mu_N\|_\varepsilon \|q_h\|_{\Omega} \quad \forall \mu_N \in \mathcal{F}_N. \quad (2.43) $$

**Remark 2.15.** Note that inequality (2.41) bounds $b_1(q_h, w_{\varepsilon,h}^{q_h})$ from below but with a residual term that corresponds to the stabilization term $s_h^p(q_h, q_h)$. Moreover (2.35) gives an estimate of $b_2(\mu_N, w_{\varepsilon,h}^{q_h})$ for any $\mu_N$ and thus “measures” the interaction between the chosen $w_{\varepsilon,h}^{q_h}$ associated with $q_h$ and the Dirichlet constraint.

**Proof of Lemma 2.14.** Let $q_h \in Q_h$ be given, since $q_h \in Q$, according to Lemma 2.4, there exists $w_{\varepsilon}^{q_h} \in V$ such that

$$ \text{div } w_{\varepsilon}^{q_h} = -q_h \quad \text{in } \Omega, \quad (2.44) $$

$$ w_{\varepsilon}^{q_h} = - \frac{1}{|\partial \omega_{\varepsilon}|} \left( \int_{\omega_{\varepsilon}} q_h \right) n^{-} \quad \text{on } \partial \omega_{\varepsilon}, \quad (2.45) $$

$$ \|w_{\varepsilon}^{q_h}\|_{1,\Omega} \lesssim \|q_h\|_{\Omega}. \quad (2.46) $$

By setting $w_{\varepsilon,h}^{q_h} = \Pi_h^{\text{SZ}} w_{\varepsilon}^{q_h} \in V_h$ and by applying (2.44), (2.46), we obtain

$$ \begin{aligned} \|q_h\|_{\Omega}^2 + (q_h, \text{div } w_{\varepsilon,h}^{q_h})_{\Omega} &= -(q_h, \text{div}(w_{\varepsilon}^{q_h} - w_{\varepsilon,h}^{q_h}))_{\Omega} = (\nabla q_h, w_{\varepsilon}^{q_h} - w_{\varepsilon,h}^{q_h})_{\Omega} \\ &\lesssim h \|\nabla q_h\|_{\Omega} \|w_{\varepsilon}^{q_h}\|_{1,\Omega} \lesssim h \|\nabla q_h\|_{\Omega} \|q_h\|_{\Omega} \\ &\le \frac{1}{4} \|q_h\|_{\Omega}^2 + Ch^2 \|\nabla q_h\|_{\Omega}^2, \end{aligned} $$

Consequently $w_{\varepsilon,h}^{q_h}$ satisfies estimate (2.41). In addition, from (2.45), we have $(\mu_N, w_{\varepsilon}^{q_h})_\varepsilon = 0$ for any $\mu_N \in \mathcal{F}_N$, since $\int_{\partial\omega_\varepsilon} \mu_N \cdot n^{-} = 0$, so that

$$ |(\mu_N, w_{\varepsilon,h}^{q_h})_\varepsilon| = |(\mu_N, w_{\varepsilon,h}^{q_h} - w_{\varepsilon}^{q_h})_\varepsilon| = |(\mu_N, \Pi_h^{\text{SZ}} w_{\varepsilon}^{q_h} - w_{\varepsilon}^{q_h})_\varepsilon| \lesssim \|\mu_N\|_\varepsilon \|\Pi_h^{\text{SZ}} w_{\varepsilon}^{q_h} - w_{\varepsilon}^{q_h}\|_\varepsilon. $$
Finally, using the same arguments as in (2.39), the scaling between the norms $\| \cdot \|_{\varepsilon}$ and $\| \cdot \|_{\partial\omega_{\varepsilon}}$ and the estimate(2.46), we get

$$ |(\boldsymbol{\mu}_N, \boldsymbol{w}_{\varepsilon,h}^{q_h})| \lesssim \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\boldsymbol{\mu}_N\|_{\varepsilon} \|\boldsymbol{w}_{\varepsilon}^{q_h}\|_{1,\Omega} \lesssim \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\boldsymbol{\mu}_N\|_{\varepsilon} \|q_h\|_{L^2(\Omega)} $$

This concludes the proof.

For the analysis of the augmented bilinear form, we introduce the following norm:

$$ \|(v_h, q_h, \boldsymbol{\mu}_N)\|_{\varepsilon} = \left( \|v_h\|_{1,\Omega}^2 + \|p_h\|_{\Omega}^2 + \|\boldsymbol{\mu}_N\|_{-\frac{1}{2},\varepsilon}^2 + \gamma_p h^2 \|\nabla p_h\|_{\Omega}^2 + \gamma_{\lambda} \frac{h}{\varepsilon} \|\boldsymbol{\mu}_N\|_{\varepsilon}^2 \right)^{\frac{1}{2}} $$

for any $(v_h, q_h, \boldsymbol{\mu}_N) \in \mathbf{V}_h \times Q_h \times \mathbf{F}_N$.

The stability of the discrete approximation can be ensured by proving an inf-sup inequality on the stabilized bilinear form $B_h$ which is addressed in the following theorem.

**Theorem 2.16.** There exists a constant $\beta > 0$ independent of $h$, $N$ and $\varepsilon$ such that

$$ \sup_{(v_h, q_h, \boldsymbol{\mu}_N) \in \mathbf{V}_h \times Q_h \times \mathbf{F}_N} \frac{B_h((\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N), (v_h, q_h, \boldsymbol{\mu}_N))}{\|(v_h, q_h, \boldsymbol{\mu}_N)\|_{\varepsilon}} \geq \beta \|(\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N)\|_{\varepsilon} \quad (2.47) $$

for any $(\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N) \in \mathbf{V}_h \times Q_h \times \mathbf{F}_N$.

*Proof of Lemma 2.16.* Let us consider $(\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N) \in \mathbf{V}_h \times Q_h \times \mathbf{F}_N$. First, by testing (2.28) with $(\boldsymbol{v}_h, q_h, \boldsymbol{\mu}_N) = (\boldsymbol{u}_h, q_h, \boldsymbol{\lambda}_N)$, we get

$$ B_h((\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N), (\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N)) = 2\mu \|\boldsymbol{\varepsilon}(\boldsymbol{u}_h)\|_{\Omega}^2 + \gamma_p h^2 \|\nabla p_h\|_{\Omega}^2 + \gamma_{\lambda} \frac{h}{\varepsilon} \|\boldsymbol{\lambda}_N\|_{\varepsilon}^2. \quad (2.48) $$

According to Lemma 2.12, there exists $\boldsymbol{v}_{\varepsilon,h}^{\lambda_N} \in \mathbf{V}_h$, and two constants $C_1, D_1 > 0$ independent of $h$, $N$ and $\varepsilon$, such that

$$ \frac{3}{4} \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon}^2 \leq - (\boldsymbol{\lambda}_N, \boldsymbol{v}_{\varepsilon,h}^{\lambda_N})_{\varepsilon} + C_1 \frac{h}{\varepsilon} \|\boldsymbol{\lambda}_N\|_{\varepsilon}^2, \quad (2.49) $$

$$ \|(p_h, \operatorname{div} \boldsymbol{v}_{\varepsilon,h}^{\lambda_N})_{\Omega}\| \le D_1 h^2 \|\nabla p_h\|_{\Omega}^2 + \frac{1}{4} \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon}^2, \quad (2.50) $$

$$ \|\boldsymbol{v}_{\varepsilon,h}^{\lambda_N}\|_{1,\Omega} \lesssim \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon}. \quad (2.51) $$

By applying the Cauchy-Schwarz and Young inequalities and using (2.51), we deduce that there exists $B_1 > 0$ independent of $h$, $N$ and $\varepsilon$, such that

$$ |2(\boldsymbol{\varepsilon}(\boldsymbol{u}_h), \boldsymbol{\varepsilon}(\boldsymbol{v}_{\varepsilon,h}^{\lambda_N}))| \le B_1 \|\boldsymbol{u}_h\|_{1,\Omega}^2 + \frac{1}{4} \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon}^2. \quad (2.52) $$

Taking $(\boldsymbol{v}_h, q_h, \boldsymbol{\mu}_N) = (\boldsymbol{v}_{\varepsilon,h}^{\lambda_N}, 0, 0)$ in (2.28) and using (2.49), (2.50) and (2.52), we get

$$ \begin{aligned} B_h((\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N), (\boldsymbol{v}_{\varepsilon,h}^{\lambda_N}, 0, 0)) &= a(\boldsymbol{u}_h, \boldsymbol{v}_{\varepsilon,h}^{\lambda_N}) - b_1(\boldsymbol{\lambda}_N, \boldsymbol{v}_{\varepsilon,h}^{\lambda_N}) - b_2(p_h, \boldsymbol{v}_{\varepsilon,h}^{\lambda_N}) \\ &\ge \frac{1}{4} \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon}^2 - B_1 \|\boldsymbol{u}_h\|_{1,\Omega}^2 - C_1 \frac{h}{\varepsilon} \|\boldsymbol{\lambda}_N\|_{\varepsilon}^2 - D_1 h^2 \|\nabla p_h\|_{\Omega}^2. \end{aligned} \quad (2.53) $$

Similarly, according to Lemma 2.14, there exists $\boldsymbol{w}_{\varepsilon,h}^{p_h} \in \mathbf{V}_h$ and $C_2, D_2 > 0$ independent of $h$, $N$ and $\varepsilon$ such that

$$ \frac{3}{4} \|p_h\|_{\Omega}^2 \leq - (p_h, \operatorname{div} \boldsymbol{w}_{\varepsilon,h}^{p_h})_{\Omega} + D_2 h^2 \|\nabla p_h\|_{\Omega}^2, \quad (2.54) $$

$$ |(\boldsymbol{\lambda}_N, \boldsymbol{w}_{\varepsilon,h}^{p_h})_{\varepsilon}| \le C_2 \frac{h}{\varepsilon} \|\boldsymbol{\lambda}_N\|_{\varepsilon}^2 + \frac{1}{4} \|p_h\|_{\Omega}^2, \quad (2.55) $$

$$ \|\boldsymbol{w}_{\varepsilon,h}^{p_h}\|_{1,\Omega} \lesssim \|p_h\|_{\Omega}. \quad (2.56) $$

From (2.56), we deduce that there exists $B_2 > 0$ independent of $h$, $N$ and $\varepsilon$ such that

$$ |2(\boldsymbol{\varepsilon}(\boldsymbol{u}_h), \boldsymbol{\varepsilon}(\boldsymbol{w}_{\varepsilon,h}^{p_h}))| \le B_2 \|\boldsymbol{u}_h\|_{1,\Omega}^2 + \frac{1}{4} \|p_h\|_{\Omega}^2. \quad (2.57) $$

Taking $(\boldsymbol{v}_h, q_h, \boldsymbol{\mu}_N) = (\boldsymbol{w}_{\varepsilon,h}^{p_h}, 0, 0)$ in (2.28) and using (2.54), (2.55) and (2.57), we get

$$ B((\boldsymbol{u}_h, p_h, \boldsymbol{\lambda}_N), (\boldsymbol{w}_{\varepsilon,h}^{p_h}, 0, 0)) \ge \frac{1}{4} \|p_h\|_{\Omega}^2 - B_2 \|\boldsymbol{u}_h\|_{1,\Omega}^2 - C_2 \frac{h}{\varepsilon} \|\boldsymbol{\lambda}_N\|_{\varepsilon}^2 - D_2 h^2 \|\nabla p_h\|_{\Omega}^2. \quad (2.58) $$
Finally, taking in (2.28) a linear combination of the previous test functions defined as

$$ (\mathbf{v}_h, q_h, \boldsymbol{\mu}_N) = (\mathbf{u}_h + \alpha_\lambda \mathbf{v}_{\varepsilon,h}^{\lambda_N} + \alpha_p \mathbf{w}_{\varepsilon,h}^{p_h}, \mathbf{p}_h, \boldsymbol{\lambda}_N), $$

with $\alpha_\lambda, \alpha_p > 0$, using Korn's inequality and collecting estimates (2.48), (2.53) and (2.58), leads to

$$ \begin{aligned} B((\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N), (\mathbf{u}_h + \alpha_\lambda \mathbf{v}_{\varepsilon,h}^{\lambda_N} + \alpha_p \mathbf{w}_{\varepsilon,h}^{p_h}, \mathbf{p}_h, \boldsymbol{\lambda}_N)) \geq & (\alpha - \alpha_\lambda B_1 - \alpha_p B_2) \|\mathbf{u}_h\|_{1,\Omega}^2 \\ & + \frac{\alpha_p}{2} \|\mathbf{p}_h\|_{\Omega}^2 + \frac{\alpha_\lambda}{4} \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon}^2 + (\gamma_\lambda - \alpha_\lambda C_1 - \alpha_p C_2) \frac{h}{\varepsilon} \|\boldsymbol{\lambda}_N\|_{\varepsilon}^2 \\ & + (\gamma_p - \alpha_\lambda D_1 - \alpha_p D_2) h^2 \|\nabla \mathbf{p}_h\|_{\Omega}^2. \end{aligned} \quad (2.59) $$

where $\alpha > 0$ denotes the coercivity constant of the bilinear form $a$, which is independent of $\varepsilon, h$ and $N$. By choosing $\alpha_p, \alpha_\lambda$ small enough, such that the following inequalities hold true

$$ \alpha - \alpha_\lambda B_1 - \alpha_p B_2 > 0, \quad \gamma_\lambda - \alpha_\lambda C_1 - \alpha_p C_2 > 0, \quad \gamma_p - \alpha_\lambda D_1 - \alpha_p D_2 > 0, $$

we get

$$ B((\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N), (\mathbf{u}_h + \alpha_\lambda \mathbf{v}_{\varepsilon,h}^{\lambda_N} + \alpha_p \mathbf{w}_{\varepsilon,h}^{p_h}, \mathbf{p}_h, \boldsymbol{\lambda}_N)) \geq C_m \|(\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N)\|_{\varepsilon}^2, \quad (2.60) $$

with $C_m > 0$ independent of $\varepsilon, h$ and $N$. On the other hand, using the triangle inequality and the a priori estimates (2.51) and (2.56), respectively, we obtain

$$ \begin{aligned} \|(\mathbf{u}_h + \alpha_\lambda \mathbf{v}_{\varepsilon,h}^{\lambda_N} + \alpha_p \mathbf{w}_{\varepsilon,h}^{p_h}, \mathbf{p}_h, \boldsymbol{\lambda}_N)\|_{\varepsilon} &\le \|(\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N)\|_{\varepsilon} + \alpha_\lambda C_3 \|\boldsymbol{\lambda}_N\|_{-\frac{1}{2},\varepsilon} + \alpha_p D_3 \|\mathbf{p}_h\|_{\Omega} \\ &\le C_M \|(\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N)\|_{\varepsilon}, \end{aligned} \quad (2.61) $$

where $C_3, D_3 > 0$ are constants independent of $h, N$ and $\varepsilon$ and

$$ C_M = 1 + \alpha_\lambda C_3 + \alpha_p D_3. $$

Finally, by combining (2.60) and (2.61) we get

$$ \frac{B_h((\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N), (\mathbf{v}_h, q_h, \boldsymbol{\mu}_N))}{\|(\mathbf{v}_h, q_h, \boldsymbol{\mu}_N)\|_{\varepsilon}} \ge \frac{C_m}{C_M} \|(\mathbf{u}_h, \mathbf{p}_h, \boldsymbol{\lambda}_N)\|_{\varepsilon}, $$

which yields (2.47) with $\beta = C_m/C_M$. This completes the proof. $\square$

## 2.5 Error analysis: convergence in $h$

To conclude the analysis of the stabilized numerical method (2.29), we provide an a priori error estimate for (2.29) as approximation of (2.3). To this purpose, we first introduce the following decomposition of the velocity and pressure errors:

$$ \begin{aligned} \mathbf{u}_{\varepsilon}^{(N)} - \mathbf{u}_{\varepsilon,h}^{(N)} &= \underbrace{\mathbf{u}_{\varepsilon}^{(N)} - \Pi_h^{SZ} \mathbf{u}_{\varepsilon}^{(N)}}_{\mathbf{e}_{\varepsilon}^{(N)}} + \underbrace{\Pi_h^{SZ} \mathbf{u}_{\varepsilon}^{(N)} - \mathbf{u}_{\varepsilon,h}^{(N)}}_{\mathbf{e}_{\varepsilon,h}^{(N)}}, \\ p_{\varepsilon}^{(N)} - p_{\varepsilon,h}^{(N)} &= \underbrace{p_{\varepsilon}^{(N)} - \bar{\Pi}_h p_{\varepsilon}^{(N)}}_{\mathbf{e}_{\varepsilon}^{(N)}} + \underbrace{\bar{\Pi}_h p_{\varepsilon}^{(N)} - p_{\varepsilon,h}^{(N)}}_{\mathbf{e}_{\varepsilon,h}^{(N)}}, \end{aligned} \quad (2.62) $$

The first step consists in providing an error estimate for the discrete counterpart of the error, viz., $(\mathbf{e}_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \boldsymbol{\lambda}_{\varepsilon}^{(N)}, \boldsymbol{\lambda}_{\varepsilon,h}^{(N)}) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_N$. Then the remaining contributions $\mathbf{e}_{\varepsilon}^{(N)}$ and $\mathbf{e}_{\varepsilon}^{(N)}$ will be simply handled via the approximation bounds (2.31) and (2.32).

To this purpose, we shall make use of the following weak consistency and continuity result.

**Lemma 2.17.** The following identity holds:

$$ B_h((\mathbf{e}_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \boldsymbol{\lambda}_{\varepsilon}^{(N)} - \boldsymbol{\lambda}_{\varepsilon,h}^{(N)}), (\mathbf{v}_h, q_h, \boldsymbol{\mu}_N)) = -a(\mathbf{e}_{\varepsilon}^{(N)}, \mathbf{v}_h) + b_1(\mathbf{e}_{\varepsilon}^{(N)}, \mathbf{v}_h) - b_1(q_h, \mathbf{e}_{\varepsilon}^{(N)}) \\ - b_2(\boldsymbol{\mu}_N, \mathbf{e}_{\varepsilon}^{(N)}) + s_h^p(q_h, \bar{\Pi}_h(p_{\varepsilon}^{(N)})) + s_h^{\lambda}(\boldsymbol{\mu}_N, \boldsymbol{\lambda}_{N,h}^{\varepsilon}) \quad (2.63) $$

for any $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_N) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_N$. Furthermore, the following bound also holds

$$ |B_h((\mathbf{e}_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \boldsymbol{\lambda}_{\varepsilon}^{(N)} - \boldsymbol{\lambda}_{\varepsilon,h}^{(N)}), (\mathbf{v}_h, q_h, \boldsymbol{\mu}_N))| \\ \lesssim \left( (1 + |\log(\varepsilon)|^{\frac{1}{2}}) \|\mathbf{e}_{\varepsilon}^{(N)}\|_{1,\Omega} + \|\mathbf{e}_{\varepsilon}^{(N)}\|_{\Omega} + h^{-1} \|\mathbf{e}_{\varepsilon}^{(N)}\|_{\Omega} + h \|\nabla \bar{\Pi}_h p_{\varepsilon}^{(N)}\|_{\Omega} + \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\boldsymbol{\lambda}_{\varepsilon}^{(N)}\|_{\varepsilon} \right) \|(\mathbf{v}_h, q_h, \boldsymbol{\mu}_N)\|_{\varepsilon} \quad (2.64) $$

for any $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_N) \in \mathbf{V}_h \times \mathbf{Q}_h \times \mathbf{F}_N$.
Proof. The relation (2.63) follows directly follows from (2.16), (2.27), (2.28) and (2.62). Finally, the continuity bound (2.64) follows from the continuity of the bilinear forms $a, b_1, b_2, s_h^\lambda$ and $s_h^p$. It is the continuity constant of $b_2$ (and Lemma 2.3) that causes a factor $|\log(\varepsilon)|^{1/2}$ to appear. $\square$

The main result of this section is stated in the next lemma.

**Lemma 2.18.** For $N \in \mathbb{N}$, let $(\mathbf{u}_\varepsilon^{(N)}, p_\varepsilon^{(N)}, \lambda_\varepsilon^{(N)}) \in \mathbf{V} \times \mathcal{Q} \times \mathbf{F}_N$ be solution of (2.3) and $(\mathbf{u}_{\varepsilon,h}^{(N)}, p_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon,h}^{(N)}) \in \mathbf{V}_h \times \mathcal{Q}_h \times \mathbf{F}_N$ be solution of (2.29). Assuming that $(\mathbf{u}_{\varepsilon,h}^{(N)}, p_{\varepsilon,h}^{(N)})$ belongs to $\mathbf{H}^{\frac{3}{2}-\delta}(\Omega) \times \mathbf{H}^{\frac{3}{2}-\delta}(\Omega)$ for $\delta > 0$ and for any $N$, there holds

$$ \| (\mathbf{u}_\varepsilon^{(N)} - \mathbf{u}_{\varepsilon,h}^{(N)}, p_\varepsilon^{(N)} - p_{\varepsilon,h}^{(N)}, \lambda_\varepsilon^{(N)} - \lambda_{\varepsilon,h}^{(N)}) \|_{\varepsilon} \lesssim (1 + |\log(\varepsilon)|^{\frac{1}{2}}) h^{\frac{1}{2}-\delta} \|\mathbf{u}_\varepsilon^{(N)}\|_{\frac{3}{2}-\delta, \Omega} \\ + h^{\frac{1}{2}-\delta} \|p_\varepsilon^{(N)}\|_{\frac{1}{2}-\delta, \Omega} + \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\lambda_\varepsilon^{(N)}\|_{\varepsilon}. \quad (2.65) $$

Proof. We start by applying Lemma 2.16 to $(e_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon,h}^{(N)} - \lambda_{\varepsilon,h}^{(N)}) \in \mathbf{V}_h \times \mathcal{Q}_h \times \mathbf{F}_N$, so that

$$ \| (e_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon}^{(N)} - \lambda_{\varepsilon,h}^{(N)}) \|_{\varepsilon} \le \frac{1}{\beta} \sup_{(v_h, q_h, \mu_N) \in \mathbf{V}_h \times \mathcal{Q}_h \times \mathbf{F}_N} \frac{B_h((e_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon}^{(N)} - \lambda_{\varepsilon,h}^{(N)}), (v_h, q_h, \mu_N))}{\|(v_h, q_h, \mu_N)\|_{\varepsilon}}. $$

Owing to (2.64), the previous estimate together with (2.32) yields

$$ \begin{aligned} \| (e_{\varepsilon,h}^{(N)}, e_{\varepsilon,h}^{(N)}, \lambda_{\varepsilon}^{(N)} - \lambda_{\varepsilon,h}^{(N)}) \|_{\varepsilon} & \lesssim (1 + |\log(\varepsilon)|^{\frac{1}{2}}) \|e_{\varepsilon}^{(N)}\|_{1,\Omega} + \|e_{\varepsilon}^{(N)}\|_{\Omega} + h^{-1} \|e_{\varepsilon}^{(N)}\|_{\Omega} + h \|\nabla \bar{\Pi}_h p_{\varepsilon}^{(N)}\|_{\Omega} \\ & \qquad + \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\lambda_{\varepsilon}^{(N)}\|_{\varepsilon} \\ & \lesssim (1 + |\log(\varepsilon)|^{\frac{1}{2}}) h^{\frac{1}{2}-\delta} \|\mathbf{u}_{\varepsilon}^{(N)}\|_{\frac{3}{2}-\delta, \Omega} + h^{\frac{1}{2}-\delta} \|p_{\varepsilon}^{(N)}\|_{\frac{1}{2}-\delta, \Omega} + \left(\frac{h}{\varepsilon}\right)^{\frac{1}{2}} \|\lambda_{\varepsilon}^{(N)}\|_{\varepsilon}. \end{aligned} \quad (2.66) $$

Estimate (2.65) then follows from the error decomposition (2.62), by using triangle's inequality and the interpolation estimates (2.31) and (2.32). This completes the proof. $\square$

These estimates show that the stabilized method is stable and convergent. Thanks to the stabilization term $s_h^\lambda$ the inf-sup constant does depend on $(1 + h/\varepsilon)^N$ (see [31] in the Poisson case) avoiding locking effect as $N$ is large. However, it suffers from the standard drawback of the fictitious domain method, for which the optimal convergence is not reached due to the lack of regularity of the solution in the whole domain $\Omega$. Moreover, we see that the convergence still depends on the ratio $h/\varepsilon$: explicitly in front of the norm of the Lagrange multiplier $\|\lambda_\varepsilon^{(N)}\|_\varepsilon$ and implicitly since $\|\mathbf{u}_\varepsilon^{(N)}\|_{\frac{3}{2}-\delta, \Omega}$ and $\|p_\varepsilon^{(N)}\|_{\frac{1}{2}-\delta, \Omega}$ may depend on $\varepsilon$. Yet observe that the rate of convergence depends only on $h$ and $\varepsilon$ and behaves in the same way for any $N$. When the size of the particle goes to zero one needs to refine the mesh, but there is no a priori impact with respect to the choice of the number of modes $N$. To recover optimal convergence rate one could enrich the discretization space with well adapted basis functions capturing the jumps as it was proposed and studied in [31] for the Poisson problem. However, note that, contrary to the Poisson case, these augmented basis functions may interact with the finite element basis due to the divergence constraint leading to additional difficulties in the stability and convergence analysis.

# 3 Numerical methods for fluid-particle interaction

In this section, we address the numerical approximation of the full reduced-order coupled model (1.17). Alternative fully discrete formulations with standard coupling conditions on the fluid-structure interface are also discussed.

## 3.1 Discrete formulation of the reduced coupled problem

We introduce in this section the time and space discretizations of the reduced problem (1.17). We denote by $\tau > 0$ the time step length and for $n \in \mathbb{N}$, we set $t^n = n\tau$. For a given function $v$ depending on $t$, we set

$$ v^n, \quad \partial_\tau v^n = \frac{1}{\tau}(v^n - v^{n-1}), \quad v^{n-\frac{1}{2}} = \frac{1}{2}(v^{n-1} + v^n) $$
for, respectively, an approximation of $v(t^n)$, the backward Euler difference and the mid point evaluation. We then consider a semi-implicit scheme in which the geometrical coupling is treated in an explicit way, as follows: knowing the approximation of the particle translation at time $(n-1)\tau$, $\mathbf{r}_\epsilon^{n-1} \in \mathbb{R}^2$ and the approximation of the particle rotation at time $(n-1)\tau$, $\theta_\epsilon^{n-1} \in \mathbb{R}$, it reads

$$ \phi^n(\hat{\mathbf{x}}) = \mathbf{r}_\epsilon^{n-1} + \mathbf{\Lambda}_\epsilon^{n-1} \hat{\mathbf{x}} \quad \forall \hat{\mathbf{x}} \in \hat{\omega}_\epsilon. $$

where $\mathbf{\Lambda}_\epsilon^{n-1}$ is obtained in terms of $\theta_\epsilon^{n-1}$ by (1.4). The approximated domain of the particle at time $n\tau$ is then defined by $\omega_\epsilon^n = \phi^n(\hat{\omega}_\epsilon)$. For the space discretization, similarly to Section 2.3, we consider a low-order stabilized formulation based on piecewise affine approximations. In particular, we set $\mathbf{V}_h = \mathbf{X}_h^1(\Omega) \cap \mathbf{V}$ and $Q_h = X_h^1(\Omega) \cap Q$ for the discretization of the velocity $\mathbf{u}_\epsilon^{(N)}$ and pressure $p_\epsilon^{(N)}$. We also introduce the discrete Lagrange multiplier space $\mathbf{F}_{N,h}^n$, discrete counterpart of $\mathbf{F}_N(t)$, as

$$ \mathbf{F}_{N,h}^n = \left\{ \hat{\mathbf{v}}_h \circ (\phi^n)^{-1} \in \mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\epsilon^n) : \hat{\mathbf{v}}_h(\epsilon)(\nu) = \mathbf{a}_0 + \sum_{k=1}^N \mathbf{a}_k \cos(k\nu) + \mathbf{b}_k \sin(k\nu), \quad \mathbf{a}_k, \mathbf{b}_k \in \mathbb{R}^2 \right\}. $$

We then consider a backward Euler scheme for the fluid and solid inertial terms. The remaining kinematic and dynamic conditions of Problem (1.17) are discretized with an implicit scheme. The resulting strongly coupled scheme is given in Algorithm 1.

<caption>**Algorithm 1** Discrete formulation of the reduced-order coupled model (1.17).</caption>

For $n \ge 1$,

**Step 1:** Update interface position:

$$ \phi^n(\hat{\mathbf{x}}) = \mathbf{r}_\epsilon^{n-1} + \mathbf{\Lambda}_\epsilon^{n-1} \hat{\mathbf{x}} \quad \forall \hat{\mathbf{x}} \in \hat{\omega}_\epsilon, \quad \omega_\epsilon^n = \phi^n(\hat{\omega}_\epsilon), $$

where $\mathbf{\Lambda}_\epsilon^{n-1}$ is obtained in terms of $\theta_\epsilon^{n-1}$ by (1.4).

**Step 2:** Find $(\mathbf{u}_{\epsilon,h}^{(N),n}, p_{\epsilon,h}^{(N),n}, \mathbf{\lambda}_{\epsilon,h}^{(N),n}, \mathbf{r}_\epsilon^n, \theta_\epsilon^n) \in \mathbf{V}_h \times Q_h^n \times \mathbf{F}_{N,h}^n \times \mathbb{R}^2 \times \mathbb{R}$ with $\dot{\mathbf{r}}_\epsilon^n = \partial_\tau \mathbf{r}_\epsilon^n$ and $\dot{\theta}_\epsilon^n = \partial_\tau \theta_\epsilon^n$ such that

$$ \left\{ \begin{array}{l} a_\Omega^f((\mathbf{u}_{\epsilon,h}^{(N),n}, p_{\epsilon,h}^{(N),n}), (\mathbf{v}_h, q_h)) - s_h^p(q_h, p_{\epsilon,h}^{(N),n}) - s_h^\lambda(\boldsymbol{\mu}_N, \boldsymbol{\lambda}_{\epsilon,h}^{(N),n}) \\ \quad + \rho_s I_b^\epsilon \partial_\tau \dot{\theta}_\epsilon^n \delta\theta + \rho_s A_b^\epsilon \partial_\tau \dot{\mathbf{r}}_\epsilon^n \cdot \delta\mathbf{r} \\ \quad - (\boldsymbol{\lambda}_{N,h}^n, \mathbf{v}_h - \delta\mathbf{r} - \delta\theta(\mathbf{x} - \mathbf{r}_\epsilon^{n-1})^\perp)_{\epsilon,n} \\ \quad + (\boldsymbol{\mu}_{N,h}, \mathbf{u}_{\epsilon,h}^{(N),n} - \dot{\mathbf{r}}_\epsilon^{n-1} - \dot{\theta}_\epsilon^{n-1}(\mathbf{x} - \mathbf{r}_\epsilon^{n-1})^\perp)_{\epsilon,n} = 0 \end{array} \right. \quad (3.1) $$

for any $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_N, \delta\mathbf{r}, \delta\theta) \in \mathbf{V}_h^n \times Q_h^n \times \mathbf{F}_{N,h}^n \times \mathbb{R}^2 \times \mathbb{R}$ and where $(\cdot, \cdot)_{e,n}$ is the discrete counterpart at time $n\tau$ of $(\cdot, \cdot)_{e,t}$ defined by (1.18).

## 3.2 Comparison with other discrete formulations

For the sake of completeness, in this section, we provide the fully discrete strategies of a fitted discretization and fictitious domain formulations with standard coupling conditions on the fluid-structure interface. These numerical methods will be used in the numerical examples of Section 4 for comparison purposes.

### 3.2.1 Coupled problem with ALE formalism in the fluid

In this section, we consider a triangulation of the reference fluid domain which is fitted to the particle boundary $\partial\omega_\epsilon^n$. We set $\mathbf{V}_h^\partial\hat{\omega}_\epsilon = \mathbf{X}_h^1(\partial\hat{\omega}_\epsilon)$ as the trace space of $\mathbf{X}_h^1(\hat{\Omega}_\epsilon)$ and $\mathbf{\Pi}_{s,h}^1$ the corresponding interpolator onto $\mathbf{V}_h^\partial\hat{\omega}_\epsilon$. The geometrical coupling is then treated explicitly. For a given displacement $\mathbf{r}_\epsilon^{n-1}$ and rotation $\theta_\epsilon^{n-1}$, we define the ALE map $\mathcal{A}_h^n$ at time $n\tau$ as

$$ \phi^n(\hat{\mathbf{x}}) = \mathbf{r}_\epsilon^{n-1} + \mathbf{\Lambda}_\epsilon^{n-1} \hat{\mathbf{x}} \quad \forall \hat{\mathbf{x}} \in \hat{\omega}_\epsilon, \quad \mathcal{A}_h^n = \mathcal{L}_h(\mathbf{\Pi}_{s,h}^1 \phi^n), $$

where $\mathbf{\Lambda}_\epsilon^{n-1}$ given in terms of $\theta_\epsilon^{n-1}$ by (1.4) and $\mathcal{L}_h$ denotes a given lifting operator (e.g., an harmonic extension operator). The one-to-one mapping $\mathcal{A}_h^n$ allows to consider, as fluid unknowns, the fluid velocity and pressure in the reference configuration defined by

$$ \hat{\mathbf{u}}_{\epsilon,h}^n = \mathbf{u}_{\epsilon,h}^n \circ \mathcal{A}_h^n, \quad \hat{p}_{\epsilon,h}^n = p_{\epsilon,h}^n \circ \mathcal{A}_h^n. $$
We then introduce the following finite element space $V_h^f = X_h^1(\hat{\Omega}_\varepsilon) \cap H_0^1(\hat{\Omega}_\varepsilon)$ and $Q_h^f = X_h^1(\hat{\Omega}_\varepsilon) \cap L^2(\hat{\Omega}_\varepsilon)$ for the approximation of the fluid unknowns. We apply a finite element approximation with an implicit treatment of the coupling conditions. The fully discrete strongly coupled scheme is given in Algorithm 2. One limitation of Algorithm 2 lies in the capability of the discrete lifting operator

**Algorithm 2** Fitted discrete formulation of the coupled problem

For $n \ge 1$,

**Step 1:** Update of the particle position and of the fluid domain

$$ egin{cases} \phi^n(\hat{\mathbf{x}}) = \mathbf{r}_\varepsilon^{n-1} + \mathbf{\Lambda}_\varepsilon^{n-1}\hat{\mathbf{x}} & \forall \hat{\mathbf{x}} \in \hat{\omega}_\varepsilon, \\ \mathbf{\Lambda}_h^n = \mathcal{L}_h(\mathbf{\Pi}_{s,h}^1\phi^n), & \Omega_\varepsilon^n = \mathbf{\Lambda}_h^n(\hat{\Omega}_\varepsilon), \end{cases} $$

with $\mathbf{\Lambda}_\varepsilon^{n-1}$ given in terms of $\theta_\varepsilon^{n-1}$ by (1.4).

**Step 2:** Find $(\hat{\mathbf{u}}_{\varepsilon,h}^n, \hat{\mathbf{p}}_{\varepsilon,h}^n, \mathbf{r}_\varepsilon^n, \theta_\varepsilon^n) \in V_h^f \times Q_h^f \times \mathbb{R}^2 \times \mathbb{R}$ with $\dot{\mathbf{r}}_\varepsilon^n = \partial_\tau \mathbf{r}_\varepsilon^n, \dot{\theta}_\varepsilon^n = \partial_\tau \theta_\varepsilon^n$ and such that

$$ a_{\omega_h^{\varepsilon,n}}^f((\mathbf{u}_{\varepsilon,h}^n, p_{\varepsilon,h}^n), (\mathbf{v}_h, q_h)) - s_h^p(q_h, p_{\varepsilon,h}^n) + \rho_b A_b^\varepsilon \partial_\tau \dot{\mathbf{r}}_\varepsilon^n \cdot \delta \mathbf{r} + \rho_b I_b^\varepsilon \partial_\tau \dot{\theta}_\varepsilon^n \delta \theta = 0 \quad (3.2) $$

for any $(\hat{\mathbf{v}}_h, \hat{\mathbf{q}}_h, \delta \mathbf{r}, \delta \theta) \in V_h^f \times Q_h^f \times \mathbb{R}^2 \times \mathbb{R}$ with $\hat{\mathbf{v}}_h = \mathbf{\Pi}_{s,h}^1[\delta \mathbf{r} + \delta \theta(\phi^n - \mathbf{r}_\varepsilon^{n-1})^\perp]$ on $\partial \hat{\omega}_\varepsilon$.

$\mathcal{L}_h$ to guarantee mesh quality when dealing with large interface deflections. This issue is even more pronounced in the context of multiple immersed particles that can move closer to each other and even get into contact. In the next section, we introduce a third algorithm, still imposing the kinematic and kinetic coupling conditions at the interface, that overcomes this issue by working with a fixed mesh and considering a fictitious domain formulation.

### 3.2.2 Coupled problem with fictitious domain method

In this section, we present the discrete formulation of the full coupled problem (1.15). We consider $V_h = X_h^1(\Omega) \cap V$, $Q_h = X_h^1(\Omega) \cap Q$ and $L_H^n = X_H^1(\partial \omega_\varepsilon^n)$ for the approximation of the velocity, pressure and Lagrange multiplier. Similarly to the stabilized discrete formulation of Problem (1.17) presented in Section (2.3), to overcome the instability of the saddle-point problem due to the choice of the Lagrange multiplier finite element space $L_H^n$, and to avoid any restriction on the ration $h/H$, we introduce the following stabilization term, see [?]:

$$ s_h^{BH} = \gamma_\lambda h (\boldsymbol{\mu}_H, \boldsymbol{\lambda}_{\varepsilon,H}^n)_{\partial \omega_\varepsilon^n}, $$

with $\gamma_\lambda > 0$ a user-defined dimensionless parameter. We apply a finite element approximation on the whole domain $\Omega$ with implicit treatment of the coupling conditions. The resulting strongly coupled scheme is given in Algorithm 3.

**Algorithm 3** Discrete formulation of the fictitious domain coupled model

For $n \ge 1$,

**Step 1:** Update interface position:

$$ \phi^n(\hat{\mathbf{x}}) = \mathbf{r}_\varepsilon^{n-1} + \mathbf{\Lambda}_\varepsilon^{n-1}\hat{\mathbf{x}} \quad \forall \hat{\mathbf{x}} \in \hat{\omega}_\varepsilon, \quad \omega_\varepsilon^n = \phi^n(\omega_\varepsilon). $$

where $\mathbf{\Lambda}_\varepsilon^{n-1}$ is obtained from $\theta_\varepsilon^{n-1}$ using nonlinear formula (1.4).

**Step 2:** Find $(\mathbf{u}_{\varepsilon,h}^n, p_{\varepsilon,h}^n, \boldsymbol{\lambda}_H^{\varepsilon,n}, \mathbf{r}_\varepsilon^n, \theta_\varepsilon^n) \in V_h \times Q_h \times L_H^n \times \mathbb{R}^2 \times \mathbb{R}$ with $\dot{\mathbf{r}}_\varepsilon^n = \partial_\tau \mathbf{r}_\varepsilon^n$ and $\dot{\theta}_\varepsilon^n = \partial_\tau \theta_\varepsilon^n$ such that

$$ \left\{ \begin{aligned} & a_{\Omega,h}^f((\mathbf{u}_{\varepsilon,h}^n, p_{\varepsilon,h}^n), (\mathbf{v}_h, q_h)) + \rho_b A_b^\varepsilon \partial_\tau \dot{\mathbf{r}}_\varepsilon^n \cdot \delta \mathbf{r} + \rho_b I_b^\varepsilon \partial_\tau \dot{\theta}_\varepsilon^n \delta \theta \\ & - (\boldsymbol{\lambda}_H^{\varepsilon,n}, \mathbf{v}_h - \delta \mathbf{r} - \delta \theta(\mathbf{x} - \mathbf{r}_\varepsilon^{n-1})^\perp)_{\partial \omega_{\varepsilon,n}} - s_h^p(q_h, p_{\varepsilon,h}^n) - s_h^{BH}(\boldsymbol{\mu}_H, \boldsymbol{\lambda}_H^{\varepsilon,n}) \\ & + (\boldsymbol{\mu}_H, \mathbf{u}_h^n - \dot{\mathbf{r}}_\varepsilon^{n-1} - \dot{\theta}_\varepsilon^{n-1}(\mathbf{x} - \mathbf{r}_\varepsilon^{n-1})^\perp)_{\partial \omega_{\varepsilon,n}} = 0 \end{aligned} \right. $$

for any $(\mathbf{v}_h, q_h, \boldsymbol{\mu}_H, \delta \mathbf{r}, \delta \theta) \in V_h \times Q_h \times L_H^n \times \mathbb{R}^2 \times \mathbb{R}$.
# 4 Numerical experiments

We first evaluate in Section 4.1 the modeling error induced by the reduced Stokes model (2.3) with respect to the full model (2.1). Numerical experiments are conducted for both a circular inclusion, for which Theorem 2.7 applies, and an elliptical inclusion, for which no formal proof has been presented in this paper. The accuracy of Algorithm 1 is illustrated in Section 4.2.1 by comparing its solution with those provided by Algorithms 2 and 3. To this purpose, we consider the case of a single particle and investigate the influence of its size $\epsilon$ and the number of Fourier modes $N$ on the accuracy solution of the reduced-order coupled model (1.17).

## 4.1 Stokes flow

In order to provide numerical evidence of the convergence result stated in Theorem 2.7, we will look at the convergence of the solution of Problem (1.17) with respect to the solution of Problem (1.15). To this purpose, we consider both the cases of a circular inclusion (aspect ratio of

<figure><img src="image_4.png" /><figcaption>Figure 2: Relative $H^1$-error for decreasing values of $\epsilon$ for a single inclusion.</figcaption></figure>

$a = 1)$ and of an elliptical one (aspect ratio of $a = 2$) centered at $(0.1, 0.2)$, and varying in size with $\epsilon \in \{0.1, 0.08, 0.06, 0.04, 0.02\}$. The fluid domain is taken as the unit square $\Omega = [0, 1]^2$ with the following boundary conditions:

$$ 
\left\{ 
\begin{array}{ll}
\mathbf{u}_{\epsilon} = \mathbf{0} & \text{on } [0, 1] \times \{0, 1\}, \\
\sigma(\mathbf{u}_{\epsilon}, p_{\epsilon})e_x = \mathbf{0} & \text{on } \{0\} \times [0, 1], \\
\sigma(\mathbf{u}_{\epsilon}, p_{\epsilon})e_x = -10e_x & \text{on } \{1\} \times [0, 1].
\end{array}
\right. 
$$

The number of modes is chosen such that $N \in \{0, 1\}$. To solve this problem, we use a conforming mesh with respect to the inclusion with the mesh size $h$ sufficiently small such that the numerical error due to the discrete approximation of both the reduced order and full order solutions can be neglected with respect to the error of the model. The $H^1$-error is reported in Figure 2. We recover the convergence rate for the two solutions which behaves as $\epsilon^{N+1}$ for the circle as derived in Section 2.2. Although no theoretical proof has been given, we observe a similar behavior for the elliptical inclusion.

## 4.2 Fluid-particle interaction

In this section,...

### 4.2.1 Particle oscillating in a Stokes flow

As the first example of fluid-structure interaction, we consider the mechanical interaction between a viscous fluid and an elastically supported circular particle in a rectangular domain $\Omega = (-1.5, 1.5) \times (-0.5, 0.5)$. The fluid satisfies the Stokes equations and the boundary conditions for the fluid are similar to the previous test case (see Figure 3):

$$ 
\left\{ 
\begin{array}{ll}
\mathbf{u}_{\epsilon} = \mathbf{0} & \text{on } [-1.5, 1.5] \times \{-0.5, 0.5\}, \\
\sigma(\mathbf{u}_{\epsilon}, p_{\epsilon})e_x = \mathbf{0} & \text{on } \{-1.5\} \times [-0.5, 0.5], \\
\sigma(\mathbf{u}_{\epsilon}, p_{\epsilon})e_x = -10e_x & \text{on } \{1.5\} \times [-0.5, 0.5].
\end{array}
\right. 
$$

We assume that the solid can oscillate along the $e_x$ direction and that its motion is governed by a spring equation with stiffness $T_x$.
<figure><img src="image_5.png" /><figcaption>Figure 3: Elastically supported particle immersed in a Stokes flow.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Time history of the displacement along the x-axis of an elastically supported particle immersed in a Stokesian flow.</figcaption></figure>

The following parameters for the particle model are considered: size $\epsilon \in \{0.06, 0.04, 0.02\}$, density $\rho_b = 1$ and string's stiffness $T_x = 10$. The particle is assumed to be initially at rest. For the fluid, we take a viscosity $\mu_f = 1$. We run the simulations over the time interval $[0, 0.1]$ with a time step $\tau = 5 \cdot 10^{-4}$ and a uniform mesh satisfying $h/\epsilon = 1/4$. In Figure 4, we plot the displacement of the particle along the x-axis obtained with Algorithm 1 and Algorithm 2. For Algorithm 1, the number of Fourier modes is taken $N = 0$. We observe that taking $N = 0$ appears to be sufficient to accurately recover the displacement of the particle, even for relatively large radius. In Figure 5, snapshots of the fluid velocity and pressure are provided at times $t = 0.03$ for $\epsilon = 0.06$. We observe as well that the pressure and velocity of the fluid are well reproduced, at least qualitatively, by the reduced order model with $N = 0$.

### 4.2.2 Elliptical particle in shear flow: Jeffery orbits

We now consider an ellipse of aspect ratio $a = 2$ centered at the origin $0$ in a rectangular domain $\Omega = (-1.5, 1.5) \times (-0.5, 0.5)$. The Stokesian flow is assumed to be driven by a linear shear flow with shear rate $G = 1$ (see Figure 6) and thus we impose the following boundary conditions:

$$ \begin{cases} \mathbf{u}_{\epsilon} = (G, 0) & \text{on } [-1.5, 1.5] \times \{0.5\}, \\ \mathbf{u}_{\epsilon} = (-G, 0) & \text{on } [-1.5, 1.5] \times \{-0.5\}, \\ \mathbf{u}_{\epsilon} \cdot \mathbf{e}_y = 0, \quad \mathbf{e}_x \cdot \boldsymbol{\sigma}(\mathbf{u}_{\epsilon}, p_{\epsilon})\mathbf{e}_x = 0 & \text{on } \{-1.5, 1.5\} \times [-0.5, 0.5]. \end{cases} $$

The physical parameters of the fluid and the solid are given by $\mu_f = 1$ and $\rho_b = 1$. The mesh size is such that $h/\epsilon = 1/4$. The behaviour of the coupled system is studied over the time interval $[0, 30]$ with time-step $\tau = 10^{-2}$. We consider the approximation provided by Algorithm 1 with $N \in \{0, 1\}$ and decreasing size $\epsilon \in \{0.1, 0.05, 0.02\}$. We compare this approximation with the exact solution reported in [37] for $L = \infty$:

$$ \begin{cases} \theta(t) = \tan^{-1} \left( a \tan \left( \frac{aGt}{a^2 + 1} \right) \right), \\ \dot{\theta}(t) = \frac{G}{a^2 + 1} (a^2 \cos^2 \theta(t) + \sin^2 \theta(t)) \end{cases} \quad (4.1) $$

for any $t \in \mathbb{R}^+$. In Figure 7, we plot the time history of the rotation angle of the ellipse. We observe that when $N = 0$, Algorithm 1 fails to accurately capture the rotational effects induced by the shear
<figure><img src="image_7.png" /><figcaption>(a) Algorithm 2, velocity snapshots</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Algorithm 1 with $N = 0$, velocity snapshots</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Algorithm 2, pressure snapshots</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Algorithm 1 with $N = 0$, pressure snapshots</figcaption></figure>

Figure 5: Snapshots of the fluid velocity magnitude (up) and pressure (down) at times $t = 0.03$ for an elastically supported particle obtained with Algorithm 1–2 for $\varepsilon = 0.06$.

<figure><img src="image_11.png" /><figcaption>Figure 6: Geometrical setting for the Jeffery orbits test case.</figcaption></figure>

flow on the ellipse. On the contrary by adding the first Fourier mode ($N = 1$), Algorithm 1 is able to capture the rotation, in particular when the size $\varepsilon$ of the ellipse is small. Note that when the

### 4.2.3 Fall of an elliptical particle in a Stokes flow

At last example, we study the free fall motion of an elliptical particle immersed in a Stokesian flow, with size $\varepsilon = 0.1$ and aspect ratio $a = 2$, subjected to a constant gravitational force $g = -e_y$. The elliptical particle is initially at rest. The fluid domain is rectangular with $\Omega = (-0.5, 0.5) \times (-1.5, 1.5)$.
<figure><img src="image_12.png" /><figcaption>Figure 7: Time history of the rotation of an elliptical particle immersed in a shear flow for different sizes $\epsilon$ and number of Fourier modes $N \in \{0, 1\}$.</figcaption></figure>

We assume the following boundary conditions for the fluid equations (see Figure 8):

$$ 
\left\{ 
\begin{array}{ll}
\mathbf{u}_{\epsilon} = \mathbf{0} & x \in \{-0.5, 0.5\}, \\
\mathbf{u}_{\epsilon} = \mathbf{0} & y = 1.5, \\
\sigma(\mathbf{u}_{\epsilon}, p^{\epsilon})e_{y} = \mathbf{0} & y = -1.5.
\end{array}
\right. 
$$

<figure><img src="image_13.png" /><figcaption>Figure 8: Free fall of an elliptical particle in a Stokesian flow.</figcaption></figure>

The physical parameters of the fluid and solid are given by $\mu_f = 1$ and $\rho_b = 2 \cdot 10^{-2}$. The coupled system is simulated over the time interval $[0, 4000]$ with time-step $\tau = 1$, and an uniform mesh of size $h = 0.025$. To investigate the impact of the number of Fourier modes on Algorithm 1, we consider three different truncation levels $N \in \{0, 1, 2\}$. We compare these approximations with the solution given by Algorithm 3 for a similar time-step $\tau = 1$, a uniform fluid mesh of size $h = 0.01$ and a discretization step for the Lagrange multiplier space of size $H = 0.01$. In Figure 9, we plot the displacement of the elliptical particle along the x- and y-axis, and its rotation. For the case $N = 0$, the ellipse exits the domain at $t = 1674$, which explains why no displacement and rotation are given for the following time steps. Regardless of the number of Fourier modes chosen, we observe that the displacement of the ellipse along the y-axis is accurately captured by Algorithm 1. This observation aligns with the conclusion drawn in Section 4.2.1, suggesting that a single mode ($N = 0$) is sufficient to capture the
translational motion of the particle. Moreover, as discussed in Section 4.2.2, our analysis indicates that at least one additional mode ($N = 1$) is necessary for an accurate representation of the rotational dynamics. Introducing a third mode ($N = 2$) further improves the approximation of the rotation and, due to the nonlinear geometrical coupling between translation and rotation, also improves the representation of the translational motion along the x-axis. In Figure 10, we provide snapshots of the fluid velocity magnitude at times {0, 500, 1000, 2000, 4000}. In addition to the previous observations, we can clearly see the improved capturing of the fluid dynamics when passing from the mode $N = 1$ to the mode $N = 2$, particularly noticeable at times $t = 2000$ and $t = 4000$.

<figure><img src="image_14.png" /><figcaption>Figure 9: Time history of the rotation and displacement of a falling ellipse in a Stokesian flow.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 10: Snapshots of the fluid velocity magnitude at time $t \in \{0, 500, 1000, 2000, 4000\}$ for a falling elliptical particle in a Stokesian flow with Algorithm 1, $N = 1$ and Algorithm 3</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 10: Snapshots of the fluid velocity magnitude at time $t \in \{0, 500, 1000, 2000, 4000\}$ for a falling elliptical particle in a Stokesian flow with Algorithm 1, $N = 1$ and Algorithm 3 (cont.)</figcaption></figure>

# 5 Conclusion

In this paper, we have introduced and analyzed a new stabilized reduced order fictitious domain method for the discretization and simulations of the motion of small rigid particles in a Stokes flow. This method has been analyzed in [31] for the Poisson problem with small inclusions and used to simulate slender elastic structure in 3D Navier-Stokes equations in [30]. The key feature of the method relies on the projection of the coupling conditions on a finite dimensional space of Fourier type. In a first part the modeling error is obtained for the Stokes equations with non-homogeneous Dirichlet boundary conditions imposed thanks to the introduction of a Lagrange multiplier: the reduced problem is well-posed and converges when $\varepsilon \to 0$ towards the problem with standard Dirichlet boundary conditions. Furthermore, the stability and convergence of the finite element stabilized scheme is derived. The crucial step is the derivation of an inf-sup condition for our problem, which writes with two constraints: the divergence constraint and the Dirichlet one. In particular the inf-sup constant is proven to be independent either of number of Fourier mode or mesh size but all of the size $\varepsilon$ of the particle. As in the Poisson case and as it is standard for fictitious domain methods, the rate of convergence is sub-optimal due on one hand to the discontinuity of pressure and velocity gradient across the immersed interface.

To overcome this issue an augmented finite element method could be proposed as in [31] and shall be the object of a forthcoming paper. Yet to avoid the coupling between the finite element basis and the additional singular functions designed to better capture the discontinuities and asymptotic singularities of the solution, one should consider stable velocity-pressure conforming finite element pairs in order to avoid interaction between these two types of basis functions.
# A Appendix

## A.1 Equivalence between continuous formulations of the full coupled problem

In this section, we justify the equivalence between formulation (1.10)–(1.12) and the one considered herein (1.14), (1.15). First we consider $(\mathbf{u}_\varepsilon, \tilde{p}_\varepsilon, \tilde{\lambda}_\varepsilon, \mathbf{r}_\varepsilon, \theta_\varepsilon)$ such that for a.e. $t \in \mathbb{R}^+$,

$$ (\mathbf{u}_\varepsilon(t), \tilde{p}_\varepsilon(t), \tilde{\lambda}_\varepsilon(t), \mathbf{r}(t), \theta(t)) \in \mathbf{V} \times \tilde{\mathcal{Q}}(t) \times \mathbf{H}^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) \times \mathbb{R}^2 \times \mathbb{R} $$

satisfies

$$ \begin{aligned} \phi_\varepsilon(t, \hat{\mathbf{x}}) = \mathbf{r}_\varepsilon + \mathbf{\Lambda}_\varepsilon \hat{\mathbf{x}} \quad \forall \hat{\mathbf{x}} \in \hat{\omega}_\varepsilon, \\ \left\{ \begin{array}{l} 2\mu a_\Omega((\mathbf{u}_\varepsilon, p_\varepsilon), (\mathbf{v}, \tilde{\mathbf{q}})) + (p_\varepsilon, \operatorname{div} \mathbf{v})_\Omega - (\tilde{\mathbf{q}}, \operatorname{div} \mathbf{u}_\varepsilon)_\Omega + \rho_b A_b^\varepsilon \ddot{\mathbf{r}}_\varepsilon \cdot \delta\mathbf{r} + \rho_b I_b^\varepsilon \ddot{\theta}_\varepsilon \delta\theta \\ - \langle \mathbf{\lambda}_\varepsilon, \mathbf{v} - \delta\mathbf{r} - \delta\theta(\mathbf{x} - \mathbf{r}_\varepsilon)^\perp \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} + \langle \tilde{\mathbf{\mu}}, \mathbf{u}_\varepsilon - \dot{\mathbf{r}}_\varepsilon - \dot{\theta}_\varepsilon(\mathbf{x} - \mathbf{r}_\varepsilon)^\perp \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} = 0, \end{array} \right. \end{aligned} $$

for any $(\mathbf{v}, \tilde{\mathbf{q}}, \tilde{\mathbf{\mu}}, \delta\mathbf{r}, \delta\theta) \in \mathbf{V} \times \tilde{\mathcal{Q}}(t) \times \mathbf{H}^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) \times \mathbb{R}^2 \times \mathbb{R}$, with

$$ \tilde{\mathcal{Q}}(t) = \{q \in L^2(\Omega), q|_{\Omega_\varepsilon(t)} \in L_0^2(\Omega_\varepsilon(t)), q|_{\omega_\varepsilon(t)} \in L_0^2(\omega_\varepsilon(t))\}. $$

We are going to build $(p_\varepsilon, \mathbf{\lambda}_\varepsilon)$ in $Q(t) \times \mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t))$ such that $(\mathbf{u}_\varepsilon, p_\varepsilon, \mathbf{\lambda}_\varepsilon, \mathbf{r}_\varepsilon, \theta_\varepsilon)$ is a solution of (1.14), (1.15) for any $(\mathbf{v}, q, \mathbf{\mu}, \delta\mathbf{r}, \delta\theta) \in \mathbf{V} \times \mathbf{Q} \times \mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) \times \mathbb{R}^2 \times \mathbb{R}$. Taking $(\mathbf{v}, \delta\mathbf{r}, \delta\theta) = \mathbf{0}$ implies, for any $(\tilde{\mathbf{\mu}}, \tilde{\mathbf{q}}) \in \mathbf{H}^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) \times \tilde{\mathcal{Q}}(t)$,

$$ \langle \tilde{\mathbf{\mu}}, \mathbf{u}_\varepsilon \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} = \langle \tilde{\mathbf{\mu}}, \dot{\mathbf{r}}_\varepsilon + \dot{\theta}_\varepsilon(\mathbf{x} - \mathbf{r}_\varepsilon)^\perp \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)}, \quad (A.1) $$

$$ (\tilde{\mathbf{q}}, \operatorname{div} \mathbf{u}_\varepsilon)_\Omega = 0. \quad (A.2) $$

Since $\mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)) \subset \mathbf{H}^{-\frac{1}{2}}(\partial\omega_\varepsilon(t))$, we have

$$ \langle \mathbf{\mu}, \mathbf{u}_\varepsilon \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} = \langle \mathbf{\mu}, \dot{\mathbf{r}}_\varepsilon + \dot{\theta}_\varepsilon(\mathbf{x} - \mathbf{r}_\varepsilon)^\perp \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)}, \quad \forall \mathbf{\mu} \in \mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t)). \quad (A.3) $$

Moreover, taking $\tilde{\mathbf{\mu}} = \mathbf{n}^-$ in (A.1) yields

$$ \langle \mathbf{n}^-, \mathbf{u}_\varepsilon \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} = 0. \quad (A.4) $$

For any domain $\mathcal{D}$ and for any $v \in L^2(\mathcal{D})$, we introduce the average operator

$$ \int_{\mathcal{D}} v = \frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} v. $$

For $q \in \mathcal{Q}$, we define $\tilde{q} \in \tilde{\mathcal{Q}}(t)$ as

$$ \tilde{q} = q - \left( \int_{\Omega_\varepsilon(t)} q \right) 1_{\Omega_\varepsilon(t)} - \left( \int_{\omega_\varepsilon(t)} q \right) 1_{\omega_\varepsilon(t)}. $$

Therefore, by integration by parts, we obtain, for any $q \in \mathcal{Q}$,

$$ \begin{aligned} (q, \operatorname{div} \mathbf{u}_\varepsilon)_\Omega &= (q, \operatorname{div} \mathbf{u}_\varepsilon)_{\Omega_\varepsilon(t)} + (q, \operatorname{div} \mathbf{u}_\varepsilon)_{\omega_\varepsilon(t)} \\ &= (\tilde{q}, \operatorname{div} \mathbf{u}_\varepsilon)_\Omega + \left( \int_{\omega_\varepsilon(t)} q - \int_{\Omega_\varepsilon(t)} q \right) \langle \mathbf{n}^-, \mathbf{u}_\varepsilon \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)}. \end{aligned} $$

Thus, taking into account (A.2) and (A.4) leads to

$$ (q, \operatorname{div} \mathbf{u}_\varepsilon)_\Omega = 0 \quad \forall q \in \mathcal{Q}. \quad (A.5) $$

Consequently the velocity $\mathbf{u}_\varepsilon$ satisfies the desired constraints. Next from the Lagrange multipliers $(\tilde{p}_\varepsilon, \tilde{\mathbf{\lambda}}_\varepsilon) \in \tilde{\mathcal{Q}}(t) \times \mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t))$ let us build $(p_\varepsilon, \mathbf{\lambda}_\varepsilon) \in \mathcal{Q} \times \mathbf{H}_0^{-\frac{1}{2}}(\partial\omega_\varepsilon(t))$. We set

$$ \mathbf{\lambda}^\varepsilon = \tilde{\mathbf{\lambda}}_\varepsilon - \left( \frac{1}{|\partial\omega_\varepsilon|} \langle \tilde{\mathbf{\lambda}}_\varepsilon, \mathbf{n}^- \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} \right) \cdot \mathbf{n}^-, $$

$$ p_\varepsilon = \tilde{p}_\varepsilon + \left( \frac{1}{|\partial\omega_\varepsilon|} \langle \tilde{\mathbf{\lambda}}_\varepsilon, \mathbf{n}^- \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_\varepsilon(t)} \right) (\alpha \mathbf{1}_{\omega_\varepsilon} - (1-\alpha) \mathbf{1}_{\Omega_\varepsilon}), $$
where $\alpha = |\Omega_{\varepsilon}|/(|\Omega_{\varepsilon}| + |\omega_{\varepsilon}|) = |\Omega_{\varepsilon}|/|\Omega|$. By definition $\langle \lambda_{\varepsilon}, n^{-} \rangle_{H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t))} = 0$, thus $\lambda_{\varepsilon} \in H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t))$ and the choice of $\alpha$ guarantees that $\int_{\Omega} p_{\varepsilon} = 0$ and thus $p_{\varepsilon} \in Q$. Hence, owing to the identity

$$ \langle n^{-}, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = \alpha (1_{\omega_{\varepsilon}}, \operatorname{div} v)_{\Omega} - (1 - \alpha) (1_{\Omega_{\varepsilon}}, \operatorname{div} v)_{\Omega} $$

that holds for any $v \in H_{0}^{1}(\Omega)$, we get

$$ (p_{\varepsilon}, \operatorname{div}(v))_{\Omega} - \langle \lambda_{\varepsilon}, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = (\tilde{p}_{\varepsilon}, \operatorname{div}(v))_{\Omega} - \langle \tilde{\lambda}_{\varepsilon}, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)}. \quad (A.6) $$

From (A.3), (A.5), (A.6), we thus have build $(u_{\varepsilon}(t), p_{\varepsilon}(t), \lambda_{\varepsilon}(t), r_{\varepsilon}(t), \theta_{\varepsilon}(t)) \in V \times Q \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R}$ that satisfies (1.15) for any test function $(v, q, \mu, \delta r, \delta\theta) \in V \times Q \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R}$. Reciprocally, let $(u_{\varepsilon}, p_{\varepsilon}, \lambda_{\varepsilon}, r_{\varepsilon}, \theta_{\varepsilon})$ such that $(u_{\varepsilon}(t), p_{\varepsilon}(t), \lambda_{\varepsilon}(t), r(t), \theta(t)) \in V \times Q \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R}$ and satisfies (1.15) for any test function $(v, q, \mu, \delta r, \delta\theta) \in V \times Q \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R}$. Let $(\tilde{q}, \tilde{\mu}) \in \tilde{Q}(t) \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t))$, defining $(q, \mu)$ by

$$ \begin{aligned} \mu &= \tilde{\mu} - \left( \frac{1}{|\partial\omega_{\varepsilon}|} \langle \tilde{\mu}, n^{-} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} \right) n^{-}, \\ q &= \tilde{q} + \left( \frac{1}{|\partial\omega_{\varepsilon}|} \langle \tilde{\mu}, n^{-} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} \right) (\alpha 1_{\omega_{\varepsilon}} - (1 - \alpha) 1_{\Omega_{\varepsilon}}), \end{aligned} $$

it satisfies, as previously, $(q, \mu) \in Q \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t))$ and

$$ (q, \operatorname{div}(u_{\varepsilon}))_{\Omega} - \langle \mu, u_{\varepsilon} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = (\tilde{q}, \operatorname{div}(u_{\varepsilon}))_{\Omega} - \langle \tilde{\mu}, u_{\varepsilon} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)}. \quad (A.7) $$

Moreover, since

$$ \langle n^{-}, \dot{r}_{\varepsilon} + \dot{\theta}_{\varepsilon}(x - r_{\varepsilon})^{\perp} \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = \int_{\partial\omega_{\varepsilon}(t)} (r_{\varepsilon} + \dot{\theta}_{\varepsilon}(x - r_{\varepsilon})^{\perp}) \cdot n^{-} = 0, $$

we obtain

$$ \langle \mu, \dot{r}_{\varepsilon} + \dot{\theta}_{\varepsilon}(x - r_{\varepsilon})^{\perp} \rangle_{H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)), H^{\frac{1}{2}}(\partial\omega_{\varepsilon}(t))} = \langle \tilde{\mu}, \dot{r}_{\varepsilon} + \dot{\theta}_{\varepsilon}(x - r_{\varepsilon})^{\perp} \rangle_{H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)), H^{\frac{1}{2}}(\partial\omega_{\varepsilon}(t))}, \quad (A.8) $$

Finally we set

$$ \begin{aligned} \tilde{\lambda}_{\varepsilon} &= \lambda_{\varepsilon} + \left( \int_{\omega_{\varepsilon}(t)} p_{\varepsilon} - \int_{\Omega_{\varepsilon}(t)} p_{\varepsilon} \right) n^{-}, \\ \tilde{p}_{\varepsilon} &= p_{\varepsilon} - \left( \int_{\Omega_{\varepsilon}(t)} p_{\varepsilon} \right) 1_{\Omega_{\varepsilon}(t)} - \left( \int_{\omega_{\varepsilon}(t)} p_{\varepsilon} \right) 1_{\omega_{\varepsilon}(t)}. \end{aligned} $$

They satisfy $(\tilde{\lambda}_{\varepsilon}, \tilde{p}_{\varepsilon}) \in H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \tilde{Q}(t)$ and from the identity

$$ \left( \left( \int_{\Omega_{\varepsilon}(t)} p_{\varepsilon} \right) 1_{\Omega_{\varepsilon}(t)} + \left( \int_{\omega_{\varepsilon}(t)} p_{\varepsilon} \right) 1_{\omega_{\varepsilon}(t)}, \operatorname{div} v \right)_{\Omega} = \left( \int_{\omega_{\varepsilon}(t)} p_{\varepsilon} - \int_{\Omega_{\varepsilon}(t)} p_{\varepsilon} \right) \langle n^{-}, v \rangle_{H^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)), H^{\frac{1}{2}}(\partial\omega_{\varepsilon}(t))}, $$

we deduce that for any $v \in H_{0}^{1}(\Omega)$,

$$ (p_{\varepsilon}, \operatorname{div}(v))_{\Omega} - \langle \lambda_{\varepsilon}, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)} = (\tilde{p}_{\varepsilon}, \operatorname{div}(v))_{\Omega} - \langle \tilde{\lambda}_{\varepsilon}, v \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega_{\varepsilon}(t)}. \quad (A.9) $$

From (A.7)-(A.9), we finally obtain that for a.e $t \in \mathbb{R}^{+}$,

$$ (u_{\varepsilon}(t), \tilde{p}_{\varepsilon}(t), \tilde{\lambda}_{\varepsilon}(t), r_{\varepsilon}(t), \theta_{\varepsilon}(t)) \in V \times \tilde{Q}(t) \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R} $$

and satisfies (1.15) for any test function $(v, \tilde{q}, \tilde{\mu}, \delta r, \delta\theta) \in V \times \tilde{Q}(t) \times H_{0}^{-\frac{1}{2}}(\partial\omega_{\varepsilon}(t)) \times \mathbb{R}^{2} \times \mathbb{R}$.

## A.2 Preliminary lemmas

In this section of the appendix, we state several results that will be used in the next sections of the appendix for the proof of some lemmas used in this paper. We recall that $\varepsilon \in (0, \varepsilon_{m})$ where $\varepsilon_{m}$ is the fixed parameter defined by $\varepsilon_{m} = \sup\{\varepsilon > 0 | \tilde{\omega}_{\varepsilon} \subset \Omega\}$. The first three lemmas and their proof are drawn from [38, Lemmas C.1–C.3], themselves based on the results presented in [39, Chapter 3].
**Lemma A.1.** For $\varphi \in H_0^{\frac{1}{2}}(\partial\Omega)$, the following problem

$$ egin{cases} 
-\Delta \mathbf{v}_{\varepsilon} + \nabla q_{\varepsilon} = \mathbf{0} & \text{in } \Omega_{\varepsilon}, \\
\operatorname{div} \mathbf{v}_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\
\mathbf{v}_{\varepsilon} = \varphi & \text{on } \partial\Omega, \\
\mathbf{v}_{\varepsilon} = \mathbf{0} & \text{on } \partial\omega_{\varepsilon}, 
\end{cases} 
\quad (A.10) $$

admits a unique weak solution in $H^1(\Omega_\varepsilon) \times L_0^2(\Omega_\varepsilon)$. Moreover, the solution satisfies the estimate

$$ \|\mathbf{v}_{\varepsilon}\|_{1,\Omega_{\varepsilon}} \lesssim \|\varphi\|_{\frac{1}{2},\partial\Omega}. $$

**Proof of Lemma A.1.** We consider $\mathbf{v}_{\varepsilon_m}$ the solution of (A.10) for $\varepsilon = \varepsilon_m$. It satisfies

$$ |\mathbf{v}_{\varepsilon_m}|_{1,\Omega_{\varepsilon_m}} = \left( \int_{\Omega_{\varepsilon_m}} |\nabla \mathbf{v}_{\varepsilon_m}|^2 \right)^{\frac{1}{2}} \lesssim \|\varphi\|_{\frac{1}{2},\partial\Omega}. \quad (A.11) $$

Now consider $\tilde{\mathbf{v}}_{\varepsilon_m} \in H^1(\Omega)$ the extension by $\mathbf{0}$ on $\Omega$. Since $\varepsilon < \varepsilon_m$, we have $\omega_\varepsilon \subset \omega_{\varepsilon_m}$ and $\Omega_{\varepsilon_m} \subset \Omega_\varepsilon$. Since $\mathbf{v}_\varepsilon$ minimizes the Dirichlet energy on

$$ \{\mathbf{v} \in H^1(\Omega_\varepsilon) \mid \mathbf{v}|_{\partial\Omega} = \varphi, \mathbf{v}|_{\partial\omega_\varepsilon} = \mathbf{0}, \operatorname{div}(\mathbf{v}) = 0\}, $$

we deduce that

$$ |\mathbf{v}_{\varepsilon}|_{1,\Omega_{\varepsilon}} \le |\tilde{\mathbf{v}}_{\varepsilon_m}|_{1,\Omega_{\varepsilon}} = |\mathbf{v}_{\varepsilon_m}|_{1,\Omega_{\varepsilon_m}}. $$

Finally, together with estimate (A.11), we obtain

$$ |\mathbf{v}_{\varepsilon}|_{1,\Omega_{\varepsilon}} \lesssim \|\varphi\|_{\frac{1}{2},\partial\Omega}. \quad (A.12) $$

Now consider $(\mathbf{v}_0, l_0) \in H^1(\Omega) \times L_0^2(\Omega)$ the solution of the system

$$ egin{cases} 
-\Delta \mathbf{v}_0 + \nabla l_0 = \mathbf{0} & \text{in } \Omega, \\
\operatorname{div} \mathbf{v}_0 = 0 & \text{in } \Omega, \\
\mathbf{v}_0 = \varphi & \text{on } \partial\Omega. 
\end{cases} $$

The unique solution of this system satisfies the energy estimate

$$ \|\mathbf{v}_0\|_{1,\Omega} \lesssim \|\varphi\|_{\frac{1}{2},\partial\Omega}. \quad (A.13) $$

Denoting by $\tilde{\mathbf{v}}_\varepsilon$ the extension by zero of $\mathbf{v}_\varepsilon$ on $\Omega$, since $\tilde{\mathbf{v}}_\varepsilon - \mathbf{v}_0 \in H_0^1(\Omega)$, we can apply the Poincaré inequality on $\Omega$ to get

$$ \|\mathbf{v}_\varepsilon\|_{1,\Omega_\varepsilon} = \|\tilde{\mathbf{v}}_\varepsilon\|_{1,\Omega} \le \|\tilde{\mathbf{v}}_\varepsilon - \mathbf{v}_0\|_{1,\Omega} + \|\mathbf{v}_0\|_{1,\Omega} \lesssim (|\mathbf{v}_\varepsilon|_{1,\Omega_\varepsilon} + |\mathbf{v}_0|_{1,\Omega}) + \|\mathbf{v}_0\|_{1,\Omega}. $$

Combining this inequality with (A.12) and (A.13), we deduce

$$ \|\mathbf{v}_\varepsilon\|_{\Omega_\varepsilon} \lesssim \|\varphi\|_{\frac{1}{2},\partial\Omega} $$

which concludes the proof of the lemma.

□

From Lemma A.1, we can deduce

**Lemma A.2.** Let $\varphi \in H^1(\Omega)$ satisfy $\operatorname{div} \varphi = 0$ in $\Omega$, and there exists $q \in L_0^2(\Omega)$ such that

$$ -\Delta \varphi + \nabla q = 0 \quad \text{in } \Omega. $$

Let $(\mathbf{v}_\varepsilon, q_\varepsilon) \in H^1(\omega_\varepsilon) \times L_0^2(\Omega_\varepsilon)$ denote the solution of the Stokes problem

$$ egin{cases} 
-\Delta \mathbf{v}_{\varepsilon} + \nabla q_{\varepsilon} = \mathbf{0} & \text{in } \Omega_{\varepsilon}, \\
\operatorname{div} \mathbf{v}_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\
\mathbf{v}_{\varepsilon} = \varphi & \text{on } \partial\Omega, \\
\mathbf{v}_{\varepsilon} = \mathbf{0} & \text{on } \partial\omega_{\varepsilon}, 
\end{cases} $$

Then $\mathbf{v}_\varepsilon$ satisfies the estimate

$$ \|\mathbf{v}_{\varepsilon}\|_{1,\Omega_{\varepsilon}} \lesssim \|\varphi\|_{\frac{1}{2},\Omega}. $$
Proof of Lemma A.2. We consider the pair $(\mathbf{w}_{\varepsilon} = \mathbf{v}_{\varepsilon} - \varphi, l_{\varepsilon} = q_{\varepsilon} - q)$ which satisfies

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{w}_{\varepsilon} + \nabla l_{\varepsilon} = \mathbf{0} & \text{in } \Omega_{\varepsilon}, \\
\operatorname{div} \mathbf{w}_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\
\mathbf{w}_{\varepsilon} = -\varphi & \text{on } \partial\Omega, \\
\mathbf{w}_{\varepsilon} = \mathbf{0} & \text{on } \partial\omega_{\varepsilon},
\end{array}
\right.
$$

It follows from Lemma A.1 that

$$
\| \mathbf{w}_{\varepsilon} \|_{1, \Omega_{\varepsilon}} \lesssim \| \varphi \|_{\frac{1}{2}, \partial \Omega}. \qquad (A.14)
$$

Since $\varphi$ is itself solution of the Stokes equations, we also have

$$
\| \varphi \|_{1, \Omega_{\varepsilon}} \le \| \varphi \|_{1, \Omega} \lesssim \| \varphi \|_{\frac{1}{2}, \partial \Omega}. \quad (A.15)
$$

Combining (A.14) and (A.15), we obtain

$$
\| \mathbf{v}_{\varepsilon} \|_{1, \Omega_{\varepsilon}} \le \| \mathbf{w}_{\varepsilon} \|_{1, \Omega_{\varepsilon}} + \| \varphi \|_{1, \Omega_{\varepsilon}} \lesssim \| \varphi \|_{\frac{1}{2}, \partial \Omega}.
$$

□

Next if the boundary data is a constant vector we can derive a better estimate. Namely,

**Lemma A.3.** For any $L \in \mathbb{R}^2$

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{v}_{\varepsilon} + \nabla q_{\varepsilon} = \mathbf{0} & \text{in } \Omega_{\varepsilon}, \\
\operatorname{div} \mathbf{v}_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\
\mathbf{v}_{\varepsilon} = \mathbf{0} & \text{on } \partial\Omega, \\
\mathbf{v}_{\varepsilon} = \mathbf{L} & \text{on } \partial\omega_{\varepsilon},
\end{array}
\right.
$$

admits a unique weak solution in $H^1(\Omega_\varepsilon) \times L_0^2(\Omega_\varepsilon)$, satisfying the following a priori estimate

$$
\| v_{\varepsilon} \|_{1, \Omega_{\varepsilon}} \lesssim |\log(\varepsilon)|^{-\frac{1}{2}} |L|.
$$

Proof of Lemma A.3. Let $(\mathbf{v}(\varepsilon), q(\varepsilon))$ denote the rescaled unknowns defined by

$$
(\mathbf{v}(\varepsilon)(\mathbf{x}), q(\varepsilon)(\mathbf{x})) = (\mathbf{v}_{\varepsilon}(\mathbf{x}/\varepsilon), q_{\varepsilon}(\mathbf{x}/\varepsilon)) \quad \forall \mathbf{x} \in \Omega.
$$

They satisfy

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{v}(\varepsilon) + \varepsilon \nabla q(\varepsilon) = \mathbf{0} & \text{in } \frac{\Omega}{\varepsilon} \setminus \bar{\omega}, \\
\operatorname{div} \mathbf{v}(\varepsilon) = 0 & \text{in } \frac{\Omega}{\varepsilon}, \\
\mathbf{v}(\varepsilon) = \mathbf{0} & \text{on } \frac{\partial \Omega}{\varepsilon}, \\
\mathbf{v}(\varepsilon) = \mathbf{L} & \text{on } \partial \omega.
\end{array}
\right.
$$

We also introduce $(\mathbf{w}(\varepsilon), l(\varepsilon))$ the solution of

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{w}(\varepsilon) + \nabla l(\varepsilon) = \mathbf{0} & \text{in } \omega_{\varepsilon m/\varepsilon} \setminus \bar{\omega}, \\
\operatorname{div} \mathbf{w}(\varepsilon) = 0 & \text{in } \omega_{\varepsilon m/\varepsilon} \setminus \bar{\omega}, \\
\mathbf{w}(\varepsilon) = \mathbf{0} & \text{on } \partial\omega_{\varepsilon m/\varepsilon}, \\
\mathbf{w}(\varepsilon) = \mathbf{L} & \text{on } \partial\omega.
\end{array}
\right.
$$

Since $\omega \subset \omega_{\varepsilon m/\varepsilon} \subset \frac{\Omega}{\varepsilon}$, we can consider $\tilde{w}(\varepsilon)$ the extension by zero of $\mathbf{w}(\varepsilon)$ in $\frac{\Omega}{\varepsilon} \setminus \omega_{\varepsilon m/\varepsilon}$ and $\mathbf{L}$ in $\omega$.
We have $\tilde{w}(\varepsilon) \in H^1(\frac{\Omega}{\varepsilon})$ and since $\mathbf{v}(\varepsilon)$ minimizes the Dirichlet energy on

$$
\{ \mathbf{v} \in H^1 \left( \frac{\Omega}{\varepsilon} \right) \mid \mathbf{v} \mid_{\frac{\partial \Omega}{\varepsilon}} = 0, \mathbf{v} \mid_{\partial \omega} = \mathbf{L}, \operatorname{div} \mathbf{v} = 0 \},
$$

we deduce

$$
|\mathbf{v}_{\varepsilon}|_{1, \Omega_{\varepsilon}} = |\mathbf{v}(\varepsilon)|_{1, \frac{\Omega}{\varepsilon} \setminus \bar{\omega}} \le |\tilde{\mathbf{w}}(\varepsilon)|_{1, \frac{\Omega}{\varepsilon} \setminus \bar{\omega}} = |\mathbf{w}(\varepsilon)|_{1, \omega_{\varepsilon m/\varepsilon} \setminus \bar{\omega}}. \quad (A.16)
$$
Let $\psi(\varepsilon) = L + 4\pi E \frac{L}{\log(\varepsilon_m/\varepsilon)}$ and $\pi(\varepsilon) = 4\pi P \cdot \frac{L}{\log(\varepsilon_m/\varepsilon)}$, with $(E, P)$ the fundamental solution of Stokes equations in $\mathbb{R}^2$ given by

$$ E(x) = \frac{1}{4\pi} \left( -\log(|x|)I_2 + \frac{xx^T}{|x|^2} \right), \quad P(x) = \frac{x}{2\pi|x|^2} \quad \forall x \in \mathbb{R}^2 \setminus \mathbf{0}. \qquad (A.17) $$

The couple $(\psi(\varepsilon), \pi(\varepsilon))$ satisfies

$$ \left\{ \begin{array}{ll} -\Delta\psi(\varepsilon) + \nabla q(\varepsilon) = \mathbf{0} & \text{in } \omega_{\varepsilon_m/\varepsilon} \setminus \bar{\omega}, \\ \operatorname{div} \psi(\varepsilon) = 0 & \text{in } \omega_{\varepsilon_m/\varepsilon} \setminus \bar{\omega}, \\ \psi(\varepsilon) = \frac{xx^T L}{|x|^2 \log(\varepsilon_m/\varepsilon)} & \text{on } \partial\omega_{\varepsilon_m/\varepsilon}, \\ \psi(\varepsilon) = L + \frac{xx^T L}{|x|^2 \log(\varepsilon_m/\varepsilon)} & \text{on } \partial\omega. \end{array} \right. $$

Moreover,

$$ \nabla\psi(\varepsilon)(x) = \frac{1}{\log(\varepsilon_m/\varepsilon)|x|^2} \left( -Lx^T + xL^T + (L \cdot x)I_2 - 2(L \cdot x)\frac{xx^T}{|x|^2} \right) $$

and a straightforward calculation yields

$$ |\psi(\varepsilon)|_{1,\omega_{\varepsilon_m/\varepsilon}\setminus\bar{\omega}} \lesssim |\log(\varepsilon)|^{-\frac{1}{2}}|L|. $$

Since the pair $(w(\varepsilon) - \psi(\varepsilon), l(\varepsilon) - q(\varepsilon))$ solves the Stokes equations with boundary conditions $-\frac{xx^T L}{|x|^2 \log(\varepsilon_m/\varepsilon)}$ on $\partial\omega \cup \partial\omega_{\varepsilon_m/\varepsilon}$, it follows from Lemma A.1 and Lemma A.2, applied to the corresponding rescaled unknowns, that

$$ |w(\varepsilon) - \psi(\varepsilon)|_{1,\omega_{\varepsilon_m/\varepsilon}\setminus\bar{\omega}} \lesssim |\log(\varepsilon)|^{-\frac{1}{2}} \left\| \frac{xx^T}{|x|^2} L \right\|_{\frac{1}{2},\partial\omega_{\varepsilon_m}} \lesssim |\log(\varepsilon)|^{-\frac{1}{2}} |L|. $$

We then deduce from (A.16), and the two previous inequalities

$$ \begin{aligned} |v_{\varepsilon}|_{1,\Omega_{\varepsilon}} &\lesssim |w(\varepsilon)|_{1,\omega_{\varepsilon_m/\varepsilon}\setminus\bar{\omega}} \lesssim |w(\varepsilon) - \psi(\varepsilon)|_{1,\omega_{\varepsilon_m/\varepsilon}\setminus\bar{\omega}} + |\psi(\varepsilon)|_{1,\omega_{\varepsilon_m/\varepsilon}\setminus\bar{\omega}}, \\ &\lesssim |\log(\varepsilon)|^{-\frac{1}{2}}|L|. \end{aligned} \qquad (A.18) $$

Next, we need to bound the $L^2$ norm of $v_{\varepsilon}$. Let $\tilde{v}_{\varepsilon}$ denote the extension of $v_{\varepsilon}$ by $L$ in $\omega$. Since this extension is in $H_0^1(\Omega)$, we can use the Poincaré inequality to get

$$ \|v_{\varepsilon}\|_{\Omega_{\varepsilon}} \lesssim \|\tilde{v}_{\varepsilon}\|_{\Omega} \lesssim |\tilde{v}_{\varepsilon}|_{1,\Omega} \lesssim |v_{\varepsilon}|_{1,\Omega_{\varepsilon}}. $$

From the estimate (A.18) of the $H^1$ semi norm, we obtain the desired result. $\square$

Next we consider a boundary data $\varphi \in H_0^1(\partial\omega_{\varepsilon})$. Once again, the following lemma and its proof are drawn from [38, Lemma B.2, Lemma 4.2], themselves based on results presented in [39, Chapter 3].

**Lemma A.4.** For any $\varphi \in H_0^1(\partial\omega_{\varepsilon})$, the problem

$$ \left\{ \begin{array}{ll} -\Delta v_{\varepsilon} + \nabla q_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\ \operatorname{div} v_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\ v_{\varepsilon} = \mathbf{0} & \text{on } \partial\Omega, \\ v_{\varepsilon} = \varphi & \text{on } \partial\omega_{\varepsilon}, \end{array} \right. $$

admits a unique weak solution $v_{\varepsilon} \in H^1(\Omega_{\varepsilon})$, which satisfies the energy estimate

$$ \|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} \lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2},\partial\omega} = \|\varphi\|_{\frac{1}{2},\varepsilon}. $$
Proof. We first consider the following exterior Stokes problem

$$
\left\{
\begin{array}{ll}
-\Delta \mathbf{v} + \nabla q = \mathbf{0} & \text{in } \mathbb{R}^2 \setminus \bar{\omega}, \\
\operatorname{div} \mathbf{v} = 0 & \text{in } \mathbb{R}^2 \setminus \bar{\omega}, \\
\mathbf{v} = \varphi(\varepsilon) & \text{on } \partial\omega.
\end{array}
\right.
\qquad (A.19)
$$

where we recall that $\varphi(\varepsilon)(\mathbf{x}) = \varphi(\varepsilon\mathbf{x})$ for any $\mathbf{x} \in \partial\omega$. According to [40, Theorem 3.4] this problem admits unique solution $(\mathbf{u}, p) \in \mathbf{H}_{-1}(\mathbb{R}^2\setminus\bar{\omega}) \times L^2(\mathbb{R}^2\setminus\bar{\omega})$ where

$$
\mathbf{H}_{-1}(\mathbb{R}^2 \setminus \bar{\omega}) = \{ \mathbf{u} \in \mathcal{D}'(\mathbb{R}^2 \setminus \bar{\omega}) \mid (\rho \log(\rho))^{-1} \mathbf{u} \in \mathcal{L}^2(\mathbb{R}^2 \setminus \bar{\omega}), \nabla \mathbf{u} \in \mathcal{L}^2(\mathbb{R}^2 \setminus \bar{\omega}), \operatorname{div} \mathbf{u} = \mathbf{0} \}
$$

with

$$
\rho(\mathbf{x}) = (1 + |\mathbf{x}|^2)^{\frac{1}{2}}.
$$

Moreover, defining

$$
\| \mathbf{u} \|_{-1,1,\mathbb{R}^2 \setminus \bar{\omega}} = \| \mathbf{u} \|_{\mathbb{R}^2 \setminus \bar{\omega}} + \|(\rho \log(\rho))^{-1} \nabla \mathbf{u} \|_{\mathbb{R}^2 \setminus \bar{\omega}}
$$

the solution satisfies the following energy estimate

$$
\|v\|_{-1,1,\mathbb{R}^2\setminus\bar{\omega}} + \|q\|_{\mathbb{R}^2\setminus\bar{\omega}} \lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2},\partial\omega} \quad (A.20)
$$

By extending $(\mathbf{v}, q)$ in $\mathbf{H}_{-1}(\mathbb{R}^2) \times L^2(\mathbb{R}^2)$ such that the extension satisfies $-\Delta \mathbf{v} + \nabla q = 0$ in $\omega$, we have

$$
-\Delta \mathbf{v} + \nabla q = [\sigma(\mathbf{v}, q)]_{\partial\omega} \delta_{\partial\omega} = \eta \delta_{\partial\omega} \text{ in } \mathbf{D}'(\mathbb{R}^2), \quad (A.21)
$$

where $[\cdot]_{\partial\omega}$ is the jump operator across the interface $\partial\omega$. We now set

$$
\mathbf{V} = \mathbf{E} * (\eta \delta_{\partial \omega}), \quad \mathbf{Q} = \mathbf{P} * (\eta \delta_{\partial \omega})
$$

where $(\mathbf{E}, \mathbf{P})$ is the fundamental solution of the Stokes equation given by (A.17) and $*$ is the convolution operator. The couple $(\mathbf{V}, \mathbf{Q})$ satisfies

$$
-\Delta \mathbf{V} + \nabla \mathbf{Q} = \eta \delta_{\partial \omega} \text{ in } \mathbf{D}'(\mathbb{R}^2)
$$

so that $(\mathbf{v} - \mathbf{V}, q - Q)$ is solution of the Stokes equations in $\mathbb{R}^2$. Since solutions of Stokes equations on $\mathbb{R}^2$ are given by polynomials [41], we can write

$$
\mathbf{v}(\mathbf{y}) = (\mathbf{E} * \eta \delta_{\partial \omega})(\mathbf{y}) + \mathbf{v}_1(\mathbf{y}) = \int_{\partial \omega} \mathbf{E}(\mathbf{y} - \mathbf{x}) \eta(\mathbf{x}) \mathrm{d}\mathbf{s}(\mathbf{x}) + \mathbf{v}_1(\mathbf{y}),
$$

$$
q(\mathbf{y}) = (\mathbf{P} * \eta \delta_{\partial \omega})(\mathbf{y}) + p_1(\mathbf{y}) = \int_{\partial \omega} \mathbf{P}(\mathbf{y} - \mathbf{x}) \cdot \eta(\mathbf{x}) \mathrm{d}\mathbf{s}(\mathbf{x}) + q_1(\mathbf{y}),
$$

where $v_1$ and $q_1$ are polynomial functions. Let $\mathbf{x} \in \partial\omega$, applying a first order Taylor expansion of
$\alpha \mapsto \mathbf{E}(\mathbf{y} - \alpha\mathbf{x})$, we get

$$
\mathbf{E}(\mathbf{y} - \mathbf{x}) = \mathbf{E}(\mathbf{y}) - \nabla \mathbf{E}(\mathbf{y}) : \mathbf{x} - \int_{0}^{1} (1-t) \nabla^2 \mathbf{E}(\mathbf{y} - t\mathbf{x}) : \mathbf{x}\mathbf{x}^{\mathrm{T}} \mathrm{d}\mathbf{t}.
$$

where for all $\mathbf{x}, \mathbf{y} \in \mathbb{R}^2 \setminus \bar{\omega}$ and $i, j \in \{1, 2\}$,

$$
(\nabla \mathbf{E}(\mathbf{y}) : \mathbf{x})_{ij} = \sum_{k=1}^{2} \partial_k E_{ij}(\mathbf{y}) x_k = \sum_{k=1}^{2} -\frac{1}{4\pi} \left[ \delta_{ij} \frac{y_k}{|\mathbf{y}|^2} - \frac{\delta_{ik} y_j + \delta_{jk} y_i}{|\mathbf{y}|^2} + 2 \frac{y_i y_j y_k}{|\mathbf{y}|^4} \right] x_k. \quad (A.22)
$$

$$
\begin{align}
(\nabla^2 \mathbf{E}(\mathbf{y}) : (\mathbf{x}\mathbf{x}^\mathrm{T}))_{ij} &= \sum_{k,l=1}^{2} \partial_l \partial_k E_{ij}(\mathbf{y}) x_k x_l \nonumber \\
&= \sum_{k,l=1}^{2} -\frac{1}{4\pi\mu} \left[ \delta_{ij} \left( \frac{\delta_{kl}}{|\mathbf{y}|^2} - 2 \frac{y_k y_l}{|\mathbf{y}|^4} \right) - \frac{\delta_{ik}\delta_{jl} + \delta_{jk}\delta_{il}}{|\mathbf{y}|^2} \right. \nonumber \\
&\qquad \left. + 2 \frac{\delta_{ik} y_j y_l + \delta_{jk} y_i y_l}{|\mathbf{y}|^4} + 2 \left( \frac{\delta_{il} y_j y_k + \delta_{jl} y_i y_k}{|\mathbf{y}|^4} + 4 \frac{y_i y_j y_k y_l}{|\mathbf{y}|^6} \right) \right] x_k x_l. \tag{A.23}
\end{align}
$$

Consequently,

$$
\int_{\partial \omega} \mathbf{E}(\mathbf{y} - \mathbf{x})\eta(\mathbf{x})\mathrm{d}\mathbf{s}(\mathbf{x}) = \mathbf{E}(\mathbf{y})\int_{\partial \omega} \eta(\mathbf{x})\mathrm{d}\mathbf{s}(\mathbf{x}) \\
- \int_{\partial \omega} \nabla \mathbf{E}(\mathbf{y}) : \mathbf{x}\eta(\mathbf{x})\mathrm{d}\mathbf{s}(\mathbf{x}) - \int_{\partial \omega} \int_{0}^{1} (1-t)\nabla^2 \mathbf{E}(\mathbf{y} - t\mathbf{x}) : \mathbf{x}\mathbf{x}^\mathrm{T}\eta(\mathbf{x})\mathrm{d}\mathbf{s}(\mathbf{x}). \quad (A.24)
$$
Let $R > 0$ such that $\omega \subset B_R$ and $\mu \in H^{\frac{1}{2}}(\partial\omega)$, there exists $v^\mu \in H^1(B_R)$ such that

$$v^\mu = \mu \text{ on } \partial\omega, \quad v^\mu = 0 \text{ on } \partial B_R \quad \text{and} \quad \|v^\mu\|_{1, B_R} \lesssim \|\mu\|_{\frac{1}{2}, \partial\omega}.$$

Moreover, for any $\mu \in H^{\frac{1}{2}}(\partial\omega)$, taking $v^\mu$ as a test in (A.21) and integrating in $B_R$

$$\langle \eta, \mu \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega} = (\nabla v, \nabla v^\mu)_{B_R} + (p, \operatorname{div}(v^\mu))_{B_R}.$$

From (A.20) and remembering that $(v, p)$ satisfies $-\Delta v + \nabla q = 0$ in $\omega$ and $v = \varphi(\varepsilon)$ on $\partial\omega$, we deduce

$$\begin{aligned} |\langle \eta, \mu \rangle_{-\frac{1}{2}, \frac{1}{2}, \partial\omega}| &\le (\|\nabla v\|_{B_R} + \|p\|_{B_R}) \|v^\mu\|_{1, B_R}, \\ &\lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega} \|\mu\|_{\frac{1}{2}, \partial\omega}, \end{aligned}$$

and thus

$$\|\eta\|_{-\frac{1}{2}, \partial\omega} \lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega}. \qquad (A.25)$$

Let us define

$$W(y) = \int_{\partial\omega} E(y-x) \eta(x) \, ds(x) - E(y) \int_{\partial\omega} \eta(x) \, ds(x), \qquad (A.26)$$

so that

$$v(y) = W(y) + E(y) \int_{\partial\omega} \eta(x) \, ds(x) + v_1(y).$$

The identity (A.24) together with (A.22), (A.23) and the estimate (A.25) imply that there exists A independent of $\varphi$ such that for any $y$ satisfying $|y| > A$,

$$|W(y)| \lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega} \frac{1}{|y|}. \qquad (A.27)$$

Now looking at the integrability of $v$ defined by (A.19), since $v \in H_{-1}(\mathbb{R}^2 \setminus \bar{\omega})$ it follows that $(\rho \log(\rho))^{-1} v \in L^2(\mathbb{R}^2 \setminus \bar{\omega})$. Moreover, from (A.27), $(\rho \log(\rho))^{-1} W \in L^2(\mathbb{R}^2 \setminus \bar{\omega})$. Since for any non zero constant vector $a$ and any non constant polynomial function $p$, $(\rho \log(\rho))^{-1} (Ea + p) \notin L^2(\mathbb{R}^2 \setminus \bar{\omega})$, we necessarily have

$$\int_{\partial\omega} \eta(x) ds(x) = 0 \quad \text{and} \quad v_1 = L,$$

where $L$ is a constant vector. Similar arguments for $q$ gives $q_1 = 0$. From (A.27), we deduce that for any $|y| > A$,

$$|L| \lesssim |v(y)| + \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega} \frac{1}{|y|}.$$

Hence

$$|L| \left( \int_{|y|>A} \frac{1}{\log(|y|)^2 |y|^2} \right)^{\frac{1}{2}} \lesssim \left( \int_{|y|>A} \frac{|v(y)|^2}{\log(|y|)^2 |y|^2} \right)^{\frac{1}{2}} + \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega} \left( \int_{|y|>A} \frac{1}{\log(|y|)^2 |y|^4} \right)^{\frac{1}{2}}.$$

Finally, since A is independent of $\varphi$ and estimate (A.20) holds true, we conclude that

$$|L| \lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega}.$$

Following the same lines as in [38, Lemma B.3], we also have

$$\left\| W \left( \frac{x}{\varepsilon} \right) \right\|_{1, \Omega_{\varepsilon}} \lesssim \| \varphi(\varepsilon) \|_{\frac{1}{2}, \partial \omega}.$$

We then define $z_\varepsilon = v_\varepsilon - W(\frac{x}{\varepsilon})$ and $p z_\varepsilon = q_\varepsilon - \frac{1}{\varepsilon} Q(\frac{x}{\varepsilon})$ such that $(z_\varepsilon, p z_\varepsilon)$ satisfies

$$\left\{ \begin{array}{ll} -\Delta z_{\varepsilon} + \nabla p z_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\ \operatorname{div} z_{\varepsilon} = 0 & \text{in } \Omega_{\varepsilon}, \\ z_{\varepsilon} = -W\left(\frac{x}{\varepsilon}\right) & \text{on } \partial\Omega, \\ z_{\varepsilon} = \varphi - W\left(\frac{x}{\varepsilon}\right) & \text{on } \partial\omega_{\varepsilon}. \end{array} \right.$$

Since $\varphi - W(\frac{x}{\varepsilon}) = v(\frac{x}{\varepsilon}) - W(\frac{x}{\varepsilon}) = L$ on $\partial\omega$, from Lemma A.1 and Lemma A.3, we get that

$$\begin{aligned} \|z_{\varepsilon}\|_{1, \Omega_{\varepsilon}} &\lesssim \|W\left(\frac{x}{\varepsilon}\right)\|_{\frac{1}{2}, \partial\Omega} + |\log(\varepsilon)|^{-\frac{1}{2}}|L|, \\ &\lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega} + |\log(\varepsilon)|^{-\frac{1}{2}}\|\varphi(\varepsilon)\|_{\frac{1}{2}, \partial\omega}. \end{aligned}$$
So finally, we obtain

$$
\begin{align*}
\|v_{\varepsilon}\|_{1,\Omega_{\varepsilon}} &\lesssim \|z_{\varepsilon}\|_{1,\Omega_{\varepsilon}} + \|W\left(\frac{x}{\varepsilon}\right)\|_{1,\Omega_{\varepsilon}}, \\
&\lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2},\partial\omega} + |\log(\varepsilon)|^{-\frac{1}{2}}\|\varphi(\varepsilon)\|_{\frac{1}{2},\partial\omega} + \|\varphi(\varepsilon)\|_{\frac{1}{2},\partial\omega} \lesssim \|\varphi(\varepsilon)\|_{\frac{1}{2},\partial\omega}.
\end{align*}
$$

□

A.3 Proof the preliminary lemmas for the well-posedness of Prob- lem (1.17)

In this section, we prove Lemma 2.3, 2.4 and 2.5 as well as Theorem 2.6 introduced in Section 2.1 for
the proof of Theorem 2.2.

Restated Lemma 2.3. for any $v \in H^1(\Omega)$,

$$
\|v\|_{\frac{1}{2}, \varepsilon} \lesssim |\log \varepsilon|^{\frac{1}{2}} \|v\|_{1, \Omega}.
$$

*Proof of Lemma* 2.3. According to [42, Section 4.1.3], the norm

$$
\langle u \rangle_{\frac{1}{2}, \partial \omega_{\varepsilon}} = \inf_{v \in H^1(\Omega_{\varepsilon}), v|_{\partial \omega_{\varepsilon}} = u} \|v\|_{1, \Omega_{\varepsilon}} \quad \forall u \in H^{\frac{1}{2}}(\partial \omega_{\varepsilon}), \quad (A.28)
$$

is equivalent, uniformly in $\varepsilon$, to the norm

$$
|\log(\varepsilon)|^{-\frac{1}{2}} \|\cdot\|_{\varepsilon} + |\cdot|_{\frac{1}{2}, \varepsilon}.
$$

By definition of the norm $\langle \cdot \rangle_{\frac{1}{2}, \partial \omega_\varepsilon}$, we have

$$
\langle v \rangle_{\frac{1}{2}, \partial \omega_{\varepsilon}} \lesssim \|v\|_{1, \Omega_{\varepsilon}} \quad \forall v \in H^1(\Omega_{\varepsilon}),
$$

and consequently

$$
|\log(\varepsilon)|^{-\frac{1}{2}} \|v\|_{\varepsilon} + |v|_{\frac{1}{2}, \varepsilon} \lesssim \|v\|_{1, \Omega_{\varepsilon}} \quad \forall v \in H^1(\Omega_{\varepsilon}).
$$

Eventually, using the definition of the norm $\| \cdot \|_{\frac{1}{2}, \varepsilon}$, we conclude

$$
\begin{align*}
\|v\|_{\frac{1}{2}, \varepsilon} &= \|v\|_{\varepsilon} + |v|_{\frac{1}{2}, \varepsilon} \lesssim |\log(\varepsilon)|^{\frac{1}{2}} (|\log(\varepsilon)|^{-\frac{1}{2}} \|v\|_{\varepsilon} + |v|_{\frac{1}{2}, \varepsilon}) \\
&\lesssim |\log(\varepsilon)|^{\frac{1}{2}} \|v\|_{1, \Omega_{\varepsilon}} \lesssim |\log(\varepsilon)|^{\frac{1}{2}} \|v\|_{1, \Omega} \quad \forall v \in H^1(\Omega).
\end{align*}
$$

□

We are now in position to derive lemmas 2.4 and 2.5.

**Restated Lemma 2.4.** Let $\eta \in H_0^\frac{1}{2}(\partial\omega_\varepsilon)$, there exists a unique pair $(v_\eta^\eta, l_\eta^\eta) \in V \times Q$ such that

$$
\left\{
\begin{array}{ll}
-\Delta \boldsymbol{v}_{\varepsilon}^{\eta} + \boldsymbol{\nabla} l_{\varepsilon}^{\eta} = 0 & \text{in } \Omega_{\varepsilon} \cup \omega_{\varepsilon}, \\
\text{div } \boldsymbol{v}_{\varepsilon}^{\eta} = 0 & \text{in } \Omega, \\
\boldsymbol{v}_{\varepsilon}^{\eta} = \boldsymbol{\eta} & \text{on } \partial\omega_{\varepsilon},
\end{array}
\right.
\qquad (A.29)
$$

Moreover, the following estimate holds:

$$
\| \boldsymbol{v}_{\varepsilon}^{\eta} \|_{1, \Omega} \lesssim \| \boldsymbol{\eta} \|_{\frac{1}{2}, \varepsilon}. \qquad (A.30)
$$

*Proof of Lemma* 2.4. Standard results on Stokes equations ensure the existence of $(v_{\epsilon}^{\eta}, l_{\epsilon}^{\eta})$ in $V \times Q$
satisfying (A.29). Now, we would like to prove the *a priori* estimate (A.30). On one hand, according
to Lemma A.4,

$$
\| \boldsymbol{v}_{\varepsilon}^{\eta} \|_{1, \Omega_{\varepsilon}} \lesssim \| \boldsymbol{\eta} \|_{\frac{1}{2}, \varepsilon}. \qquad (A.31)
$$

On the other hand, by defining $(v_{\epsilon}^{\eta}(\epsilon), l_{\epsilon}^{\eta}(\epsilon))$ such that $(v_{\epsilon}^{\eta}(\epsilon)(\boldsymbol{x}), l_{\epsilon}^{\eta}(\epsilon)(\boldsymbol{x})) = (v_{\epsilon}^{\eta}(\epsilon\boldsymbol{x}), l_{\epsilon}^{\eta}(\epsilon\boldsymbol{x}))$ for any $\boldsymbol{x} \in \omega$, the pair $(v_{\epsilon}^{\eta}(\epsilon), l_{\epsilon}^{\eta}(\epsilon))$ satisfies

$$
\left\{
\begin{array}{ll}
-\Delta \boldsymbol{v}_{\epsilon}^{\eta}(\epsilon) + \epsilon \nabla l_{\epsilon}^{\eta}(\epsilon) = 0 & \text{in } \omega, \\
\text{div } \boldsymbol{v}_{\epsilon}^{\eta}(\epsilon) = 0 & \text{in } \omega, \\
\boldsymbol{v}_{\epsilon}^{\eta}(\epsilon) = \boldsymbol{\eta}(\epsilon) & \text{on } \partial\omega,
\end{array}
\right.
$$
together with the a priori estimate

$$ \|\mathbf{v}_{\varepsilon}^{\eta}(\varepsilon)\|_{1,\omega} \lesssim \|\eta(\varepsilon)\|_{\frac{1}{2},\partial\omega}. $$

By rescaling, this yields

$$ \varepsilon^{-2} \|\mathbf{v}_{\varepsilon}^{\eta}\|_{\omega_{\varepsilon}}^{2} + \|\nabla \mathbf{v}_{\varepsilon}^{\eta}\|_{\omega_{\varepsilon}}^{2} \lesssim \|\eta(\varepsilon)\|_{\frac{1}{2}, \partial\omega}. $$

Consequently, we obtain

$$ \|\mathbf{v}_{\varepsilon}^{\eta}\|_{1,\omega_{\varepsilon}} \lesssim \|\eta\|_{\frac{1}{2},\varepsilon}. \qquad (A.32) $$

Gathering (A.31) and (A.32) together, we obtain inequality (A.30). □

**Restated Lemma 2.5.** Let $q \in Q$, there exists $\mathbf{w}_{\varepsilon}^{q} \in H_{0}^{1}(\Omega)$ such that

$$ \begin{cases} \operatorname{div} \mathbf{w}_{\varepsilon}^{q} = q & \text{in } \Omega, \\ \mathbf{w}_{\varepsilon}^{q} = \frac{1}{|\partial \omega_{\varepsilon}|} \left( \int_{\omega_{\varepsilon}} q \right) \mathbf{n}^{-} & \text{on } \partial \omega_{\varepsilon}, \end{cases} $$

and

$$ \|\mathbf{w}_{\varepsilon}^{q}\|_{1,\Omega} \lesssim \|q\|_{\Omega}. $$

**Proof of Lemma 2.5.** Let $\mathbf{v}_{\varepsilon}^{q} \in H_{0}^{1}(\Omega)$ be the unique solution of

$$ \begin{cases} -\Delta \mathbf{v}_{\varepsilon}^{q} = 0 & \text{in } \Omega_{\varepsilon} \cup \omega_{\varepsilon}, \\ \mathbf{v}_{\varepsilon}^{q} = \frac{1}{|\partial \omega_{\varepsilon}|} \left( \int_{\omega_{\varepsilon}} q \right) \mathbf{n}^{-} & \text{on } \partial \omega_{\varepsilon}. \end{cases} $$

According to [31, Lemma 3.3], it satisfies the estimate

$$ \|\mathbf{v}_{\varepsilon}^{q}\|_{1,\Omega} \lesssim \left| \frac{1}{|\partial \omega_{\varepsilon}|} \left( \int_{\omega_{\varepsilon}} q \right) \right| \|\mathbf{n}^{-}\|_{\frac{1}{2},\varepsilon} \lesssim \|q\|_{0,\Omega}. \qquad (A.33) $$

Since

$$ \int_{\Omega_{\varepsilon}} (q - \operatorname{div}(\mathbf{v}_{\varepsilon}^{q})) = \int_{\Omega_{\varepsilon}} q + \int_{\partial \omega_{\varepsilon}} \mathbf{v}_{\varepsilon}^{q} \cdot \mathbf{n}^{-} = \int_{\Omega_{\varepsilon}} q + \int_{\omega_{\varepsilon}} q = 0, $$

according to [43, Lemma 2.2.4], there exists $\mathbf{w}_{\varepsilon}^{1} \in H_{0}^{1}(\Omega_{\varepsilon})$ such that

$$ \operatorname{div}(\mathbf{w}_{\varepsilon}^{1}) = q - \operatorname{div}(\mathbf{v}_{\varepsilon}^{q}) \quad \text{in } \Omega_{\varepsilon}, $$

and

$$ \|\mathbf{w}_{\varepsilon}^{1}\|_{1,\Omega_{\varepsilon}} \lesssim \|q + \operatorname{div}(\mathbf{v}_{\varepsilon}^{q})\|_{\Omega_{\varepsilon}} \lesssim \|q\|_{\Omega}. $$

Moreover, in the same way $q - \operatorname{div}(\mathbf{v}_{\varepsilon}^{q}) \in L_{0}^{2}(\omega_{\varepsilon})$. Using the standard Bogovskii operator, there exists $\mathbf{w}_{\varepsilon}^{2}(\varepsilon) \in H_{0}^{1}(\omega)$ such that

$$ \operatorname{div}(\mathbf{w}_{\varepsilon}^{2}(\varepsilon)) = \varepsilon q(\varepsilon) - \operatorname{div}(\mathbf{v}_{\varepsilon}^{q}(\varepsilon)) \quad \text{in } \omega, $$

and

$$ \|\mathbf{w}_{\varepsilon}^{2}(\varepsilon)\|_{1,\omega} \lesssim \|\varepsilon q(\varepsilon) - \operatorname{div}(\mathbf{v}_{\varepsilon}^{q}(\varepsilon))\|_{\omega} \lesssim \varepsilon \|q(\varepsilon)\|_{\omega} + \|\nabla \mathbf{v}_{\varepsilon}(\varepsilon)\|_{1,\omega}. $$

After rescaling, and thanks to (A.33), we then obtain

$$ \|\mathbf{w}_{\varepsilon}^{2}\|_{1,\omega_{\varepsilon}} \lesssim \|q\|_{\omega_{\varepsilon}} + \|\nabla \mathbf{v}_{\varepsilon}^{q}\|_{1,\omega_{\varepsilon}} \lesssim \|q\|_{0,\Omega}, $$

Setting $\tilde{\mathbf{w}}_{\varepsilon}^{q} = \mathbf{w}_{\varepsilon}^{1} \mathbf{1}_{\Omega_{\varepsilon}} + \mathbf{w}_{\varepsilon}^{2} \mathbf{1}_{\omega_{\varepsilon}}$ and $\mathbf{w}_{\varepsilon}^{q} = \mathbf{v}_{\varepsilon}^{q} + \tilde{\mathbf{w}}_{\varepsilon}^{q}$ gives the desired result. □

**Restated Theorem 2.6.** Let $Q_1$ and $Q_2$ be two reflexive Banach spaces, $a: X \times X \to \mathbb{R}$, $b_1: Q_1 \times X \to \mathbb{R}$, $b_2: Q_2 \times X \to \mathbb{R}$ three continuous bilinear forms, $c_1: Q_1 \to \mathbb{R}$, $c_2: Q_2 \to \mathbb{R}$, $d: X \to \mathbb{R}$ three continuous linear forms, we consider the twofold saddle point problem: find $(u, \lambda_1, \lambda_2) \in X \times Q_1 \times Q_2$ such that

$$ \left\{ \begin{array}{ll} a(u, v) - b_1(\lambda_1, v) - b_2(\lambda_2, v) = d(v), & \forall v \in X, \\ b_1(\mu_1, u) = c_1(\mu_1), & \forall \mu_1 \in Q_1, \\ b_2(\mu_2, u) = c_2(\mu_2), & \forall \mu_2 \in Q_2. \end{array} \right. \qquad (A.34) $$

Let

$$ Z_{b_i} = \{v \in X \mid b_i(\mu_i, v) = 0 \quad \forall \mu_i \in Q_i\} \subset X, \quad i = 1, 2. \qquad (A.35) $$
We suppose that there exists $\beta_1 > 0$, $\beta_2 > 0$ such that

$$ \sup_{v \in Z_{b_2}} \frac{b_1(\lambda_1, v)}{\|v\|_X} \ge \beta_1 \|\lambda_1\|_{Q_1} \quad \forall \lambda_1 \in Q_1, \qquad (A.36) $$

$$ \sup_{v \in Z_{b_1}} \frac{b_2(\lambda_2, v)}{\|v\|_X} \ge \beta_2 \|\lambda_2\|_{Q_2} \quad \forall \lambda_2 \in Q_2, \qquad (A.37) $$

and that there exists $\alpha > 0$

$$ a(v, v) \ge \alpha \|v\|_X^2 \quad \forall v \in X. \qquad (A.38) $$

Then Problem (A.34) is well posed and we have the following estimates on $u$, $\lambda_1$ and $\lambda_2$:

$$ \|u\|_X \le \alpha^{-1} \|d\| + \beta_1^{-1} (1 + \alpha^{-1} \|a\|) \|c_1\| + \beta_2^{-1} (1 + \alpha^{-1} \|a\|) \|c_2\|, \qquad (A.39) $$

and

$$ \|\lambda_1\|_{Q_1} \le \beta_1^{-1} (\|d\| + \|a\| \|u\|_X), \quad \|\lambda_2\|_{Q_2} \le \beta_2^{-1} (\|d\| + \|a\| \|u\|_X). \qquad (A.40) $$

Proof of Theorem 2.6. Let $Q = Q_1 \times Q_2$ and $b: Q \times X \to \mathbb{R}$ the continuous bilinear form defined for $\lambda = (\lambda_1, \lambda_2) \in Q$ and $v \in X$ by

$$ b(\lambda, v) = b_1(\lambda_1, v) + b_2(\lambda_2, v). $$

First we prove that conditions (A.36) and (A.37) implies that the bilinear form $b$ satisfies an inf-sup condition. Let $\lambda = (\lambda_1, \lambda_2) \in Q$, assumptions (A.36) and (A.37) imply that there exists $(v_1, v_2) \in Z_{b_2} \times Z_{b_1}$ such that $\|v_1\|_X = \|v_2\|_X = 1$ and

$$ \frac{b_1(\lambda_1, v_1)}{\|v_1\|_X} \ge \frac{\beta_1}{2} \|\lambda_1\|_{Q_1}, \quad \frac{b_2(\lambda_2, v_2)}{\|v_2\|_X} \ge \frac{\beta_2}{2} \|\lambda_2\|_{Q_2}, $$

Defining $v = v_1 + v_2$, it satisfies $\|v\|_X \le 2$ and

$$ \frac{b(\lambda, v)}{\|v\|_X} \ge \frac{1}{4} (\beta_1 \|\lambda_1\|_{Q_1} + \beta_2 \|\lambda_2\|_{Q_2}) \ge \frac{1}{4} \min\{\beta_1, \beta_2\} \|\lambda\|_Q. $$

In particular, there exists $\beta > 0$ such that

$$ \sup_{v \in X} \frac{b(\lambda, v)}{\|v\|_X} \ge \beta \|\lambda\|_Q. $$

We then deduce that Problem (A.34) admits a unique solution $(u, (\lambda_1, \lambda_2))$ in $X \times Q$. Let us now prove the a priori estimates (A.39) and (A.40). First the conditions (A.36) and (A.37) imply that there exists $(w_1, w_2) \in Z_{b_1} \times Z_{b_2}$ such that

$$ \begin{aligned} b_1(\mu_1, w_1) &= c_1(\mu_1) \quad \forall \mu_1 \in Q_1, \quad \|w_1\|_X \le \beta_1^{-1} \|c_1\|, \\ b_2(\mu_2, w_2) &= c_2(\mu_2) \quad \forall \mu_2 \in Q_2, \quad \|w_2\|_X \le \beta_2^{-1} \|c_2\|. \end{aligned} $$

Let us define

$$ \kappa = u - w \qquad (A.41) $$

with $w = w_1 + w_2$. We have, for any $(\mu_1, \mu_2) \in Q_1 \times Q_2$,

$$ \begin{aligned} b_1(\mu_1, \kappa) &= b_1(\mu_1, w_1) + b_1(\mu_1, w_2) - b_1(\mu_1, u) = c_1(\mu_1) - c_1(\mu_1) = 0, \\ b_2(\mu_2, \kappa) &= b_2(\mu_2, w_1) + b_2(\mu_2, w_2) - b_2(\mu_2, u) = c_2(\mu_2) - c_2(\mu_2) = 0. \end{aligned} $$

We then deduce that $\kappa \in Z_{b_1} \cap Z_{b_2}$. Furthermore, we have, for any $v \in Z_{b_1} \cap Z_{b_2}$,

$$ a(\kappa, v) = a(u, v) - a(w, v) = d(v) - a(w, v). $$

The continuity assumption on $a$ and $d$, as well as the coercivity assumption (A.38) on $a$, imply

$$ \alpha \|\kappa\|_X \le \frac{a(\kappa, \kappa)}{\|\kappa\|_X} = \frac{d(\kappa) - a(w, \kappa)}{\|\kappa\|_X} \le \|d\| + \|a\| \|w\|_X. $$

So, we get

$$ \|\kappa\|_X \le \alpha^{-1} (\|d\| + \|a\| \|w\|_X). $$

From the definition (A.41) and the estimates satisfied by $w_1$ and $w_2$, we obtain

$$ \|u\|_X \le \|\kappa\|_X + \|w\|_X \le \alpha^{-1} \|d\| + \beta_1^{-1} (\alpha^{-1} \|a\| + 1) \|c_1\| + \beta_2^{-1} (\alpha^{-1} \|a\| + 1) \|c_2\|. $$
which corresponds to (A.39). To obtain the estimates on the Lagrange multipliers we take $v \in Z_{b_1}$ in the first equation of system (A.34), to obtain

$$b_2(\lambda_2, v) = d(v) - a(u, v).$$

By continuity of $a$ and $d$ and the inf-sup condition (A.37) on $b_2$, we get

$$\beta_2 \|\lambda_2\|_{Q_2} \le \sup_{v \in Z_{b_1}, v \ne 0} \frac{b_2(\lambda_2, v)}{\|v\|_X} = \sup_{v \in Z_{b_1}, v \ne 0} \frac{d(v) - a(u, v)}{\|v\|_X} \le \|d\| + \|a\| \|u\|_X.$$

We thus obtain

$$\|\lambda_2\|_{Q_2} \le \beta_2^{-1}(\|d\| + \|a\|\|u\|_X).$$

Finally, we have in a similar way

$$\|\lambda_1\|_{Q_1} \le \beta_1^{-1}(\|d\| + \|a\|\|u\|_X).$$

□

## A.4 Preliminary lemmas for the convergence of Problem (1.17)

**Restated Lemma 2.8.** Let $(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)}, \lambda_{\varepsilon}^{(N)}) \in V \times Q \times F_N$ be the solution of Problem (2.16). We have in $H^{-\frac{1}{2}}(\partial\omega_{\varepsilon})$

$$\lambda_{\varepsilon}^{(N)} = \varepsilon(\sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{+} n^{+} + \sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{-} n^{-}), \quad (A.42)$$

where $n^+$ and $n^-$ are the exterior normal vectors to $\partial\Omega_{\varepsilon}$ and $\partial\omega_{\varepsilon}$ respectively, and $\sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^+$ and $\sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^-$ are the restrictions of $\sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})$ to $\Omega_{\varepsilon}$ and $\omega_{\varepsilon}$, respectively.

*Proof of Lemma 2.8.* On the one hand, integrating by parts, we have for any $\phi \in H_0^1(\Omega)$,

$$\begin{align*}
2 \int_{\Omega} \varepsilon(u_{\varepsilon}^{(N)}) : \varepsilon(\phi) - \int_{\Omega} p_{N}^{\varepsilon} \operatorname{div} \phi \\
&= 2 \int_{\Omega_{\varepsilon}} \varepsilon(u_{\varepsilon}^{(N)}) : \varepsilon(\phi) + 2 \int_{\omega_{\varepsilon}} \varepsilon(u_{\varepsilon}^{(N)}) : \varepsilon(\phi) - \int_{\Omega_{\varepsilon}} p_{\varepsilon}^{(N)} \operatorname{div} \phi - \int_{\omega_{\varepsilon}} p_{\varepsilon}^{(N)} \operatorname{div} \phi, \\
&= - \int_{\Omega_{\varepsilon}} \Delta u_{\varepsilon}^{(N)} \cdot \phi - \int_{\omega_{\varepsilon}} \Delta u_{\varepsilon}^{(N)} \cdot \phi + \int_{\Omega_{\varepsilon}} \nabla p_{\varepsilon}^{(N)} \cdot \phi + \int_{\omega_{\varepsilon}} \nabla p_{\varepsilon}^{(N)} \cdot \phi \\
&\quad + \langle (\sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{+} n^{+} + \sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{-} n^{-}, \phi) \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}}.
\end{align*}$$

Since

$$-\Delta u_{\varepsilon}^{(N)} + \nabla p_{\varepsilon}^{(N)} = 0 \text{ in } \Omega_{\varepsilon}, \quad -\Delta u_{\varepsilon}^{(N)} + \nabla p_{\varepsilon}^{(N)} = 0 \text{ in } \omega_{\varepsilon},$$

we deduce that

$$2 \int_{\Omega} \varepsilon(u_{\varepsilon}^{(N)}) : \varepsilon(\phi) - \int_{\Omega} p_{\varepsilon}^{(N)} \operatorname{div} \phi = \langle \sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{+} n^{+} + \sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{-} n^{-}, \phi \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}}$$

On the other hand, we have

$$2 \int_{\Omega} \varepsilon(u_{\varepsilon}^{(N)}) : \varepsilon(\phi) - \int_{\Omega} p_{\varepsilon}^{(N)} \operatorname{div} \phi = (\lambda_{\varepsilon}^{(N)}, \phi)_{\varepsilon} \quad \forall \phi \in H_{0}^{1}(\Omega).$$

Identifying these two relations yields

$$\begin{align*}
(\lambda_{\varepsilon}^{(N)}, \phi)_{\varepsilon} &= \varepsilon^{-1}(\lambda_{\varepsilon}^{(N)}, \phi)_{\partial\omega_{\varepsilon}} \\
&= \langle \sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{+} n^{+} + \sigma(u_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})^{-} n^{-}, \phi \rangle_{-\frac{1}{2}, \partial\omega_{\varepsilon}}
\end{align*}$$

for any $\phi \in H_0^1(\Omega)$. Finally, since the trace operator mapping $H_0^1(\Omega)$ to $H^{\frac{1}{2}}(\partial\omega_{\varepsilon})$ is surjective, we get the desired relation (A.42). □

**Restated Lemma 2.9.** Let $u_{\varepsilon}^{(N)} \in V$ be the solution of Problem (2.16) and let $u_b = v_b + \varepsilon w_b e_{\theta}$ with $v_b \in \mathbb{R}^3$ and $w_b \in \mathbb{R}$. There exists $a_{k,r}^{\varepsilon}, a_{k,\theta}^{\varepsilon}, b_{k,r}^{\varepsilon}, b_{k,\theta}^{\varepsilon} \in \mathbb{R}$ such that for $\varepsilon \in (0, \rho)$,

$$\begin{align}
(u_{\varepsilon}^{(N)}(\varepsilon)(\theta) - u_b(\varepsilon)(\theta)) \cdot e_r &= \left(\frac{\varepsilon}{\rho}\right)^{N+1} \sum_{k=N}^{\infty} (a_{k,r}^{\varepsilon} \cos(k\theta) + b_{k,r}^{\varepsilon} \sin(k\theta)), \tag{A.43} \\
(u_{\varepsilon}^{(N)}(\varepsilon)(\theta) - u_b(\varepsilon)(\theta)) \cdot e_{\theta} &= \left(\frac{\varepsilon}{\rho}\right)^{N+1} \sum_{k=N}^{\infty} (a_{k,\theta}^{\varepsilon} \cos(k\theta) + b_{k,\theta}^{\varepsilon} \sin(k\theta)). \nonumber
\end{align}$$
Moreover, there exists $0 < \tilde{\rho} < \rho$ such that for $\varepsilon \in (0, \tilde{\rho})$, the following estimate holds:

$$ \left( \sum_{k=N}^{\infty} (1+k)(|a_{k,r}^{\varepsilon}|^2 + |a_{k,\theta}^{\varepsilon}|^2 + |b_{k,r}^{\varepsilon}|^2 + |b_{k,\theta}^{\varepsilon}|^2) \right)^{\frac{1}{2}} \lesssim (1+N)(|v_b| + \varepsilon|w_b|). \quad (A.44) $$

**Proof of Lemma 2.9.** Since div $\mathbf{u}_{\varepsilon}^{(N)} = 0$ in $\Omega$, there exists a scalar function $\psi_{\varepsilon}^{(N)}$, called the stream function, such that

$$ \mathbf{u}_{\varepsilon}^{(N)} = \nabla^{\perp} \psi_{\varepsilon}^{(N)} \stackrel{\text{def}}{=} \left( \frac{1}{r} \partial_{\theta} \psi_{\varepsilon}^{(N)} \right) \mathbf{e}_r - \left( \partial_r \psi_{\varepsilon}^{(N)} \right) \mathbf{e}_{\theta}. \quad (A.45) $$

In each subdomain $\omega_{\rho} \setminus \bar{\omega}_{\varepsilon}$ and $\omega_{\varepsilon}$, the pair $(\mathbf{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})$ satisfies the homogeneous Stokes equations, so the stream function satisfies the biharmonic equation

$$ -\Delta^2 \psi_{\varepsilon}^{(N)} = 0 \text{ in } \omega_{\rho} \setminus \bar{\omega}_{\varepsilon}, \quad -\Delta^2 \psi_{\varepsilon}^{(N)} = 0 \text{ in } \omega_{\varepsilon}. $$

In polar coordinates, solutions of the biharmonic equations can be constructed by separation of variables. This yields the general expansion

$$ \psi_{\varepsilon}^{(N)}(r, \theta) = \Phi_0(r) + \sum_{k=1}^{\infty} \Phi_k^c(r) \cos(k\theta) + \Phi_k^s(r) \sin(k\theta), \quad (A.46) $$

where the restriction of the radial functions $\Phi_k^c$ and $\Phi_k^s$ to each subdomain $\omega_{\varepsilon}$ and $\omega_{\rho} \setminus \bar{\omega}_{\varepsilon}$ satisfy, for all $k \ge 2$,

$$ \left\{ \begin{array}{l} \Phi_0 \in \operatorname{span}\{r^2, \ln(r), r, r^2 \ln(r)\}, \quad \Phi_1^c, \Phi_1^s \in \operatorname{span}\{r^3, r, r \ln(r), r^{-1}\}, \\ \Phi_k^c, \Phi_k^s \in \operatorname{span}\{r^{k+2}, r^k, r^{-k+2}, r^{-k}\}. \end{array} \right. $$

We refer to [44] for the derivation of (A.46). Substituting (A.46) in (A.45) and differentiating term by term gives the following expansions for $\mathbf{u}_{\varepsilon}^{(N)}$:

$$ \begin{aligned} \mathbf{u}_{\varepsilon}^{(N)}(r, \theta) \cdot \mathbf{e}_r = & -(c_{1,1}^{\oplus} r^2 + c_{2,1}^{\oplus} \sin(\theta) + (s_{1,1}^{\oplus} r^2 + s_{2,1}^{\oplus} \cos(\theta) \\ & + \sum_{k=2}^{\infty} -k(c_{1,k}^{\oplus} r^{k+1} + c_{2,k}^{\oplus} r^{k-1}) \sin(k\theta) + \sum_{k=2}^{\infty} k(s_{1,k}^{\oplus} r^{k+1} + s_{2,k}^{\oplus} r^{k-1}) \cos(k\theta) \end{aligned} \quad (A.47) $$

$$ \begin{aligned} \mathbf{u}_{\varepsilon}^{(N)}(r, \theta) \cdot \mathbf{e}_{\theta} = & -2c_{1,0}^{\oplus}r - c_{3,0}^{\oplus}r(1 + 2\ln(r)) - (3c_{1,1}^{\oplus}r^2 + c_{2,1}^{\oplus}\cos(\theta) - (3s_{1,1}^{\oplus}r^2 + s_{2,1}^{\oplus}\sin(\theta) \\ & - \sum_{k=2}^{\infty} ((k+2)c_{1,k}^{\oplus}r^{k+1} + kc_{2,k}^{\oplus}r^{k-1})\cos(k\theta) - \sum_{k=2}^{\infty} ((k+2)s_{1,k}^{\oplus}r^{k+1} + ks_{2,k}^{\oplus}r^{k-1})\sin(k\theta) \end{aligned} \quad (A.48) $$

for any $(r, \theta) \in [0, \varepsilon] \times (0, 2\pi)$ and

$$ \begin{aligned} \mathbf{u}_{\varepsilon}^{(N)}(r, \theta) \cdot \mathbf{e}_r = & -(c_{1,1}^{\oplus} r^2 + c_{2,1}^{\oplus} + c_{3,1}^{\oplus} \ln(r) + c_{4,1}^{\oplus} r^{-2}) \sin(\theta) \\ & + (s_{1,1}^{\oplus} r^2 + s_{2,1}^{\oplus} + s_{3,1}^{\oplus} \ln(r) + s_{4,1}^{\oplus} r^{-2}) \cos(\theta) \\ & - \sum_{k=2}^{\infty} k(c_{1,k}^{\oplus} r^{k+1} + c_{2,k}^{\oplus} r^{k-1} + c_{3,k}^{\oplus} r^{-k+1} + c_{4,k}^{\oplus} r^{-k-1}) \sin(k\theta) \\ & + \sum_{k=2}^{\infty} k(s_{1,k}^{\oplus} r^{k+1} + s_{2,k}^{\oplus} r^{k-1} + s_{3,k}^{\oplus} r^{-k+1} + s_{4,k}^{\oplus} r^{-k-1}) \cos(k\theta), \end{aligned} \quad (A.49) $$

$$ \begin{aligned} \mathbf{u}_{\varepsilon}^{(N)}(r, \theta) \cdot \mathbf{e}_{\theta} = & -(2c_{1,0}^{\oplus}r + c_{2,0}^{\oplus}r^{-1} + c_{3,0}^{\oplus}r(1 + 2\ln(r))) \\ & - (3c_{1,1}^{\oplus}r^2 + c_{2,1}^{\oplus} + c_{3,1}^{\oplus}(\ln(r) + 1) - c_{4,1}^{\oplus}r^{-2}) \cos(\theta) \\ & - (3s_{1,1}^{\oplus}r^2 + s_{2,1}^{\oplus} + s_{3,1}^{\oplus}(\ln(r) + 1) - s_{4,1}^{\oplus}r^{-2}) \sin(\theta) \\ & - \sum_{k=2}^{\infty} ((k+2)c_{1,k}^{\oplus}r^{k+1} + kc_{2,k}^{\oplus}r^{k-1}) \cos(k\theta) \\ & - \sum_{k=2}^{\infty} ((-k+2)c_{3,k}^{\oplus}r^{-k+1} - kc_{4,k}^{\oplus}r^{-k-1}) \cos(k\theta) \\ & - \sum_{k=2}^{\infty} ((k+2)s_{1,k}^{\oplus}r^{k+1} + ks_{2,k}^{\oplus}r^{k-1}) \sin(k\theta) \\ & - \sum_{k=2}^{\infty} ((-k+2)s_{3,k}^{\oplus}r^{-k+1} - ks_{4,k}^{\oplus}r^{-k-1}) \sin(k\theta) \end{aligned} \quad (A.50) $$

for any $(r, \theta) \in [\varepsilon, \rho] \times (0, 2\pi)$.
In what follows, we focus on the coefficients $c_{j,k}^{\ominus}, c_{j,k}^{\oplus}$; the coefficients $s_{j,k}^{\ominus}, s_{j,k}^{\oplus}$ can be treated in an analogous way due to symmetry. The proof proceeds as follows. First, we derive a set of equations for the coefficients $c_{j,k}^{\ominus}, c_{j,k}^{\oplus}$ from the interface conditions on $\partial\omega_{\varepsilon}$,

$$ (\boldsymbol{\mu}_N, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = (\boldsymbol{\mu}_N, \boldsymbol{v}_b + \varepsilon w_b \boldsymbol{e}_{\theta})_{\varepsilon} \quad \forall \boldsymbol{\mu}_N \in \mathbf{F}_N, \qquad (\text{A.51}) $$

$$ [\sigma(\boldsymbol{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})]_{\partial\omega_{\varepsilon}} = \boldsymbol{\lambda}_{N}, \quad \boldsymbol{\lambda}_{N} \in \mathbf{F}_{N}. \qquad (\text{A.52}) $$

More precisely, for Fourier modes belonging to $\mathbf{F}_N$, we use the kinematic condition (A.51), whereas for modes orthogonal to $\mathbf{F}_N$, we apply (A.52), which implies continuity of the stress at the interface for those modes. Next, we use the uniqueness of the Fourier expansion on $\partial\omega_{\varepsilon}$. Denoting by + and – the restrictions of $\boldsymbol{u}_{\varepsilon}^{(N)}$ to $\omega_{\rho}\setminus\bar{\omega}_{\varepsilon}$ and $\omega_{\varepsilon}$, respectively, we have

$$ \mathcal{F}_k^j((\boldsymbol{u}_{\varepsilon}^{(N)}(\varepsilon))^{-} \cdot \boldsymbol{e}_r) = \mathcal{F}_k^j((\boldsymbol{u}_{\varepsilon}^{(N)}(\varepsilon))^{+} \cdot \boldsymbol{e}_r), \quad \mathcal{F}_k^j((\boldsymbol{u}_{\varepsilon}^{(N)}(\varepsilon))^{-} \cdot \boldsymbol{e}_{\theta}) = \mathcal{F}_k^j((\boldsymbol{u}_{\varepsilon}^{(N)}(\varepsilon))^{+} \cdot \boldsymbol{e}_{\theta}), \quad (\text{A.53}) $$

for all $k \ge 0$ and $j \in \{c, s\}$, where the operators $\mathcal{F}_k^j: (0, L) \to \mathbb{R}$ are given in (2.5). Finally, to close the system, we impose the boundary conditions on $\partial\omega_{\rho}$. For any $k \ge 0$ and $j \in \{c, s\}$, these are expressed as

$$ \mathcal{F}_k^j(\boldsymbol{u}_{\varepsilon}^{(N)}(\rho) \cdot \boldsymbol{e}_r) = \xi_{k,r}^{\varepsilon,j}, \quad \mathcal{F}_k^j(\boldsymbol{u}_{\varepsilon}^{(N)}(\rho) \cdot \boldsymbol{e}_{\theta}) = \xi_{k,\theta}^{\varepsilon,j} \qquad (\text{A.54}) $$

where $\xi_{k,r}^{\varepsilon,j}, \xi_{k,\theta}^{\varepsilon,j}$ are given scalars. In particular, the coefficients $c_{j,k}^{\ominus}, c_{j,k}^{\oplus}$ will be expressed in terms of the data $\boldsymbol{v}_b, w_b$ and the Fourier coefficients $\xi_{k,r}^{\varepsilon,j}, \xi_{k,\theta}^{\varepsilon,j}$. We finally use the a priori estimates on $\boldsymbol{u}_{\varepsilon}^{(N)}$ provided by Theorem 2.2 to establish the desired bounds (A.44) in terms of the given data $\boldsymbol{v}_b$ and $w_b$.

We start the analysis considering the reduced kinematic conditions (A.51). By writing these conditions for each basis element of the finite-dimensional space $\mathbf{F}_N$ defined in (1.19), we obtain the following set of equations:

$$ \star k = 0: \qquad (\boldsymbol{e}_x, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = \boldsymbol{v}_b \cdot \boldsymbol{e}_x, \quad (\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = \boldsymbol{v}_b \cdot \boldsymbol{e}_y, \qquad (\text{A.55}) $$

$$ \star k = 1: \qquad \left\{ \begin{array}{l} (\sin(\theta)\boldsymbol{e}_x - \cos(\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = -\varepsilon w_b, \\ (\sin(\theta)\boldsymbol{e}_x + \cos(\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = (\cos(\theta)\boldsymbol{e}_x - \sin(\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0, \end{array} \right. \qquad (\text{A.56}) $$

$$ \star k \in \{2, \dots, N\}: \qquad \left\{ \begin{array}{l} (\cos(k\theta)\boldsymbol{e}_x, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = (\cos(k\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0 \\ (\sin(k\theta)\boldsymbol{e}_x, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = (\sin(k\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0. \end{array} \right. \qquad (\text{A.57}) $$

Since the expressions (A.47)-(A.48) and (A.49)-(A.50) of $\boldsymbol{u}_N^\varepsilon$ are given in polar coordinates, we rewrite the system (A.57) in polar coordinates. To do so, we notice that for $k \in \{2, \dots, N\}$, system (A.57) is equivalent to

$$ \left\{ \begin{array}{l} (\cos(k\theta)\boldsymbol{e}_x + \sin(k\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0, \\ (\sin(k\theta)\boldsymbol{e}_x - \cos(k\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0, \\ (\sin(k\theta)\boldsymbol{e}_x + \cos(k\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0, \\ (\cos(k\theta)\boldsymbol{e}_x - \sin(k\theta)\boldsymbol{e}_y, \boldsymbol{u}_{\varepsilon}^{(N)})_{\varepsilon} = 0, \end{array} \right. $$

such that the kinematic coupling condition (A.51) writes equivalently in polar coordinates

$$ \star k = 0: \qquad \left\{ \begin{array}{l} (\cos(\theta)\boldsymbol{e}_r - \sin(\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = \boldsymbol{v}_b \cdot \boldsymbol{e}_x, \\ (\sin(\theta)\boldsymbol{e}_r + \cos(\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = \boldsymbol{v}_b \cdot \boldsymbol{e}_y, \end{array} \right. \qquad (\text{A.58}) $$

$$ \star k = 1: \qquad \left\{ \begin{array}{l} (-\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = -\varepsilon w_b, \\ (\sin(2\theta)\boldsymbol{e}_r + \cos(2\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = (\cos(2\theta)\boldsymbol{e}_r - \sin(2\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = 0, \end{array} \right. \qquad (\text{A.59}) $$

$$ \star k \in \{2, \dots, N\}: \qquad \left\{ \begin{array}{l} (\cos((k-1)\theta)\boldsymbol{e}_r + \sin((k-1)\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = 0, \\ (\sin((k-1)\theta)\boldsymbol{e}_r - \cos((k-1)\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = 0, \\ (\sin((k+1)\theta)\boldsymbol{e}_r + \cos((k+1)\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = 0, \\ (\cos((k+1)\theta)\boldsymbol{e}_r - \sin((k+1)\theta)\boldsymbol{e}_\theta, \boldsymbol{u}_\varepsilon^{(N)})_\varepsilon = 0. \end{array} \right. \qquad (\text{A.60}) $$

Writing the interface conditions in terms of the expansions (A.47)-(A.48) then yields the following
equations on the coefficients $c_{j,k}^{\ominus}$:

$$ \bullet \quad N = 0: \quad -4\epsilon^2 c_{1,1}^{\ominus} - 2c_{2,1}^{\ominus} = v_b \cdot e_y, \qquad (A.61) $$

$$ \bullet \quad N = 1: \quad -2c_{1,0}^{\ominus} - c_{3,0}^{\ominus} (1 + 2 \ln(\epsilon)) = w_b, \quad -4\epsilon^2 c_{1,1}^{\ominus} - 2c_{2,1}^{\ominus} = v_b \cdot e_y, \qquad (A.62) $$

$$ \bullet \quad N \ge 2: \quad \left\{ \begin{array}{ll} \star k = 0: & -2c_{1,0}^{\ominus} - c_{3,0}^{\ominus} (1 + 2 \ln(\epsilon)) = w_b \\ \star k = 1: & c_{1,1}^{\ominus} = 0, \quad -2c_{2,1}^{\ominus} = v_b \cdot e_y, \\ \star k \in \{2, \dots, N-1\}: & c_{1,k}^{\ominus} = c_{2,k}^{\ominus} = 0, \\ \star k \in \{N, N+1\}: & -(k+1)\epsilon^2 c_{1,k}^{\ominus} - kc_{2,k}^{\ominus} = 0. \end{array} \right. \qquad (A.63) $$

Equations (A.63) can be understood as follows. For $k \in \{2, \dots, N-1\}$, from (A.60) it follows that $\mathbf{u}_{\epsilon}^{N}$ satisfies

$$ (\cos(k\theta)e_r, \mathbf{u}_{\epsilon}^{(N)})_{\epsilon} = (\cos(k\theta)e_{\theta}, \mathbf{u}_{\epsilon}^{(N)})_{\epsilon} = (\sin(k\theta)e_r, \mathbf{u}_{\epsilon}^{(N)})_{\epsilon} = (\sin(k\theta)e_{\theta}, \mathbf{u}_{\epsilon}^{(N)})_{\epsilon} = 0. $$

This directly implies $c_{1,k}^{\ominus} = c_{2,k}^{\ominus} = 0$. On the other hand, for $\{k \in N, N+1\}$, only the two conditions

$$ (\sin((k\theta)e_r + \cos(k\theta)e_{\theta}, \mathbf{u}_{\epsilon}^{(N)})_{\epsilon} = (\cos(k\theta)e_r - \sin(k\theta)e_{\theta}, \mathbf{u}_{\epsilon}^{(N)})_{\epsilon} = 0, $$

apply, leading to a linear relation between $c_{1,k}^{\ominus}$ and $c_{2,k}^{\ominus}$. We proceed by applying the relation (A.52), which follows from Lemma 2.8. More precisely, since the stress jump belongs to $F_N$ by (A.52), its projection onto orthogonal modes to $F_N$ vanishes. To express the stress jump only in terms of the coefficients $c_{k,j}^{\ominus}, c_{k,j}^{\oplus}$ (and $s_{k,j}^{\ominus}, s_{k,j}^{\oplus}$), we use the relations

$$ \Delta \mathbf{u}_{\epsilon}^{(N)} = \nabla p_{\epsilon}^{(N)} \text{ in } \omega_{\rho} \setminus \bar{\omega}_{\epsilon}, \quad \Delta \mathbf{u}_{\epsilon}^{(N)} = \nabla p_{\epsilon}^{(N)} \text{ in } \omega_{\epsilon}. \qquad (A.64) $$

Taking the divergence of (A.64) shows that $p_{\epsilon}^{(N)}$ is harmonic in both $\omega_{\epsilon}$ and $\omega_{\rho} \setminus \bar{\omega}_{\epsilon}$. Using the method of separation of variables, we then obtain the general form of $p_{\epsilon}^{(N)}$:

$$ p_{\epsilon}^{(N)}(r, \theta) = \gamma_{1,k}^{\ominus} + \sum_{k=1}^{\infty} \gamma_{1,k}^{\ominus} r^k \cos(k\theta) + \sigma_{1,k}^{\ominus} r^k \sin(k\theta) \qquad (A.65) $$

for all $(r, \theta) \in [0, \epsilon) \times (0, 2\pi)$, and

$$ p_{\epsilon}^{(N)}(r, \theta) = \gamma_{1,0}^{\oplus} + \gamma_{2,0}^{\oplus} \ln(r) + \sum_{k=1}^{\infty} (\gamma_{1,k}^{\oplus} r^k + \gamma_{2,k}^{\oplus} r^{-k}) \cos(k\theta) + (\sigma_{1,k}^{\oplus} r^k + \sigma_{2,k}^{\oplus} r^{-k}) \sin(k\theta) \quad (A.66) $$

for all $(r, \theta) \in (\epsilon, \rho] \times (0, 2\pi)$. From (A.64), the coefficients $\sigma_{j,k}^{\ominus}$ (resp. $\sigma_{j,k}^{\oplus}$) can be expressed in terms of $c_{j,k}^{\ominus}$ (resp. $c_{j,k}^{\oplus}$) as follows:

$$ \begin{aligned} c_{3,0}^{\ominus} &= 0, & c_{3,0}^{\oplus} &= 0, & \gamma_{2,0}^{\ominus} &= 0, & \\ \sigma_{1,1}^{\ominus} &= -8c_{1,1}^{\ominus}, & \sigma_{1,1}^{\oplus} &= -8c_{1,1}^{\oplus}, & \sigma_{2,1}^{\ominus} &= 2c_{3,1}^{\oplus}, & \sigma_{1,k}^{\ominus} &= -4(k+1)c_{1,k}^{\ominus}, & \sigma_{1,k}^{\oplus} &= -4(k+1)c_{1,k}^{\oplus}, & \sigma_{2,k}^{\ominus} &= -4(k-1)c_{3,k}^{\oplus}, & \sigma_{2,k}^{\oplus} &= -4(k-1)c_{3,k}^{\ominus}, & \forall k \ge 2. \end{aligned} \qquad (A.67) $$

Substituting (A.67) into (A.65) and (A.66) provides an explicit expression for the pressure jump across the interface $\partial\omega_{\epsilon}$ in terms of the coefficients $c_{j,k}^{\ominus}$ and $c_{j,k}^{\oplus}$. Similarly, differentiating (A.65) and (A.66) yields an expression for the gradient jump across the interface. Projecting (A.52) onto modes orthogonal to $F_N$ gives the following relations:

$$ \bullet \quad \{N = 0, k = 0\}: \qquad c_{1,0}^{\ominus} = c_{1,0}^{\oplus}, \qquad (A.68) $$

$$ \bullet \quad \{N \in \{0, 1\}, k = 1\}: \qquad 2\epsilon^2 c_{1,1}^{\ominus} = 2\epsilon^2 c_{1,1}^{\oplus} - c_{3,1}^{\ominus} + 2\epsilon^{-2} c_{4,1}^{\oplus} \qquad (A.69) $$

$$ \bullet \quad \{N \ge 1, k \ge 2 \cap k \in \{N, N+1\}\}: \qquad 4(k+1)\epsilon^{k+1}c_{1,k}^{\oplus} + 4(k-1)(k+1)\epsilon^{-k+1}c_{3,k}^{\oplus} + 4k(k+1)\epsilon^{-k-1}c_{4,k}^{\oplus} \\ = 4(k+1)\epsilon^{k+1}c_{1,k}^{\ominus}. \qquad (A.70) $$
• {$N \ge 0, k \ge N + 2$}:

$$
\left\{
\begin{array}{l}
(k+1)(k-2)\varepsilon^{k+1} c_{1,k}^{\oplus} + k(k-1)\varepsilon^{k-1} c_{2,k}^{\oplus} \\
\quad - (k+2)(k-1)c_{3,k}^{\oplus}\varepsilon^{-k+1} - k(k+1)c_{4,k}^{\oplus}\varepsilon^{-k-1} \\
\quad = (k+1)(k-2)\varepsilon^{k+1} c_{1,k}^{\ominus} + k(k-1)\varepsilon^{k-1} c_{2,k}^{\ominus} \\
k(k+1)\varepsilon^{k+1} c_{1,k}^{\oplus} + k(k-1)\varepsilon^{k-1} c_{2,k}^{\oplus} \\
\quad + k(k-1)c_{3,k}^{\oplus}\varepsilon^{-k+1} + k(k+1)c_{4,k}^{\oplus}\varepsilon^{-k-1} \\
\quad = k(k+1)\varepsilon^{k+1} c_{1,k}^{\ominus} + k(k-1)\varepsilon^{k-1} c_{2,k}^{\ominus}
\end{array}
\right.
\qquad (A.71)
$$

Notice that, similarly to what explained for equation (A.63), for $k \in \{N, N+1\}$, the projection of the jump stress on the orthogonal mode of $F_N$ only implies

$$
(\sin(k\theta)e_r - \cos(k\theta)e_\theta, [\boldsymbol{\sigma}(\mathbf{u}_\varepsilon^{(N)}, p_\varepsilon^{(N)})])_{\partial\omega_\varepsilon} = (\cos(k\theta)e_r + \sin(k\theta)e_\theta, [\boldsymbol{\sigma}(\mathbf{u}_\varepsilon^{(N)}, p_\varepsilon^{(N)})])_\varepsilon = 0,
$$

instead of

$$
(\cos(k\theta)e_r, [\boldsymbol{\sigma}(\mathbf{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})])_{\partial\omega_{\varepsilon}} = (\cos(k\theta)e_{\theta}, [\boldsymbol{\sigma}(\mathbf{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})])_{\varepsilon} = 0,
$$

$$
(\sin(k\theta)e_r, [\boldsymbol{\sigma}(\mathbf{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})])_{\partial\omega_{\varepsilon}} = (\sin(k\theta)e_{\theta}, [\boldsymbol{\sigma}(\mathbf{u}_{\varepsilon}^{(N)}, p_{\varepsilon}^{(N)})])_{\varepsilon} = 0
$$

for $k \ge N + 2$.

We continue by expressing the boundary conditions on $\partial\omega_\rho$ given in (A.54). By evaluating equa-
tions (A.49)-(A.50) for $r = \rho$ and subsequently projecting the resulting expressions onto the different
Fourier modes, we obtain the following result.

$$
\star \quad k = 0 : - (2c_{1,0}^{\oplus} \rho + c_{2,0}^{\oplus} \rho^{-1} + c_{3,0}^{\oplus} \rho (1 + 2 \ln(\rho))) = \xi_{0,\theta}^{\varepsilon,c}, \quad (A.72)
$$

$$
\begin{align}
\star \quad k = 1: & \left\{ -(c_{1,1}^{\oplus} \rho^2 + c_{2,1}^{\oplus} + c_{3,1}^{\oplus} \ln(\rho) + c_{4,1}^{\oplus} \rho^{-2}) = \xi_{1,r}^{\varepsilon,s} \right. \nonumber \\
& \left. -(3c_{1,1}^{\oplus} \rho^2 + c_{2,1}^{\oplus} + c_{3,1}^{\oplus} (\ln(\rho) + 1) - c_{4,1}^{\oplus} \rho^{-2}) = \xi_{1,\theta}^{\varepsilon,c}, \right. \tag{A.73}
\end{align}
$$

$$
\begin{equation}
\begin{aligned}
\star \quad k \ge 2: \quad & -k(\rho^{k+1} c_{1,k}^{\oplus} + \rho^{k-1} c_{2,k}^{\oplus}) \\
& \quad -k(\rho^{-k+1} c_{3,k}^{\oplus} + \rho^{-k-1} c_{4,k}^{\oplus}) = \xi_{k,r}^{\varepsilon,s}, \\
& \quad -((k+2)\rho^{k+1} c_{1,k}^{\oplus} + k\rho^{k-1} c_{2,k}^{\oplus}) \\
& \quad -((-k+2)\rho^{-k+1} c_{3,k}^{\oplus} - k\rho^{-k-1} c_{4,k}^{\oplus}) = \xi_{k,\theta}^{\varepsilon,c}.
\end{aligned}
\tag{A.74}
\end{equation}
$$

We conclude the derivation of the systems for the coefficients $c_{j,k}^{\ominus}$, $c_{j,k}^{\oplus}$ by considering the $k^{\text{th}}$ Fourier modes of $\mathbf{u}_{\varepsilon}^{(N)}$ on $\partial\omega_{\varepsilon}$. Namely, taking $r = \varepsilon$ both in (A.47)-(A.48) and (A.49)-(A.50) we get

$$
\begin{equation}
\star \quad k = 0 : - (2\varepsilon c_{1,0}^{\oplus} + \varepsilon^{-1} c_{2,0}^{\oplus} + c_{3,0}^{\oplus} \varepsilon (1 + 2 \ln(\varepsilon))) = -2\varepsilon c_{1,0}^{\ominus} - c_{3,0}^{\ominus} \varepsilon (1 + 2 \ln(\varepsilon)), \tag{A.75}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\star \quad k = 1: \quad & - (\varepsilon^2 c_{1,1}^{\oplus} + c_{2,1}^{\oplus} + \ln(\varepsilon) c_{3,1}^{\oplus} + \varepsilon^{-2} c_{4,1}^{\oplus}) = - (\varepsilon^2 c_{1,1}^{\ominus} + c_{2,1}^{\ominus}), \\
& - (3\varepsilon^2 c_{1,1}^{\oplus} + c_{2,1}^{\oplus} + (\ln(\varepsilon) + 1) c_{3,1}^{\oplus} - \varepsilon^{-2} c_{4,1}^{\oplus}) = - (3\varepsilon^2 c_{1,1}^{\ominus} + c_{2,1}^{\ominus}),
\end{aligned}
\tag{A.76}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
\star \quad k \ge 2: \quad & -k(\varepsilon^{k+1} c_{1,k}^{\oplus} + \varepsilon^{k-1} c_{2,k}^{\oplus}) - k(\varepsilon^{-k+1} c_{3,k}^{\oplus} + \varepsilon^{-k-1} c_{4,k}^{\oplus}) = \\
& \qquad -k(\varepsilon^{k+1} c_{1,k}^{\ominus} + \varepsilon^{k-1} c_{2,k}^{\ominus}), \\
& \quad -((k+2)\varepsilon^{k+1} c_{1,k}^{\oplus} + k\varepsilon^{k-1} c_{2,k}^{\oplus}) \\
& \qquad -((-k+2)\varepsilon^{-k} c_{3,k}^{\oplus} - k\varepsilon^{-k-1} c_{4,k}^{\oplus}) = \\
& \qquad -((k+2)\varepsilon^{k+1} c_{1,k}^{\ominus} + k\varepsilon^{k-1} c_{2,k}^{\ominus}).
\end{aligned}
\tag{A.77}
\end{equation}
$$

Thanks to the system of equations obtained above, we can express the coefficients $c_{k,j}^{\ominus}, c_{k,j}^{\oplus}$ in
terms of $\xi_{k,r}^{\varepsilon,j}, \xi_{k,\theta}^{\varepsilon,j}$ defined in (A.54). However, to derive (A.43) and the estimate (A.44), only the
expressions for $c_1^{\ominus}$ and $c_2^{\ominus}$ are needed. In details, we have

• For $k=0$ and $N=0$, gathering (A.68), (A.72) and (A.75),

$$
c_{1,0}^{\ominus} = - \frac{\xi_{0,\theta}^{\varepsilon,s}}{2\rho}. \qquad (A.78)
$$

• For $k = 0$ and $N \ge 1$, gathering (A.62), (A.63), (A.72) and (A.75),

$$
c_{1,0}^{\ominus} = - \frac{w_b}{2}.
$$
• For $k=1$ and $N \in \{0, 1\}$, gathering (A.61), (A.69), (A.73) and (A.76),

$$c_{1,1}^{\ominus} = \frac{(-1 + \bar{\epsilon}^2)(\xi_{1,r}^{\epsilon,s} + \xi_{1,\theta}^{\epsilon,c} + v_b \cdot e_y) + 2\log(\bar{\epsilon})(\xi_{1,r}^{\epsilon,s} - \xi_{1,\theta}^{\epsilon,c})}{4(\log(\bar{\epsilon}) + 1 - 2\bar{\epsilon}^2 + \bar{\epsilon}^4)}, \quad c_{2,1}^{\ominus} = -2\epsilon^2 c_{1,1}^{\ominus} - \frac{1}{2}v_b \cdot e_y. \qquad (A.79)$$

• For $k=1$ and $N \ge 2$, gathering (A.63), (A.73) and (A.76),

$$c_{1,1}^{\ominus} = 0, \quad c_{2,1}^{\ominus} = -\frac{1}{2}v_b \cdot e_y$$

• For $k \in \{2, \dots, N-1\}$ and $N \ge 2$, gathering (A.63), (A.74) and (A.77),

$$c_{1,k}^{\ominus} = c_{2,k}^{\ominus} = 0.$$

• For $k \ge 2$ and $k \in \{N, N+1\}$, gathering (A.63), (A.70), (A.74) and (A.77),

$$c_{1,k}^{\ominus} = \frac{(\xi_{k,r}^{\epsilon,c} - \xi_{k,\theta}^{\epsilon,s}) + \bar{\epsilon}^{2(k-1)}((k-2)\xi_{k,r}^{\epsilon,s} + k\xi_{k,\theta}^{\epsilon,s}) - (k-1)\bar{\epsilon}^{2k}(\xi_{k,r}^{\epsilon,c} + \xi_{k,\theta}^{\epsilon,s})}{2\rho^{k+1}(1 - k^2\bar{\epsilon}^2(2^{k-1}) + 2(k^2-1)\bar{\epsilon}^{2k} - (k^2-1)\bar{\epsilon}^{2(k+1)})}, \quad c_{2,k}^{\ominus} = -\bar{\epsilon}^2 \frac{k+1}{k} c_{1,k}^{\ominus}. \qquad (A.80)$$

• For $N=0$ and $k \ge 2$, or $k \ge N+2$, gathering (A.74), (A.77), and (A.71),

$$c_{1,k}^{\ominus} = \frac{\xi_{k,r}^{\epsilon,s} - \xi_{k,\theta}^{\epsilon,c}}{2\rho^{k+1}}, \quad c_{2,k}^{\ominus} = -\frac{(k+2)\xi_{k,r}^{\epsilon,s} - k\xi_{k,\theta}^{\epsilon,c}}{2k\rho^{k-1}}.$$

We deduce that $u_N^\epsilon - u_b$ satisfies (2.24) with

• $N=0$

$$\star k=0: \quad a_{\theta,0}^{\epsilon} = \xi_{0,\theta}^{\epsilon,c}$$

$$\star k=1: \quad b_{r,1}^{\epsilon} = \bar{\epsilon}\rho^2 c_{1,1}^{\ominus}, \quad a_{\theta,1}^{\epsilon} = -\bar{\epsilon}\rho^2 c_{1,1}^{\ominus},$$

$$\star k \ge 2: \quad \left\{ \begin{array}{l} b_{r,k}^{\epsilon} = -k \frac{\xi_{k,r}^{\epsilon,s} - \xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N} + \frac{(k+2)\xi_{k,r}^{\epsilon,s} - k\xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N-2}, \\ a_{\theta,k}^{\epsilon} = -(k+2) \frac{\xi_{k,r}^{\epsilon,s} - \xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N} + \frac{(k+2)\xi_{k,r}^{\epsilon,s} - k\xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N-2}. \end{array} \right.$$

• $N \ge 1$,

$$\star k \in \{N, N+1\}: \quad b_{r,k}^{\epsilon} = \bar{\epsilon}^{k-N} \rho^{k+1} c_{1,k}^{\ominus}, \quad a_{\theta,k}^{\epsilon} = -\bar{\epsilon}^{k-N} \rho^{k+1} c_{1,k}^{\ominus}$$

$$\star k \ge N+2: \quad \left\{ \begin{array}{l} b_{r,k}^{\epsilon} = -k \frac{\xi_{k,r}^{\epsilon,s} - \xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N} + \frac{(k+2)\xi_{k,r}^{\epsilon,s} - k\xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N-2}, \\ a_{\theta,k}^{\epsilon} = -(k+2) \frac{\xi_{k,r}^{\epsilon,s} - \xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N} + \frac{(k+2)\xi_{k,r}^{\epsilon,s} - k\xi_{k,\theta}^{\epsilon,c}}{2} \bar{\epsilon}^{k-N-2}. \end{array} \right.$$

To complete the proof of the lemma, it only remains to estimate the coefficients $c_{1,k}^{\oplus}$ for $k \in \{N, N+1\}$ with explicit bounds in terms of $k$ and $\bar{\epsilon}$.

• For $N=0$ and $k=0$, the coefficient $c_{1,0}^{\ominus}$ is given directly by (A.78), and therefore

$$|c_{1,0}^{\ominus}| = \frac{1}{2\rho} |\xi_{0,\theta}^{\epsilon,s}|.$$

For $N=0$ and $k=1$, the expression of $c_{1,1}^{\ominus}$ is given in (A.79). If $\xi_{1,r}^{\epsilon,s} - \xi_{1,\theta}^{\epsilon,c} \ne 0$, we have

$$\rho^2 c_{1,1}^{\oplus} \xrightarrow{\bar{\epsilon} \to 0} \frac{1}{2} (\xi_{1,r}^{\epsilon,s} - \xi_{1,\theta}^{\epsilon,c}).$$

Hence, there exists $\rho_1 > 0$ such that for $\epsilon \in (0, \rho_1)$,

$$\rho^2 |c_{1,1}^{\oplus}| \lesssim |\xi_{1,r}^{\epsilon,s}| + |\xi_{1,\theta}^{\epsilon,c}|.$$

Otherwise,

$$\rho^2 c_{1,1}^{\oplus} \xrightarrow{\bar{\epsilon} \to 0} 0,$$

and it follows that there exists $\rho_1 > 0$ such that for $\epsilon \in (0, \rho_1)$,

$$\rho^2 |c_{1,1}^{\oplus}| \lesssim 1.$$
• For $k \ge 2$ and $k \in \{N, N+1\}$, the coefficient $c_{1,k}^{\ominus}$ is given in (A.80). If $\xi_{k,r}^{\varepsilon,s} - \xi_{k,\theta}^{\varepsilon,c} \ne 0$, there exists $\bar{\rho} > 0$ such that

$$ \frac{2\rho^{k+1} c_{1,k}^{\ominus}}{\xi_{k,r}^{\varepsilon,s} - \xi_{k,\theta}^{\varepsilon,c}} \xrightarrow{k \to \infty} 1 \quad \text{uniformly on } (0, \bar{\rho}). $$

We deduce that there exists $k_0 \in \mathbb{N}$ such that for $k \ge k_0$ and any $\varepsilon \in (0, \bar{\rho})$,

$$ 2\rho^{k+1} c_{1,k}^{\ominus} \le c (|\xi_{k,r}^{\varepsilon,s}| + |\xi_{k,\theta}^{\varepsilon,c}|) $$

For the finite number of indices $k$ inside $\{2, \dots, k_0\}$, we have

$$ \rho^{k+1} c_{1,k}^{\ominus} \xrightarrow{\varepsilon \to 0} \frac{1}{2} (\xi_{k,r}^{\varepsilon,s} - \xi_{k,\theta}^{\varepsilon,c}) $$

so that there exists $\bar{\rho}_k > 0$ such that for $\varepsilon \in (0, \bar{\rho}_k)$,

$$ \rho^{k+1} c_{1,k}^{\ominus} \le c(k)(|\xi_{k,r}^{\varepsilon,s}| + |\xi_{k,\theta}^{\varepsilon,c}|). $$

Taking

$$ \rho_2 \stackrel{\text{def}}{=} \min_{k \in \{2, \dots, k_0\}} \{\bar{\rho}_M, \bar{\rho}_k\} $$

and $\varepsilon \in (0, \rho_2)$, we obtain that

$$ \rho^{k+1} c_{1,k}^{\ominus} \lesssim |\xi_{k,r}^{\varepsilon,s}| + |\xi_{k,\theta}^{\varepsilon,c}|, $$

with an implicit constant independent of $k$ and $\varepsilon$.

• For $k \ge N+2$, we directly have

$$ |b_{r,k}^{\varepsilon}| + |a_{\theta,k}^{\varepsilon}| \lesssim (1+k)\varepsilon^{N-k-2} $$

Finally, by taking $\tilde{\rho} = \min(\rho_1, \rho_2)$ and setting

$$ \xi_k = (|\xi_{k,r}^{\varepsilon,s}|^2 + |\xi_{k,\theta}^{\varepsilon,c}|^2 + |\xi_{k,r}^{\varepsilon,c}|^2)^{\frac{1}{2}} \quad (A.81) $$

we deduce that, for any $\varepsilon \in (0, \bar{\rho})$,

$$ \star \quad k \in \{N, N+1\} : |b_{r,k}^{\varepsilon}|^2 + |a_{\theta,k}|^2 \lesssim \left(\frac{\varepsilon}{\rho}\right)^{2(N-k)} \xi_k^2, \quad (A.82) $$

$$ \star \quad k \ge N+2 : |b_{r,k}^{\varepsilon}|^2 + |a_{\theta,k}^{\varepsilon}|^2 \lesssim (1+k)^2 \left(\frac{\varepsilon}{\rho}\right)^{2(N-k-2)} \xi_k^2. \quad (A.83) $$

We can proceed in a similar way to obtain a priori estimates for $b_{\theta,k}^{\varepsilon}$ and $a_{r,k}^{\varepsilon}$. We now decompose $\|\mathbf{u}_{\varepsilon} - \mathbf{u}_{\varepsilon}^{(N)}\|_{\frac{1}{2},\varepsilon}^{\mathcal{F}}$ as

$$ \|\mathbf{u}_{\varepsilon} - \mathbf{u}_{\varepsilon}^{(N)}\|_{\frac{1}{2},\varepsilon}^{\mathcal{F}} = \sum_{k=N}^{N+1} (1+k)(|a_{r,k}^{\varepsilon}|^2 + |a_{\theta,k}^{\varepsilon}|^2 + |b_{r,k}^{\varepsilon}|^2 + |b_{\theta,k}^{\varepsilon}|^2) \\ + \sum_{k=N+2}^{\infty} (1+k)(|a_{r,k}^{\varepsilon}|^2 + |a_{\theta,k}^{\varepsilon}|^2 + |b_{r,k}^{\varepsilon}|^2 + |b_{\theta,k}^{\varepsilon}|^2). \quad (A.84) $$

Focusing on the first term of (A.84), we deduce from (A.83) that, for any $\varepsilon \in (0, \bar{\rho})$,

$$ \sum_{k=N+2}^{\infty} (1+k)(|a_{r,k}^{\varepsilon}|^2 + |a_{\theta,k}^{\varepsilon}|^2 + |b_{r,k}^{\varepsilon}|^2 + |b_{\theta,k}^{\varepsilon}|^2) \\ \lesssim \sum_{k=N+2}^{\infty} \left(\frac{\bar{\rho}}{\rho}\right)^{2(N-k-2)} (1+k)^3 \xi_k^2 = \sum_{k=0}^{\infty} \left(\frac{\bar{\rho}}{\rho}\right)^{2k} (1+N+k+2)^3 \xi_{k+N+2}^2 $$

Using the inequality comes from $(1+k+N+2)^2 \lesssim (1+k)^2(1+N)^2$. and the boundedness of the sequence $(\frac{\bar{\rho}}{\rho})^{2k}$ $(1+k)^2$, we deduce that

$$ \sum_{k=N+2}^{\infty} (1+k)(|a_{r,k}^{\varepsilon}|^2 + |a_{\theta,k}^{\varepsilon}|^2 + |b_{r,k}^{\varepsilon}|^2 + |b_{\theta,k}^{\varepsilon}|^2) \\ \lesssim (1+N)^2 \sum_{k=0}^{\infty} (1+k+N+2) \xi_{k+N+2}^2 = (1+N)^2 \sum_{k=N+2}^{\infty} (1+k) \xi_k^2 $$