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

HAL
open science

# A new class of uniformly stable time-domain Foldy-Lax models for scattering by small particles. Acoustic sound-soft scattering by circles. Extended version

Maryna Kachanovska

► To cite this version:

Maryna Kachanovska. A new class of uniformly stable time-domain Foldy-Lax models for scattering by small particles. Acoustic sound-soft scattering by circles. Extended version. Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal, In press. hal-03664569v2

HAL Id: hal-03664569

https://hal.science/hal-03664569v2

Submitted on 8 Jun 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A NEW CLASS OF UNIFORMLY STABLE TIME-DOMAIN FOLDY-LAX MODELS FOR SCATTERING BY SMALL PARTICLES. ACOUSTIC SOUND-SOFT SCATTERING BY CIRCLES (EXTENDED VERSION)

MARYNA KACHANOVSKA*

**Abstract.** In this work we study time-domain sound-soft scattering by small circles. Our goal is to derive an asymptotic model for this problem valid when the size of the particles tends to zero. We present a systematic approach to constructing such models, based on a well-chosen Galerkin discretization of a boundary integral equation. The convergence of the method is achieved by decreasing the asymptotic parameter rather than increasing the number of basis functions. We prove the second-order convergence of the field error with respect to the particle size. Our findings are illustrated with numerical experiments.

**1. Introduction.** The problem of wave scattering by many particles has been a subject of active research since more than a century, cf. the monograph by Martin [32] for an overview of the various approaches to multiple scattering, as well as various multiple scattering problem settings, or the monograph by Mishchenko [36] for studies of electromagnetic multiple scattering. This type of wave propagation problems often appears in applications, for example in non-destructive testing (elastic or sound wave scattering by small defects) and atmospheric optics (scattering of light by atmospheric particles). Depending on the nature of the phenomena, such multiple wave scattering can be studied either from deterministic or probabilistic viewpoint, see [41]. E.g. the original seminal articles by Foldy [15] and Lax [25] were concerned with scattering by randomly positioned obstacles. In this work we will concentrate on the time-domain sound scattering by small particles in a deterministic regime.

In the asymptotic regime when the size of the particles tends to zero, and the distance between them is fixed (or decreases 'slowly'), it is possible to obtain simpler (and, in particular, easier for the computational treatment) models. To our knowledge, there exist two principal approaches to do so: either by using matching (near-field and far-field) asymptotic expansions, where one usually works with the original PDE, or by deriving asymptotics from the integral equation representations. Somewhat apart stands the original method of Foldy [15], developed further by Lax [25] (Foldy-Lax methods), cf. [32, Chapter 8.3] for its detailed description.

The literature on the frequency-domain asymptotic models for various types of wave propagation problems is quite rich; a non-exhaustive list of works exploiting either of the above approaches includes [6, 7, 24, 10, 38, 11, 9, 26, 27]. It seems that there exist fewer results in the time domain (see the recent monograph by Martin [33]). Asymptotic models for acoustic and electromagnetic wave propagation were obtained in [34, 4] (3D scattering by a single obstacle), [31], [19] and in [20] (transmission problems) by using matched expansion method. The integral equation approach was applied to the 3D wave scattering by multiple obstacles in the recent work [40]. However, the stability of the model of [40] was proven under some geometric conditions relating an asymptotic parameter, the number of particles and the minimal distance between them. In general, time-domain asymptotic models may exhibit instabilities, unlike the original model. We will show that this is the case for the time-domain counterpart of the Foldy-Lax model for the acoustic scattering by circles of [9].

The goal of this work is to derive an asymptotic Foldy-Lax model for the time domain scattering by small particles, which would be stable for arbitrary geometries (i.e. provided that particles do not touch each other). One of the novelties is the way such a model will be obtained: we start with a well-chosen boundary integral formulation, which possesses some coercivity properties, and semi-discretize it in space with the help of a very coarse Galerkin method. The resulting model is then stable due to the coercivity of the underlying boundary integral formulation. The convergence of the method does not rely on increasing the number of the basis functions, but on decreasing the asymptotic parameter. Let us remark that this idea is not entirely new: it

<footnote>*POEMS, CNRS, ENSTA Paris, INRIA, Institut Polytechnique de Paris, Palaiseau, France maryna.kachanovska@inria.fr</footnote>
appeared in the form of numerical experiments in the PhD thesis [23] for frequency-domain electromagnetic scattering by spheres.

We proceed as follows. Section 2 is an introduction: we present the problem, introduce the notation, and discuss whether the scattered field can be approximated by zero. We finish Section 2 by a motivation to the present work: we examine the time-domain counterpart of the Foldy-Lax model for the acoustic scattering by circles [9] and demonstrate that it is unstable for some geometric configurations. Section 3 is dedicated to the introduction of the new Galerkin Foldy-Lax model. In Section 4 we present convergence analysis of the new model. A short Section 5 extends the results of the article to the case when particles densify. In Section 6 we present numerical experiments, and conclude with a discussion of open questions in Section 7.

## 2. Problem setting and motivation.

**2.1. Problem setting.** Let $\Omega_j, j \in \mathcal{N} := \{1, \dots, N\}$, be a collection of $N > 0$ circular domains (particles), with the $j$th particle being centered at $c_j \in \mathbb{R}^2$ and having a radius $R_j > 0$, i.e. $\Omega_j = B(c_j, R_j)$. We assume that $\Omega_j \cap \Omega_i = \emptyset$ if $i \neq j$. For $0 < \varepsilon \le 1$, we introduce the rescaled domains $\Omega_j^\varepsilon = B(c_j, r_j^\varepsilon)$, where $r_j^\varepsilon = \varepsilon R_j$. Their interior, exterior and the boundary are denoted by

$$ \Omega^\varepsilon = \cup_j \Omega_j^\varepsilon, \quad \Omega^{\varepsilon,c} = \mathbb{R}^2 \setminus \bar{\Omega}^\varepsilon, \quad \Gamma_j^\varepsilon = \partial \Omega_j^\varepsilon, \quad \Gamma^\varepsilon = \cup_j \Gamma_j^\varepsilon = \partial \Omega^\varepsilon. $$

For the particular case $\varepsilon = 1$, we use a simplified notation $\Omega = \Omega^1$, $\Omega^c = \Omega^{1,c}$ and $\Gamma = \Gamma^1$. The trace operator on $\Gamma^\varepsilon$ is denoted by $\gamma_0$. Additionally, we define

$$ d_*^\varepsilon = \min_{i \neq j} \mathrm{dist}(\Omega_i^\varepsilon, \Omega_j^\varepsilon), \quad d_* = d_*^1, \quad R_* = \min_j R_j, \quad R^* = \max_j R_j. $$

Evidently, $d_* > d_*^\varepsilon$ for all $0 < \varepsilon < 1$.

We look for the solution of the sound-soft scattering problem. Provided sufficiently regular initial data $u_0, u_1: \Omega \to \mathbb{R}$, s.t. $\mathrm{supp} u_0 \cap \bar{\Omega} = \mathrm{supp} u_1 \cap \bar{\Omega} = \emptyset$, we are given the solution $u^{\mathrm{inc}}: \mathbb{R}_+ \times \mathbb{R}^2 \to \mathbb{R}$ of the free-space wave equation:

$$ (2.1) \qquad \begin{aligned} \partial_t^2 u^{\mathrm{inc}}(t, \boldsymbol{x}) - \Delta u^{\mathrm{inc}}(t, \boldsymbol{x}) &= 0, \quad (t, \boldsymbol{x}) \in \mathbb{R}_+^* \times \mathbb{R}^2, \\ u^{\mathrm{inc}}(0, \boldsymbol{x}) &= u_0(\boldsymbol{x}), \quad \partial_t u^{\mathrm{inc}}(0, \boldsymbol{x}) = u_1(\boldsymbol{x}). \end{aligned} $$

Presence of the obstacles $\Omega^\varepsilon$ alters the field; the new field $u_{\mathrm{tot}}^\varepsilon$ satisfies the following boundary-value problem (BVP):

$$ \begin{aligned} & \partial_t^2 u_{\mathrm{tot}}^\varepsilon - \Delta u_{\mathrm{tot}}^\varepsilon = 0 \quad \text{in} \quad \mathbb{R}_+^* \times \Omega^{\varepsilon,c}, \\ & \gamma_0 u_{\mathrm{tot}}^\varepsilon(t) = 0, \quad t \ge 0, \\ & u_{\mathrm{tot}}^\varepsilon(0) = u_0, \quad \partial_t u_{\mathrm{tot}}^\varepsilon(0) = u_1, \quad \text{in} \quad \Omega^{\varepsilon,c}. \end{aligned} $$

The auxiliary scattered field $u^\varepsilon := u_{\mathrm{tot}}^\varepsilon - u^{\mathrm{inc}}$ solves the inhomogeneous exterior BVP

$$ (2.2) \qquad \begin{aligned} & \partial_t^2 u^\varepsilon - \Delta u^\varepsilon = 0, \quad \text{in} \quad \mathbb{R}_+^* \times \Omega^{\varepsilon,c}, \\ & \gamma_0 u^\varepsilon(t) = g^\varepsilon(t), \quad t \ge 0, \quad \text{where } g^\varepsilon(t) := -\gamma_0 u^{\mathrm{inc}}(t), \\ & u^\varepsilon(0) = \partial_t u^\varepsilon(0) = 0 \quad \text{in} \quad \Omega^{\varepsilon,c}. \end{aligned} $$

**2.2. Goal of the present article.** We are interested in solving the above problem efficiently in the asymptotic regime, when the data $g^\varepsilon$ is approximately frequency-bandlimited to the frequency $\omega_{\mathrm{max}}$, and $\varepsilon$ is small compared to the smallest wavelength present in the system, i.e. $\varepsilon\omega_{\mathrm{max}} \ll 1$.
For each fixed $t > 0$, for a compact $K \subset \Omega^c$, $\lim_{\epsilon \to 0} \|u^\epsilon(t)\|_{L^2(K)} = 0$. In 2D, we can expect that (cf. Section 2.4.1) $\|u^\epsilon\|_{L^\infty(\mathbb{R}_+; L^2(K))} \ge c|\log^{-1}\epsilon|$, with some $c > 0$ independent of $\epsilon$. Our goal is to find an approximation $\tilde{u}^\epsilon$ of $u^\epsilon$ which would satisfy the following properties (provided that the initial data is sufficiently regular):

* (convergence) this approximation is more accurate than approximating $u^\epsilon$ by zero, i.e., for all $t > 0$, compact $K \subset \Omega^c$ and regular initial data $u_0, u_1$,

$$ (2.3) \qquad \|u^\epsilon(t) - \tilde{u}^\epsilon(t)\|_{L^2(K)} = o(\log^{-1} \epsilon). $$

* (uniform w.r.t. $\epsilon$ stability) the approximation is uniformly stable, i.e. there exist $C > 0$, $\ell_0, \ell_1 \in \mathbb{N}$, $p \in \mathbb{N}$, s.t., for all $0 < \epsilon \le 1$, and all $T > 0$,

$$ (2.4) \qquad \|\tilde{u}^\epsilon\|_{L^\infty(0,T;L^2(\Omega^{\epsilon,c}))} \le C(1+T)^p (\|u_0\|_{H^{\ell_0}(\mathbb{R}^2)} + \|u_1\|_{H^{\ell_1}(\mathbb{R}^2)}). $$

The uniform stability requirement may seem too strong: indeed, in general it is not necessary for the convergence. One could have authorized a super-polynomial growth (e.g. $e^{aT}$, with $a > 0$ potentially depending on $\epsilon$) of the solution. Such a model would likely be unsuitable for long time simulations. Alternatively, one could think of a CFL like condition on $\epsilon$ which would be dependent on the geometry and which would guarantee the bound (2.4) (for example, this is the case for the model [40]). However, in practical situations such a condition may be difficult to ensure. Thus we prefer to have models possessing a uniform stability property.

**2.3. Notation. Fourier-Laplace transform.** We will use the following convention of the Fourier-Laplace transform; for $v \in L^1(\mathbb{R}; X)$ (with $X$ being a Banach space), s.t. $v = 0$ on $\mathbb{R}_-$ (i.e. a causal function), we define

$$ \hat{v}(\omega) := \int_{0}^{+\infty} e^{i\omega t} v(t) dt, \quad \omega \in \mathbb{C}. $$

The above integral converges in particular for $\omega \in \mathbb{C}^+ := \{\omega \in \mathbb{C} : \operatorname{Im}\omega > 0\}$.

In the frequency domain (2.2) reads, for $\omega \in \mathbb{C}^+$,

$$ (2.5) \qquad \begin{aligned} & -\omega^2 \hat{u}^\epsilon - \Delta \hat{u}^\epsilon = 0 \quad \text{in } \Omega^{\epsilon,c}, \\ & \gamma_0 \hat{u}^\epsilon = \hat{g}^\epsilon, \quad \text{where } \hat{g}^\epsilon = -\gamma_0 \hat{u}^{\text{inc}}. \end{aligned} $$

Let us introduce the fundamental solution for $-\Delta - \omega^2$:

$$ G_\omega(r) = \frac{i}{4} H_0^{(1)}(\omega r), \quad r > 0, $$

where $H_0^{(1)}$ is a Hankel function, cf. [13, Chapter 10.2]. It is a Fourier-Laplace transform of the Green function for the 2D wave equation

$$ \mathcal{G}(t, r) = \frac{1_{\{t>r\}}}{2\pi\sqrt{t^2 - r^2}}. $$

*Sobolev spaces on the boundary of a circle.* Let $r > 0$, and let $C_r := \partial B(0, r)$. Given a function $u \in L^2(C_r)$, we define its Fourier decomposition (with $\hat{s}_\theta = (\cos \theta, \sin \theta)$):

$$ u(\hat{s}_\theta) = \sum_{m \in \mathbb{Z}} u_m e^{im\theta}, \quad u_m = \frac{1}{2\pi} \int_0^{2\pi} e^{-im\theta} u(\hat{s}_\theta) d\theta. $$
The fractional Sobolev spaces on the boundary of the circle are defined as follows: $u \in H^s(C_r)$, $s \ge 0$, iff the following norm is finite:

$$ (2.6) \qquad \|u\|_{H^s(C_r)}^2 = 2\pi r \sum_{m \in \mathbb{Z}} \left(1 + \frac{m^2}{r^2}\right)^s |u_m|^2. $$

For $s < 0$, we define $H^s(C_r)$ as a completion of $L^2(C_r)$ for the norm (2.6).

We will further make use of the following observation: by [18, proof of Lemma 4.2.5], the Sobolev-Slobodeckij seminorm

$$ (2.7) \qquad |v|_{H^{1/2}(C_r)}^2 = \iint_{C_r \times C_r} \frac{|v(x) - v(y)|^2}{|x - y|^2} dx dy $$

is equivalent to $\sum_{m \neq 0} m |v_m|^2$, with the equivalence constants independent of $r$. Sobolev spaces on $\Gamma^\epsilon$. We define

$$ H^s(\Gamma^\epsilon) := \prod_{k \in \mathcal{N}} H^s(\Gamma_k^\epsilon), \quad s \in \mathbb{R}. $$

For a function $v \in H^s(\Gamma^\epsilon)$, let us set $v_m := v|_{\Gamma_m^\epsilon}$. We equip $H^s(\Gamma^\epsilon)$ with the norm

$$ \|v\|_{H^s(\Gamma^\epsilon)}^2 := \sum_{k \in \mathcal{N}} \|v_k\|_{H^s(\Gamma_k^\epsilon)}^2. $$

Sometimes we will use a shortened notation $\|v\|_s$ for $\|v\|_{H^s(\Gamma^\epsilon)}$.

By $(\rho_j, \theta_j)$ we will denote the polar coordinates centered in $c_j$. The coefficients of the Fourier expansion of $v_j$ on $\Gamma_j^\epsilon$ will be denoted by $v_{j,m}$, $m \in \mathbb{Z}$.

The (real) duality bracket in $H^{-1/2}(\Gamma^\epsilon)$, $H^{1/2}(\Gamma^\epsilon)$ will be denoted by $\langle \cdot, \cdot \rangle$. More precisely, for $u, v \in L^2(\Gamma^\epsilon)$, $\langle u, v \rangle = \int_{\Gamma^\epsilon} u v d\Gamma$.

An energy norm. Given $a > 0$ and a domain $\mathcal{O}$, we denote

$$ \|v\|_{a, \mathcal{O}}^2 = a^2 \|v\|_{L^2(\mathcal{O})}^2 + \|\nabla v\|_{L^2(\mathcal{O})}^2. $$

A notation for a complement domain. Given a domain $\mathcal{O}$, we denote by $\mathcal{O}^c = \mathbb{R}^3 \setminus \bar{\mathcal{O}}$.

Trace operators. Later in the article it will be necessary to distinguish between the exterior and the interior traces.

We define a normal $n$ to $\Omega^\epsilon$ so that it points into the exterior of $\Omega^\epsilon$. Let us consider a function $\phi \in H_\Delta^1(\mathbb{R}^2 \setminus \Gamma^\epsilon)$, where

$$ H_\Delta^1(\mathbb{R}^2 \setminus \Gamma^\epsilon) = \{v \in L^2(\mathbb{R}^2) : v|_{\Omega^\epsilon} \in H_\Delta^1(\Omega^\epsilon), v|_{\Omega^{\epsilon,c}} \in H_\Delta^1(\Omega^{\epsilon,c})\}, \quad \text{with} \\ H_\Delta^1(\mathcal{O}) = \{v \in H^1(\mathcal{O}) : \Delta v \in L^2(\mathcal{O})\}. $$

We then denote by $\gamma_0^+ \phi$ and $\gamma_1^+ \phi$ its exterior trace and its exterior normal trace, and by $\gamma_0^- \phi$ and $\gamma_1^- \phi$ its interior traces. More precisely, if $\phi$ is sufficiently regular, we have that $\gamma_1^- \phi = (\gamma_0^- \nabla \phi) \cdot n$. By $[\gamma_0 \phi]$, $[\gamma_1 \phi]$ we denote the respective jumps:

$$ [\gamma_0 \phi] = \gamma_0^- \phi - \gamma_0^+ \phi, \quad [\gamma_1 \phi] = \gamma_1^- \phi - \gamma_1^+ \phi. $$

Inequality notation. We will use $a \lesssim b$ (resp., $a \gtrsim b$) to indicate that $a \le Cb$ (resp., $a \ge Cb$), for a generic constant $C$ independent of $\epsilon$, $N$, geometry, final time, data (i.e. any parameters of the problem).

Notation for $\min(1, a)$. Where convenient, we will use the notation $a = \min(1, a)$.
**2.4. Motivation for designing a new asymptotic model.** In this section we argue that constructing a non-trivial asymptotic model for (2.2) is necessary, since the exact solution $u^\epsilon$ decays to zero as $\epsilon \to 0$ at a fairly slow rate $O(\log^{-1} \epsilon)$, cf. Section 2.4.1. On the other hand, simple rewriting of an available frequency-domain Foldy-Lax model in the time domain by applying the inverse Fourier-Laplace transform may yield a potentially unstable in the time domain model, see Section 2.4.2.

<figure><img src="image_3.png" /><figcaption>FIG. 2.1. Illustration to Section 2.4.1. Left: dependence of $\|u^\epsilon(\cdot, \mathbf{x}_0)\|_{L^\infty(0,T)}$ at $\mathbf{x}_0 = (0.2, 0.2)$, $T = 8$, on $\epsilon$. Right: dependence of $u^\epsilon(t, \mathbf{x}_0)$ on time for different values of $\epsilon$.</figcaption></figure>

**2.4.1. Slow convergence of the solution $u^\epsilon$ to zero.** First of all, let us consider scattering by a single particle in the frequency domain. In the result below we provide an asymptotic expansion of the solution with respect to $\epsilon > 0$.

**PROPOSITION 2.1.** Let $\omega > 0$, $N = 1$, $c_1 = 0$, $R_1 = 1$ and $\rho > 0$ be fixed. Additionally, let $\hat{u}^{\text{inc}} \in C^1(\mathbb{R}^2)$. Then, as $\epsilon \to 0$, the following holds true uniformly in $\theta \in [0, 2\pi)$:

$$ \hat{u}^{\epsilon}(\omega, \rho\hat{s}_{\theta}) = \frac{i\pi}{2} H_0^{(1)}(\omega\rho) \hat{u}^{\text{inc}}(\omega, 0) \log^{-1} \epsilon + o(\log^{-1} \epsilon). $$

*Proof.* See Appendix A. □

This result can be translated into the time domain. Since this is technical, we rather illustrate the time-domain counterpart with numerical experiments. Let us finally remark that a related (but more general) question was studied in [4].

To obtain a numerical time-domain counterpart of Proposition 2.1, we compute the field scattered by a single obstacle ($N = 1$, $R = 1$, $c_1 = 0$) on the interval $(0, 8)$ for different $\epsilon > 0$. The incident field is given by $u^{\text{inc}}(t, \mathbf{x}) = -e^{-20(t-d\cdot\mathbf{x}-1)^2}$, with $d = (0, 1)$. The solution is computed using the time-domain BEM, semi-discretized in time using the trapezoid rule Convolution Quadrature (CQ) method [30, 28, 29], and in space using the Galerkin method with the basis functions $\{e^{im\theta}\}_{m=-N_s}^{N_s}$. The results are given in Figure 2.1. We observe a good agreement of the results with the lower bound $\|u^\epsilon(\cdot, \mathbf{x}_0)\|_{L^\infty(0,T)} \ge c_T \log^{-1} \epsilon$.

**2.4.2. Potential instability of frequency-domain models recast into the time domain.** Let us consider one of frequency-domain models, namely the Foldy-Lax model analyzed in [9], which is $O(\epsilon/\log\epsilon)$-accurate in the frequency domain. We will show that there exist geometric configurations for which this model is unstable.
Foldy-Lax model of [9]. In [9] it was suggested to approximate the field $\hat{u}^{\varepsilon}(\mathbf{x})$ solving (2.5) with $\omega > 0$ by the following linear combination:

$$ (2.8) \qquad \hat{u}^{\varepsilon}(\mathbf{x}) \approx \hat{u}_{FL}^{\varepsilon}(\mathbf{x}) = \sum_{k \in \mathcal{N}} \frac{G_{\omega}(\|\mathbf{x} - \mathbf{c}_k\|)}{G_{\omega}(r_k^{\varepsilon})} \hat{\lambda}_{FL,k}^{\varepsilon}, \quad \text{where } \hat{\lambda}_{FL}^{\varepsilon} \in \mathbb{C}^N \text{ solves} $$

$$ (2.9) \qquad \hat{\lambda}_{FL,n}^{\varepsilon} + \sum_{k \in \mathcal{N} \setminus \{n\}} \frac{G_{\omega}(\|\mathbf{c}_n - \mathbf{c}_k\|)}{G_{\omega}(r_k^{\varepsilon})} \hat{\lambda}_{FL,k}^{\varepsilon} = -\hat{u}^{\text{inc}}(\mathbf{c}_n), \quad n = 1, \dots, N. $$

To rewrite the above in the time domain, first of all we remark that $\hat{\lambda}_{FL}^{\varepsilon}$ is frequency-dependent, and can be shown to be a Fourier-Laplace transform of a causal time-dependent function. With $\hat{\mu}_{FL,n}^{\varepsilon} = (G_{\omega}(r_n^{\varepsilon}))^{-1}\hat{\lambda}_{FL,n}^{\varepsilon}$, the time-domain Foldy-Lax approximation of the scattered field reads:

$$ (\text{FL1}) \qquad \begin{aligned} u^{\varepsilon}(t, \mathbf{x}) &\approx u_{FL}^{\varepsilon}(t, \mathbf{x}) = \sum_{k \in \mathcal{N}} \mathcal{G}(t, \|\mathbf{x} - \mathbf{c}_k\|) *_t \mu_{FL,k}^{\varepsilon} \\ &= \frac{1}{2\pi} \sum_{k \in \mathcal{N}} \int_0^t \frac{\mathbb{1}_{|t-\tau| > \|\mathbf{x}-\mathbf{c}_k\|} \|\mathbf{x}-\mathbf{c}_k\|}{\sqrt{(t-\tau)^2 - \|\mathbf{x}-\mathbf{c}_k\|^2}} \mu_{FL,k}^{\varepsilon}(\tau) d\tau, \end{aligned} $$

where the functions $\mu_{FL}^{\varepsilon}: \mathbb{R}_{+} \to \mathbb{R}^N$ satisfy, for all $n = 1, \dots, N$,

$$ (\text{FL2}) \qquad \mathcal{G}(t, r_n^{\varepsilon}) *_t \mu_{FL,n}^{\varepsilon} + \sum_{k \in \mathcal{N} \setminus \{n\}} \mathcal{G}(t, \|\mathbf{c}_n - \mathbf{c}_k\|) *_t \mu_{FL,k}^{\varepsilon} = -u^{\text{inc}}(t, \mathbf{c}_n). $$

It is possible to show that, for $N \le 2$, the model (FL1-FL2) is uniformly stable, i.e. (2.4) holds for $\tilde{u}^{\varepsilon} = u_{FL}^{\varepsilon}$. However, when $N > 2$, this is no longer the case; e.g. (FL1-FL2) can exhibit instabilities when the particles are very close to each other.

**Instability of the model (FL1-FL2).** It appears that for a particular geometry consisting of three particles that are located in the vertices of an equilateral triangle, s.t. the distance between them is smaller than their diameter, the Foldy-Lax model (FL1-FL2) is unstable. This is formalized in the statement below.

**PROPOSITION 2.2.** Let $N = 3$ and $R_i = r$ for all $i$. Let $\|\mathbf{c}_i - \mathbf{c}_j\| = c > 0$ for all $i \ne j$ (in other words, the centers of the particles are located in the vertices of an equilateral triangle of side length $c$). Assume that $c/r < 4$.

Then there exist $u_0, u_1 \in C_0^{\infty}(\Omega^c)$, a compact $K \subset \Omega^c$ and $A, C > 0$ s.t.

$$ \limsup_{t \to +\infty} (e^{-At} \|u_{FL}^1(t, \cdot)\|_{L^2(K)}) \ge C. $$

The proof of this proposition can be found in Appendix B. Let us explain its main idea. Examining (FL1), (FL2) shows that the temporal behaviour of $u_{FL}^1$ is closely related to that one of $\lambda_{FL}^1$. In order to show that $t \mapsto \|\lambda_{FL}^1(t)\|$ admits exponential growth, we prove that its Fourier-Laplace transform has poles in $\mathbb{C}^+$. Its (infinitely many) poles $\{\omega_n\}_{n \in \mathbb{Z}}$ satisfy (using the notation of Proposition 2.2)

$$ (2.10) \qquad \lim_{n \to \pm\infty} \operatorname{Im} \omega_n = \frac{1}{r(\eta-1)} \log \frac{2}{\sqrt{\eta}} > 0, \quad \text{where } \eta = cr^{-1} < 4. $$

One can then expect the rate $A$ in Proposition 2.2 to be close to $\sup_n \operatorname{Im} \omega_n$.
A numerical illustration to Proposition 2.2. Let us now illustrate the statement of Proposition 2.2. We computed the Foldy-Lax solution $u_{FL}^{\epsilon}$ for the configuration of Proposition 2.2, with $r = 0.1$, $c_i$ located in the vertices of the equilateral triangle with $c_1 = (-0.105, 0)$, $c_2 = (0.105, 0)$ and $c_3 \approx (0, 0.1819)$ (so that $\eta = 2.1$). Moreover, the incident field and the final simulation time are given by

$$u^{\text{inc}}(t, \mathbf{x}) = e^{-20(t-x_2-2)^2} \sin(30(t-x_2-2)), \quad T = 4.4.$$

An illustration to this experiment is given in Figure 2.2. As discussed after the statement of Proposition 2.2,

<figure><img src="image_4.png" /><figcaption>FIG. 2.2. The dependence of the absolute value of the solution $u_{FL}^{\epsilon}(t, \mathbf{x}_0)$, $\mathbf{x}_0 = (0.2, 0.2)$ on time, computed for the experiment of Section 2.4.2.</figcaption></figure>

we expect the rate of the exponential growth of the solution to be close to $\sup_n \text{Im}\,\omega_n$. The asymptotic rate given by (2.10) is approximately 2.93, which, as we see in Figure 2.2, predicts fairly well the behaviour of the solution.

REMARK 1. Performing a direct numerical simulation of the problem in Proposition 2.2 is delicate, due to the time semi-discretization. The CQ method, cf. the end of Section 2.4.1, can be shown to converge for this problem, see the analysis in [30], but the issue lies in its implementation, namely computation of the convolution weights. This is usually done in the Laplace domain using contour integration techniques [29], with the contour that depends on the final simulation time and the time step $\Delta t$. This procedure may fail because of the presence of the resonances in $\mathbb{C}^+$. Hence it is difficult to ensure convergence of the algorithm in this case, though the exponential blow-up is of course seen numerically. Because of such convergence problems, we computed the Foldy-Lax solution only on a fairly short time interval, where we could ensure the validity of the method used for approximating the convolution weights.

Discussion. Proposition 2.2 states that for some geometric configurations, the Foldy-Lax model (FL1-FL2) exhibits instabilities. However, as $\epsilon \to 0$, one can prove that the three-particle configuration of Proposition 2.2 with $r_j^{\epsilon} = \epsilon r$ becomes again stable. For the moment, we do not know whether for any system of particles there exists $\epsilon_0 > 0$ depending on $N$, $d_*$, s.t. (FL1-FL2) is stable for all $\epsilon < \epsilon_0$ (we conjecture that this is the case).

In this paper we develop an alternative Foldy-Lax model, which is uniformly stable and $O(\epsilon^2)$-convergent. This model it is stable for any geometric configuration where the particles do not touch each other. It can be viewed as an $O(\epsilon^2)$-perturbation of the Foldy-Lax model (FL1-FL2), see Section 3.4.3. More experiments comparing the two models can be found in Section 6.2.

### 3. Galerkin Foldy-Lax model.
**3.1. Derivation.** The solution to (2.2) can be represented as the single-layer potential of an unknown density $\boldsymbol{\mu}^{\varepsilon}$:

$$ (3.1) \qquad u^{\varepsilon}(t, \boldsymbol{x}) = S^{\varepsilon}(\partial_t)\phi(t, \boldsymbol{x}) = \int_0^t \int_{\Gamma^{\varepsilon}} \mathcal{G}(t - \tau, \boldsymbol{x} - \boldsymbol{y})\boldsymbol{\mu}^{\varepsilon}(\tau, \boldsymbol{y})d\Gamma_{\boldsymbol{y}}d\tau, \quad \boldsymbol{x} \in \Omega^{\varepsilon,c}. $$

The trace of the single-layer potential is continuous across $\Gamma^{\varepsilon}$, and thus the density $\boldsymbol{\mu}^{\varepsilon}$ can be found from the following boundary integral equation:

$$ (3.2) \qquad \begin{aligned} g^{\varepsilon}(t, \boldsymbol{x}) &= \gamma_0 \boldsymbol{u}^{\varepsilon}(t, \boldsymbol{x}) = S^{\varepsilon}(\partial_t)\boldsymbol{\mu}^{\varepsilon}(t, \boldsymbol{x}), \\ S^{\varepsilon}(\partial_t)\boldsymbol{\mu}^{\varepsilon}(t, \boldsymbol{x}) &= \gamma_0 S^{\varepsilon}(\partial_t)\boldsymbol{\mu}^{\varepsilon}(t, \boldsymbol{x}) = \int_0^t \int_{\Gamma^{\varepsilon}} \mathcal{G}(t - \tau, \boldsymbol{x} - \boldsymbol{y})\boldsymbol{\mu}^{\varepsilon}(\tau, \boldsymbol{y})d\Gamma_{\boldsymbol{y}}d\tau, \quad \boldsymbol{x} \in \Gamma^{\varepsilon}. \end{aligned} $$

We will not make precise the spaces in which the single-layer operator $S^{\varepsilon}$ is invertible; those can be found in [2]. Nonetheless, since we are going to construct the asymptotic model as its Galerkin semi-discretization in space, we remark that the operator $S^{\varepsilon}$ is an inverse Fourier-Laplace transform of the Helmholtz single-layer boundary integral operator $\hat{S}^{\varepsilon}$. The latter is defined for a sufficiently regular $\phi$ as follows:

$$ (3.3) \qquad \hat{S}^{\varepsilon} \phi := \int_{\Gamma^{\varepsilon}} G_{\omega}(\|\boldsymbol{x} - \boldsymbol{y}\|) \phi(\boldsymbol{y}) d\Gamma_{\boldsymbol{y}}, \quad \boldsymbol{x} \in \Gamma^{\varepsilon}, $$

and is continuous in appropriate spaces: $\hat{S}^{\varepsilon} \in \mathcal{L}(H^{-1/2}(\Gamma^{\varepsilon}), H^{1/2}(\Gamma^{\varepsilon}))$. Therefore, a spatial variational formulation for (3.2) can be written as follows:

$$ (3.4) \qquad \begin{aligned} \text{find } \boldsymbol{\mu}^{\varepsilon} : \mathbb{R}_{+} &\rightarrow H^{-1/2}(\Gamma^{\varepsilon}), \text{ s.t.} \\ \langle \boldsymbol{v}, \boldsymbol{g}^{\varepsilon}(t, \cdot) \rangle &= \langle \boldsymbol{v}, S^{\varepsilon}(\partial_t)\boldsymbol{\mu}^{\varepsilon}(t, \cdot) \rangle, \quad \text{for all } \boldsymbol{v} \in H^{-1/2}(\Gamma^{\varepsilon}). \end{aligned} $$

Let us now introduce the following coarse Galerkin space:

$$ (3.5) \qquad \mathcal{V}_0^{\varepsilon} \subset H^{-1/2}(\Gamma^{\varepsilon}), \quad \mathcal{V}_0^{\varepsilon} = \prod_{n=1}^{N} \mathcal{V}_0(\Gamma_n^{\varepsilon}), \quad \mathcal{V}_0(\Gamma_n^{\varepsilon}) = \operatorname{span}\{1, \boldsymbol{x} \in \Gamma_n^{\varepsilon}\}. $$

Evidently, $\mathcal{V}_0^{\varepsilon}$ is an $N$-dimensional space, with the basis $\{\boldsymbol{e}^{\varepsilon,n}(\boldsymbol{x})\}_{n=1}^N$ defined as

$$ \boldsymbol{e}^{\varepsilon,n}(\boldsymbol{x}) = \begin{cases} 1, & \boldsymbol{x} \in \Gamma_n^{\varepsilon}, \\ 0, & \text{otherwise.} \end{cases} $$

The main idea is to discretize (3.4) by using the Galerkin method with the trial and test space $\mathcal{V}_0^{\varepsilon}$. Such a discretization rewrites:

$$ (3.6) \qquad \begin{aligned} \text{find } \boldsymbol{\mu}_G^{\varepsilon} : \mathbb{R}_{+} &\rightarrow \mathcal{V}_0^{\varepsilon}, \text{ s.t.} \\ \langle \boldsymbol{e}^{\varepsilon,n}, \boldsymbol{g}^{\varepsilon}(t, \cdot) \rangle &= \langle \boldsymbol{e}^{\varepsilon,n}, S^{\varepsilon}(\partial_t) \boldsymbol{\mu}_G^{\varepsilon}(t, \cdot) \rangle, \quad \text{for all } n \in \mathcal{N}. \end{aligned} $$

Let us introduce

$$ (3.7) \qquad \mathcal{G}_{nk}^{\varepsilon}(t) = \iint_{\Gamma_n^{\varepsilon} \times \Gamma_k^{\varepsilon}} \mathcal{G}(t, \|\boldsymbol{x} - \boldsymbol{y}\|) d\boldsymbol{x} d\boldsymbol{y}, $$
and rewrite $\mu_G^\epsilon(t, \mathbf{x}) = \sum_{n \in \mathcal{N}} \mu_{G,n}^\epsilon(t) \mathbf{e}^{\epsilon,n}(\mathbf{x})$. Then (3.6) can be rewritten as a convolutional system of equations:

$$ (GFL1) \qquad \int_{\Gamma_n^\epsilon} g_n^\epsilon(t, \mathbf{x}) d\Gamma_n = \int_0^t \mathcal{G}_{nn}^\epsilon(t-\tau) \mu_{G,n}^\epsilon(\tau) d\tau \\ \qquad \qquad + \sum_{k \in \mathcal{N} \setminus \{n\}} \int_0^t \mathcal{G}_{nk}^\epsilon(t-\tau) \mu_{G,k}^\epsilon(\tau) d\tau, \quad n \in \mathcal{N}. $$

Knowing the approximate density $\mu_G^\epsilon$ allows to compute the approximate solution according to (3.1):

$$ (GFL2) \qquad u_G^\epsilon(t, \mathbf{x}) = \sum_{k \in \mathcal{N}} \int_0^t \mathcal{G}_k^\epsilon(t-\tau, \mathbf{x}) \mu_{G,k}^\epsilon(\tau) d\tau, \quad \mathcal{G}_k^\epsilon(t, \mathbf{x}) = \int_{\Gamma_k^\epsilon} \mathcal{G}(t, \mathbf{x}-\mathbf{y}) d\Gamma_y. $$

The problem (GFL1) together with (GFL2) is the desired asymptotic model.

**REMARK 2.** The above model is a priori stable for each $0 < \epsilon \le 1$. This is formalized in Theorem 4.1 and in Proposition 4.3. The reason for this is that the system of equations in (GFL1) stems from the Galerkin semidiscretization of the single-layer boundary integral operator $\mathcal{S}^\epsilon(\partial_t)$. By applying the Fourier-Laplace transform to (3.6) we obtain the following identity for each $\omega \in \mathbb{C}^+$:

$$ (3.8) \qquad \langle \mathbf{e}^{\epsilon,n}, \hat{\mathbf{g}}^{\epsilon} \rangle = \langle \mathbf{e}^{\epsilon,n}, \hat{\mathcal{S}}^{\epsilon} \hat{\mu}_{G}^{\epsilon} \rangle, \quad \text{for all } n \in \mathcal{N}. $$

The above problem is well-posed because the frequency-domain single layer boundary integral operator is coercive for all $\omega \in \mathbb{C}^+$: $\text{Im}(-\langle v, \omega \hat{\mathcal{S}}^{\epsilon} v \rangle) \gtrsim c_1(\epsilon, \omega) \|v\|_{H^{-1/2}(\Gamma^{\epsilon})}^2$. Combined with some analyticity arguments, this argument allows to show the well-posedness of (3.6). Moreover, the frequency-domain bound $\|\hat{\mu}_G^\epsilon\|_{-1/2} \lesssim c_2(\epsilon, \omega) \|\hat{\mathbf{g}}^\epsilon\|_{1/2}$ can be translated into a time-domain stability bound. See the proof of Proposition 4.3 for a more detailed description of this approach.

**3.2. Numerical illustration of the stability of the Galerkin Foldy-Lax model for the configuration of Proposition 2.2.** Let us illustrate the stability of the Galerkin Foldy-Lax model for the configuration described in Proposition 2.2. For this we use the same data as in Section 2.4.2. The field $u_G^\epsilon(t, \mathbf{x}_0)$ is shown in Figure 3.1. It remains bounded, unlike $u_{\tilde{F}_L}^\epsilon$, cf. Figure 2.2, although we see that the solution obtained with the help of the asymptotic model is quite far from the reference solution, computed with the help of high-order BEM.

**3.3. Galerkin Foldy-Lax model in the frequency domain.** For the analysis purposes, we will need to work with the Fourier-Laplace transformed Galerkin Foldy-Lax model. We fix $\omega \in \mathbb{C}^+$. Then the exact solution to (2.5) can be written as the single layer potential of the density $\hat{\mu}^\epsilon(\mathbf{x})$:

$$ (3.9) \qquad \hat{u}^\epsilon(\mathbf{x}) = \hat{\mathcal{S}}^\epsilon \hat{\mu}^\epsilon(\mathbf{x}) = \int_{\Gamma^\epsilon} G_\omega(\|\mathbf{x} - \mathbf{y}\|) \hat{\mu}^\epsilon(\mathbf{y}) d\Gamma_y, \quad \mathbf{x} \in \Omega^{\epsilon,c}. $$

The above operator satisfies: $\hat{\mathcal{S}}^\epsilon \in \mathcal{L}(H^{-1/2}(\Gamma^\epsilon), H^1(\mathbb{R}^d \setminus \Gamma^\epsilon))$. The density can be found from the frequency-domain counterpart of (3.2),

$$ (3.10) \qquad \hat{\mathbf{g}}^{\epsilon}(\mathbf{x}) = \hat{\mathcal{S}}^{\epsilon} \hat{\mu}^{\epsilon} = \int_{\Gamma^{\epsilon}} G_{\omega}(\|\mathbf{x} - \mathbf{y}\|) \hat{\mu}^{\epsilon}(\mathbf{y}) d\Gamma_{\mathbf{y}}, \quad \mathbf{x} \in \Gamma^{\epsilon}, $$
<figure><img src="image_5.png" /><figcaption>FIG. 3.1. An illustration to the numerical experiment of Section 3.2. The dependence of the approximate solution $u_G^\epsilon(t, \mathbf{x}_0)$ and the reference solution $u^\epsilon(t, \mathbf{x}_0)$ at $\mathbf{x}_0 = (0.2, 0.2)$ on time.</figcaption></figure>

see also (3.3). Introducing

$$ \hat{S}_{nk}^{\epsilon} : H^{-1/2}(\Gamma_k^{\epsilon}) \rightarrow H^{1/2}(\Gamma_n^{\epsilon}), \quad n, k \in \mathcal{N}, $$

defined for a sufficiently regular $\phi$ by

$$ (3.11) \qquad (\hat{S}_{nk}^{\epsilon}\phi)(\mathbf{x}) = \int_{\Gamma_k^{\epsilon}} G_{\omega}(\|\mathbf{x} - \mathbf{y}\|)\phi(\mathbf{y})d\Gamma_{\mathbf{y}}, \quad \mathbf{x} \in \Gamma_n^{\epsilon}, $$

we can further rewrite the equation (3.10) as follows:

$$ (3.12) \qquad \hat{g}_n^{\epsilon}(\mathbf{x}) = \hat{S}_{nn}^{\epsilon} \hat{\mu}_n^{\epsilon} + \sum_{k \in \mathcal{N} \setminus \{n\}} \hat{S}_{nk}^{\epsilon} \hat{\mu}_k^{\epsilon}, \quad n \in \mathcal{N}. $$

The frequency-domain counterpart of the Foldy-Lax model (GFL1) then reads:

$$ (3.13) \qquad \int_{\Gamma_n^{\epsilon}} \hat{g}_n^{\epsilon}(\mathbf{x}) d\Gamma_{\mathbf{x}} = \hat{\mu}_{G,n}^{\epsilon} \iint_{\Gamma_n^{\epsilon} \times \Gamma_n^{\epsilon}} G_{\omega}(\|\mathbf{x} - \mathbf{y}\|) d\Gamma_{\mathbf{y}} d\Gamma_{\mathbf{x}} \\ \qquad + \sum_{k \in \mathcal{N} \setminus \{n\}} \hat{\mu}_{G,k}^{\epsilon} \iint_{\Gamma_n^{\epsilon} \times \Gamma_k^{\epsilon}} G_{\omega}(\|\mathbf{x} - \mathbf{y}\|) d\Gamma_{\mathbf{x}} d\Gamma_{\mathbf{y}}, \quad n \in \mathcal{N}. $$

An approximation $\hat{u}_G^\epsilon$ to the exact solution $\hat{u}^\epsilon$ can be computed according to (3.9):

$$ (3.14) \qquad \hat{u}_G^{\epsilon}(\mathbf{x}) = \sum_{k=1}^{N} \hat{\mu}_{G,k}^{\epsilon} \int_{\Gamma_k^{\epsilon}} G_{\omega}(\|\mathbf{x} - \mathbf{y}\|) d\Gamma_{\mathbf{y}}, \quad \mathbf{x} \in \Omega^{\epsilon,c}. $$

### 3.4. Explicit Galerkin Foldy-Lax model for circles.
Practical implementation of the Galerkin Foldy-Lax model relies on the evaluation of the quantities (3.7) (computing the respective double integrals). In the case of the circular scatterers, this can be done analytically, exploiting the Graf addition theorem, on one hand, and the fact that constant functions are eigenfunctions of the single-layer boundary integral operator on a circle, on the other hand.
<figure><img src="image_6.png" /><figcaption>FIG. 3.2. An illustration to the two applications of the Graf addition theorem.</figcaption></figure>

### 3.4.1. Galerkin Foldy-Lax model in the frequency domain.

Derivation of the system of equations for unknown coefficients. We start with the variational formulation (3.8), which we rewrite in accordance with (3.12):

$$ (3.15) \qquad \langle \mathbf{e}^{\varepsilon,n}, \hat{\mathbf{g}}^{\varepsilon} \rangle = \langle \mathbf{e}^{\varepsilon,n}, \hat{S}_{nn}^{\varepsilon} \mathbf{e}^{\varepsilon,n} \rangle \hat{\mu}_{G,n}^{\varepsilon} + \sum_{k \in \mathcal{N} \setminus \{n\}} \langle \mathbf{e}^{\varepsilon,n}, \hat{S}_{nk}^{\varepsilon} \mathbf{e}^{\varepsilon,k} \rangle \hat{\mu}_{G,k}^{\varepsilon}. $$

Let us rewrite the right-hand side of the above identity. The operator $\hat{S}_{nn}^{\varepsilon}$ is normal, as proven in [8]; moreover, its (non-normalized) eigenfunctions are $\{e^{ik\theta_n(x)}\}_{k \in \mathbb{Z}}$, and eigenvalues are given by [22, 21],

$$ (3.16) \qquad \hat{S}_{nn}^{\varepsilon} e^{ik\theta_n} = \frac{i\pi r_n^{\varepsilon}}{2} H_{|k|}^{(1)}(\omega r_n^{\varepsilon}) J_{|k|}(\omega r_n^{\varepsilon}) e^{ik\theta_n}, \quad \text{therefore} \\ \langle \mathbf{e}^{\varepsilon,n}, \hat{S}_{nn}^{\varepsilon} \mathbf{e}^{\varepsilon,n} \rangle = i\pi^2 (r_n^{\varepsilon})^2 H_0^{(1)}(\omega r_n^{\varepsilon}) J_0(\omega r_n^{\varepsilon}). $$

The operator $\hat{S}_{nk}^{\varepsilon}$ can also be rewritten in a more convenient form, by using the Graf addition theorem [1, 9.1.79]. According to it, for $z, p \in \mathbb{C}$, s.t. $p = |p|e^{i\theta^p}$, $z = |z|e^{i\theta^z}$, with $|p| > |z|$, and $w = p - z = |w|e^{i\theta^w}$, $\chi = \theta^w - \theta^p$, it holds that

$$ (3.17) \qquad H_m^{(1)}(\omega|p|e^{i\theta^p} - |z|e^{i\theta^z}) = e^{-im\chi} \sum_{k \in \mathbb{Z}} H_{m+k}^{(1)}(\omega|p|) J_k(\omega|z|) e^{ik(\theta^p - \theta^z)}. $$

Given $\mathbf{x} \in \Gamma_k^{\varepsilon}$, $\mathbf{y} \in \Gamma_n^{\varepsilon}$, first we apply the Graf theorem to $H_0^{(1)}(\omega\|\mathbf{x} - \mathbf{y}\|)$, with $p = \|\mathbf{x} - \mathbf{c}_n\|e^{i\theta_n(\mathbf{x})}$ and $z = \|\mathbf{y} - \mathbf{c}_n\|e^{i\theta_n(\mathbf{y})}$. We obtain its expansion into the series of $H_m^{(1)}(\omega\|\mathbf{x} - \mathbf{c}_m\|)$, and next we re-apply the Graf theorem to each of the terms of the series, with $p = \|\mathbf{c}_n - \mathbf{c}_k\|e^{i\theta_k(\mathbf{c}_n)}$ and $z = \|\mathbf{x} - \mathbf{c}_k\|e^{i\theta_k(\mathbf{x})}$ (in this latter case $w = p - z = \|\mathbf{x} - \mathbf{c}_n\|e^{i\theta_n(\mathbf{x})}$, $\theta_n(\mathbf{c}_k) = \pi + \theta_k(\mathbf{c}_n)$ and $\chi = \theta_n(\mathbf{x}) - \theta_n(\mathbf{c}_k)$). These two applications are illustrated in Figure 3.2. We finally obtain the following expansion:

$$ H_0^{(1)}(\omega\|\mathbf{x} - \mathbf{y}\|) = \sum_{m=-\infty}^{\infty} J_m(\omega r_n^{\varepsilon}) e^{-im\theta_n(\mathbf{y})} \\ \times \sum_{j=-\infty}^{\infty} H_{m+j}^{(1)}(\omega\|\mathbf{c}_n - \mathbf{c}_k\|) e^{i(m\theta_n(\mathbf{c}_k) + j\theta_k(\mathbf{c}_n))} J_j(\omega r_k^{\varepsilon}) e^{ij\theta_k(\mathbf{x})}. $$

Replacing $H_0^{(1)}$ in the definition of $\hat{S}_{nk}^{\varepsilon}$ by the above expression yields

$$ (3.18) \qquad \langle \mathbf{e}^{n,\varepsilon}, \hat{S}_{nk}^{\varepsilon} \mathbf{e}^{k,\varepsilon} \rangle = i\pi^2 r_k^{\varepsilon} r_n^{\varepsilon} J_0(\omega r_k^{\varepsilon}) J_0(\omega r_n^{\varepsilon}) H_0^{(1)}(\omega\|\mathbf{c}_k - \mathbf{c}_n\|). $$
Therefore, the Foldy-Lax model (3.15) then can be rewritten in the following form:

$$ (3.19) \qquad \int_{\Gamma_n^\varepsilon} \hat{g}^\varepsilon(\mathbf{x}) d\Gamma_\mathbf{x} = i\pi^2 (r_n^\varepsilon)^2 \mathbf{H}_0^{(1)}(\omega r_n^\varepsilon) \mathbf{J}_0(\omega r_n^\varepsilon) \hat{\mu}_{G,n}^\varepsilon \\ \qquad \qquad + i\pi^2 \sum_{k \in \mathcal{N} \setminus \{n\}} r_k^\varepsilon r_n^\varepsilon \mathbf{J}_0(\omega r_k^\varepsilon) \mathbf{J}_0(\omega r_n^\varepsilon) \mathbf{H}_0^{(1)}(\omega \|\mathbf{c}_k - \mathbf{c}_n\|) \hat{\mu}_{G,k}^\varepsilon. $$

*Approximation of the field.* Once the coefficients $\hat{\mu}_G^\varepsilon$ are found by solving the system of equations (3.19), it remains to approximate the field $\hat{u}^\varepsilon$. We start with (3.14), which we further rewrite by using the Graf addition theorem (3.17):

$$ (3.20) \qquad \begin{aligned} \hat{u}_G^\varepsilon &= \frac{i}{4} \sum_{k \in \mathcal{N}} \hat{\mu}_{G,k}^\varepsilon \int_{\Gamma_k^\varepsilon} \mathbf{H}_0^{(1)}(\omega \|\mathbf{x} - \mathbf{y}\|) d\Gamma_\mathbf{y} \\ &= \frac{i}{4} \sum_{k \in \mathcal{N}} \hat{\mu}_{G,k}^\varepsilon \mathbf{H}_0^{(1)}(\omega \|\mathbf{x} - \mathbf{c}_k\|) \int_{\Gamma_k^\varepsilon} \mathbf{J}_0(\omega \|\mathbf{y} - \mathbf{c}_k\|) d\Gamma_\mathbf{y} \\ &= \frac{i\pi}{2} \sum_{k \in \mathcal{N}} \hat{\mu}_{G,k}^\varepsilon r_k^\varepsilon \mathbf{H}_0^{(1)}(\omega \|\mathbf{x} - \mathbf{c}_k\|) \mathbf{J}_0(\omega r_k^\varepsilon). \end{aligned} $$

**3.4.2. Galerkin Foldy-Lax model in the time domain.** To rewrite the Galerkin Foldy-Lax model for circular obstacles in the time domain, let us introduce the following causal distributions (with 'sc' standing for 'scaled', compared to the analogous definitions in Section 3.1):

$$ \begin{aligned} \mathcal{G}_{nn}^{sc,\varepsilon} &= \mathcal{F}^{-1} \frac{i\pi r_n^\varepsilon}{2} \mathbf{J}_0(\omega r_n^\varepsilon) \mathbf{H}_0^{(1)}(\omega r_n^\varepsilon), \\ \mathcal{G}_{nk}^{sc,\varepsilon} &= \mathcal{F}^{-1} \frac{i\pi r_k^\varepsilon}{2} \mathbf{J}_0(\omega r_k^\varepsilon) \mathbf{J}_0(\omega r_n^\varepsilon) \mathbf{H}_0^{(1)}(\omega \|\mathbf{c}_k - \mathbf{c}_n\|), \quad k \neq n, \\ \mathcal{G}_k^{sc,\varepsilon}(\mathbf{x}) &= \mathcal{F}^{-1} \frac{i\pi r_k^\varepsilon}{2} \mathbf{H}_0^{(1)}(\omega \|\mathbf{x} - \mathbf{c}_k\|) \mathbf{J}_0(\omega r_k^\varepsilon). \end{aligned} $$

Then the Galerkin Foldy-Lax model (3.19), (3.20) in the time domain rewrites as

$$ (3.21) \qquad \begin{aligned} & -\frac{1}{2\pi} \int_0^{2\pi} u^{\text{inc}}(t, \mathbf{c}_n + r_n^\varepsilon \hat{s}_\theta) d\theta = \mathcal{G}_{nn}^{sc,\varepsilon} *_t \mu_{G,n}^\varepsilon + \sum_{k \in \mathcal{N} \setminus \{n\}} \mathcal{G}_{nk}^{sc,\varepsilon} *_t \mu_{G,k}^\varepsilon, \quad n \in \mathcal{N}, \\ & \qquad u_G^\varepsilon(t, \mathbf{x}) = \sum_{k \in \mathcal{N}} \mathcal{G}_k^{sc,\varepsilon}(\mathbf{x}) *_t \mu_{G,k}^\varepsilon. \end{aligned} $$

**3.4.3. Connection to the Foldy-Lax model (2.9).** The goal of this section is to show the relation between the original Foldy-Lax model and the new Galerkin one; we will show it in the frequency domain, i.e. for (2.9) and (3.19). We start by rewriting (3.19) in a form similar to (2.9). Let us consider the left-hand side of (2.9). First, $\hat{u}^{\text{inc}}(\mathbf{x})$ satisfies the homogeneous Helmholtz equation inside $B(\mathbf{c}_n, r_n^\varepsilon)$, therefore

$$ \hat{u}^{\text{inc}}(\mathbf{c}_n + \rho \hat{s}_\theta) = - \sum_{m \in \mathbb{Z}} \hat{g}_{n,m}^\varepsilon \frac{J_m(\omega \rho)}{J_m(\omega r_n^\varepsilon)} e^{-im\theta}, \quad \rho < r_n^\varepsilon. $$

Because $J_m(0) = \delta_{0m}$,

$$ (3.22) \qquad -\mathbf{J}_0(\omega r_n^\varepsilon) \hat{u}^{\text{inc}}(\mathbf{c}_n) = \hat{g}_{n,0}^\varepsilon = (2\pi r_n^\varepsilon)^{-1} \int_{\Gamma_n^\varepsilon} \hat{g}^\varepsilon(\mathbf{x}) d\Gamma_\mathbf{x}. $$
Next, let us rewrite the right-hand side of (2.9). Let

$$ \hat{\lambda}_{G,n}^{\varepsilon} = \frac{i\pi}{2} r_n^{\varepsilon} H_0^{(1)}(\omega r_n^{\varepsilon}) \hat{\mu}_{G,n}^{\varepsilon}. $$

Replacing the left-hand side of (3.19) by (3.22) and $\hat{\mu}_G^\varepsilon$ by its expression of $\hat{\lambda}_G^\varepsilon$ yields

$$ (3.23) \qquad -\hat{u}^{\text{inc}}(\mathbf{c}_n) = \hat{\lambda}_{G,n}^{\varepsilon} + \sum_{k \in \mathcal{N} \setminus \{n\}} J_0(\omega r_k^{\varepsilon}) \frac{H_0^{(1)}(\omega \|\mathbf{c}_k - \mathbf{c}_n\|)}{H_0^{(1)}(\omega r_k^{\varepsilon})} \hat{\lambda}_{G,k}^{\varepsilon}, \quad n \in \mathcal{N}. $$

And the field approximation then reads

$$ (3.24) \qquad -\hat{u}_G^{\varepsilon}(\mathbf{x}) = \sum_{k \in \mathcal{N}} J_0(\omega r_k^{\varepsilon}) \frac{H_0^{(1)}(\omega \|\mathbf{x} - \mathbf{c}_k\|)}{H_0^{(1)}(\omega r_k^{\varepsilon})} \hat{\lambda}_{G,k}^{\varepsilon}. $$

We recognize in (3.23) the system of equations (2.9), and in (3.24) the approximation (2.8), modulo the term $J_0(\omega r_k^\varepsilon)$. For a fixed frequency $\omega$, as $r_n^\varepsilon \to 0$, $J_0(\omega r_n^\varepsilon) = 1 + O(\varepsilon^2)$, and thus approximating $J_0(\omega r_n^\varepsilon)$ by 1 does not seem to affect the convergence rates of the model, at least in the frequency domain [9]. However, in the time domain, the absence of this term is accountable for a potential lack of stability.

**4. Convergence Analysis.** We aim at obtaining an estimate of the type

$$ \|u^{\varepsilon}(t) - u_{G}^{\varepsilon}(t)\|_{L^2(\mathcal{K})} < C \times \text{error}(\varepsilon), $$

where we will specify the type of dependence (e.g. polynomial) of the constant $C$ on the number of particles $N$, the minimal distance between particles $d_*$, simulation time, the largest and smallest radii $R^*$ and $R_*$, and make explicit the dependence of the error on a certain norm of the data $u^{\text{inc}}$. In the analysis, we never use an assumption about $N$ being a fixed number with respect to $\varepsilon > 0$, or $d_*^\varepsilon > d_* = \text{const}$. This enables us to extend our results to the case when $N$ and $d_*^\varepsilon$ depend on the asymptotic parameter $\varepsilon$ (Section 5).

**4.1. Main results.** Let us first introduce a technical assumption.

**ASSUMPTION 1.** We assume that the incident field satisfies

$$ u^{\text{inc}} \in H^8(0, T; W^{1,\infty}(\mathbb{R}^2)). $$

Such a regularity assumption holds true notably when $u_0$, $u_1 \in C_0^\infty(\mathbb{R}^2)$.

We will additionally need the weighted norm

$$ \|v\|_{H_{\varepsilon}^m(0,T;X)} := \sum_{k=0}^{m} \varepsilon^{2k} \|\partial_t^k v\|_{L^2(0,T;X)}^2. $$

**REMARK 3.** We would like to use such an $\varepsilon$-dependent norm in the convergence estimates, where suitable, since it takes into account the interplay between the wavelength and the asymptotic parameter $\varepsilon$: for $v(t) = e^{i\omega t}$, with $w > 0$, one has that $\|v\|_{H_{\varepsilon}^m(0,T)} \lesssim T(\sum_{k=0}^m (\varepsilon w)^{2k})^{1/2}$.

The principal result of this section is summarized below.
**THEOREM 4.1** (Convergence of $u_G^\epsilon$). Let Assumption 1 hold true. Let $\delta = \delta_\epsilon > 0$. Then, for all $0 < \epsilon \le 1$, $0 < t < T$, for all $\mathbf{x} \in \Omega^{\epsilon+\delta_\epsilon, c}$, the field $u_G^\epsilon$ satisfies the following error bound, with a constant $C_G$ depending polynomially on $R^*, R_*^{-1}$, $N, d_*^{-1}$ and $T$:

$$e^\epsilon(t, \mathbf{x}) = |u_G^\epsilon(t, \mathbf{x}) - u^\epsilon(t, \mathbf{x})| \le \epsilon^2 \times \max(1, \delta_\epsilon^{-1}) C_G \|\partial_t^5 u^{\text{inc}}\|_{H_\epsilon^2(0,T; W^{1,\infty}(\mathbb{R}^2))}.$$

In particular,

* when $\delta_\epsilon \ge \delta = \text{const} > 0$ ('far-field'), then, for all $\mathbf{x} \in \Omega^{\delta, c}$, as $\epsilon \to 0$,

$$e^\epsilon(t, \mathbf{x}) = O(\epsilon^2).$$

* when $\delta_\epsilon = \epsilon^\alpha$, with $\alpha \in (0, 2)$ ('near-field'), then, for $\mathbf{x}_\epsilon$ s.t. $\text{dist}(\mathbf{x}_\epsilon, \Omega^\epsilon) = \delta_\epsilon$, it holds that

$$e^\epsilon(t, \mathbf{x}_\epsilon) = O(\epsilon^{2-\alpha}), \text{ as } \epsilon \to 0.$$

**REMARK 4.** The error bound in the statement of the theorem can be refined to:

$$e^\epsilon(t, \mathbf{x}) \le \epsilon^2 \times \max(1, \delta_\epsilon^{-1}) C_G (\|\partial_t^5 u^{\text{inc}}\|_{H_\epsilon^3(0,T; L^\infty(\mathbb{R}^2))} + \|\partial_t^4 u^{\text{inc}}\|_{H_\epsilon^2(0,T; W^{1,\infty}(\mathbb{R}^2))}).$$

The proof of this theorem can be found in Section 4.6. The above result shows that the Galerkin Foldy-Lax model yields a second-order accurate approximation of the far field and a first-order approximation at the distance $\epsilon$ to the particle, provided that $\epsilon$ is sufficiently small with respect to the wavelength.

**REMARK 5.** Let us quantify the term 'sufficiently small'. Let the incident field be given by $u^{\text{inc}}(t, \mathbf{x}) = \Phi(t - \mathbf{x} \cdot \mathbf{d})$, with $\|\mathbf{d}\| = 1$ and $\Phi: \mathbb{R} \to \mathbb{R}$ being approximately bandlimited to the frequency $\omega_{\text{max}} > 1$ (a typical example being a modulated Gaussian $\Phi(p) = \sin(\omega_0 p)e^{-\sigma p^2}\chi(p)$, with $\chi(p)$ being a smooth truncation function that vanishes for sufficiently large $|p|$). In this case, cf. Remark 4, informally, we have

$$\|\partial_t^5 u^{\text{inc}}\|_{H_\epsilon^3(0,T; L^\infty(\mathbb{R}^2))} \lesssim \omega_{\text{max}}^5 \sum_{j=0}^3 (\omega_{\text{max}}\epsilon)^j, \\ \|\partial_t^4 u^{\text{inc}}\|_{H_\epsilon^2(0,T; W^{1,\infty}(\mathbb{R}^2))} \lesssim \omega_{\text{max}}^5 \sum_{j=0}^2 (\omega_{\text{max}}\epsilon)^j.$$

This shows that the Galerkin Foldy-Lax model starts converging in the far-field as soon as $\epsilon^2\omega_{\text{max}}^5 \lesssim 1$. We believe this condition to be non-optimal (probably $\epsilon\omega_{\text{max}} \lesssim 1$ is sufficient). It is due to the technique chosen to analyze the problem, see Section 4.7.

With Theorem 4.1, we can show that the Galerkin Foldy-Lax model has the desired properties listed in Section 2.2.

**COROLLARY 4.2.** Under Assumption 1, the Galerkin Foldy-Lax model is convergent and uniformly stable, as defined in Section 2.2.

*Proof.* See Appendix C.

Theorem 4.1 relies on the following result about convergence of the density $\mu^\epsilon$.

**PROPOSITION 4.3** (Convergence of the density). Let Assumption 1 hold true. Let $T > 0$. For all $0 < \epsilon \le 1$, the density $\mu_G^\epsilon \in L^\infty(0, T; V_0^\epsilon)$ satisfies the following bound, with the constant $C_{\mu 0}$ depending polynomially on $R^*, R_*^{-1}$, $N, d_*^{-1}$ and $T$:

$$\|\mu_G^\epsilon\|_{L^\infty(0,T; H^{-1/2}(\Gamma^\epsilon))} \le \epsilon^{-1/2} \times C_{\mu 0} \|\partial_t^2 u^{\text{inc}}\|_{H_\epsilon^1(0,T; L^\infty(\mathbb{R}^2))}.$$
The error $e^\epsilon = \mu_G^\epsilon - \mu^\epsilon$ satisfies, with the constant $C_{\mu,G}$ depending polynomially on $R^*$, $R_*^{-1}$, $N$, $d_*^{-1}$ and $T$:

$$ \|e^\epsilon\|_{L^\infty(0,T;H^{-1/2}(\Gamma^\epsilon))} \le \epsilon \times C_{\mu,G} \|\partial_t^5 u^{\text{inc}}\|_{H^3(0,T;W^{1,\infty}(\mathbb{R}^2))}. $$

The proof of this proposition is presented in Section 4.5. It states two results:

* the $H^{-1/2}(\Gamma^\epsilon)$-norm of the Galerkin density $\mu_G^\epsilon(t)$ grows at most as $O(\epsilon^{-1/2})$; this implies in particular that $\|\mu_G^\epsilon(t)\|_{L^\infty(\Gamma^\epsilon)} \lesssim \epsilon^{-1}$.

* the Galerkin Foldy-Lax method approximates $\mu_G^\epsilon$ with the error $O(\epsilon)$, i.e. we can expect the relative error to be of order $O(\epsilon^{3/2})$.

The growth of the density $\mu_G^\epsilon$ w.r.t. $\epsilon$ is due to the choice of the Foldy-Lax unknown $\mu^\epsilon$ as a solution to the single-layer integral equation. It is possible to rescale the formulation in order to avoid this problem, by taking as the new unknown

$$ \lambda_k^\epsilon = \mathbb{S}_{kk}^\epsilon(\partial_t)\mu_k^\epsilon = \int_0^t \int_{\Gamma_k^\epsilon} \mathcal{G}(t - \tau, \cdot, -\mathbf{y})\mu_k^\epsilon(\tau, \mathbf{y})d\tau d\Gamma_{\mathbf{y}}, \quad k = 1, \dots, N. $$

It can be shown that the Galerkin method applied to the integral formulation with this new unknown will yield a problem equivalent to (3.6), because constant functions are eigenfunctions of $\hat{\mathbb{S}}_{kk}^\epsilon$. The analysis of the respective model is out of scope of the present article, since, conceptually, it will not differ much from the analysis presented here, however, requires some extra technical results.

**REMARK 6.** Comparing Proposition 4.3 and Theorem 4.1, one may find it surprising that the error in the density of $O(\epsilon)$ can produce an error of $O(\epsilon^2)$ in the field. See Section 4.6 for a detailed explanation to such a 'super-convergence' phenomenon.

The sections that follow are dedicated to the proof of Theorem 4.1 and Proposition 4.3. We will perform the analysis in the Fourier-Laplace domain, see the seminal work by Bamberger and Ha-Duong [2] and the recent monograph by Sayas [39].

For this we first study the convergence of the asymptotic method, reformulated in the frequency domain, for the frequencies with absorption ($\omega \in \mathbb{C}^+$). The main goal is then to obtain convergence estimates that would be explicit in the asymptotic parameter $\epsilon$ and the frequency $\omega$. Next they need to be translated into the time domain based on the Plancherel theorem.

The rest of this section is organized as follows. In Section 4.2 we introduce the basis of our analysis. We derive an expression for the Galerkin error of the density $\hat{\mu}^\epsilon$, which would involve in particular the restrictions of the operator $\hat{\mathbb{S}}^\epsilon$ to different subspaces of $H^{-1/2}(\Gamma^\epsilon)$, as well as $\epsilon$-dependent $H^{1/2}(\Gamma^\epsilon)$-norms of the data. Then the bound on the error can be derived by obtaining explicit in $\epsilon$ estimates on the operator norms and on the data.

Section 4.3 is dedicated to obtaining the estimates on the operators in the frequency domain. For our purposes, it will be necessary to obtain the estimates on the data directly in the time domain, which is done in Section 4.4. In Section 4.5, we prove Proposition 4.3, by combining the results of the previous section. In Section 4.6 we prove Theorem 4.1. Finally, Section 4.7 is dedicated to the discussion of the optimality of the results, see also Remark 5.

**4.2. Error expression in the frequency domain.** To obtain an expression for the error that will be convenient to analyze, we will work with the Foldy-Lax model (FL1), (FL2) in the frequency domain, namely with (3.13) and (3.14). In this section we will obtain an expression for the error $\|\hat{\mu}^\epsilon - \hat{\mu}_G^\epsilon\|_{-1/2}$ between the densities (cf. Proposition 4.3). It is organized as follows. In Section 4.2.1, we introduce a decomposition of the spaces that will allow to consider separately the constant component of the error (i.e. the one that belongs to the Galerkin space) and the remainder; these errors will have different scalings with respect to the asymptotic
parameter $\epsilon$. Section 4.2.2 is dedicated to deriving an expression for the error that will serve as a basis for our analysis.

**4.2.1. Decomposition of spaces.** Let us decompose the space $H^{1/2}(\Gamma^\epsilon)$ into a direct sum of the spaces (orthogonal with respect to $L^2(\Gamma^\epsilon)$-scalar product)

$$H^{1/2}(\Gamma^\epsilon) = \mathcal{V}_0^\epsilon + H_*^{1/2}(\Gamma^\epsilon),$$

where $\mathcal{V}_0^\epsilon$ was defined in (3.5) as follows:

$$\mathcal{V}_0^\epsilon = \prod_{n \in \mathcal{N}} \mathcal{V}_0(\Gamma_n^\epsilon), \quad \mathcal{V}_0(\Gamma_n^\epsilon) = \operatorname{span}\{1, \boldsymbol{x} \in \Gamma_n^\epsilon\}, \quad \text{and} \\ H_*^{1/2}(\Gamma^\epsilon) = \{\boldsymbol{v} \in H^{1/2}(\Gamma^\epsilon) : (\boldsymbol{e}^{\epsilon,n}, \boldsymbol{v})_{L^2(\Gamma^\epsilon)} = 0, \text{ for all } n \in \mathcal{N}\}.$$

The respective $L^2(\Gamma^\epsilon)$-orthogonal projector on the space $\mathcal{V}_0^\epsilon$ is denoted by $\mathbb{P}_0^\epsilon$, and on the space $H_*^{1/2}(\Gamma^\epsilon)$ by $\mathbb{P}_\perp^\epsilon$. For any $\boldsymbol{v} \in H^{1/2}(\Gamma^\epsilon)$, we introduce the decomposition $\boldsymbol{v} = \boldsymbol{v}_0 + \boldsymbol{v}_\perp$, with $\boldsymbol{v}_0 \in \mathcal{V}_0^\epsilon$, $\boldsymbol{v}_\perp \in H_*^{1/2}(\Gamma^\epsilon)$. In a similar manner we decompose

$$H^{-1/2}(\Gamma^\epsilon) = \mathcal{V}_0^\epsilon + H_*^{-1/2}(\Gamma^\epsilon), \quad \text{where} \\ H_*^{-1/2}(\Gamma^\epsilon) = \{\boldsymbol{v} \in H^{-1/2}(\Gamma^\epsilon) : \langle \boldsymbol{v}, \boldsymbol{e}^{\epsilon,n} \rangle = 0, \text{ for all } n \in \mathcal{N}\}.$$

The respective adjoint projectors will be denoted by

$$\mathbb{P}_0^{\epsilon,*}: H^{-1/2}(\Gamma^\epsilon) \mapsto \mathcal{V}_0^\epsilon, \quad \mathbb{P}_\perp^{\epsilon,*}: H^{-1/2}(\Gamma^\epsilon) \mapsto H_*^{-1/2}(\Gamma^\epsilon).$$

They are defined as follows. For $\boldsymbol{v} \in H^{-1/2}(\Gamma^\epsilon)$, $\boldsymbol{\lambda} \in H^{1/2}(\Gamma^\epsilon)$:

$$\langle \mathbb{P}_0^{\epsilon,*} \boldsymbol{v}, \bar{\boldsymbol{\lambda}} \rangle = \langle \boldsymbol{v}, \mathbb{P}_0^{\epsilon} \bar{\boldsymbol{\lambda}} \rangle, \quad \langle \mathbb{P}_\perp^{\epsilon,*} \boldsymbol{v}, \bar{\boldsymbol{\lambda}} \rangle = \langle \boldsymbol{v}, \mathbb{P}_\perp^{\epsilon} \bar{\boldsymbol{\lambda}} \rangle.$$

Any $\boldsymbol{v} \in H^{-1/2}(\Gamma^\epsilon)$ can be written as $\boldsymbol{v} = \boldsymbol{v}_0 + \boldsymbol{v}_\perp$, with $\boldsymbol{v}_0 \in \mathcal{V}_0^\epsilon$, $\boldsymbol{v}_\perp \in H_*^{-1/2}(\Gamma^\epsilon)$.

It remains to introduce the following operators:

$$(4.1) \qquad \hat{S}_{00}^{\epsilon} := \mathbb{P}_{0}^{\epsilon} \hat{S}^{\epsilon} \mathbb{P}_{0}^{\epsilon,*}, \quad \hat{S}_{0\perp}^{\epsilon} = \mathbb{P}_{0}^{\epsilon} \hat{S}^{\epsilon} \mathbb{P}_{\perp}^{\epsilon,*}, \quad \hat{S}_{\perp 0}^{\epsilon} = \mathbb{P}_{\perp}^{\epsilon} \hat{S}^{\epsilon} \mathbb{P}_{0}^{\epsilon,*}, \quad \hat{S}_{\perp\perp}^{\epsilon} = \mathbb{P}_{\perp}^{\epsilon} \hat{S}^{\epsilon} \mathbb{P}_{\perp}^{\epsilon,*}.$$

The main idea of the error analysis is to exploit different scaling of the norms of these operators with respect to $\epsilon \to 0$.

**REMARK 7.** *Remark that there should be no confusion between $\hat{S}_{00}^{\epsilon}$ defined above and the operators $\hat{S}_{kk}^{\epsilon}$, $k = 1, \dots, N$, corresponding to the single layer boundary integral operators for a single obstacle, cf. (3.11).*

**REMARK 8.** *Where convenient, we will view the operators (4.1) as operators acting on the subspaces $\mathcal{V}_0^\epsilon$, $H_*^{\pm 1/2}(\Gamma^\epsilon)$, without specifying this explicitly. E.g.*

$$\hat{S}_{00}^{\epsilon} : (\mathcal{V}_0^{\epsilon}, \|\cdot\|_{-1/2}) \to (\mathcal{V}_0^{\epsilon}, \|\cdot\|_{1/2}), \\ \hat{S}_{0\perp}^{\epsilon} : H_*^{-1/2}(\Gamma^{\epsilon}) \to (\mathcal{V}_0^{\epsilon}, \|\cdot\|_{1/2}), \text{ and so on.}$$

In what follows, we will also omit the indices in the operator norms $\|\cdot\|$ when working with the above operators.
### 4.2.2. Error expression in the frequency domain.
With the decomposition of the spaces introduced in the previous section, we rewrite the Galerkin Foldy-Lax model as: find $\hat{\mu}_G^\epsilon \in V_0^\epsilon$, s.t.

$$ (4.2) \qquad \hat{g}_0^\epsilon = \hat{S}_{00}^\epsilon \hat{\mu}_G^\epsilon. $$

The exact solution $\hat{\mu}^\epsilon$ satisfies the system

$$ \begin{pmatrix} \hat{S}_{00}^\epsilon & \hat{S}_{0\perp}^\epsilon \\ \hat{S}_{\perp 0}^\epsilon & \hat{S}_{\perp \perp}^\epsilon \end{pmatrix} \begin{pmatrix} \hat{\mu}_0^\epsilon \\ \hat{\mu}_\perp^\epsilon \end{pmatrix} = \begin{pmatrix} \hat{g}_0^\epsilon \\ \hat{g}_\perp^\epsilon \end{pmatrix}. $$

From the above it is straightforward to obtain an equation for the error $\hat{e}^\epsilon = \hat{\mu}^\epsilon - \hat{\mu}_G^\epsilon$:

$$ \begin{pmatrix} \hat{S}_{00}^\epsilon & \hat{S}_{0\perp}^\epsilon \\ \hat{S}_{\perp 0}^\epsilon & \hat{S}_{\perp \perp}^\epsilon \end{pmatrix} \begin{pmatrix} \hat{e}_0^\epsilon \\ \hat{e}_\perp^\epsilon \end{pmatrix} = \begin{pmatrix} 0 \\ \hat{g}_\perp^\epsilon - \hat{S}_{\perp 0}^\epsilon \hat{\mu}_G^\epsilon \end{pmatrix}. $$

The constant component of the error can be obtained by the usual Schur complement:

$$ (4.3) \qquad \hat{e}_0^\epsilon = - (\hat{S}_{00}^\epsilon)^{-1} \hat{S}_{0\perp}^\epsilon \hat{e}_\perp^\epsilon, $$

and the non-constant component is nothing else than

$$ \hat{e}_\perp^\epsilon = \mathbb{P}_\perp^{\epsilon,*}(\hat{S}^\epsilon)^{-1} \mathbb{P}_\perp^\epsilon (\hat{g}_\perp^\epsilon - \hat{S}_{\perp 0}^\epsilon \hat{\mu}_G^\epsilon) = \mathbb{P}_\perp^{\epsilon,*}(\hat{S}^\epsilon)^{-1} \mathbb{P}_\perp^\epsilon (\hat{g}_\perp^\epsilon - \hat{S}_{\perp 0}^\epsilon (\hat{S}_{00}^\epsilon)^{-1} \hat{g}_0^\epsilon). $$

Next, we estimate

$$ (4.4) \qquad \begin{aligned} \|\hat{e}_0^\epsilon\|_{-1/2} &\le \|(\hat{S}_{00}^\epsilon)^{-1}\| \|\hat{S}_{0\perp}^\epsilon\| \|\hat{e}_\perp^\epsilon\|_{-1/2}, \\ \|\hat{e}_\perp^\epsilon\|_{-1/2} &\le \|\mathbb{P}_\perp^{\epsilon,*}(\hat{S}^\epsilon)^{-1}\| \|\hat{g}_\perp^\epsilon\|_{1/2} + \|\hat{S}_{\perp 0}^\epsilon\| \|(\hat{S}_{00}^\epsilon)^{-1}\| \|\hat{g}_0^\epsilon\|_{1/2}. \end{aligned} $$

Thus, the convergence of the Galerkin Foldy-Lax model relies on the estimates of the data $\hat{g}^\epsilon$ and on the estimates on the operator norms of the (inverses of) restrictions of $\hat{S}^\epsilon$ to different subspaces, as $\epsilon \to 0$.

Let us anticipate the results that we will prove in the sections that follow. For the data, we will see that as $\epsilon \to 0$, it is a constant component of $\hat{g}^\epsilon$ that provides the most significant contribution in the norm of $\|\hat{g}^\epsilon\|_{H^{1/2}(\Gamma^\epsilon)}$:

$$ (4.5) \qquad \|\hat{g}_0^\epsilon\|_{1/2} = O(\epsilon^{1/2}), \quad \|\hat{g}_\perp^\epsilon\|_{1/2} = O(\epsilon). $$

For the operators, this is the case as well: the most significant contribution to $(\hat{S}^\epsilon)^{-1}$ is provided by the inverse of $\hat{S}_{00}^\epsilon$, which scales as $O(\epsilon^{-1})$. Moreover, as $\epsilon \to 0$, the operator matrix $\hat{S}^\epsilon$ approaches a diagonal operator matrix. This is reflected below:

$$ (4.6) \qquad \|(\hat{S}_{00}^\epsilon)^{-1}\| = O(\epsilon^{-1}), \quad \|\mathbb{P}_\perp^{\epsilon,*}(\hat{S}^\epsilon)^{-1}\| = O(1), \quad \|\hat{S}_{0\perp}^\epsilon\| = \|\hat{S}_{\perp 0}^\epsilon\| = O(\epsilon^{3/2}). $$

These estimates, combined with (4.4), will imply the following error behaviour:

$$ \|\hat{e}_\perp^\epsilon\|_{-1/2} = O(\epsilon), \quad \|\hat{e}_0^\epsilon\|_{-1/2} = O(\epsilon^{3/2}). $$

We proceed as follows. Because our final goal is to obtain time-domain estimates, we

* first obtain frequency-domain estimates on the operators, explicit in the frequency and in the asymptotic parameter;

* next obtain estimates on the data, directly in the time domain;

* translate all the estimates into the time domain.
**4.3. Estimates on the operators in the frequency domain.** The goal of this section is to prove the bounds (4.6), explicit in the frequency.

Our proof of the estimates on the operator norms $\|(\hat{S}_{00}^{\varepsilon})^{-1}\|$ and $\|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^{\varepsilon})^{-1}\mathbb{P}_{\perp}^{\varepsilon}\|$ will follow the path of [2]. It is based on three elements: 1) coercivity of the underlying operators; 2) their connection to solutions of the Helmholtz transmission/boundary-value problems (BVPs); 3) a so-called *lifting lemma*. More precisely, with the help of the first two ingredients, we will obtain the following bounds, with $C_0, C_{\perp} > 0$,

$$ (4.7) \qquad \begin{aligned} & \|(\hat{S}_{00}^{\varepsilon})^{-1}\| \le C_0 \sup_{\lambda \in V_0^{\varepsilon} \setminus \{0\}} L_{|\omega|}^{\varepsilon}(\lambda), \\ & \|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^{\varepsilon})^{-1}\mathbb{P}_{\perp}^{\varepsilon}\| \le C_{\perp} \sup_{\lambda \in H_*^{1/2}(\Gamma^{\varepsilon}) \setminus \{0\}} L_{|\omega|}^{\varepsilon}(\lambda), \\ & \text{where } L_{|\omega|}^{\varepsilon}(\lambda) = \inf_{\substack{\Lambda \in H^1(\mathbb{R}^2): \\ \gamma_0 \Lambda = \lambda}} \frac{\|\Lambda\|_{|\omega|, \mathbb{R}^2}}{\|\lambda\|_{1/2}}. \end{aligned} $$

Remark the difference in the sets over which supremum is taken in the bounds in the operator norms. The constants in the above bounds depend on frequency only, and all the information about the geometry and asymptotic parameter is contained in $L_{|\omega|}^{\varepsilon}(\lambda)$. Appropriate bounds on $L_{|\omega|}^{\varepsilon}(\lambda)$ are then given by the lifting lemma.

To prove the remaining estimates on $\|\hat{S}_{0\perp}^{\varepsilon}\|, \|\hat{S}_{00}^{\varepsilon}\|$, we will make use of the explicit representation of the operator $\hat{S}^{\varepsilon}$.

This section is organized as follows. In Section 4.3.1, we derive a lifting lemma (Proposition 4.4). This result is then used in Section 4.3.2 to obtain an estimate on $\|(\hat{S}_{00}^{\varepsilon})^{-1}\|$ and in Section 4.3.3 to obtain an estimate on $\|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^{\varepsilon})^{-1}\mathbb{P}_{\perp}^{\varepsilon}\|$. Finally, Section 4.3.4 is dedicated to the proof of the bounds on $\|\hat{S}_{0\perp}^{\varepsilon}\| \|\hat{S}_{00}^{\varepsilon}\|$.

**4.3.1. Lifting lemmas.** The goal of this section is to provide an estimate on the quantity (4.7), for the case when $\lambda \in V_0^{\varepsilon}$ or $H_*^{1/2}(\Gamma^{\varepsilon})$. Finding an optimal upper bound on (4.7) amounts to finding a lifting $\Lambda^{\varepsilon}$ of $\lambda$ that would minimize the norm $\|\cdot\|_{|\omega|, \mathbb{R}^2}$. Such a lifting is given by the solution of the boundary-value problem defined in the statement of the proposition below, with $a = |\omega|$.

**PROPOSITION 4.4 (Lifting lemma for multiple scatterers).** Let $a > 0$, $\lambda \in H^{1/2}(\Gamma^{\varepsilon})$, and $\Lambda^{\varepsilon} \in H^1(\mathbb{R}^2)$ satisfy the following two boundary-value problems:

$$ \begin{aligned} & -\Delta\Lambda^{\varepsilon} + a^2\Lambda^{\varepsilon} = 0 \text{ in } \Omega^{\varepsilon}, \quad \gamma_0^{-}\Lambda^{\varepsilon} = \lambda, \\ & -\Delta\Lambda^{\varepsilon} + a^2\Lambda^{\varepsilon} = 0 \text{ in } \Omega^{\varepsilon,c}, \quad \gamma_0^{+}\Lambda^{\varepsilon} = \lambda. \end{aligned} $$

Let us define

$$ (4.8) \qquad C_0(\varepsilon, \underline{d}_{\ast}^{\varepsilon}) := 1 + \frac{\varepsilon}{\underline{d}_{\ast}^{\varepsilon}}, \quad \underline{d}_{\ast}^{\varepsilon} = \min(1, \underline{d}_{\ast}^{\varepsilon}). $$

If $\lambda \in V_0^{\varepsilon}$, then, with a constant $c_{m,0}$ depending polynomially on $R_{\ast}^{-1}$ and $R^{*}$,

$$ (4.9) \qquad \|\Lambda^{\varepsilon}\|_{a, \mathbb{R}^2}^2 \le c_{m,0} \varepsilon^{-1} \max(1, \varepsilon a) C_0(\varepsilon, \underline{d}_{\ast}^{\varepsilon}) \|\lambda\|_{H^{1/2}(\Gamma^{\varepsilon})}^2. $$

If $\lambda \in H_*^{1/2}(\Gamma^{\varepsilon})$, then, with a constant $c_{m,\perp}$ depending polynomially on $R_{\ast}^{-1}$ and $R^{*}$,

$$ (4.10) \qquad \|\Lambda^{\varepsilon}\|_{a, \mathbb{R}^2}^2 \le c_{m,\perp} \max(1, \varepsilon a) C_0(\varepsilon, \underline{d}_{\ast}^{\varepsilon}) \|\lambda\|_{H^{1/2}(\Gamma^{\varepsilon})}^2. $$
### REMARK 9. The index 'm' in $c_{m,0}$, $c_{m,\perp}$ stands for 'multiple'.

In the above result we distinguish between the data from different subspaces of $H^{1/2}(\Gamma^\epsilon)$. This is important for the analysis, in particular since two obtained bounds differ by a factor of $\epsilon^{-1}$, see the discussion in the beginning of Section 4.3.

Remark as well that the constants in (4.9) and (4.10) are independent of the number of the particles $N$.

The rest of the section is dedicated to the proof of Proposition 4.4. The principal difficulty in the proof of the above result is to derive estimates explicit in the small parameter $\epsilon > 0$. To do so, we will base our considerations on the idea from the recent work by Hassan and Stamm [17], which addresses the three-dimensional electrostatic problem in the exterior of multiple spheres. They adapt the lifting lemma to their geometry, and we will make use of their construction, with some modifications due to the nature of our model.

We first obtain the lifting bounds for the case of a single scatterer, see Lemmas 4.6, 4.7. Next, we argue how these results can be extended to multiple scatterers.

Lifting lemmas for a single obstacle. We start by recalling the lifting lemma of [2], in the form presented in the monograph by Sayas [39].

PROPOSITION 4.5 (Proposition 2.5.1 in [39]). Let $\mathcal{O}$ be a Lipschitz domain. Then there exists $C_\mathcal{O} > 0$, s.t. for all $\xi \in H^{1/2}(\partial\mathcal{O})$ and all $a > 0$, the solution $v \in H^1(\mathcal{O})$ to the Dirichlet problem

$$-\Delta v + a^2 v = 0 \text{ in } \mathcal{O}, \quad \gamma_0 v = \xi,$$

is bounded as follows:

$$\|v\|_{a,\mathcal{O}} = (\|\nabla v\|_{L^2(\mathcal{O})}^2 + a^2 \|v\|_{L^2(\mathcal{O})}^2)^{1/2} \le C_{\mathcal{O}} \max(1, a^{1/2}) \|\xi\|_{H^{1/2}(\mathcal{O})}.$$

Our goal is to prove lifting lemmas analogous to the above proposition, when $\mathcal{O} = B(0, r)$ and where the constant $C_\mathcal{O}$ depends explicitly on $r$. We will also need to consider separately the two cases: when the data belongs to the space $H_*^{1/2}(C_r)$ and the data is constant on $C_r$, cf. Proposition 4.4. We start with the latter case.

In what follows, we use the notation

$$B_{r,r+d} = B(0, r+d) \setminus \overline{B(0, r)}.$$

LEMMA 4.6 (Lifting lemma for constants). Let $a, r > 0$, and $V^r \in H^1(\mathbb{R}^2)$ satisfy the interior boundary value problem

$$-\Delta V^r + a^2 V^r = 0 \text{ in } B(0, r), \quad \gamma_0^- V^r = 1,$$

and the exterior boundary value problem

$$-\Delta V^r + a^2 V^r = 0 \text{ in } B^c(0, r), \quad \gamma_0^+ V^r = 1.$$

Then, for all $d > 0$, $V^r$ satisfies the following bounds

$$(4.11) \qquad \|V^r\|_{a,\mathbb{R}^2}^2 \lesssim \max(1, ar),$$

$$(4.12) \qquad \|V^r\|_{L^2(B_{r,r+d})}^2 \lesssim rd + d^2.$$

*Proof.* The solution $V^r$ is defined via the modified Bessel functions:

$$(4.13) \qquad V^r(\rho) = \begin{cases} \frac{I_0(a\rho)}{I_0(ar)} & \rho < r, \\ \frac{K_0(a\rho)}{K_0(ar)} & \rho > r. \end{cases}$$
Proof of the bound (4.12). Using Lemma D.2 to bound $K_0(a\rho)/K_0(ar)$,

$$ \|V^r\|_{L^2(B_{r,r+d})}^2 = 2\pi \int_r^{r+d} \left| \frac{K_0(a\rho)}{K_0(ar)} \right|^2 \rho d\rho \lesssim d(r+d). $$

Proof of the bound (4.11). To prove this result, one could have used an explicit expression (4.13) and next proved some bounds on the modified Bessel functions. We have tried this approach, however, the final bound did not improve significantly over the bound obtained from a simple scaling argument, which we will present below.

In what follows, we use the fact that $V^r$ is rotation invariant, and we make its dependence on $a$ explicit by writing $V^r(\rho, a)$ instead of $V^r(\rho)$. By a scaling argument, one can verify that $V^r(\rho, a) = V^1(\rho, ar)$. Moreover (where we use $V^r \in H^1(\mathbb{R}^2)$):

$$ (4.14) \qquad \|V^r(\cdot, a)\|_{a, \mathbb{R}^2}^2 = \int_0^{+\infty} (|\partial_\rho V^r(\rho, a)|^2 + a^2|V^r(\rho, a)|^2) \rho d\rho = \|V^1(\cdot, ar)\|_{ar, \mathbb{R}^2}^2. $$

To obtain (4.11), it suffices to bound

$$ \|V^1(\cdot, ar)\|_{ar, \mathbb{R}^2}^2 = \|V^1(\cdot, ar)\|_{ar, B(0,1)}^2 + \|V^1(\cdot, ar)\|_{ar, B^c(0,1)}^2, $$

which is done by applying Proposition 4.5 twice, first with $O = B(0,1)$ and next with $O = B^c(0,1)$. $\square$

A similar scaling approach can be used to prove a lifting lemma for the case when the boundary data is orthogonal to constants.

LEMMA 4.7 (Lifting lemma for elements of $H_*^{1/2}(C_r)$). Given $a, r > 0$, $g^r \in H_*^{1/2}(C_r)$, let $V^r \in H^1(\mathbb{R}^2)$ satisfy the interior boundary value problem

$$ -\Delta V^r + a^2 V^r = 0 \text{ in } B(0, r), \quad \gamma_0^- V^r = g^r, $$

and the exterior boundary value problem

$$ -\Delta V^r + a^2 V^r = 0 \text{ in } B^c(0, r), \quad \gamma_0^+ V^r = g^r. $$

Then, for all $d > 0$, $V^r$ satisfies the following bounds:

$$ (4.15) \qquad \|V^r\|_{a, \mathbb{R}^2}^2 \lesssim \max(1, ar) \|g^r\|_{H^{1/2}(C_r)}^2, $$

$$ (4.16) \qquad \|V^r\|_{L^2(B_{r,r+d})}^2 \lesssim rd \|g^r\|_{H^{1/2}(C_r)}^2. $$

Proof. Provided that $g^r = \sum_{m \in \mathbb{Z} \setminus \{0\}} g_m^r e^{im\theta}$, the solution to the above BVPs can be written explicitly, as a series of modified Bessel functions:

$$ (4.17) \qquad V^r(\rho, \theta) = \sum_{m \in \mathbb{Z} \setminus \{0\}} v_m^r(\rho) e^{im\theta}, \quad v_m^r(\rho) = \begin{cases} \frac{I_{|m|}(a\rho)}{I_{|m|}(ar)} g_m^r e^{im\theta}, & \rho < r, \\ \frac{K_{|m|}(a\rho)}{K_{|m|}(ar)} g_m^r e^{im\theta}, & \rho > r. \end{cases} $$

Proof of the bound (4.16). We start by bounding $\|v_m^r\|_{L^2(B_{r,r+d})}^2$. With Lemma D.3,

$$ \int_r^{r+d} \left| \frac{K_m(a\rho)}{K_m(ar)} \right|^2 \rho d\rho \lesssim r^{2m} \int_r^{r+d} \rho^{-2m+1} d\rho \lesssim r \int_r^{r+d} d\rho \lesssim rd, \quad \forall m \ge 1. $$
From the definition of $V^r(\rho, \theta)$, it is straightforward to see that

$$ \|V^r\|_{L^2(B_{r,r+d})}^2 \lesssim rd \sum_{m \in \mathbb{Z} \setminus \{0\}} |g_m^r|^2 \lesssim rd \|g^r\|_{H^{1/2}(C^r)}^2. $$

**Proof of the bound (4.15).** For the same reason as in the proof of Lemma 4.6, we will use the scaling argument to prove (4.15).

First of all, let $g^1(\mathbf{x}) := g^r(r\mathbf{x})$, for $\mathbf{x} \in C_1$. Let $V^1(\rho, b)$ satisfy:

$$ - \Delta V^1 + b^2 V^1 = 0 \text{ in } B(0, 1), \quad \gamma_0^{-} V^1 = g^1, $$

$$ - \Delta V^1 + b^2 V^1 = 0 \text{ in } B^c(0, 1), \quad \gamma_0^{+} V^1 = g^1. $$

Again, by the scaling argument, one can verify that $V^r(\rho, a) = V^1(\rho, ar)$, and one has the norm identity (4.14). It remains to use Proposition 4.5 which yields

$$ (4.18) \qquad \|V^r(\cdot, a)\|_{a, \mathbb{R}^2}^2 \le C \max(1, ar) \|g^1\|_{H^{1/2}(C_1)}^2. $$

Here the constant $C$ is independent of $r$, but encodes some information about the circular shape of the scatterer. The functions $g^1, g^r$ are defined by their Fourier series, for $\mathbf{x} = (\cos \theta, \sin \theta)$,

$$ g^1(\mathbf{x}) = g^r(r\mathbf{x}) = \sum_{m \in \mathbb{Z} \setminus \{0\}} g_m^1 e^{im\theta} = \sum_{m \in \mathbb{Z} \setminus \{0\}} g_m^r e^{im\theta}. $$

We remark that

$$ \begin{aligned} \|g^1\|_{H^{1/2}(C_1)}^2 &= \sum_{m \in \mathbb{Z} \setminus \{0\}} (1+m^2)^{1/2} |g_m^1|^2 \le 2 \sum_{m \in \mathbb{Z} \setminus \{0\}} |m| |g_m^1|^2 \\ &\le 2r \sum_{m \in \mathbb{Z}} (1+\frac{m^2}{r^2})^{1/2} |g_m^r|^2 = 2 \|g^r\|_{H^{1/2}(C_r)}^2, \end{aligned} $$

which, with (4.18), yields the desired inequality. $\square$

**REMARK 10.** Comparing the results of Lemmas 4.6 and 4.7, it may seem that the bounds on $\|V^r\|_{a, \mathbb{R}^2}$ are identical with respect to the norm of the boundary data $\|g\|_{H^{1/2}(C_r)}$ and the scaling parameter $r$ (with $g$ being equal to 1 or $g^r$). This is however not the case: the bound (4.11) can alternatively be rewritten as $\|V^r\|_{a, \mathbb{R}^2}^2 \lesssim \max(1, ar)r^{-1}\|g\|_{H^{1/2}(C_r)}^2$.

To prove the analogous results for multiple circles (Proposition 4.4), let us introduce an auxiliary truncation function, following [17].

*An auxiliary truncation function.* Let

$$ (4.19) \qquad \chi_{r,d}(\rho) = \begin{cases} 1, & \rho < r + d/8, \\ 0, & \rho > r + d/4, \\ p_{r,d}(\rho), & \text{otherwise,} \end{cases} $$

with $p_{r,d} = \sum_{k=0}^7 a_k(\rho - r - d/8)^k$, with the coefficients $a_k$ chosen so that $\chi_{r,d} \in C^3(\mathbb{R}_{\ge 0})$. One can verify that $a_k = c_k/d^k$ (with some $c_k$ independent of $d$ and $r$), for $k \ge 4$, and $a_3 = a_2 = a_1 = a_0 = 0$. We then introduce

$$ (4.20) \qquad \chi_n^\epsilon(\mathbf{x}) := \chi_{r_n^\epsilon, d_n^\epsilon}(\|\mathbf{x} - \mathbf{c}_n\|), \quad \chi^\epsilon(\mathbf{x}) := \sum_{n \in \mathcal{N}} \chi_n^\epsilon(\mathbf{x}). $$
With this definition,

$$ (4.21) \qquad \operatorname{supp} \chi_n^\varepsilon \cap \operatorname{supp} \chi_k^\varepsilon = \emptyset, \quad \text{if } n \neq k, $$

$$ (4.22) \qquad \operatorname{supp} \chi_n^\varepsilon \subseteq \bar{B}(c_n, r_n^\varepsilon + \underline{d}_*^\varepsilon / 2). $$

It is then straightforward to see that

$$ (4.23) \qquad \|\chi^\varepsilon\|_{L^\infty(\mathbb{R}^2)} \lesssim 1, \quad \|\chi^\varepsilon\|_{W^{s,\infty}(\mathbb{R}^2)} \lesssim (\underline{d}_*^\varepsilon)^{-s}, \quad s \ge 0. $$

**Proof of Proposition 4.4.** Remark that $\Lambda^\varepsilon = \arg \min_{v \in H^1(\mathbb{R}^2)}: \gamma_0^\pm v = \lambda \|v\|_{a,\mathbb{R}^2}$. Thus, to bound $\|\Lambda^\varepsilon\|_{a,\mathbb{R}^2}$, we can construct a lifting $v \in H^1(\mathbb{R}^2)$ of $\lambda$, obtain an explicit stability bound $\|v\|_{a,\mathbb{R}^2} \le C_\varepsilon \|\lambda\|_{H^{1/2}(\Gamma^\varepsilon)}$, and conclude that $\|\Lambda^\varepsilon\|_{a,\mathbb{R}^2} \le C_\varepsilon \|\lambda\|_{H^{1/2}(\Gamma^\varepsilon)}$. Let us define the following lifting of $\lambda$:

$$ \Lambda_\chi^\varepsilon := \sum_{n \in \mathcal{N}} \Lambda_n^\varepsilon \chi_n^\varepsilon, $$

where $\chi_n^\varepsilon$ are as in (4.20) and each $\Lambda_n^\varepsilon \in H^1(\mathbb{R}^2)$ solves the following exterior and interior boundary value problems, with $\lambda_n \in H^{1/2}(\Gamma_n^\varepsilon)$,

$$ (4.24) \qquad \begin{aligned} & -\Delta \Lambda_n^\varepsilon + a^2 \Lambda_n^\varepsilon = 0 \quad \text{in } \Omega_n^\varepsilon, \quad \gamma_0^- \Lambda_n^\varepsilon = \lambda_n, \\ & -\Delta \Lambda_n^\varepsilon + a^2 \Lambda_n^\varepsilon = 0 \quad \text{in } \Omega_n^{\varepsilon,c}, \quad \gamma_0^+ \Lambda_n^\varepsilon = \lambda_n. \end{aligned} $$

We next relate $\|\Lambda_\chi^\varepsilon\|_{a,\mathbb{R}^2}$ to the norms $\|\Lambda_n^\varepsilon\|_{a,\mathbb{R}^2}$, which we will bound based on Lemmas 4.6, 4.7. By definition,

$$ \|\Lambda_\chi^\varepsilon\|_{a,\mathbb{R}^2}^2 = a^2 \left\| \sum_{n \in \mathcal{N}} \chi_n^\varepsilon \Lambda_n^\varepsilon \right\|^2 + \left\| \sum_{n \in \mathcal{N}} \chi_n^\varepsilon \nabla \Lambda_n^\varepsilon \right\|^2 + \left\| \sum_{n \in \mathcal{N}} \nabla \chi_n^\varepsilon \Lambda_n^\varepsilon \right\|^2. $$

By the property (4.21), we have that

$$ (4.25) \qquad \begin{aligned} & \|\Lambda_\chi^\varepsilon\|_{a,\mathbb{R}^2}^2 = a^2 \sum_{n \in \mathcal{N}} \|\chi_n^\varepsilon \Lambda_n^\varepsilon\|^2 + \sum_{n \in \mathcal{N}} \|\chi_n^\varepsilon \nabla \Lambda_n^\varepsilon\|^2 + \sum_{n \in \mathcal{N}} \|\nabla \chi_n^\varepsilon \Lambda_n^\varepsilon\|^2 \\ & \le \sum_{n \in \mathcal{N}} \|\chi_n^\varepsilon\|_{L^\infty}^2 \|\Lambda_n^\varepsilon\|_{a,\mathbb{R}^2}^2 + \sum_{n \in \mathcal{N}} \|\nabla \chi_n^\varepsilon\|_{L^\infty}^2 \|\Lambda_n^\varepsilon\|_{\text{supp } \nabla \chi_n^\varepsilon}^2 \\ & \stackrel{(4.21)}{=} \|\chi^\varepsilon\|_{L^\infty}^2 \sum_{n \in \mathcal{N}} \|\Lambda_n^\varepsilon\|_{a,\mathbb{R}^2}^2 + \|\nabla \chi^\varepsilon\|_{L^\infty}^2 \sum_{n \in \mathcal{N}} \|\Lambda_n^\varepsilon\|_{\text{supp } \nabla \chi_n^\varepsilon}^2 \\ & \stackrel{(4.23)}{\lesssim} \sum_{n \in \mathcal{N}} \|\Lambda_n^\varepsilon\|_{a,\mathbb{R}^2}^2 + (\underline{d}_*^\varepsilon)^{-2} \sum_{n \in \mathcal{N}} \|\Lambda_n^\varepsilon\|_{\text{supp } \nabla \chi_n^\varepsilon}^2. \end{aligned} $$

It remains to bound all the above quantities, by recalling that, see (4.22), $\operatorname{supp} \nabla \chi_n^\varepsilon \subset \overline{B(c_n, r_n^\varepsilon + \underline{d}_*^\varepsilon / 2)} \setminus B(c_n, r_n^\varepsilon)$. We consider two cases: $\lambda \in \mathcal{V}_0^\varepsilon$ and $\lambda \in H_*^{1/2}(\Gamma^\varepsilon)$.

**Case 1.** Bound (4.9) for $\lambda \in \mathcal{V}_0^\varepsilon$. Lemma 4.6 yields

$$ \|\Lambda_\chi^\varepsilon\|_{a,\mathbb{R}^2}^2 \lesssim \sum_{n \in \mathcal{N}} \max(1, ar_n^\varepsilon) |\lambda_{n,0}|^2 + (\underline{d}_*^\varepsilon)^{-2} \sum_{n \in \mathcal{N}} \underline{d}_*^\varepsilon (r_n^\varepsilon + \underline{d}_*^\varepsilon) |\lambda_{n,0}|^2. $$
Because $\|\lambda\|_{H^{1/2}(\Gamma^\epsilon)}^2 = 2\pi \sum_{n \in \mathcal{N}} r_n^\epsilon |\lambda_{n,0}|^2$, we have the following bound, with a constant $C_1 > 0$ depending on $R_*^{-1}$ and $R^*$ polynomially:

$$ \|\Lambda_{\chi}^{\epsilon}\|_{a, \mathbb{R}^2}^2 \le C_1 \epsilon^{-1} \left( \max(1, \epsilon a) + \frac{\epsilon}{\underline{d}_{*}^{\epsilon}} + 1 \right) \|\lambda\|_{H^{1/2}(\Gamma^{\epsilon})}^2. $$

With the use of $1 \le \max(1, a\epsilon)$, the bound of the proposition follows straightforwardly from the above and $\|\Lambda^\epsilon\|_{a, \mathbb{R}^2} \le \|\Lambda_\chi^\epsilon\|_{a, \mathbb{R}^2}$.

Case 2. Bound (4.10) for $\lambda \in H_*^{1/2}(\Gamma^\epsilon)$. An application of bounds of Lemma 4.7 to bound the right-hand side of the inequality (4.25) yields the following bound, with a constant $C_2$ depending on on $R_*^{-1}$ and $R^*$ polynomially:

$$ \|\Lambda_{\chi}^{\epsilon}\|_{a, \mathbb{R}^2}^2 \le C_2 \left( \max(1, \epsilon a) + \frac{\epsilon}{\underline{d}_{*}^{\epsilon}} \right) \|\lambda\|_{H^{1/2}(\Gamma^{\epsilon})}^2. $$

The bound in the statement of the proposition is then obtained like in Case 1. $\square$

**4.3.2. An estimate on** $\| (\hat{S}_{00}^{\epsilon})^{-1} \|$. Our goal is to prove that $\| (\hat{S}_{00}^{\epsilon})^{-1} \|$ behaves as $O(\epsilon^{-1})$ as $\epsilon \to 0$. In the proof of this result, we will make use of the connection between the solution to a certain transmission problem and the potential $\hat{S}^\epsilon$.

Given $\omega \in \mathbb{C}^+$, $\eta \in H^{1/2}(\Gamma^\epsilon)$, we denote $v_\eta^\epsilon := \hat{S}^\epsilon \eta \in H^1(\mathbb{R}^2)$. We remark that $v_\eta^\epsilon$ solves the following transmission problem:

$$ \begin{gathered} -\Delta v_{\eta}^{\epsilon} - \omega^2 v_{\eta}^{\epsilon} = 0 \text{ in } \mathbb{R}^2 \setminus \Gamma^{\epsilon}, \\ [\gamma_0 v_{\eta}^{\epsilon}] = 0, \quad [\gamma_1 v_{\eta}^{\epsilon}] = \eta, \end{gathered} \quad (4.26) $$

and, moreover, $\gamma_0 v_\eta^\epsilon = \hat{S}^\epsilon \eta$.

**THEOREM 4.8.** Let $\omega \in \mathbb{C}^+$. The operator $\hat{S}_{00}^\epsilon : \mathcal{V}_0^\epsilon \to \mathcal{V}_0^\epsilon$ satisfies, with $c_{m,0}$, $C_0(\epsilon, \underline{d}_*)$ as in Proposition 4.4,

$$ \| (\hat{S}_{00}^{\epsilon})^{-1} \| \le c_{m,0} C_0(\epsilon, \underline{d}_*) \epsilon^{-1} \max(1, \epsilon|\omega|) \frac{|\omega|}{\operatorname{Im} \omega}. $$

*Proof.* As usual in the theory of time-domain boundary integral equations, we will prove a coercivity bound on $\hat{S}_{00}^\epsilon$, cf. [2]. Let $\eta \in \mathcal{V}_0^\epsilon$. Then

$$ \mathbb{P}_0^\epsilon \gamma_0 v_\eta^\epsilon = \mathbb{P}_0^\epsilon \hat{S}^\epsilon \mathbb{P}_0^{\epsilon,*} \eta = \hat{S}_{00}^\epsilon \eta. $$

Therefore, since $\eta \in \mathcal{V}_0^\epsilon$,

$$ \langle \eta, \overline{\hat{S}^\epsilon \eta} \rangle = \langle \eta, \overline{\hat{S}_{00}^\epsilon \eta} \rangle. $$

Hence,

$$ \begin{align} (4.27) \qquad & - \operatorname{Im} \langle \eta, \overline{\hat{S}_{00}^{\epsilon} \eta} \rangle = - \operatorname{Im} \langle \eta, \overline{\hat{S}^{\epsilon} \eta} \rangle = - \operatorname{Im} \langle [\gamma_1 v_{\eta}^{\epsilon}], \overline{\omega \gamma_0 v_{\eta}^{\epsilon}} \rangle \nonumber \\ & = - \operatorname{Im} (\overline{\omega} \|\nabla v_{\eta}^{\epsilon}\|^2 - \omega |\omega|^2 \|v_{\eta}^{\epsilon}\|^2) = \operatorname{Im} \omega \|v_{\eta}^{\epsilon}\|_{\omega, \mathbb{R}^2}^2, \end{align} $$
where we used the Green formula and (4.26). Next, we need to relate the right-hand side of the above to $\eta$, i.e. to prove the bound $\|\eta\|_{H^{-1/2}(\Gamma^\varepsilon)} \le C_\omega \|v_\eta^\varepsilon\|_{\omega|,\mathbb{R}^2}$ with an explicit dependence on $\omega$ in the constant $C_\omega$. This is done in a classical way:

$$ \|\eta\|_{H^{-1/2}(\Gamma^\varepsilon)} = \sup_{\lambda \in V_0^\varepsilon \setminus \{0\}} \frac{\langle \eta, \bar{\lambda} \rangle}{\|\lambda\|_{H^{1/2}(\Gamma^\varepsilon)}} = \sup_{\lambda \in V_0^\varepsilon \setminus \{0\}} \inf_{\Lambda \in H^1(\mathbb{R}^2): \gamma_0 \Lambda = \lambda} \frac{\langle \eta, \overline{\gamma_0 \Lambda} \rangle}{\|\lambda\|_{H^{1/2}(\Gamma^\varepsilon)}}. $$

Again, by the Green identity,

$$ \langle \eta, \overline{\gamma_0 \Lambda} \rangle = \int_{\mathbb{R}^2} \nabla v_\eta^\varepsilon \overline{\nabla \Lambda} dx - \omega^2 \int_{\mathbb{R}^2} v_\eta^\varepsilon \overline{\Lambda} dx. $$

Then

$$ (4.28) \qquad \|\eta\|_{H^{-1/2}(\Gamma^\varepsilon)} \le \|v_\eta^\varepsilon\|_{\omega|,\mathbb{R}^2} \sup_{\lambda \in V_0^\varepsilon \setminus \{0\}} \inf_{\Lambda \in H^1(\mathbb{R}^2): \gamma_0 \Lambda = \lambda} \frac{\|\Lambda\|_{\omega|,\mathbb{R}^2}}{\|\lambda\|_{H^{1/2}(\Gamma^\varepsilon)}}. $$

The $\inf_{\Lambda \in H^1(\mathbb{R}^2): \gamma_0 \Lambda = \lambda} \|\Lambda\|_{\omega|,\mathbb{R}^2}$ is realized for $\Lambda = \Lambda^\varepsilon$, where $\Lambda^\varepsilon$ is from Proposition 4.4, with $a = |\omega|$. With the notation of Proposition 4.4, we obtain

$$ \|\eta\|_{H^{-1/2}(\Gamma^\varepsilon)}^2 \le c_{m,0} \varepsilon^{-1} \max(1, |\varepsilon|) C_0(\varepsilon, d_*^\varepsilon) \|v_\eta^\varepsilon\|_{\omega|,\mathbb{R}^2}^2. $$

Therefore, from (4.27) and the above, we obtain the coercivity bound

$$ - \operatorname{Im} \langle \eta, \omega \overline{\hat{S}_{00}^\varepsilon \eta} \rangle \ge \varepsilon c_{m,0}^{-1} \min(1, (\varepsilon|\omega|)^{-1}) C_0^{-1}(\varepsilon, d_*^\varepsilon) \operatorname{Im} \omega \|\eta\|_{H^{-1/2}(\Gamma^\varepsilon)}^2, $$

and the bound for the inverse $(\hat{S}_{00}^\varepsilon)^{-1}$ stated in the theorem follows from the above. $\square$

**4.3.3. An estimate on** $\|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^\varepsilon)^{-1}\mathbb{P}_{\perp}^\varepsilon\|$. In the theorem below we show that $\|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^\varepsilon)^{-1}\mathbb{P}_{\perp}^\varepsilon\| = O(1)$, and thus this term does not have an effect on the asymptotic error estimate.

**THEOREM 4.9.** Let $\omega \in \mathbb{C}^+$. The operator $\hat{S}^\varepsilon: H^{-1/2}(\Gamma^\varepsilon) \to H^{1/2}(\Gamma^\varepsilon)$ satisfies, with $c_{m,\perp}$ like in Proposition 4.4,

$$ \|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^\varepsilon)^{-1}\mathbb{P}_{\perp}^\varepsilon\| \le c_{m,\perp} \frac{|\omega|}{\operatorname{Im} \omega} C_0(\varepsilon, d_*^\varepsilon) \max(1, |\varepsilon|\omega). $$

*Proof.* We have that

$$ \|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^\varepsilon)^{-1}\mathbb{P}_{\perp}^\varepsilon\|_{H^{1/2}(\Gamma^\varepsilon) \to H^{-1/2}(\Gamma^\varepsilon)} = \|\mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^\varepsilon)^{-1}\|_{H_*^{1/2}(\Gamma^\varepsilon) \to H^{-1/2}(\Gamma^\varepsilon)}. $$

Let $\lambda \in H_*^{1/2}(\Gamma^\varepsilon)$, and $\eta = (\hat{S}^\varepsilon)^{-1}\lambda \equiv (\hat{S}^\varepsilon)^{-1}\mathbb{P}_{\perp}^\varepsilon\lambda$. Our goal is to find an estimate on $\eta_\perp = \mathbb{P}_{\perp}^{\varepsilon,*}(\hat{S}^\varepsilon)^{-1}\mathbb{P}_{\perp}^\varepsilon\lambda$ by $\lambda$. For this we remark that

$$ \langle \eta, \overline{\hat{S}^\varepsilon \eta} \rangle = \langle \eta, \overline{\lambda} \rangle = \langle \eta_\perp, \overline{\lambda} \rangle, $$

because $\lambda \in H_*^{1/2}(\Gamma^\varepsilon)$. We proceed now like in the proof of Theorem 4.8. We define $v_\eta^\varepsilon := \hat{S}^\varepsilon \eta \in H^1(\mathbb{R}^2)$ the solution to the transmission problem (4.26). The Green identity yields

$$ - \operatorname{Im} \langle \eta_\perp, \omega \overline{\lambda} \rangle = \operatorname{Im} \omega \|v_\eta^\varepsilon\|_{\omega|,\mathbb{R}^2}^2. $$
It remains to estimate $\eta_{\perp}$ by $\|v_{\eta}^{\epsilon}\|_{\omega, \mathbb{R}^2}$, like in (4.28):

$$ \|\eta_{\perp}\| \le \|v_{\eta}^{\epsilon}\|_{\omega, H^1(\mathbb{R}^2)} \sup_{\lambda \in H_*^{1/2}(\Gamma^{\epsilon}) \setminus \{0\}} \inf_{\substack{\Lambda \in H^1(\mathbb{R}^2) \\ \gamma_0 \Lambda = \lambda}} \frac{\|\Lambda\|_{\omega, \mathbb{R}^2}}{\|\lambda\|_{H^{1/2}(\Gamma^{\epsilon})}}. $$

The desired estimate again follows from Proposition 4.4, cf. proof of Theorem 4.8. We have

$$ - \operatorname{Im} \langle \eta_{\perp}, \omega \lambda \rangle \ge c_{m, \perp}^{-1} \operatorname{Im} \omega \min(1, (\varepsilon|\omega|)^{-1}) C_0^{-1} (\varepsilon, d_*)^{\varepsilon} \|\eta_{\perp}\|_{H^{-1/2}(\Gamma^{\varepsilon})}^2, $$

which yields the estimate in the statement of the theorem.

**4.3.4. An estimate on $\|\hat{S}_{0\perp}^{\epsilon}\|$ and $\|\hat{S}_{\perp 0}^{\epsilon}\|$**. The upper bounds on the norms of these operators will be obtained in a different manner compared to the proofs of Theorems 4.8 and 4.9. The first result shows that it is sufficient to obtain a bound on either of the norms $\|\hat{S}_{0\perp}^{\epsilon}\|$ or $\|\hat{S}_{\perp 0}^{\epsilon}\|$.

$$ \textbf{PROPOSITION 4.10.} \quad \|\hat{S}_{0\perp}^{\epsilon}\| = \|\hat{S}_{\perp 0}^{\epsilon}\|. $$

*Proof.* We start with

$$ \|\hat{S}_{0\perp}^{\epsilon}\| = \sup_{\mu \in V_0^{\epsilon}} \sup_{\lambda \in H_*^{-1/2}(\Gamma^{\epsilon})} \frac{\langle \mu, \hat{S}^{\epsilon} \lambda \rangle}{\|\lambda\|_{H^{-1/2}(\Gamma^{\epsilon})} \|\mu\|_{H^{-1/2}(\Gamma^{\epsilon})}}, $$

$$ \|\hat{S}_{\perp 0}^{\epsilon}\| = \sup_{\mu \in V_0^{\epsilon}} \sup_{\lambda \in H_*^{-1/2}(\Gamma^{\epsilon})} \frac{\langle \lambda, \hat{S}^{\epsilon} \mu \rangle}{\|\lambda\|_{H^{-1/2}(\Gamma^{\epsilon})} \|\mu\|_{H^{-1/2}(\Gamma^{\epsilon})}}. $$

With the notation like in (4.26), third Green's identity yields

$$ \langle \boldsymbol{\mu}, \hat{S}^{\epsilon} \boldsymbol{\lambda} \rangle - \langle \boldsymbol{\lambda}, \hat{S}^{\epsilon} \boldsymbol{\mu} \rangle = \int_{\mathbb{R}^2 \setminus \Gamma^{\epsilon}} (\Delta v_{\mu}^{\epsilon} + \omega^2 v_{\mu}^{\epsilon}) v_{\lambda}^{\epsilon} - \int_{\mathbb{R}^2 \setminus \Gamma^{\epsilon}} (\Delta v_{\lambda}^{\epsilon} + \omega^2 v_{\lambda}^{\epsilon}) v_{\mu}^{\epsilon} = 0, $$

and hence the desired result.

Before proving the main theorem of this section, let us state an auxiliary lemma, which we will use also later.

**LEMMA 4.11.** Let $\omega \in \mathbb{C}^+$, $r > 0$, and the point $x \in B^c(0, r)$ be such that $\mathrm{dist}(x, B(0, r)) = d$. Let $h: C_r \to \mathbb{C}$ be defined by $h(y) = G_\omega(y - x)$. Then the following bounds hold true:

$$ \|h\|_{L^2(C_r)} \lesssim r^{1/2} \min \left( \max \left( 1, \log \frac{1}{d \operatorname{Im} \omega} \right), |\omega d|^{-1/2} \right), $$

$$ \|h_{\perp}\|_{H^{1/2}(C_r)} \lesssim r \max(1, r^{1/2}) d^{-1} \max(1, |\omega d|^{1/2}). $$

*Proof.* Evidently, $\|h\|_{L^2(C_r)}^2 \lesssim \|h\|_{L^\infty(C_r)}^2 r$. For all $y \in C_r$, the function $\omega \mapsto G_\omega(y)$ is analytic in $\omega \in \mathbb{C}^+$ and its restriction to $\mathbb{C}^+$ is continuous up to the set $\mathbb{R} \setminus \{0\}$. Therefore, we make use of the following bound that stems from the asymptotic behaviour of the Hankel functions for $\omega \in \mathbb{C}^+$, when $\omega \to 0$ and when $|\omega| \to +\infty$, cf. [13, Section 10.8, Section 10.17(iv)] (remark that the asymptotic expansions are uniform in Arg $\omega$)):

$$ |h(y)| \lesssim \min \left( \max \left( 1, \log \frac{1}{|\omega| \|x - y\|} \right), \frac{1}{|\omega|^{1/2} \|x - y\|^{1/2}} \right). $$

To obtain the desired bound we use $\operatorname{Im} \omega < |\omega|$.
To compute $\|h_{\perp}\|_{H^{1/2}(C_r)}$, we proceed as follows. First of all, we compute the $L^2$-norm:

$$ \begin{aligned} \|h_{\perp}\|_{L^2(C_r)}^2 &= \inf_{c \in \mathbb{C}} \|h-c\|^2 \le \|h(\cdot)-h(0)\|_{L^2(C_r)}^2 \\ &\lesssim \|\nabla h\|_{L^{\infty}(B(0,r))}^2 r^3 = \|\nabla G_{\omega}(\cdot, \mathbf{x})\|_{L^{\infty}(B(0,r))}^2 r^3. \end{aligned} $$

Because $\frac{d}{dz} H_0^{(0)}(z) = -H_1^{(0)}(z)$,

$$ |\nabla_y G_{\omega}(\mathbf{y} - \mathbf{x})| \lesssim |\omega| |H_1^{(1)}(\omega \|\mathbf{y} - \mathbf{x}\|)| \lesssim |\omega| \max \left( \frac{1}{|\omega| \|\mathbf{y} - \mathbf{x}\|}, \frac{1}{\sqrt{|\omega| \|\mathbf{y} - \mathbf{x}\|}} \right), $$

and we again used the fact that the restriction $z \mapsto H_1^{(1)}(z)$ to $\mathbb{C}^+$ is analytic, and continuous up to $\mathbb{R} \setminus \{0\}$, as well as the asymptotics of the Hankel functions, cf. [13, Section 10.8, Section 10.17(iv)]. This finally yields

$$ \|h_{\perp}\|_{L^2(C_r)}^2 \lesssim r^3 d^{-2} \max(1, |\omega|d). $$

Next, to obtain a complete bound for $\|h_{\perp}\|_{H^{1/2}(C_r)}$, it remains to compute its Sobolev-Slobodeckii seminorm, where we again use the Lipschitz regularity of $h$:

$$ \begin{aligned} \|h_{\perp}\|_{H^{1/2}(C_r)}^2 &= \int_{C_r} \int_{C_r} \frac{|h(\mathbf{y}) - h(\mathbf{y}')|^2}{\|\mathbf{y} - \mathbf{y}'\|^2} d\Gamma_{\mathbf{y}} d\Gamma_{\mathbf{y}'} \lesssim \|\nabla h\|_{L^{\infty}(B(0,r))}^2 \int_{C_r} d\Gamma_{\mathbf{y}} d\Gamma_{\mathbf{y}'} \\ &\lesssim r^2 d^{-2} \max(1, |\omega|d). \end{aligned} \quad \square $$

Finally we have all the necessary results to formulate and prove the principal result of this section.

**THEOREM 4.12.** Let $\omega \in \mathbb{C}^+$. Then the following bound holds for all $\varepsilon > 0$, with a constant $C_{0\perp}$ depending on $R_*^{-1}$ and $R^*$ polynomially:

$$ \|\hat{S}_{\perp 0}^{\varepsilon}\| = \|\hat{S}_{0\perp}^{\varepsilon}\| \le C_{0\perp} \varepsilon^{3/2} N(d_{\perp}^{\varepsilon})^{-1} |\omega|^{1/2} \max(1, (\mathrm{Im} \, \omega)^{-1/2}). $$

*Proof.* By Proposition 4.10, it suffices to prove the bound on $\hat{S}_{\perp 0}^{\varepsilon}$. This will be done using its explicit representation.

Let $\eta = \hat{S}^{\varepsilon}\lambda$, with $\lambda \in H_*^{-1/2}(\Gamma^{\varepsilon})$. Our goal is to bound the constant component $\eta_0$ by $\lambda$; by the density argument it suffices to obtain the bound for $\lambda$ sufficiently regular (e.g. $\lambda \in C^0(\Gamma^{\varepsilon})$). On $\Gamma_{\ell}^{\varepsilon}$, it holds that

$$ \eta_{\ell} = \hat{S}_{\ell\ell}^{\varepsilon} \lambda_{\ell} + \sum_{n \in \mathbb{N} \setminus \{\ell\}} \hat{S}_{\ell n}^{\varepsilon} \lambda_n, $$

see (3.11) for the definition of the above operators. Let us consider the first term in the above. Because the eigenfunctions of $\hat{S}_{\ell}$ are given by $\{e^{im\phi}\}_{m \in \mathbb{Z}}$, see e.g. [8] and references therein, and $\lambda \in H_*^{-1/2}(\Gamma^{\varepsilon})$, it follows that $\hat{S}_{\ell\ell}^{\varepsilon} \lambda_{\ell} \in H_*^{1/2}(\Gamma_{\ell}^{\varepsilon})$. Therefore, we have, after integrating both sides of the above over $\Gamma_{\ell}^{\varepsilon}$,

$$ r_{\ell}^{\varepsilon} \eta_{\ell,0} = \frac{1}{2\pi} \sum_{n \in \mathbb{N} \setminus \{\ell\}} \int_{\Gamma_{\ell}^{\varepsilon}} \int_{\Gamma_{n}^{\varepsilon}} G_{\omega}(\mathbf{x} - \mathbf{y}) \lambda_{n}(\mathbf{y}) d\Gamma_{\mathbf{y}} d\Gamma_{\mathbf{x}}. $$
It remains to use the definition of $\|\eta_0\|_{L^2(\Gamma^\epsilon)}^2$:

$$ (4.29) \qquad \begin{aligned} \|\eta_0\|_{L^2(\Gamma^\epsilon)}^2 & \lesssim \sum_{\ell \in \mathcal{N}} (r_\ell^\epsilon)^{-1} \left| \sum_{n \in \mathcal{N} \setminus \{\ell\}} \iint_{\Gamma_\ell^\epsilon \times \Gamma_n^\epsilon} G_\omega(\mathbf{x} - \mathbf{y}) \lambda_n(\mathbf{y}) d\Gamma_\mathbf{y} d\Gamma_\mathbf{x} \right|^2 \\ & \lesssim R_*^{-1} \epsilon^{-1} \sum_{\ell \in \mathcal{N}} \left| \sum_{n \in \mathcal{N} \setminus \{\ell\}} r_\ell^\epsilon \sup_{\mathbf{x} \in \Gamma_\ell^\epsilon} \left| \int_{\Gamma_n^\epsilon} H_0^{(1)}(\omega \|\mathbf{x} - \mathbf{y}\|) \lambda_n(\mathbf{y}) d\Gamma_\mathbf{y} \right| \right|^2. \end{aligned} $$

We then use the bound of Lemma 4.11 to obtain, with $C_1$ depending polynomially on $R^*$ and $R_*^{-1}$, for $\mathbf{x} \in \Gamma_\ell^\epsilon$:

$$ \left| \int_{\Gamma_n^\epsilon} H_0^{(1)}(\omega \|\mathbf{x} - \mathbf{y}\|) \lambda_n(\mathbf{y}) d\Gamma_\mathbf{y} \right| \lesssim \|G_\omega(\mathbf{x} - \cdot)\|_{H^{1/2}(\Gamma_n^\epsilon)} \|\lambda_n\|_{H_*^{-1/2}(\Gamma_n^\epsilon)} \\ \le C_1 \epsilon / d_*^\epsilon \max(1, |\omega d_*^\epsilon|^{1/2}) \|\lambda_n\|_{H_*^{-1/2}(\Gamma_n^\epsilon)} $$

The rest follows by the Cauchy-Schwarz inequality and by remarking that

$$ (d_*^\epsilon)^{-1} \max(1, |\omega d_*^\epsilon|^{1/2}) \le (d_*^\epsilon)^{-1} \max(1, (d_*^\epsilon)^{1/2}) \max(1, |\omega|^{1/2}) \\ \le \max(1, (d_*^\epsilon)^{-1}) |\omega|^{1/2} \max(1, (\operatorname{Im} \omega)^{-1/2}). $$

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

**4.4. Estimates on the data in the time domain.** The purpose of this section is to prove (4.5), more precisely

$$ \|\hat{g}_0^\epsilon\|_{1/2} = O(\epsilon^{1/2}), \quad \|\hat{g}_\perp^\epsilon\|_{1/2} = O(\epsilon). $$

For technical reasons, see the discussion in the end of Step 1.2 in Section 4.5.2, we provide these bounds in the time domain directly; they of course can be translated to the frequency domain.

The two main results of this section are Propositions 4.15 and 4.16, which relate the $\epsilon$-dependent norms $\|g_\perp^\epsilon\|_{1/2}$ and $\|g_0^\epsilon\|_{1/2}$ to $\epsilon$-independent norms of the field $u^{\text{inc}}$.

The proof of these propositions relies on two auxiliary lemmas, formulated for the case of a single circle of a fixed radius $r > 0$. We will make use of the same decomposition of the spaces, as the one introduced in Section 4.2.1, however, for the case of the single boundary $C_r$. The projection operators corresponding to $\mathbb{P}_0^\epsilon$ and $\mathbb{P}_\perp^\epsilon$ will be denoted by $P_0$ and $P_\perp$. We start with estimating $P_0\gamma_0 v$.

**LEMMA 4.13.** There exists $C > 0$ s.t., for all $r > 0$, $v \in L^\infty(\mathbb{R}^2)$,

$$ \|P_0\gamma_0 v\|_{L^2(C_r)} \le Cr^{1/2} \|v\|_{L^\infty(\mathbb{R}^2)}. $$

*Proof.* We have that $P_0\gamma_0 v = \frac{1}{2\pi r} \int_{C_r} v d\Gamma$, thus, $\|P_0\gamma_0 v\|_{L^2(C_r)} \le (2\pi r)^{1/2} \|v\|_{L^\infty(C_r)}$. □

Next, let us estimate $P_\perp\gamma_0 v$.

**LEMMA 4.14.** There exists $C > 0$, s.t., for all $r > 0$, $v \in W^{1,\infty}(\mathbb{R}^2)$,

$$ \|P_\perp\gamma_0 v\|_{H^{1/2}(C_r)} \le Cr \max(1, r^{1/2}) \|v\|_{W^{1,\infty}(\mathbb{R}^2)}. $$
*Proof.* Decomposing $\gamma_0 v$ into its Fourier series (2.6) yields

$$ \|P_{\perp} \gamma_0 v\|_{H^{1/2}(C_r)}^2 = r \sum_{m \in \mathbb{Z} \setminus \{0\}} \left(1 + \frac{m^2}{r^2}\right)^{1/2} |v_m|^2 \le 2 \max(1, r) \sum_{m \in \mathbb{Z} \setminus \{0\}} m |v_m|^2. $$

To estimate the above, we use the seminorm (2.7). Defining $\hat{s}_{\theta} = (\cos \theta, \sin \theta)$, we rewrite it in the following way:

$$ \|P_{\perp} \gamma_0 v\|_{H^{1/2}(C_r)}^2 = \int_0^{2\pi} \int_0^{2\pi} \frac{|v(r\hat{s}_{\theta}) - v(r\hat{s}_{\psi})|^2}{|\hat{s}_{\theta} - \hat{s}_{\psi}|^2} d\theta d\psi \lesssim r^2 \sup_{\boldsymbol{x} \in B(0,r)} |\nabla v(\boldsymbol{x})|^2. \quad \square $$

The above two lemmas then allow us to estimate $g_0^\epsilon$ and $g_\perp^\epsilon$.

**PROPOSITION 4.15.** Let $\ell \in \mathbb{N}$, and let for all $t \ge 0$, $\partial_t^\ell u^{\text{inc}}(t) \in L^\infty(\mathbb{R}^2)$. There exists $c_0 > 0$, depending polynomially on $R^*$, s.t.

$$ \|\partial_t^\ell g_0^\epsilon(t)\|_{L^2(\Gamma^\epsilon)} \le c_0 \epsilon^{1/2} N^{1/2} \|\partial_t^\ell u^{\text{inc}}(t)\|_{L^\infty(\mathbb{R}^2)}, \quad t \ge 0. $$

*Proof.* We prove the result for $\ell = 0$; for $\ell > 0$ the proof is a straightforward extension. By Lemma 4.13, we have the bound

$$ \sum_{n \in \mathcal{N}} \|g_0^\epsilon(t)\|_{L^2(\Gamma_n^\epsilon)}^2 \lesssim \sum_{n \in \mathcal{N}} r_n^\epsilon \|u^{\text{inc}}(t)\|_{L^\infty(\mathbb{R}^2)}^2 \le R^* N \epsilon \|u^{\text{inc}}(t)\|_{L^\infty(\mathbb{R}^2)}^2. \quad \square $$

The counterpart of Lemma 4.14, which relates $g_\perp^\epsilon$ to $u^{\text{inc}}$, is given below.

**PROPOSITION 4.16.** Let $\ell \in \mathbb{N}$, and let, for all $t \ge 0$, $\partial_t^\ell u^{\text{inc}}(t) \in W^{1,\infty}(\mathbb{R}^2)$. There exists $c_\perp > 0$, depending polynomially on $R^*$, s.t.

$$ \|\partial_t^\ell g_\perp^\epsilon(t)\|_{H^{1/2}(\Gamma^\epsilon)} \le c_\perp \epsilon N^{1/2} \|\partial_t^\ell u^{\text{inc}}(t)\|_{W^{1,\infty}(\mathbb{R}^2)}, \quad t \ge 0. $$

*Proof.* The proof mimics almost verbatim the proof Proposition 4.15. □

We have now finished with the proofs of auxiliary results, namely, the bounds on the data and the bounds on the operators. In the section that follows we will use these ingredients to prove Proposition 4.3 about convergence of the density, and use this result in order to prove the principal statement of Section 4, i.e. Theorem 4.1.

**4.5. Proof of Proposition 4.3.** The proof of Proposition 4.3 is based on the Plancherel theorem. We split the proof into two parts: first we summarize the bounds in the frequency domain (Section 4.5.1), and next show how they can be translated to the time domain in Section 4.5.2.

**4.5.1. Intermediate frequency domain bounds.** Here we will summarize the frequency-domain bounds corresponding to the quantities in the statement of Proposition 4.3, which will serve as a base for deriving the time-domain bounds.

The constants in the bounds in this and the section that follows depend on $R_*^{-1}$, $R^*$, $d_*^{-1}$, $N$ polynomially.

First of all, the frequency-domain bound for the density $\hat{\mu}_G^\epsilon = (\hat{S}_{00}^\epsilon)^{-1} \hat{g}_0^\epsilon$ follows from Theorem 4.8. With $C_1 > 0$, we have that

$$ (4.30) \qquad \|\hat{\mu}_G^\epsilon\|_{-1/2} \le C_1 \epsilon^{-1} \max(1, \epsilon|\omega|) \frac{|\omega|}{\text{Im} \omega} \|\hat{g}_0^\epsilon\|_{1/2} \le C_1 \epsilon^{-1} \frac{|\omega|}{\text{Im} \omega} \|(1 - i\epsilon\omega)\hat{g}_0^\epsilon\|_{1/2}, $$
where we used $\max(a, b) \le a+b$. We bound $\max(1, \varepsilon|\omega|)\|\hat{g}_0^\varepsilon\|_{1/2}$ by $\|(1 - i\varepsilon\omega)\hat{g}_0^\varepsilon\|_{1/2}$, since the latter expression is more convenient when passing to the time domain. Let

$$ (4.31) \qquad \hat{P}_\omega^\varepsilon = (1 - i\varepsilon\omega). $$

The error for the orthogonal component of the density can be bounded using the expression (4.4), based on Theorems 4.8, 4.12, 4.9:

$$ (4.32) \qquad \|\hat{e}_\perp^\varepsilon\|_{-1/2} \le C_2 \frac{|\omega|}{\operatorname{Im}\omega} \left( \|\hat{P}_\omega^\varepsilon \hat{g}_\perp^\varepsilon\|_{1/2} + \varepsilon^{1/2} |\omega|^{3/2} \max(1, (\operatorname{Im}\omega)^{-3/2}) \|(\hat{P}_\omega^\varepsilon)^2 \hat{g}_0^\varepsilon\|_{1/2} \right) $$

$$ (4.33) \qquad \le C_2 \frac{|\omega|}{\operatorname{Im}\omega} \left( \|\hat{P}_\omega^\varepsilon \hat{g}_\perp^\varepsilon\|_{1/2} + \varepsilon^{1/2} |\omega|^2 \max(1, (\operatorname{Im}\omega)^{-2}) \|(\hat{P}_\omega^\varepsilon)^2 \hat{g}_0^\varepsilon\|_{1/2} \right), $$

where to obtain the last expression we used the bound $|\omega|^{1/2} \max(1, (\operatorname{Im}\omega)^{-1/2}) \le |\omega| \max(1, (\operatorname{Im}\omega)^{-1})$. The formula (4.33), compared to (4.32), is easier to use to pass to the time domain.

Finally, (4.4) can be used to bound the constant component of the density error, combined with (4.32), Theorems 4.8 and 4.12:

$$ (4.34) \qquad \|\hat{e}_0^\varepsilon\|_{-1/2} \le C_3 \varepsilon^{1/2} \frac{|\omega|^2}{(\operatorname{Im}\omega)^2} \left( \begin{aligned} & |\omega| \max(1, (\operatorname{Im}\omega)^{-1}) \|(\hat{P}_\omega^\varepsilon)^2 \hat{g}_\perp^\varepsilon\|_{1/2} + \\ & + \varepsilon^{1/2} |\omega|^2 \max(1, (\operatorname{Im}\omega)^{-2}) \|(\hat{P}_\omega^\varepsilon)^3 \hat{g}_0^\varepsilon\|_{1/2} \end{aligned} \right), $$

where we used the bound $|\omega|^{1/2} \max(1, (\operatorname{Im}\omega)^{-1/2}) \le |\omega| \max(1, (\operatorname{Im}\omega)^{-1})$.

**4.5.2. Translating the bounds into the time domain.** Step 1. *Stability bound for the density.* To obtain a time-domain bound on $\mu_G^\varepsilon$, we apply the Plancherel identity

$$ (4.35) \qquad \int_{-\infty}^{+\infty} e^{-2\eta t} \|f(t)\|^2 dt = \frac{1}{2\pi} \int_{-\infty+i\eta}^{+\infty+i\eta} \|\hat{f}(\omega + i\eta)\|^2 d\omega, \quad \eta > 0, $$

to (4.30). This will make appear the terms of the type $\|\partial_t^k g_0^\varepsilon\|_{H^{1/2}(\Gamma^\varepsilon)}$, which we will further bound by a certain $\varepsilon$-independent norm of $u^{\text{inc}}$, as stated in Section 4.4.

*Step 1.1. Application of the Plancherel identity and the causality argument.* Rewriting the bound (4.30) in the time domain is quite classical, and, e.g. can be found in the monograph [39], or [3], [5], see the references therein. For the convenience of the reader, we will outline the corresponding procedure below.

Let us first of all remark that in the time-domain, $t \mapsto \mu_G^\varepsilon(t)$ vanishes with all its derivatives in the vicinity of 0 as a convolution of a causal operator-valued distribution, whose symbol is given by $(\hat{S}_{00}^\varepsilon)^{-1}$, and $g_0^\varepsilon$, which itself vanishes in the vicinity of 0 (see Section 2.1).

Let $T > 0$. First, we use the Plancherel identity (4.35) and the causality of $\mu_G^\varepsilon$ to obtain the following inequality, valid for any $\eta > 0$,

$$ \int_{\mathbb{R}_+} e^{-2\eta t} \|\mu_G^\varepsilon(t)\|_{-1/2}^2 dt = \frac{1}{2\pi} \int_{\mathbb{R}} \|\hat{\mu}_G^\varepsilon(\omega + i\eta)\|_{-1/2}^2 d\omega \lesssim C_1 \varepsilon^{-2} \eta^{-2} \\ \times \int_{\mathbb{R}} |\omega + i\eta|^2 \|\hat{P}_{\omega+i\eta}^\varepsilon \hat{g}_0^\varepsilon(\omega + i\eta)\|_{1/2}^2 d\omega. $$

We next apply the Plancherel identity again to the right-hand-side of the above inequality, in order to obtain a bound on $\mu_G^\varepsilon$ by $g_0^\varepsilon$, rather than its Fourier-Laplace transform. In particular, recalling that $g_0^\varepsilon$ vanishes in the
vicinity of 0, and choosing $\eta = \frac{1}{T}$ yields the following bound (recall the definition (4.31) of $\hat{P}_{\omega}^{\epsilon} = (1 - i\omega\epsilon)$)

$$ (4.36) \qquad \int_{\mathbb{R}_+} e^{-2t/T} \|\boldsymbol{\mu}_G^{\epsilon}(t)\|_{-1/2}^2 dt \lesssim C_1 \epsilon^{-2} T^2 \int_{\mathbb{R}_+} e^{-2t/T} \|(1 + \epsilon \partial_t) \partial_t \boldsymbol{g}_0^{\epsilon}(t)\|_{1/2}^2 dt. $$

Our goal is now to replace in the above the integrals over $\mathbb{R}_+$ by the integrals over $(0, T)$. This is trivial for the integral on the left; however, to do so for the right-hand side we should make use of the causality argument. More precisely, as argued before, $\boldsymbol{\mu}_G^{\epsilon}$ is a convolution of two causal distributions. This implies that the density $\boldsymbol{\mu}_G^{\epsilon}(t)$ depends only on $\boldsymbol{g}_0^{\epsilon}(\tau)$ for $\tau \le t$. Hence, if $\boldsymbol{g}_{T,0}^{\epsilon}$ is s.t.

$$ \boldsymbol{g}_{T,0}^{\epsilon} = \boldsymbol{g}_0^{\epsilon} \text{ on } (-\infty, T), $$

and

$$ \boldsymbol{\mu}_{G,T}^{\epsilon} = (\mathbb{S}_{00}^{\epsilon}(\partial_t))^{-1} \boldsymbol{g}_{T,0}^{\epsilon}, $$

then we have that $\boldsymbol{\mu}_G^{\epsilon} = \boldsymbol{\mu}_{G,T}^{\epsilon}$ on $(-\infty, T)$. Moreover, if $\boldsymbol{g}_{T,0}^{\epsilon}$ is sufficiently regular, the bound (4.36) holds with $\boldsymbol{\mu}_G^{\epsilon}, \boldsymbol{g}_0^{\epsilon}$ replaced respectively by $\boldsymbol{\mu}_{G,T}^{\epsilon}, \boldsymbol{g}_{T,0}^{\epsilon}$. We choose $\boldsymbol{g}_{T,0}^{\epsilon}$ so that it solves the following Cauchy problem (cf. the proof of Proposition 3.2.2 in [39]), cf. (4.36),

$$ (1 + \epsilon \partial_t) \partial_t \boldsymbol{g}_{T,0}^{\epsilon}(t) = 0 \text{ on } (T, \infty), \quad \partial_t^{\ell} \boldsymbol{g}_{T,0}^{\epsilon}(T) = \partial_t^{\ell} \boldsymbol{g}_0^{\epsilon}(T), \text{ for } \ell = 0, 1. $$

With this choice, (4.36) becomes

$$ (4.37) \qquad \int_{\mathbb{R}_+} e^{-2t/T} \|\boldsymbol{\mu}_{G,T}^{\epsilon}(t)\|_{-1/2}^2 dt \lesssim C_1 \epsilon^{-2} T^2 \int_0^T e^{-2t/T} \|(1 + \epsilon \partial_t) \partial_t \boldsymbol{g}_0^{\epsilon}(t)\|_{1/2}^2 dt. $$

We use $\boldsymbol{\mu}_{G,T}^{\epsilon}(t) = \boldsymbol{\mu}_G^{\epsilon}(t)$ on $(0, T)$ and $1 > e^{-t/T} > e^{-1}$ on $(0, T)$ to obtain further

$$ (4.38) \qquad \int_0^T \|\boldsymbol{\mu}_G^{\epsilon}(t)\|_{-1/2}^2 dt \lesssim C_1 \epsilon^{-2} T^2 \int_0^T \|(1 + \epsilon \partial_t) \partial_t \boldsymbol{g}_0^{\epsilon}(t)\|_{1/2}^2 dt. $$

Thus, we have derived a time-domain bound on $\boldsymbol{\mu}_G^{\epsilon}$ in terms of $\boldsymbol{g}_0^{\epsilon}$.

To get a stronger $L^{\infty}(0, T; H^{-1/2}(\Gamma^{\epsilon}))$-bound on $\boldsymbol{\mu}_G^{\epsilon}$, we use the following bound (see the explanation below on how it can be obtained)

$$ (4.39) \qquad \|\boldsymbol{\mu}_G^{\epsilon}\|_{L^{\infty}(0,T;\nu_0^{\epsilon})} \le T^{1/2} \|\partial_t \boldsymbol{\mu}_G^{\epsilon}\|_{L^2(0,T;\nu_0^{\epsilon})}, $$

and the bound on $\|\partial_t \boldsymbol{\mu}_G^{\epsilon}\|_{L^2(0,T;H^{-1/2}(\Gamma^{\epsilon}))}$ can be obtained like (4.38), with $\boldsymbol{g}_0^{\epsilon}$ replaced by $\partial_t \boldsymbol{g}_0^{\epsilon}$. Finally, to derive (4.39), we use the following identity, valid for sufficiently regular $t \mapsto v(t)$, with $v$ taking values in a Banach space $X$ and s.t. $v(0) = 0$:

$$ (4.40) \qquad v(t) = \int_0^t \partial_{\tau} v(\tau) d\tau \implies \|v(t)\|_X \le t^{1/2} \|\partial_t v\|_{L^2(0,t;X)}, $$

where the last bound follows from the Cauchy-Schwarz inequality.

Altogether, we obtain

$$ (4.41) \qquad \|\boldsymbol{\mu}_G^{\epsilon}\|_{L^{\infty}(0,T;\nu_0^{\epsilon})} \lesssim C_1^{1/2} \epsilon^{-1} T^{3/2} \left( \int_0^T \|(1 + \epsilon \partial_t) \partial_t \boldsymbol{g}_0^{\epsilon}(t)\|_{1/2}^2 dt \right)^{1/2}. $$
Step 1.2. Bounding the right-hand side with respect to $u^{inc}$. Finally, to obtain the bound in the statement of the proposition, we use in (4.41) the estimate of Proposition 4.15. It yields the following bound, where the constant $C'_{\mu 0}$ depends on $R^*$, $R_*^{-1}$, $N$, $d_*^{-1}$ polynomially:

$$ \|\boldsymbol{\mu}_G^\varepsilon\|_{L^\infty(0,T;\mathcal{V}_0^\varepsilon)} \le C'_{\mu 0} \varepsilon^{-1/2} T^{3/2} \|\partial_t^2 u^{inc}\|_{H^1_\varepsilon(0,T;L^\infty(\mathbb{R}^2))}^2. $$

The above yields the desired stability bound in the statement of the proposition.

Let us remark that such a bound would have been difficult to obtain, had we used immediately the frequency-domain counterpart of Proposition 4.15 in the frequency-domain bound (4.30), and next applied the Plancherel theorem. The reason for this is that when applying the Plancherel theorem to $-i\omega\hat{u}^{inc}$, we need to take the non-vanishing initial conditions for $u^{inc}$ when using $\mathcal{F}^{-1}(-i\omega\hat{u}^{inc}) = \partial_t u^{inc}$.

Step 2. Error bounds. The bounds on the errors of different components of the density can be computed like in Step 1, by making use of the inequalities (4.33) and (4.34), the Plancherel theorem, the causality argument and Propositions 4.15, 4.16.

In order to avoid technicalities, we will not specify explicit dependence of the constants in the error estimates on $T$, since it is clear that the resulting bounds will depend on it polynomially. In what follows, the constants $C_{i,G,\perp}$, $C_{i,G,0}$, $C_{i,G}$, $i \in \mathbb{N}$, depend on polynomially on $T$, $N$, $d_*^{-1}$, $R^*$, $R_*^{-1}$.

Step 2.1. An error bound for the orthogonal component of the density. Our next goal is to estimate the error $\boldsymbol{e}_\perp^\varepsilon$, starting with the bound (4.33). Applying the same arguments as in Step 1 yields the following bound:

$$ (4.42) \qquad \|\boldsymbol{e}_\perp^\varepsilon\|_{L^\infty(0,T;H^{-1/2}(\Gamma^\varepsilon))} \le C_{1,G,\perp} (\|\partial_t^2 \boldsymbol{g}_\perp^\varepsilon\|_{H^1_\varepsilon(0,T;H^{1/2}(\Gamma^\varepsilon))} + \varepsilon^{1/2} \|\partial_t^3 \boldsymbol{g}_0^\varepsilon\|_{H^2_\varepsilon(0,T;\mathcal{V}_0^\varepsilon)}). $$

Step 2.2. An error bound for the constant component of the density. We start with (4.34), which in the time domain yields:

$$ (4.43) \qquad \|\boldsymbol{e}_0^\varepsilon\|_{L^\infty(0,T;H^{-1/2}(\Gamma^\varepsilon))} \le \varepsilon^{1/2} \\ \qquad \times C_{1,G,0} (\|\partial_t^4 \boldsymbol{g}_\perp^\varepsilon\|_{H^2_\varepsilon(0,T;H^{1/2}(\Gamma^\varepsilon))} + \varepsilon^{1/2} \|\partial_t^5 \boldsymbol{g}_0^\varepsilon\|_{H^3_\varepsilon(0,T;H^{1/2}(\Gamma^\varepsilon))}), $$

Step 2.3. Obtaining final error bounds. It remains to sum up (4.42), (4.43):

$$ \|\boldsymbol{e}^\varepsilon\|_{L^\infty(0,T;H^{-1/2}(\Gamma^\varepsilon))} \lesssim C_{1,G} (\|\partial_t^4 \boldsymbol{g}_\perp^\varepsilon\|_{H^2_\varepsilon(0,T;H^{1/2}(\Gamma^\varepsilon))} + \varepsilon^{1/2} \|\partial_t^5 \boldsymbol{g}_0^\varepsilon\|_{H^3_\varepsilon(0,T;H^{1/2}(\Gamma^\varepsilon))}), $$

where we used the bound (4.40). Applying Propositions 4.15 and 4.16 yields the final bound in the statement of the proposition:

$$ \|\boldsymbol{e}^\varepsilon\|_{L^\infty(0,T;H^{-1/2}(\Gamma^\varepsilon))} \lesssim C_{2,G} \varepsilon (\|\partial_t^4 u^{inc}\|_{H^2_\varepsilon(0,T;W^{1,\infty}(\mathbb{R}^2))} + \|\partial_t^5 u^{inc}\|_{H^3_\varepsilon(0,T;L^\infty(\mathbb{R}^2))}). $$

Let us remark that with the bounds of Propositions 4.15 and 4.16, we have

$$ (4.44) \qquad \|\hat{\boldsymbol{e}}_\perp^\varepsilon\|_{-1/2} = O(\varepsilon) \quad \text{and} \quad \|\hat{\boldsymbol{e}}_0^\varepsilon\|_{-1/2} = O(\varepsilon^{3/2}). $$

This fact is not reflected in the statement of the proposition, but will be of use later.

**4.6. Proof of Theorem 4.1.** Let $\delta_\varepsilon > 0$, and let $\boldsymbol{x} \in \Omega^{\varepsilon+\delta_\varepsilon,c}$. Remark that $\mathrm{dist}(\boldsymbol{x}, \Omega^{\varepsilon,c}) \ge \delta_\varepsilon$. First, we will obtain an explicit in frequency bound for the error

$$ \hat{v}^\varepsilon := \hat{u}^\varepsilon(\omega, \boldsymbol{x}) - \hat{u}_G^\varepsilon(\omega, \boldsymbol{x}), $$
and next translate it into the time domain like in the proof of Proposition 4.3.

We use an explicit representation of $\hat{v}^{\epsilon}$ and split

$$ \hat{v}^{\epsilon} = \hat{v}_{0}^{\epsilon} + \hat{v}_{\perp}^{\epsilon}, $$

$$ \hat{v}_{0}^{\epsilon} := \int_{\Gamma^{\epsilon}} G_{\omega}(\mathbf{x} - \mathbf{y}) \hat{e}_{0}^{\epsilon}(\mathbf{y}) d\Gamma_{\mathbf{y}}, \quad \hat{v}_{\perp}^{\epsilon} := \int_{\Gamma^{\epsilon}} G_{\omega}(\mathbf{x} - \mathbf{y}) \hat{e}_{\perp}^{\epsilon}(\mathbf{y}) d\Gamma_{\mathbf{y}}. $$

We would like to show in particular that the far-field error behaves as $|\hat{v}^{\epsilon}| = O(\epsilon^2)$, while the density error loses an order of $\epsilon$: $\|\hat{e}^{\epsilon}\|_{H^{-1/2}(\Gamma^{\epsilon})} = O(\epsilon)$, cf. Proposition 4.3. In order to obtain such a super-convergence result, we will treat the terms $\hat{v}_{0}^{\epsilon}$ and $\hat{v}_{\perp}^{\epsilon}$ separately, in order to make use of the different scalings of the error components (4.44). We will explain the interest in this approach afterwards.

In what follows, by $C_{i,G}$, $i \in \mathbb{N}$, we will denote a constant depending polynomially on $R^*$, $R_*^{-1}$, $N$, $d_*^{-1}$, and, where applicable, on $T$.

We start with the error generated by the constant component of the density error:

$$ \begin{align} |\hat{v}_0^\epsilon|^2 &\le \left| \sum_{n \in \mathcal{N}} \int_{\Gamma_n^\epsilon} G_\omega(\mathbf{x}-\mathbf{y}) \hat{e}_{n,0}^\epsilon(\mathbf{y}) d\Gamma_\mathbf{y} \right|^2 \tag{4.45} \\ &\le N \sup_n \|G_\omega(\mathbf{x}-\cdot)\|_{L^2(\Gamma_n^\epsilon)}^2 \|\hat{e}_0^\epsilon\|_{L^2(\Gamma^\epsilon)}^2. \nonumber \end{align} $$

The desired bound for $\hat{v}_0^\epsilon$ then follows from Lemma 4.11 (used in a simplified form):

$$ \begin{align} |\hat{v}_0^\epsilon| &\le C_{1,G} \epsilon^{1/2} \max(1, \log \frac{1}{\delta_\epsilon \operatorname{Im} \omega}) \|\hat{e}_0^\epsilon\|_{L^2(\Gamma^\epsilon)} \tag{4.46} \\ &\lesssim C_{1,G} \epsilon^{1/2} \max(1, \log \frac{1}{\delta_\epsilon}) \max(1, \log \frac{1}{\operatorname{Im} \omega}) \|\hat{e}_0^\epsilon\|_{L^2(\Gamma^\epsilon)}, \nonumber \end{align} $$

where to obtain the latter inequality we used $\max(1, a+b) \le \max(1, a) + \max(1, b) \le \max(1, a)(1+\max(1, b)) \le 2\max(1, a)\max(1, b)$.

In a similar manner, we can bound the field error generated by the error in the orthogonal component of the density:

$$ (4.47) \qquad |\hat{v}_{\perp}^{\epsilon}|^2 \le N \sup_n \|G_{\omega}(\mathbf{x} - \cdot)\|_{H^{1/2}(\Gamma_n^{\epsilon})}^2 \|\hat{e}_{\perp}^{\epsilon}\|_{H^{-1/2}(\Gamma^{\epsilon})}^2. $$

Using the result of Lemma 4.11,

$$ \begin{align} |\hat{v}_{\perp}^{\epsilon}| &\le C_{2,G} \epsilon \delta_{\epsilon}^{-1} \max(1, |\omega \delta_{\epsilon}|^{1/2}) \|\hat{e}_{\perp}^{\epsilon}\|_{-1/2} \tag{4.48} \\ &\le C_{2,G} \epsilon \delta_{\epsilon}^{-1} \max(1, \delta_{\epsilon}^{1/2}) |\omega| \max(1, (\operatorname{Im} \omega)^{-1}) \|\hat{e}_{\perp}^{\epsilon}\|_{-1/2}, \nonumber \end{align} $$

where to get the last bound we used in particular $\max(1, |\omega|^{1/2}) < \max(1, |\omega|) < |\omega| \max(1, (\operatorname{Im} \omega)^{-1})$.

At this point we can see the interest in splitting $\hat{v}^{\epsilon}$ into two components. When $\delta_{\epsilon} = \text{const}$, both $|\hat{v}_{0}^{\epsilon}|$ and $|\hat{v}_{\perp}^{\epsilon}|$ have the same scaling with respect to $\epsilon$: they behave as $O(\epsilon^2)$, cf. (4.44). Had we used from the beginning the same argument as we used to obtain (4.48), we would have the bound

$$ |\hat{v}^{\epsilon}| \lesssim \sup_{n} \|G_{\omega}(\mathbf{x}_{\epsilon} - \cdot)\|_{H^{1/2}(\Gamma_{n}^{\epsilon})} \|\hat{e}^{\epsilon}\|_{H^{-1/2}(\Gamma^{\epsilon})}, $$

which, with Lemma 4.11 and Proposition 4.3 would yield, for $\delta_{\epsilon} = \text{const}$, $|\hat{v}^{\epsilon}| = O(\epsilon^{3/2})$.
With bounds (4.46) and (4.48), and remarking that both $\max(1, \log \frac{1}{\delta_{\varepsilon}})$ and $\delta_{\varepsilon}^{-1} \max(1, \delta_{\varepsilon}^{1/2})$ are bounded by $\max(1, \delta_{\varepsilon}^{-1})$, we arrive at the following expression:

$$|\hat{\nu}^{\varepsilon}| \le C_{3,G} \max(1, \delta_{\varepsilon}^{-1}) \max(1, (\mathrm{Im} \, \omega)^{-1}) \left( \varepsilon |\omega| \|\hat{e}_{\perp}^{\varepsilon}\|_{-1/2} + \varepsilon^{1/2} \|\hat{e}_{0}^{\varepsilon}\|_{-1/2} \right).$$

It remains to replace in the above bound $\varepsilon|\omega|\|\hat{e}_{\perp}^{\varepsilon}\|_{-1/2} + \varepsilon^{1/2}\|\hat{e}_{0}^{\varepsilon}\|_{-1/2}$ by their bounds (4.33), (4.34), and next proceed like in the proof of Proposition 4.3.

**4.7. Discussion of the obtained results and the techniques used in the proofs.** As one can notice already in the proof of Proposition 4.3, the powers of $\omega$ in the frequency-domain bounds are translated into the derivatives in the time domain, while powers of $(\mathrm{Im} \, \omega)^{-1}$ yield the constants growing polynomially in time. As discussed in Remark 5, obtaining optimal error bounds in $\omega$ and $\mathrm{Im} \, \omega$ is necessary to derive optimal convergence results, and to understand the interplay between the asymptotic parameter $\varepsilon$ and the wavelength. The approach that we use here most likely yields non-optimal results. First of all, we bound norms of products of operators by products of their norms, cf. (4.4), see the related discussions on p. 87 of [39]. Some information is also lost when passing from the Fourier-Laplace domain into the time domain, see the bottom of p. 44 in [39]. Moreover, we obtain bounds without distinguishing the cases of the frequencies belonging to a bounded range and asymptotic $|\omega| \to +\infty$, cf. the proofs of Theorems 4.8, 4.9.

**5. Extension of the results: a cluster of particles.** The goal of this section is to show how the obtained convergence bounds can be adapted to the case of sound-soft scattering problem in the case when the centres of the particles and distances between them are no longer fixed with $\varepsilon$. We will assume that $N = O(\varepsilon^{-\gamma})$, and $d_{*}^{\varepsilon} = O(\varepsilon^{\alpha})$, with $\gamma, \alpha \ge 0$. Let us additionally assume that $\alpha \le 1$. We will limit our discussion to the case when the particles fill a fixed volume, which implies that $N \lesssim \min(\varepsilon^{-2}, \varepsilon^{-2\alpha})$, thus

$$(5.1) \qquad \gamma \le 2\alpha.$$

Our goal is to extend the (far-field) result of Theorem 4.1 to this case. In other words, we would like to bound the following quantity:

$$|\nu^{\varepsilon}(t)| = |u_{G}^{\varepsilon}(t, \boldsymbol{x}) - u^{\varepsilon}(t, \boldsymbol{x})|, \quad \boldsymbol{x} \in \Omega^{c}, \quad \mathrm{dist}(\boldsymbol{x}, \Omega^{c, c}) = \delta > 0.$$

It is straightforward to verify that the only assumption used in our proofs was an assumption on the circles not touching each other. Therefore, as seen from the proof of Proposition 4.3 and Theorem 4.1, it suffices to (carefully) keep track of the dependence of the error bounds in the frequency domain on $N$, $d_{*}^{\varepsilon}$ and $\varepsilon$. Like before, we will rely on the inequalities (4.4). For convenience, let us rewrite the bounds of Theorems 4.8, 4.9, 4.12, as well as of Propositions 4.15, 4.16. First, according to (4.8),

$$C_0(\varepsilon, \frac{d_*^\varepsilon}{\varepsilon}) = 1 + \frac{\varepsilon}{d_*^\varepsilon} = O(1), \quad \varepsilon \to 0.$$

We then have the following bounds for the operator norms:

$$\| (\hat{S}_{00}^{\varepsilon})^{-1} \| = O(\varepsilon^{-1}), \quad \| \mathbb{P}_{\perp}^{\varepsilon, *} (\hat{S}^{\varepsilon})^{-1} \mathbb{P}_{\perp}^{\varepsilon} \| = O(1), \quad \| \hat{S}_{0\perp}^{\varepsilon} \| = \| \hat{S}_{\perp 0}^{\varepsilon} \| = O(\varepsilon^{3/2-\gamma-\alpha}).$$

Moreover,

$$\| \hat{g}_0^\varepsilon \|_{1/2} = O(\varepsilon^{1/2-\gamma/2}), \quad \| \hat{g}_\perp^\varepsilon \|_{1/2} = O(\varepsilon^{1-\gamma/2}).$$
It remains to rewrite the inequalities (4.4):

$$ \|\hat{e}_{\perp}^{\epsilon}\|_{-1/2} = O(\epsilon^{1-\frac{\gamma}{2}}) + O(\epsilon^{1-3\gamma/2-\alpha}) = O(\epsilon^{1-3\gamma/2-\alpha}). $$

The constant component of the density error, cf. (4.4), satisfies:

$$ \|\hat{e}_0^\epsilon\| = O(\epsilon^{3/2-5\gamma/2-2\alpha}). $$

Finally, to obtain the bound for the field, we proceed like in the proof of Theorem 4.1. In particular, by (4.45), (4.46), as well as (4.47), (4.48) it follows that

$$ |\hat{v}_0^\epsilon| = O(\epsilon^{-\gamma/2+1/2} \|\hat{e}_0^\epsilon\|_{-1/2}) = O(\epsilon^{2-3\gamma-2\alpha}) \quad \text{and} $$

$$ |\hat{v}_{\perp}^{\epsilon}| = O(\epsilon^{-\gamma/2+1} \|\hat{e}_{\perp}^{\epsilon}\|_{-1/2}) = O(\epsilon^{2-2\gamma-\alpha}). $$

Therefore,

$$ |\hat{v}^{\epsilon}| = O(\epsilon^{2-3\gamma-2\alpha}). $$

Evidently, we cannot expect convergence in the case when $\alpha = 1$, independently of the number of particles. For $\alpha < 1$, in the particular case $\gamma = 2\alpha$, we can expect the convergence rate of $O(\epsilon^{2(1-4\alpha)})$ for $\alpha < 1/4$.

**6. Numerical experiments.** In all the experiments we use the trapezoid convolution quadrature method for the time semi-discretization. The reference solutions are computed with the use of the spatial Galerkin BEM with $2N_s + 1$ spectral Galerkin basis functions $\{e^{im\theta}\}_{m=-N_s}^{N_s}$ on each of the obstacles.

**6.1. Numerical validation of Theorem 4.1 and Proposition 4.3.** We study a configuration of $N = 20$ particles with centers located on the boundary of a unit circle, see Figure 6.1, left. We start with the $R_i = R = 0.1$ (in this case $d_* \approx 0.11$). The incident field $u^{\text{inc}}(t, \mathbf{x}) = e^{-100(t - \mathbf{d} \cdot \mathbf{x} - 2)^2}$, $\mathbf{d} = (0, 1)$, is approximately band-limited, with the smallest wavelength $\lambda_{\text{min}} \approx 0.12$. The simulations are performed on the time interval $(0, T)$, $T = 8$. We change $\epsilon$ and compare $\mu_G^\epsilon(t)$ and $u_G^\epsilon(t, \mathbf{x}_0)$ in the centre of the circle $\mathbf{x}_0 = 0$ to the reference density $\mu^\epsilon(t, \mathbf{x})$ and solution $u^\epsilon(t, \mathbf{x}_0)$ by computing

$$ (6.1) \qquad e_{\mu}^{\epsilon} = \|e^{\epsilon}\|_{L^{\infty}(0,T;H^{-1/2}(\Gamma^{\epsilon}))}, \quad e_{u}^{\epsilon} = \|u_{G}^{\epsilon}(., \mathbf{x}_{0}) - u^{\epsilon}(., \mathbf{x}_{0})\|_{L^{\infty}(0,T)}. $$

The dependence of the solutions $u^\epsilon(t, \mathbf{x}_0)$ and $u_G^\epsilon(t, \mathbf{x}_0)$ on $t$ is depicted in Figure 6.1, right. The dependence of $e_\mu^\epsilon$ and $e_u^\epsilon$ on $\epsilon$ is shown in Figure 6.2, left. The numerical experiments support the claims of Proposition 4.3 and Theorem 4.1: as expected, $e_\mu^\epsilon = O(\epsilon)$ and $e_u^\epsilon = O(\epsilon^2)$.

**6.2. Comparison with the Foldy-Lax model (FL1-FL2).** The goal of this section is to compare the Foldy-Lax model (FL1-FL2) to the new Galerkin Foldy-Lax model. In the case when the Foldy-Lax model converges, we expect it to have the same convergence order as the Galerkin Foldy-Lax model, see Section 3.4.3 for the related discussion, as well as Proposition 4.3. The reason for this is that in the frequency domain, for a fixed frequency, the difference between the coefficients of the two models is $O(\epsilon^2)$, which is of convergence order of the solution.

We present two experiments: the first one in which the Foldy-Lax model (FL1-FL2) has a performance comparable to the Galerkin Foldy-Lax model, and the second one where the Galerkin Foldy-Lax model outperforms the Foldy-Lax model.

The first experiment is performed with the same data as in the previous section. In Figure 6.2, right, we depict $e_u^\epsilon$ defined for the Galerkin Foldy-Lax model in (6.1) and the analogous quantity $e_{FL,u}^\epsilon$ for the Foldy-Lax model (FL1-FL2). We were not able to make the first experiment with $\epsilon R = 0.1$ for the original Foldy-Lax model.
<figure><img src="image_8.png" /><figcaption>FIG. 6.1. Left: the geometric configuration for the numerical experiment of Section 6 for $\epsilon = 1$. Right: the dependence of the solutions $u^\epsilon(t, \mathbf{x}_0)$ in $\mathbf{x}_0 = 0$ on time $t$ obtained with the help of the Galerkin Foldy-Lax method and the high-order BEM for $\epsilon R = 0.05$ (top) and $\epsilon R = 0.01$ (bottom).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6.2. Left: an illustration to the experiment of Section 6.1. Convergence of the errors $e_u^\epsilon$ and $e_\mu^\epsilon$. Right: an illustration to the experiment of Section 6.2. Convergence of the error $e_{FL,u}^\epsilon$.</figcaption></figure>

model, since we observed an exponential blow-up of the solution on the given time interval. Nonetheless, for small $\epsilon$, we see that the errors of the Galerkin Foldy-Lax model and the model (FL1-FL2) almost coincide.

The second experiment is performed on longer times. We take $N = 16$ particles; the centers of the 15 particles of radius $R_i = 0.01$, $1 \le i \le 15$, are equidistant, and are located on the boundary of a circle of radius
<figure><img src="image_10.png" /><figcaption>FIG. 6.3. An illustration to the numerical experiment of Section 6.2. Left: The dependence of the solutions $u^1(t, \mathbf{x}_0)$ and $u_G^1(t, \mathbf{x}_0)$ on $t$. Right: the dependence of the condition number of the Foldy-Lax matrix (FL1) on $\omega$ plotted in the axis $(\text{Re}\omega, \text{Im}\omega)$.</figcaption></figure>

0.1 centered in the origin, and the remaining particle of radius $R_{16} = 0.02$ is again centered in the origin. We set

$$u^{\text{inc}}(t, \mathbf{x}) = -\sin(8(t - s(\mathbf{x})))\mathcal{H}(t - s(\mathbf{x})), \quad s(\mathbf{x}) = x_2 + 2, \quad \mathcal{H}(t) = \frac{1}{1 + e^{-20t}}.$$

Let us remark that the initial data for the Cauchy problem satisfied by the incident field (2.1) is no longer compactly supported, but the analysis of the article can be extended to this case as well.

We perform the simulations on the time interval $(0, T)$ with $T = 20$, and measure the field in the point $\mathbf{x}_0 = (0.07, 0)$. We show the reference solution versus the Galerkin Foldy-Lax approximation in Figure 6.3, left. The relative $L^\infty(0, T)$-error between the two solutions does not exceed 3%. At the same time, the solution obtained with the Foldy-Lax model (FL1-FL2) explodes, at least for the chosen time step and the simulation time.<fnref n="1" /> Numerically, we observe that the matrix in the right-hand side of (2.9) is not invertible for $\omega \approx 39.47 + 0.03i$, see Figure 6.3, right, and this probably accounts for the exponential growth of the computed solution.

Remark that as the source is harmonic in time, we expect the limiting amplitude principle [14] to hold for the exact solution; this is the case for the approximation of the solution computed with the help of the Galerkin Foldy-Lax model, as shown in Figure 6.3, left.

## 7. Conclusions and open questions.
In this work we have constructed an asymptotic model for scattering by small particles as a Galerkin discretization of the single-layer boundary integral formulation. This procedure yields a model a priori stable in the time-domain. For a particular case of scattering by circles, we have shown the second order convergence of the method, confirmed by the numerical experiments. Nonetheless, many questions still remain open, and we list them below:

1. optimality of the obtained estimates in regularity requirements;

2. construction of a (higher order) Galerkin Foldy-Lax model for particles of arbitrary shapes;

3. efficient numerical methods for such models, and comparison of their performance to the finite element method based simulations;

4. advantages/disadvantages in using other integral formulations (e.g. time-domain CFIE or direct integral formulation);

5. performance of Galerkin Foldy-Lax models purely in the frequency domain (in particular for the problems without absorption).

<footnote n="1">¹Remarkably, one still may have stability for larger time steps, see Section 3.2 for a related discussion.</footnote>
We plan to address some of these questions in future works, as well as extend the analysis to the 3D electro-magnetic case.

**Acknowledgements.** I am grateful to Maxence Cassier (CNRS, Institut Fresnel), Patrick Joly (INRIA, IP Paris) for fruitful discussions. In particular, I thank Maxence for explaining the approach of [9], and Patrick for suggesting the name 'Galerkin Foldy-Lax'.

I am also thankful to the reviewers for their remarks that allowed to improve the quality of the presentation significantly.

**Appendix A. Proof of Proposition 2.1.** The outgoing solution to (2.5) is given by (we write $\hat{g}_n^\epsilon$ for $\hat{g}_{1,n}^\epsilon$) the series of Hankel functions:

$$ \hat{u}^{\epsilon}(\omega, \rho\hat{s}_{\theta}) = \sum_{n \in \mathbb{Z}} \frac{H_{n}^{(1)}(\omega\rho)}{H_{n}^{(1)}(\omega\epsilon)} \hat{g}_{n}^{\epsilon} e^{in\theta}, \quad \hat{s}_{\theta} = (\cos\theta, \sin\theta). $$

We further split $\hat{u}^\epsilon$ into two terms that have different scalings with respect to $\epsilon$:

$$ \hat{u}^{\epsilon} = \hat{u}_{0}^{\epsilon} + \hat{u}_{\perp}^{\epsilon}, \quad \hat{u}_{0}^{\epsilon}(\omega, \rho\hat{s}_{\theta}) = \frac{H_{0}^{(1)}(\omega\rho)}{H_{0}^{(1)}(\omega\epsilon)} \hat{g}_{0}^{\epsilon}, \quad \hat{u}_{\perp}^{\epsilon}(\omega, \rho\hat{s}_{\theta}) = \sum_{n \in \mathbb{Z} \setminus \{0\}} \frac{H_{n}^{(1)}(\omega\rho)}{H_{n}^{(1)}(\omega\epsilon)} \hat{g}_{n}^{\epsilon} e^{in\theta}. $$

Asymptotic behaviour of Hankel functions [13, Chapter 10.8] yields

$$ (A.1) \qquad \hat{u}_0^\epsilon(\omega, \rho\hat{s}_\theta) = -\frac{i\pi}{2} H_0^{(1)}(\omega\rho) \hat{g}_0^\epsilon \log^{-1} \epsilon + o(\log^{-1} \epsilon), \quad \text{as } \epsilon \to 0, $$

cf. the statement of Proposition 2.1. The term $\hat{u}_{\perp}^{\epsilon}(\omega, \rho\hat{s}_{\theta})$ is bounded with the Cauchy-Schwarz inequality and Lemma D.1:

$$ (A.2) \qquad |\hat{u}_{\perp}^{\epsilon}(\omega, \rho\hat{s}_{\theta})|^2 \le \sum_{n \in \mathbb{Z} \setminus \{0\}} \left| \frac{H_{n}^{(1)}(\omega\rho)}{H_{n}^{(1)}(\omega\epsilon)} \right|^2 \sum_{n \in \mathbb{Z} \setminus \{0\}} |\hat{g}_{n}^{\epsilon}|^2 \le C\epsilon^2 \sum_{n \in \mathbb{Z} \setminus \{0\}} |\hat{g}_{n}^{\epsilon}|^2. $$

Let us now estimate the quantities related to the data. For sufficiently regular $\hat{u}^{\text{inc}}$,

$$ \begin{align*} \hat{g}_0^\epsilon &= -\frac{1}{2\pi\epsilon} \int_{\Gamma^\epsilon} \hat{u}^{\text{inc}}(\omega, \mathbf{x}) d\Gamma_\mathbf{x} = -\hat{u}^{\text{inc}}(\omega, 0) + O(\epsilon), \\ \sum_{n \in \mathbb{Z} \setminus \{0\}} |\hat{g}_n^\epsilon|^2 &= (2\pi\epsilon)^{-1} \|\hat{g}^\epsilon - \hat{g}_0^\epsilon\|_{L^2(\Gamma^\epsilon)}^2 = (2\pi\epsilon)^{-1} \inf_{c \in \mathcal{C}} \|\hat{g}^\epsilon - c\|_{L^2(\Gamma^\epsilon)}^2 \\ &\lesssim \epsilon^{-1} \int_{\Gamma^\epsilon} |\hat{u}^{\text{inc}}(\omega, \mathbf{y}) - \hat{u}^{\text{inc}}(\omega, 0)|^2 d\Gamma_\mathbf{y} = O(\epsilon^2). \end{align*} $$

Inserting the above expressions into (A.2), (A.1) yields the desired result in the statement of the proposition.

**Appendix B. Proof of Proposition 2.2.** Proof of Proposition 2.2 relies on several auxiliary lemmas. First of all, let us examine (FL1), (FL2). Evidently, the temporal behaviour of $u_{FL}^1$ is defined by the behaviour of $\lambda_{FL}^1$. In order to show that $t \mapsto \|\lambda_{FL}^1(t)\|$ admits an exponential growth, we will prove that its Fourier-Laplace transform has poles with a strictly positive imaginary part. The poles of $\hat{\lambda}_{FL}^1$ are closely related to
the points $\omega$ where the matrix in the left-hand side of (2.9) is not invertible. Let us denote this matrix for the geometry described in Proposition 2.2 by $\Lambda(\omega)$:

$$ \Lambda_{nk}(\omega) = \delta_{nk} + (1 - \delta_{nk})P(\omega), \quad P(\omega) = \frac{H_0^{(1)}(\omega c)}{H_0^{(1)}(\omega r)}, \quad k, n = 1, \dots, N. $$

We then have the following result.

**LEMMA B.1.** Under assumptions of Proposition 2.2, with $\eta = cr^{-1} < 4$,

(i) The matrix $\Lambda = \Lambda(\omega)$ is not invertible for a countable set of frequencies $\omega_n \in \mathbb{C}^+$, $n \in \mathbb{Z}$.

(ii) These frequencies satisfy

$$ (B.1) \qquad \lim_{n \to \pm\infty} \operatorname{Im} \omega_n = \frac{1}{r(\eta - 1)} \log \frac{2}{\sqrt{\eta}} > 0. $$

(iii) The respective eigenvalues are simple, moreover $\operatorname{Ker} \Lambda(\omega_n) = \operatorname{span}\{(1, 1, 1)\}$.

We present the complete proof later, while here we will briefly outline how (i) and (ii) are proven. One can show that the points $\omega \in \mathbb{C}^+$ where the matrix $\Lambda$ is not invertible are in particular the roots $\omega \in \mathbb{C}^+$ of the equation $P(\omega) = -1/2$. From the known asymptotics of the Hankel functions it follows that

$$ P(\omega) = \eta^{-\frac{1}{2}} e^{i(\eta-1)\omega r} + o(1), \quad |\omega| \to +\infty. $$

The inverse function theorem can then be used to prove that the roots of $P(\omega) = -1/2$ are close to the roots of $\eta^{-\frac{1}{2}}e^{i(\eta-1)\omega r} = -1/2$, which, in turn, are given by

$$ (B.2) \qquad \omega_{\infty,n} = \frac{1}{r(\eta-1)} \left( (2n+1)\pi + i \log \frac{2}{\sqrt{\eta}} \right), \quad n \in \mathbb{Z}. $$

Evidently, $\operatorname{Im} \omega_{\infty,n} > 0$ for $\eta < 4$.

To illustrate the statement of Lemma B.1, we plot $E(\omega) = |P(\omega) + 1/2|$ in Figure B.1 (see the proof of the lemma). We compare the numerically computed roots of $E(\omega)$ (zeros of $\det \Lambda(\omega)$ in $\mathbb{C}^+$) to the values $\omega_{\infty,n}$ cf. (B.2). As seen in Figure B.1, for the chosen parameters, the values $\omega_{\infty,n}$ provide a good approximation to the roots of $E(\omega)$.

<figure><img src="image_11.png" /><figcaption>FIG. B.1. $|P(\omega) + 1/2|$ in the coordinates $(\operatorname{Re}\omega, \operatorname{Im}\omega)$, $\eta = 2.1$ and $r = 1$. The values $\omega_{\infty,n}$ are shown in white.</figcaption></figure>

To prove Lemma B.1, we will rely on the following result.

**LEMMA B.2.** [12, Chapter XII, Lemma 1.2] Let $f$ be analytic on $B(0, R)$, $f(0) = 0$, $|f'(0)| = \mu > 0$, and $|f(z)| \le M$ for all $z \in B(0, R)$. Then

$$ f(B(0; R)) \supset B \left( 0; \frac{R^2 \mu^2}{6M} \right). $$
We will also need an asymptotic behaviour of Hankel functions [37, pp. 266-267]. For all $z \in \mathbb{C}^+$ (where for $x < 0$ we write $H_m^{(1)}(x) = \lim_{\nu \to 0+} H_m^{(1)}(x + i\nu)$):

$$ (B.3) \qquad H_0^{(1)}(z) = \left(\frac{2}{\pi z}\right)^{\frac{1}{2}} e^{i(z-\frac{\pi}{4})} (1 + O(|z|^{-1})), \quad |z| \to +\infty, $$

$$ (B.4) \qquad H_1^{(1)}(z) = \left(\frac{2}{\pi z}\right)^{\frac{1}{2}} e^{i(z-\frac{3\pi}{4})} (1 + O(|z|^{-1})), \quad |z| \to +\infty. $$

The error constants in the above expansions can be bounded uniformly in $Arg z \in [0, \pi]$.

*Proof of Lemma B.1.* The matrix $\Lambda$ with $\det \Lambda = 1+2P^3-3P^2$ is not invertible for $\omega$ s.t. $P(\omega) \in \{-1/2, 1\}$. Let us look for $\omega \in \mathbb{C}^+$, s.t.

$$ (B.5) \qquad P(\omega) = -1/2. $$

Let us remark that in general it is possible to show that for $\omega \in \mathbb{C}^+$, $|P(\omega)| < 1$, and therefore, $\Lambda(\omega)$ is not invertible in $\omega \in \mathbb{C}^+$ if and only if (B.5) holds true.

*Proof of (i) and (ii).* The set of the roots of (B.5) is at most countable, since $\omega \mapsto P(\omega)$ is analytic in $\mathbb{C}^+$. To show that there are infinitely many roots, we will rely on

* an explicit asymptotic behaviour of $P(\omega)$, which will allow us to obtain an ansatz for the location of the solutions of (B.5) in $\mathbb{C}^+$;

* Lemma B.2 which will that some solutions of (B.5) are close to this ansatz.

*Step 1. Ansatz for the roots.* By (B.3), for all $|\omega r|$ sufficiently large, we have that

$$ (B.6) \qquad P(\omega) = \frac{H_0^{(1)}(\omega\eta r)}{H_0^{(1)}(\omega r)} = \eta^{-\frac{1}{2}} e^{i(\eta-1)\omega r} + O(|\omega|^{-1}). $$

The solutions of

$$ (B.7) \qquad \eta^{-1/2} e^{i(\eta-1)\omega r} = -1/2 $$

are given by

$$ (B.8) \qquad \omega_{\infty,n} = \frac{1}{r(\eta-1)} \left( (2n+1)\pi + i \log \frac{2}{\sqrt{\eta}} \right), \quad n \in \mathbb{Z}. $$

Because $1 < \sqrt{\eta} < 2$ by the assumption of the proposition, all these roots satisfy $\omega_{\infty,n} \in \mathbb{C}^+$, with $\mathrm{Im}\,\omega_{\infty,n}$ independent of $n$. As $n \to \infty$, we also have that $|\omega_{\infty,n}| = O(n)$. Moreover, $|\omega_{\infty,n+1} - \omega_{\infty,n}| = \frac{2\pi}{r(\eta-1)}$. Without loss of generality, in the proof we will consider the case $n > 0$.

*Step 2. Roots of (B.5).* Based on the above, let us show that countably many solutions of (B.5) lie in $\mathbb{C}^+$. Let us consider the following function:

$$ p_n(\omega) := P(\omega) - P(\omega_{\infty,n}). $$

Remark that, by (B.6), we have in particular

$$ (B.9) \qquad |P(\omega_{\infty,n}) + 1/2| \le C_P n^{-1}, \quad C_P > 0. $$
We will show that for all $n$ sufficiently large, there exists $\delta_n > 0$, s.t. $B(\omega_{\infty,n}, \delta_n) \subset \mathbb{C}^+$ and, with some $c_p > 0$,

$$ (B.10) \qquad p_n(B(\omega_{\infty,n}, \delta_n)) \supset B(0, c_p n^{-1/2}). $$

For the function

$$ P(\omega) + 1/2 = p_n(\omega) + z_n, \quad z_n = P(\omega_{\infty,n}) + 1/2, $$

this will imply that for all $z \in B(z_n, c_p n^{-1/2})$, there exists $\omega_z \in B(\omega_{\infty,n}, \delta_n)$, s.t.

$$ P(\omega_z) + 1/2 = z. $$

In particular, by the bound (B.9), $0 \in B(z_n, c_p n^{-1/2})$, and hence the desired result. We will show a stronger result, namely, we will prove that in the above one can choose $\delta_n$ s.t. $\lim_{n \to +\infty} \delta_n = 0$, which will yield the statement (ii).

It remains to prove (B.10). For this we will use Lemma B.2, applied to $f(\cdot) = \tilde{p}_n(\cdot) = p_n(\cdot + \omega_{\infty,n})$. A priori, we will choose $\delta_n$ so that

$$ (B.11) \qquad \delta_n < \frac{1}{4} \min(\operatorname{Im} \omega_{\infty,n}, \sup_k |\omega_{\infty,k+1} - \omega_{\infty,k}|) \le \frac{1}{4r(\eta-1)} \log \frac{2}{\sqrt{\eta}}, $$

so that $B(\omega_{\infty,n}, \delta_n) \subset \mathbb{C}^+$, and $B(\omega_{\infty,n}, \delta_n) \cap B(\omega_{\infty,n+1}, \delta_{n+1}) = \emptyset$.

Let us now prove (B.10). We need a lower bound on $p'_n(\omega_{\infty,n})$ and an upper bound on $p_n$ inside $B(\omega_{\infty,n}, \delta)$. Step 2.1. A lower bound on $p'_n(\omega_{\infty,n})$. Because $(H_0^{(1)}(z))' = -H_1^{(1)}(z)$ ([1, 9.1.3, 9.1.28]), we have that

$$ (B.12) \qquad p'_n(\omega_{\infty,n}) = -c \frac{H_1^{(1)}(\omega_{\infty,n}c)}{H_0^{(1)}(\omega_{\infty,n}r)} + r \frac{H_1^{(1)}(\omega_{\infty,n}r)P(\omega_{\infty,n})}{H_0^{(1)}(\omega_{\infty,n}r)}. $$

We consider (B.12) for $n \to +\infty$, i.e. $|\omega_{\infty,n}| \to +\infty$. Using the asymptotic expansions (B.4), (B.6) and (B.9) in (B.12) results in the following expansion for $p'_n$:

$$ p'_n(\omega_{\infty,n}) = i \frac{c}{\sqrt{\eta}} e^{i(\eta-1)\omega_{\infty,n}r} + \frac{ir}{2} + O(n^{-1}), \quad n \to +\infty. $$

Because $\omega_{\infty,n}$ solves (B.7), we have

$$ p'_n(\omega_{\infty,n}) = -\frac{i(c-r)}{2} + O(n^{-1}), $$

hence for $n$ sufficiently large it holds that

$$ (B.13) \qquad |p'_n(\omega_{\infty,n})| \ge \frac{r}{4}(\eta - 1). $$

Step 2.2. An upper bound on $p_n$. Again, we are interested in the case when $n \to +\infty$, and $\omega \in B(\omega_{\infty,n}, \delta_n)$.

From the definition of $p_n(\omega)$, the bound (B.9) and (B.6) it follows that, for all $n$ sufficiently large, with some $C > 0$ independent of $n$,

$$ |p_n(\omega)| \le \eta^{-1/2} e^{-(\eta-1) \operatorname{Im} \omega r} + \frac{1}{2} + C n^{-1}. $$
For $\omega \in B(\omega_{\infty,n}, \delta_n) \subset \mathbb{C}^+$, and using $\eta > 1$, one sees that there exists $C_\eta > 0$ s.t.

$$ (B.14) \qquad |p_n(\omega)| \le C_\eta, \text{ for all } n \text{ sufficiently large.} $$

**Step 2.3.** *Proof of (B.10).* Using (B.13) and (B.14), and applying Lemma B.2 to $\tilde{p}_n(\cdot) := p_n(\cdot + \omega_{\infty,n})$, we conclude that, for all $n$ sufficiently large

$$ p_n(B(\omega_{\infty,n}, \delta_n)) \supset B \left( 0; \frac{\delta_n^2 r^2 (\eta - 1)^2}{96C_\eta} \right). $$

To obtain the desired result, it remains to choose $\delta_n$ e.g. so that $|\delta_n| < n^{-1/4}$ (for all $n$ sufficiently large (B.11) will be automatically satisfied).

To summarize, it follows that for all $n$ sufficiently large, there exists $\omega_{k(n)} \in B(\omega_{\infty,n}; \delta_n) \subset \mathbb{C}^+$ s.t. (B.5) is satisfied with $\omega = \omega_{k(n)}$. This proves (i).

Because $\omega_{k(n)} \in B(\omega_{\infty,n}; \delta_n)$, $\delta_n \to 0$ as $n \to +\infty$, and using explicit expressions of $\omega_{\infty,n}$ in (B.8), we see the validity of (ii).

Finally, (iii) follows by considering an explicit expression of $\Lambda$, when $P(\omega) = -\frac{1}{2}$. $\square$

We now have all necessary ingredients to prove Proposition 2.2.

*Proof of Proposition 2.2.* For brevity, we will omit the index 1 in $\lambda_{FL}^1$ and $u_{FL}^1$. First of all, let us prove the result for $\lambda_{FL}(t)$, and next discuss how to extend the reasoning to $u_{FL}$. More precisely, we will show that there exist $u_0$, $u_1 \in C_0^\infty(\Omega^c)$, $c, B > 0$ s.t., for some sequence $t_n \to +\infty$, one has that $\|\lambda_{FL}(t_n)\| \ge c e^{Bt_n}$.

Let us assume the opposite, i.e. that for all $u_0, u_1 \in C_0^\infty(\Omega^c)$, for all $B > 0$,

$$ (B.15) \qquad \limsup_{t \to +\infty} (e^{-Bt} \|\lambda_{FL}(t)\|) = 0. $$

The assumption implies that for all $\delta > 0$, there exists $C_\delta > 0$, s.t.

$$ \|\lambda_{FL}(t)\| \le C_\delta e^{\delta t}, \quad \text{for all } t > 0. $$

Then $\omega \mapsto \hat{\lambda}_{FL}$ is $\mathbb{C}^3$-valued analytic function in $\mathbb{C}^+$. We will arrive at contradiction by choosing $u^{\text{inc}}$ so that $\hat{\lambda}_{FL}$ has a pole in $\mathbb{C}^+$. Let us introduce

$$ \mathbf{a}(\omega) = - (\hat{u}^{\text{inc}}(\omega, \mathbf{c}_k))_{k=1}^3 \in \mathbb{C}^3, \quad \text{so that } \hat{\lambda}_{FL}(\omega) = \Lambda(\omega)^{-1} \mathbf{a}(\omega). $$

In the vicinity of $\omega_n$ defined as in Lemma B.1, one has that $\Lambda(\omega)^{-1} = (\omega - \omega_n)^{-1} P_1 + O(1)$, where $P_1$ is an orthogonal projector on $\text{Span}\{\mathbf{1} = (1, 1, 1)^t\}$. To show that $\hat{\lambda}_{FL}$ has a pole in $\omega = \omega_n$, it suffices to find $u_0, u_1$ so that for some $n \in \mathbb{Z}$, $P_1 \mathbf{a}(\omega_n) \ne 0$.

Let us fix $n = n_*$, The choice of $u_0, u_1$ is then as follows. We fix $u_0 = 0$. Let a priori $\text{supp } u_1 = K_0$, where $K_0$ is a compact in $\Omega^c$. Without loss of generality, let us assume that $K_0$ includes a small ball in the vicinity of the midpoint $\mathbf{x}^*$ of the equilateral triangle formed by the points $\mathbf{c}_k$, $k = 1, \dots, 3$. Then

$$ \hat{u}^{\text{inc}}(\omega, \mathbf{x}) = - \int_{K_0} G_{\omega}(\|\mathbf{x} - \mathbf{y}\|) u_1(\mathbf{y}) d\mathbf{y}, $$

and we need to choose $u_1$ so that

$$ (B.16) \qquad \alpha = (\mathbf{1}, \mathbf{a}(\omega_{n_*})) = - \int_{K_0} \sum_{k=1}^{3} G_{\omega_{n_*}}(\|\mathbf{c}_k - \mathbf{y}\|) u_1(\mathbf{y}) d\mathbf{y} \ne 0. $$
Evidently, such a function $u_1 \in C_0^\infty(K_0)$ exists, since otherwise one would have had that $\sum_{k=1}^3 G_{\omega_{n_*}}(\|c_k - y\|) = 0$ for all $y \in K_0$, and, in particular, $G_{\omega_{n_*}}(\|c_k - x^*\|) = 0$; this is impossible since $H_0^{(1)}(z)$ does not vanish in $\mathbb{C}^+$.

This construction ensures that $\hat{\lambda}(\omega)$ has a pole in $\omega_{n_*} \in \mathbb{C}^+$, thus a contradiction to (B.15).

Let us now argue that with the above choice of the initial data, $\omega \mapsto \hat{u}_{FL}(\omega)$ is not $L^2$-analytic in $\mathbb{C}^+$; this will prove the statement of the proposition. By (2.8),

$$ \hat{u}_{FL}(\omega, \mathbf{x}) = \sum_{k=1}^{3} \frac{G_{\omega}(\|\mathbf{x} - \mathbf{c}_k\|)}{G_{\omega}(r)} \hat{\lambda}_{FL,k}. $$

In the vicinity of $\omega_{n_*}$, $\hat{\lambda}_{FL}(\omega) = (\omega - \omega_{n_*})^{-1}\alpha\mathbf{1} + O(1)$, cf. (B.16), and thus $\hat{u}_{FL}$ has the following expansion:

$$ \hat{u}_{FL}(\omega, \mathbf{x}) = (\omega - \omega_{n_*})^{-1} (\mathbf{f}(\omega_{n_*}, \mathbf{x}), \mathbf{1}) + O(1), \quad f_k(\omega, \mathbf{x}) = \alpha \frac{G_{\omega}(\|\mathbf{x} - \mathbf{c}_k\|)}{G_{\omega}(r)}. $$

It remains to prove that there exists a compact $K$, s.t. $\hat{u}_{FL}: \mathbb{C}^+ \to L^2(K)$ is not analytic. For this it is sufficient to show that $(\mathbf{f}(\omega_{n_*}, \mathbf{x}), \mathbf{1}) \neq 0$ for $\mathbf{x}$ lying in a subset of $K$ of a non-zero Lebesgue measure. In particular, we will choose $K$ as a small ball around $\mathbf{x}^*$, where $\mathbf{x}^*$ is as defined above. We remark that

$$ \mathbf{f}(\omega_{n_*}, \mathbf{x}^*) = \alpha \frac{G_{\omega_{n_*}}(\|\mathbf{x}^* - \mathbf{c}_k\|)}{G_{\omega_{n_*}}(r)} \mathbf{1} \neq 0, $$

because $z \mapsto H_0^{(1)}(z)$ does not vanish in $\mathbb{C}^+$. Moreover, $\mathbf{x} \mapsto \mathbf{f}(\omega_{n_*}, \mathbf{x})$ is continuous in the vicinity of $\mathbf{x}^*$, which proves the desired result. □

**Appendix C. Proof of Corollary 4.2.** The convergence property is immediate from Theorem 4.1 and (2.3).

The uniform stability property relies on the triangle inequality

$$ |u_G^\varepsilon(t, \mathbf{x}_\varepsilon)| \leq \text{error} + |u^\varepsilon(t, \mathbf{x}_\varepsilon)|. $$

It remains to bound $|u^\varepsilon(t, \mathbf{x}_\varepsilon)|$ and $\|\partial_t^5 u^{\text{inc}}\|_{H_\varepsilon^3(0,T; W^{1,\infty}(\mathbb{R}^2))}$ by some Sobolev norms of $u_0$, $u_1$. Let us sketch how this can be done.

To obtain a bound on $u^\varepsilon$ we start by remarking that, by the Sobolev embedding theorem [16, Theorem 1.4.4.1], $|u^\varepsilon(t, \mathbf{x}_0)| \lesssim \|u^\varepsilon(t)\|_{H^2(B(\mathbf{x}_0, r))}$, where $r > 0$ is s.t. $B(\mathbf{x}_0, r) \subset \Omega^c$. To obtain a bound on $\|u^\varepsilon(t)\|_{H^2(B(\mathbf{x}_0, r))}$, we use the definition $u^\varepsilon = u_{\text{tot}}^\varepsilon - u^{\text{inc}}$ and the energy argument applied to $u_{\text{tot}}^\varepsilon$ and $u^{\text{inc}}$. More precisely, by the triangle inequality, we have that

$$ \|\partial_t u^\varepsilon(t)\|_{\Omega^c} + |u^\varepsilon(t)|_{H^1(\Omega^c)} \leq 2(\|u_0\| + \|u_1\|). $$

This immediately yields a bound on $|u^\varepsilon(t)|_{H^1(B(\mathbf{x}_0, r))}$.

To control $\|u^\varepsilon(t)\|_{B(\mathbf{x}_0, r)}$, we use that for $v(0) = 0$, $|v(t)| \leq t\|v'\|_{L^2(0,t)}$, thus

$$ \|u^\varepsilon(t)\|_{\Omega^c} \leq 2t(\|u_0\| + \|u_1\|). $$

Finally, by using elliptic regularity results, cf. [35, Theorem 4.16], to bound $\|D^2 u^\varepsilon(t)\|_{B(\mathbf{x}_0, r)}$, it suffices to bound $\|\Delta u^\varepsilon(t)\|_{B(\mathbf{x}_0, r+\delta)}$, where $\delta > 0$ is s.t. $B(\mathbf{x}_0, r+\delta) \subset \Omega^c$. We have that $\|\Delta u^\varepsilon(t)\|_{B(\mathbf{x}_0, r+\delta)} = \|\partial_t^2 u^\varepsilon(t)\|_{B(\mathbf{x}_0, r+\delta)}$. The bound on $\|\partial_t^2 u^\varepsilon(t)\|_{B(\mathbf{x}_0, r+\delta)}$ can then be obtained by the triangle inequality and the
energy argument applied to $\partial_t u_{tot}^{\varepsilon}$, $\partial_t u^{inc}$, which in their turn, solve the respective BVP and free space problem with the initial data $(u_1, \Delta u_0)$. This yields

$$ \|D^2 u^{\varepsilon}(t)\|_{B(x_0,r)} \lesssim \|\partial_t^2 u^{\varepsilon}(t)\|_{B(x_0,r+\delta)} \le 2(\|\Delta u_0\| + \|u_1\|). $$

To bound $\|\partial_t^5 u^{inc}\|_{H^3(0,T;W^{1,\infty}(\mathbb{R}^2))}$, we use the Sobolev embedding theorem [16], to show that $\|\partial_t^5 u^{inc}\|_{H^3(0,T;W^{1,\infty}(\mathbb{R}^2))} \le \|\partial_t^5 u^{inc}\|_{H^3(0,T;H^3(\mathbb{R}^2))}$. Next we remark that $\partial_t^k u^{inc}$ solves the wave equation, and, like before, make use of the initial data regularity and energy arguments. We leave the remaining details to the reader.

## Appendix D. Estimates on Hankel functions.

LEMMA D.1. Let $\omega, \rho > 0$. Then

$$ \sum_{n \in \mathbb{Z} \setminus \{0\}} \left| \frac{H_n^{(1)}(\omega \rho)}{H_n^{(1)}(\omega \varepsilon)} \right|^2 = O(\varepsilon^2), \quad \text{as } \varepsilon \to 0. $$

*Proof.* Since it holds that $H_n^{(1)}(z) = (-1)^n H_{-n}^{(1)}(z)$, see [13, §10.4], it suffices to derive a bound on $\sum_{n \in \mathbb{N}_+^*} \left| \frac{H_n^{(1)}(\omega \rho)}{H_n^{(1)}(\omega \varepsilon)} \right|^2$.

First of all, there exist $\varepsilon_1 > 0$, $c_1 > 0$, s.t. for all $\varepsilon < \varepsilon_1$ and $n \ge 1$,

$$ (D.1) \qquad \frac{1}{|H_n^{(1)}(\omega\varepsilon)|} \le c_1 \frac{(\omega\varepsilon)^n}{n!}. $$

For $n=1$, the above bound follows from the expansion of $H_1^{(1)}$ into series, cf. [13, §10.8 and 10.2.2]. For $n \ge 2$, by Lemma 7 in [9], we have that

$$ \frac{1}{|H_n^{(1)}(\omega\varepsilon)|} \le \frac{2(\omega\varepsilon)^{n-2}}{n!|H_2^{(1)}(\omega\varepsilon)|}, \quad \text{for } 0 < \omega\varepsilon < 1 \text{ and } n \ge 2. $$

Then the bound (D.1) can be obtained again by considering the series expansion of $H_2^{(1)}$, cf. [13, §10.8 and 10.2.2].

Next, by the asymptotic expansions of Hankel functions as $n \to +\infty$ [13, §10.19], it follows that there exist $N_{\omega\rho}, c_2 > 0$, s.t. for all $n \ge N_{\omega\rho}$,

$$ (D.2) \qquad |H_n^{(1)}(\omega\rho)| \le c_2 \sqrt{\frac{2}{\pi n}} \left(\frac{2n}{\varepsilon\omega\rho}\right)^n \le c_2(n-1)! \left(\frac{2}{\rho\omega}\right)^n. $$

Finally, we split the sum in the statement of the lemma into two terms and bound them using (D.1), (D.2), assuming in particular that $\varepsilon < \rho/2$:

$$ \begin{aligned} \sum_{n \in \mathbb{N}_+^*} \left| \frac{H_n^{(1)}(\omega \rho)}{H_n^{(1)}(\omega \varepsilon)} \right|^2 &\le \sup_{0 < n \le N_{\omega \rho}} |H_n^{(1)}(\omega \rho)|^2 \sum_{n=1}^{N_{\omega \rho}-1} |H_n^{(1)}(\omega \varepsilon)|^{-2} + \sum_{n=N_{\omega \rho}}^{+\infty} \left| \frac{H_n^{(1)}(\omega \rho)}{H_n^{(1)}(\omega \varepsilon)} \right|^2 \\ &\le C_{\omega \rho} \varepsilon^2. \end{aligned} \quad \square $$

LEMMA D.2. For $\mu > 1$, all $x > 0$, $K_0(\mu x) < K_0(x)$.
Proof. By [13, 10.32.9], the function

$$ (D.3) \qquad x \mapsto K_0(x) = \int_0^{+\infty} e^{-x \cosh t} dt, $$

is strictly monotonically decreasing.

LEMMA D.3. For $\mu \ge 1$, $m \ge 1$, $x \in \mathbb{R}^+$, it holds that $K_m(\mu x) < \mu^{-m}K_m(x)$.

Proof. By [13, 10.32.11], for all $\mu > 0$,

$$ (D.4) \qquad K_m(\mu x) = \frac{\Gamma(m + 1/2)(2x)^m}{\sqrt{\pi}\mu^m} \int_0^{+\infty} \frac{\cos \mu t}{(t^2 + x^2)^{m+1/2}} dt. $$

A straightforward computation yields

$$ (D.5) \qquad \operatorname{sign}(\mu^m K_m(\mu x) - K_m(x)) = \operatorname{sign} G(\mu, x), \quad G(\mu, x) = \int_0^{+\infty} \frac{\cos \mu t - \cos t}{(t^2 + x^2)^{m+1/2}} dt. $$

It remains to prove that for $\mu > 1$ the quantity in the right hand side is negative. Let $x > 0$ be fixed. We have $G(1, x) = 0$ and, integrating by parts (with $m \ge 1$),

$$ \partial_{\mu} G(\mu, x) = - \int_{0}^{+\infty} \frac{t \sin \mu t}{(t^2 + x^2)^{m+1/2}} dt = - \frac{\mu}{(2m-1)} \int_{0}^{+\infty} \frac{\cos \mu t}{(t^2 + x^2)^{m-1/2}} dt. $$

By (D.4), $\operatorname{sign} \partial_{\mu} G = - \operatorname{sign} K_{m-1}(\mu x) < 0$, since $K_{\nu}(x) > 0$ for $x > 0$ (compare [13, 10.32.8]). Thus $G(\mu, x) < G(1, x) = 0$, and the desired result follows from (D.5). □

REFERENCES

[1] M. ABRAMOWITZ AND I. A. STEGUN, Handbook of mathematical functions with formulas, graphs, and mathematical tables, vol. 55 of National Bureau of Standards Applied Mathematics Series, 1964.

[2] A. BAMBERGER AND T. H. DUONG, Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d'une onde acoustique. I, Math. Methods Appl. Sci., 8 (1986), pp. 405-435.

[3] L. BANJAI, C. LUBICH, AND F.-J. SAYAS, Stable numerical coupling of exterior and interior problems for the wave equation, Numerische Mathematik, 129 (2015), pp. 611-646.

[4] H. BARUCQ, J. DIAZ, V. MATTESI, AND S. TORDEUX, Asymptotic behavior of acoustic waves scattered by very small obstacles, ESAIM Math. Model. Numer. Anal., 55 (2021), pp. S705-S731.

[5] E. BÉCACHE AND M. KACHANOVSKA, Stability and convergence analysis of time-domain perfectly matched layers for the wave equation in waveguides, SIAM J. Numer. Anal., 59 (2021), pp. 2004-2039.

[6] A. BENDALI, P.-H. COCQUET, AND S. TORDEUX, Scattering of a scalar time-harmonic wave by n small spheres by the method of matched asymptotic expansions, Numerical Analysis and Applications, 5 (2012), pp. 116-123.

[7] A. BENDALI, P.-H. COCQUET, AND S. TORDEUX, Approximation by multipoles of the multiple acoustic scattering by small obstacles in three dimensions and application to the Foldy theory of isotropic scattering, Arch. Ration. Mech. Anal., 219 (2016), pp. 1017-1059.

[8] T. BETCKE, J. PHILLIPS, AND E. A. SPENCE, Spectral decompositions and nonnormality of boundary integral operators in acoustic scattering, IMA J. Numer. Anal., 34 (2014), pp. 700-731.

[9] M. CASSIER AND C. HAZARD, Multiple scattering of acoustic waves by small sound-soft obstacles in two dimensions: mathematical justification of the Foldy-Lax model, Wave Motion, 50 (2013), pp. 18-28.

[10] D. P. CHALLA AND M. SINI, On the justification of the Foldy-Lax approximation for the acoustic scattering by small rigid bodies of arbitrary shapes, Multiscale Model. Simul., 12 (2014), pp. 55-108.

[11] ——, The Foldy-Lax approximation of the scattered waves by many small bodies for the Lamé system, Math. Nachr., 288 (2015), pp. 1834-1872.

[12] J. B. CONWAY, Functions of one complex variable, Springer-Verlag, New York-Heidelberg, 1973. Graduate Texts in Mathematics, 11.
[13] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.1 of 2021-03-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.

[14] D. M. ĔDUS, The principle of limiting amplitude, Uspehi Mat. Nauk, 24 (1969), pp. 91–156.

[15] L. L. FOLDY, *The multiple scattering of waves. I. General theory of isotropic scattering by randomly distributed scatterers*, Phys. Rev. (2), 67 (1945), pp. 107–119.

[16] P. GRISVARD, Elliptic problems in nonsmooth domains, vol. 69 of Classics in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011. Reprint of the 1985 original [MR0775683], With a foreword by Susanne C. Brenner.

[17] M. HASSAN AND B. STAMM, *An integral equation formulation of the $N$-body dielectric spheres problem. Part I: Numerical analysis*, ESAIM Math. Model. Numer. Anal., 55 (2021), pp. S65–S102.

[18] G. C. HSIAO AND W. L. WENDLAND, *Boundary integral equations*, vol. 164 of Applied Mathematical Sciences, Springer-Verlag, Berlin, 2008.

[19] D. V. KORIKOV, Asymptotic behavior of solutions of the wave equation in a domain with a small hole, Algebra i Analiz, 26 (2014), pp. 164–199.

[20] D. V. KORIKOV AND B. A. PLAMENEVSKIİ, Asymptotics of solutions for stationary and nonstationary Maxwell systems in a domain with small holes, Algebra i Analiz, 28 (2016), pp. 102–170.

[21] R. KRESS, *Minimizing the condition number of boundary integral operators in acoustic and electromagnetic scattering*, Quart. J. Mech. Appl. Math., 38 (1985), pp. 323–341.

[22] R. KRESS AND W. T. SPASSOV, *On the condition number of boundary integral operators for the exterior Dirichlet problem for the Helmholtz equation*, Numer. Math., 42 (1983), pp. 77–95.

[23] J. LABAT, *Modélisation multi-échelle de la diffraction des ondes électromagnétiques par de petits obstacles*, PhD thesis, 2019. Thèse de doctorat dirigée par Tordeux, Sébastien et Péron, Victor Mathématiques Pau 2019.

[24] J. LABAT, V. PÉRON, AND S. TORDEUX, *Equivalent multipolar point-source modeling of small spheres for fast and accurate electromagnetic wave scattering computations*, Wave Motion, 92 (2020), pp. 102409, 21.

[25] M. LAX, Multiple scattering of waves, Rev. Modern Physics, 23 (1951), pp. 287–310.

[26] J. LIAO, Foldy-Lax approximation on multiple scattering by many small scatterers, Appl. Anal., 92 (2013), pp. 2547–2560.

[27] J. LIAO AND C. JI, Extended Foldy-Lax approximation on multiple scattering, Math. Model. Anal., 19 (2014), pp. 85–98.

[28] C. LUBICH, Convolution quadrature and discretized operational calculus. I, Numer. Math., 52 (1988), pp. 129–145.

[29] ——, Convolution quadrature and discretized operational calculus. II, Numer. Math., 52 (1988), pp. 413–425.

[30] C. LUBICH, On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations, Numer. Math., 67 (1994), pp. 365–389.

[31] S. MARMORAT, Modèles asymptotiques et simulation numérique pour la diffraction d'ondes par des petites hétérogénéités, PhD thesis, Université Paris-Saclay, 11 2015.

[32] P. A. MARTIN, Multiple scattering, vol. 107 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, Cambridge, 2006. Interaction of time-harmonic waves with $N$ obstacles.

[33] ——, *Time-domain scattering*, vol. 180 of Encyclopedia of Mathematics and its Applications, Cambridge University Press, Cambridge, 2021.

[34] V. MATTESI, Propagation des ondes dans un domaine comportant des petites hétérogénéités : modélisation asymptotique et calcul numérique, theses, Université de Pau et des pays de l'Adour, Dec. 2014.

[35] W. MCLEAN, Strongly elliptic systems and boundary integral equations, Cambridge University Press, Cambridge, 2000.

[36] M. I. MISHCHENKO, *Electromagnetic scattering by particles and particle groups: an introduction*, Cambridge University Press, 2014.

[37] F. W. J. OLVER, *Asymptotics and special functions*, AKP Classics, A K Peters, Ltd., Wellesley, MA, 1997. Reprint of the 1974 original [Academic Press, New York; MR0435697 (55 #8655)].

[38] A. G. RAMM, Wave scattering by small bodies, Rep. Math. Phys., 21 (1985), pp. 69–77.

[39] F.-J. SAYAS, Retarded potentials and time domain boundary integral equations, vol. 50 of Springer Series in Computational Mathematics, Springer, [Cham], 2016. A road map.

[40] M. SINI, H. WANG, AND Q. YAO, *Analysis of the acoustic waves reflected by a cluster of small holes in the time-domain and the equivalent mass density, Multiscale Model. Simul.*, 19 (2021), pp. 1083–1114.

[41] A. TOURIN, M. FINK, AND A. DERODE, *Multiple scattering of sound*, Waves in random media, 10 (2000), p. R31.