# Discretization in multilayered media with high contrasts: is it all about the boundaries?

Camille Carvalho, Stéphanie Chaillat, Chrysoula Tsogka, Elsie A Cortes

## ► To cite this version:

Camille Carvalho, Stéphanie Chaillat, Chrysoula Tsogka, Elsie A Cortes. Discretization in multilayered media with high contrasts: is it all about the boundaries?. 2026. hal-05520237

HAL Id: hal-05520237

https://hal.science/hal-05520237v1

Preprint submitted on 20 Feb 2026

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

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

<figure><img src="image_1.png" /></figure>
# DISCRETIZATION IN MULTILAYERED MEDIA WITH HIGH CONTRASTS: IS IT ALL ABOUT THE BOUNDARIES?

CAMILLE CARVALHO, STÉPHANIE CHAILLAT, ELSIE CORTES, AND CHRYSOULA TSOGKA

**ABSTRACT.** Wave propagation in multilayered media with high material contrasts poses significant numerical challenges, as large variations in wavenumbers lead to strong reflections and complex transmission of the incoming wave field. To address these difficulties, we employ a boundary integral formulation thereby avoiding volumetric discretization. In this framework, the accuracy of the numerical solution depends strongly on how the material interfaces are discretized. In this work, we demonstrate that standard meshing strategies based on resolving the maximum wavenumber across the domain become computationally inefficient in multilayered configurations, where high wavenumbers are confined to localized subdomains. Through a systematic study of multilayer transmission problems, we show that no simple discretization rule based on the maximum wavenumber or material contrasts emerges. Instead, the wavenumber of the background (exterior) medium plays a dominant role in determining the optimal boundary resolution. Building on these insights, we propose an adaptive approach that achieves uniform accuracy and efficient computation across multiple layers. Numerical experiments for a range of multilayer configurations demonstrate the scalability and robustness of the proposed approach.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. The Multilayer Transmission Problem</td><td>3</td></tr><tr><td>3. Boundary Mesh Optimization for BEM</td><td>5</td></tr><tr><td>4. Adaptive Mesh Refinement</td><td>11</td></tr><tr><td>5. Conclusions and Perspectives</td><td>20</td></tr><tr><td>Appendix A. Analytic Solution for Multilayer Transmission Problems</td><td>21</td></tr><tr><td>Appendix B. Analytic Solutions for the Sound-Hard Boundary Condition</td><td>22</td></tr><tr><td>References</td><td>22</td></tr></table>

## 1. INTRODUCTION

We consider wave propagation in 2D multilayered media modeled by the Helmholtz equation. The domain is modeled as a sequence of $N+1$ concentric layers, separated by smooth interfaces across which transmission conditions are imposed. Such multilayer models are important in various applications such as noninvasive medical imaging [1][2], non-destructive testing [3][4], and especially invisibility cloaking [5][6]. Our long-term motivation is related to optical cloaking, which can be obtained with metamaterials (artificially structured materials with deliberately engineered variations in material properties that enable complex wave behavior). Simulating metamaterials and optimizing their design are of particular interest as an insight into their underlying physics, bridging the gap between theoretical understanding and practical implementation. One of the challenges in designing metamaterials is dealing with large contrasts in wavenumbers $k_i$, $i \in [0, N]$ across the interfaces (e.g., $\frac{k_{i+1}}{k_i} \ll 1$ or $\frac{k_{i+1}}{k_i} \gg 1$), since they lead to strong reflection or complex transmission of the incident wave field.

*Date:* February 20, 2026.

*Key words and phrases.* Multilayered Scattering, Automatic Mesh Adaptation, Wave Propagation, 2D Helmholtz, Boundary Integral Equations.
Configurations with high material contrasts between layers inherently give rise to both high and low wavenumber regimes. Mesh-based numerical approaches encounter significant challenges in the high-$k_i$ regions, where the shorter wavelength $\lambda_i$ of the wave field demands finer discretization and greater computational effort. A common approach to discretizing wave problems is to maintain mesh size $h$ such that $kh \ll 1$. This conservative choice prioritizes accuracy by ensuring sufficient resolution of the solution, even at the cost of overestimating the number of points per wavelength [7]. For multiple layers with different wavenumbers, the discretization is often governed by the maximum wavenumber $k_{max}$ across all layers, i.e., according to the most constraining domain. This naive approach constrains the mesh to uniformly resolve the highest wavenumber although some parts of the domain correspond to lower wavenumbers. Furthermore, the pollution effect implies that the number of degrees of freedom must grow faster than the magnitude of $k_i$ to maintain accuracy [8]. This effect implies that increasingly fine volumetric meshes are necessary to maintain accuracy. Consequently, for multilayer problems, volumetric methods suffer from rapidly increasing computational costs since achieving accurate resolution requires mesh refinement even in regions where the wave field remains relatively simple. While domain decomposition approaches combined with adaptive mesh refinement in each subdomain can provide effective solutions to this problem (see [9] and references therein), we follow a different approach here. This deterioration in computational scalability motivates our use of a boundary integral formulation, where discretization is confined to material interfaces instead of the entire volume.

Boundary element methods (BEM) offer an efficient framework for approximating solutions to partial differential equations (PDE). They rely on boundary integral equations (BIE), which reformulate the original PDE problem as an equation defined solely on the boundaries or interfaces. BEM have been applied to a wide range of applications e.g. [10]) including wave propagation problems, as they naturally enforce radiation conditions and largely avoid the pollution effect [11, 12]. A key advantage of this approach is the reduction of dimensionality: only boundary data must be discretized, avoiding the need to mesh the entire computational domain, unlike Finite Element Methods (FEM) or Finite-Difference (FD) methods. This makes BEM especially attractive for layered media [13, 14], although multilayer configurations lead to densely populated matrices [15]. Moreover, the singular or nearly-singular kernels appearing in the BIE require careful treatment [16]. While techniques for weakly [11] and nearly-singular kernels [17, 18] are well-established for single-domain problems, multilayered domains require a boundary integral framework configured to accurately handle multiple interfaces in relation to each other. We have previously derived a modified boundary equation method to find accurate solutions to the multilayered transmission problem when interfaces are close to each other [19]. By reducing the influence of nearly singular kernels on the boundary solution, this method avoids the need for fine interface refinement to obtain a precise solution.

One of the critical factors in BEM is the discretization of the interfaces, which strongly influences both the accuracy of the solution and the efficiency of the computation. Coarse meshes or meshes that fail to capture the boundary geometry or key solution features can produce significant errors, while overly refined meshes increase computational cost. In the current work, we tackle this challenge with adaptive mesh refinement. We propose an effective, systematic way to optimize the discretization by refining the mesh only where it is needed, improving both the quality of our solutions and the computational performance. In general, mesh adaptation in the context of BEM remains an open area of research where several theoretical and practical questions are being explored [20, 21]. While no universal method exists in the BEM community, the situation is even more challenging for multilayered domains. In this work, we specifically investigate the use of metric-based mesh refinement for this problem [22, 23].

By focusing on the discretization of the **boundaries** rather than the entire domain, this work presents a new perspective on discretization of multilayer media for wave propagation problems. In Section 2, we begin by formulating the multilayer transmission problem with boundary integral equations and present high-order quadrature rules for their discretization.
In Section 3, we explore boundary mesh optimization, using reference analytic solutions for single-material domains to construct a posteriori estimates that provide insight into solution accuracy for multilayer configurations with increasing contrasts. We observe that the accuracy of the numerical solutions is not driven by the highest wavenumber. Building on this, in Section 4, we extend our focus to adaptive mesh generation, automatically refining boundary meshes in configurations with multiple boundaries and more complex geometries. Section 5 concludes with a discussion of possible directions for future research.

## 2. THE MULTILAYER TRANSMISSION PROBLEM

### 2.1. Problem Formulation for Concentric Regions.
Consider a domain composed of $N+1$ concentric regions in $\mathbb{R}^2$, separated by closed, smooth interfaces $\Gamma_j$ for $j = 0, \dots, N-1$ as illustrated in Figure 1. Let $\Omega_j \subset \mathbb{R}^2$ for $j = 0, \dots, N$ denote the concentric subdomains, such that $\Omega_0$ is the outermost region, or the background medium, and $\Omega_N$ is the innermost region. The interfaces $\Gamma_j$ separate $\Omega_j$ and $\Omega_{j+1}$ for $j = 0, \dots, N-1$. Each region is characterized by a wavenumber $k_j \in \mathbb{R}_+$ (determined by its material properties). We define the **contrast** between the material properties of adjacent regions $\Omega_j$ and $\Omega_{j+1}$ as $\beta_j := k_{j+1}^2/k_j^2$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Schematic of the concentric layered transmission problem. Each interface $\Gamma_j$ separates regions with wavenumbers $k_j$ and $k_{j+1}$. The exterior field $u_0$ includes an incoming incident wave $u_0^{in}$.</figcaption></figure>

Working in the time-harmonic regime, the total field solution in $\Omega_0$ is the sum of the incident and scattered fields, and the incident field is a plane wave of the form $u_0^{in} = e^{ik_0 \mathbf{a} \cdot \mathbf{x}}$ ($^1$). In each layer $\Omega_j$ the total field $u_j$ satisfies the Helmholtz equation [24]:

$$ (1) \qquad \Delta u_j + k_j^2 u_j = 0 \quad \text{in } \Omega_j. $$

Furthermore, across each interface we impose the following transmission conditions:

$$ (2) \qquad \begin{array}{ll} u_{j+1} = u_j & \text{on } \Gamma_j \quad j = 0, \dots, N-1, \\ \frac{\partial u_{j+1}}{\partial n_j} = \beta_j \frac{\partial u_j}{\partial n_j} & \text{on } \Gamma_j \quad j = 0, \dots, N-1, \end{array} $$

where $n_j$ denotes the outward unit normal on $\Gamma_j$. Additionally, we assume the scattered solution in $\Omega_0$, $u_0^{sc} = u_0 - u_0^{in}$, satisfies the Sommerfeld radiation condition:

$$ (3) \qquad \lim_{r \to \infty} \sqrt{r} \left( \frac{\partial u_0^{sc}}{\partial r} - ik_0 u_0^{sc} \right) = 0. $$

This condition ensures that the outgoing wave solution decays as $r$ tends to infinity.

<footnote n="1">¹In this work, we choose $\mathbf{a} = (0, 1)$ so that the incident field $u^{in} = e^{ik_0 y}$ for $\mathbf{x} := (x, y)$ is an incoming plane wave propagating in the positive y-direction. But all results in this paper hold for any $\mathbf{a}$.</footnote>
With the multilayer transmission model in place, we now turn to finding the boundary integral formulation for the BEM.

## 2.2. Boundary Integral Equation Solutions.
BIE, as the name implies, are representations of the solutions in terms of boundary integrals. For the total field satisfying the Helmholtz equation in each layer, a boundary integral formulation may be obtained using a combination of single- and double-layer boundary integral operators $S_{i,j}$ and $D_{i,j}$ respectively [11]. They are defined as

$$ (4) \qquad S_{i,j}[\mu](\mathbf{x}) = \int_{\Gamma_i} \Phi_j(\mathbf{x}, \mathbf{y})\mu(\mathbf{y})d\Gamma_{i,\mathbf{y}}, \qquad \mathbf{x} \in \mathring{\Omega}_j, $$

$$ (5) \qquad D_{i,j}[\mu](\mathbf{x}) = \int_{\Gamma_i} \frac{\partial \Phi_j}{\partial n_i(\mathbf{y})}(\mathbf{x}, \mathbf{y})\mu(\mathbf{y})d\Gamma_{i,\mathbf{y}}, \qquad \mathbf{x} \in \mathring{\Omega}_j, $$

where $\Gamma_i$ (for $i = j - 1$ or $j$) denotes the interface on which they are defined, and $\Phi_j$ is the fundamental solution to the Helmholtz equation in $\Omega_j$. In 2D the fundamental solution is given by:

$$ (6) \qquad \Phi_j(\mathbf{x}, \mathbf{y}) = \frac{i}{4} H_0^{(1)}(k_j |\mathbf{x} - \mathbf{y}|) $$

where $H_0^{(1)}$ denotes the Hankel function of first kind. There are multiple boundary integral formulations, but for the one we choose the unknown densities $\mu \in C^2(\Gamma_i)$ in $S_{i,j}$ and $D_{i,j}$ define the boundary data. More precisely, in the double-layer integral operator, $\mu$ corresponds to the trace $u_j$, and in the single-layer integral operator it corresponds to the normal trace $\frac{\partial u_j}{\partial n_i}$. For simplicity, we denote the traces $u_j(\mathbf{y})$ and $\frac{\partial u_j}{\partial n_i}(\mathbf{y})$ respectively by

$$ u_{i,j} := u_j(\mathbf{x}), \quad \text{and} \quad \partial_n u_{i,j} := \frac{\partial u_j}{\partial n_i}(\mathbf{x}), \qquad \mathbf{x} \in \Gamma_i. $$

Using the representation formula [11], and utilizing the layer potentials, we construct the following boundary integral representations for the solutions for $j = 1, \dots, N - 1$:

$$ (7) \qquad \begin{aligned} u_0(\mathbf{x}) &= u_0^{\text{in}}(\mathbf{x}) + \mathcal{D}_{0,0}[u_{0,0}](\mathbf{x}) + \mathcal{S}_{0,0}[\partial_n u_{0,0}](\mathbf{x}) & \qquad \text{in } \mathring{\Omega}_0, \\ u_j(\mathbf{x}) &= -\mathcal{D}_{j-1,j}[u_{j-1,j}](\mathbf{x}) + \mathcal{S}_{j-1,j}[\partial_n u_{j-1,j}](\mathbf{x}) \\ & \qquad + \mathcal{D}_{j,j}[u_{j,j}](\mathbf{x}) - \mathcal{S}_{j,j}[\partial_n u_{j,j}](\mathbf{x}) & \qquad \text{in } \mathring{\Omega}_j, \\ u_N(\mathbf{x}) &= -\mathcal{D}_{N-1,N}[u_{N-1,N}](\mathbf{x}) + \mathcal{S}_{N-1,N}[\partial_n u_{N-1,N}](\mathbf{x}) & \qquad \text{in } \mathring{\Omega}_N \end{aligned} $$

The solution in each domain is determined from the traces of $u$ and $\partial u$. To obtain these traces $u_{i,j}$ and $\partial u_{i,j}$ we construct a system of representation equations defined on the interfaces $\Gamma_j$ in terms of layer potentials. The layer potentials can be obtained by taking the limits of the single- and double-layer layer integral operators on each side of the boundaries. The single-layer integral operator is a continuous operator such that

$$ S_{i,j}[\mu](\mathbf{x}) \underset{\mathbf{x} \to \mathbf{x}^{b, \pm}}{\longrightarrow} S_{i,j}[\mu](\mathbf{x}^b), $$

but the double-layer operator exhibits a jump such that

$$ \mathcal{D}_{i,j}[\mu](\mathbf{x}) \underset{\mathbf{x} \to \mathbf{x}^{b, \pm}}{\longrightarrow} \pm \frac{1}{2}\mu(\mathbf{x}^b) + D_{i,j}[\mu](\mathbf{x}^b) $$

where $\mathbf{x}^b$ is on the boundary.

By introducing the transmission conditions (2) in the boundary integral equation system (8), we obtain a solvable system of $2N$ equations with $2N$ unknowns, namely the traces $u_{i,j}$ and normal traces $\partial u_{i,j}$ respectively that we can directly solve for these quantities. For example, if
we substitute $u_{j-1,j} = u_{j-1,j-1}$ and $\partial_n u_{j-1,j} = \beta_{j-1} \partial_n u_{j-1,j-1}$, we obtain the following boundary integral equations:

$$ \begin{array}{ll} \frac{1}{2} u_0(\mathbf{x}^b) = u_0^{in}(\mathbf{x}^b) + D_{0,0}[u_{0,0}](\mathbf{x}^b) + S_{0,0}[\partial_n u_{0,0}](\mathbf{x}^b) & \text{on } \Gamma_0, \\ \\ \frac{1}{2} u_j(\mathbf{x}^b) = -D_{j-1,j}[u_{j-1,j-1}](\mathbf{x}^b) + \beta_{j-1} S_{j-1,j}[\partial_n u_{j-1,j-1}](\mathbf{x}^b) & \\ \qquad + \mathcal{D}_{j,j}[u_{j,j}](\mathbf{x}^b) - \mathcal{S}_{j,j}[\partial_n u_{j,j}](\mathbf{x}^b) & \text{on } \Gamma_{j-1}, \\ \\ \frac{1}{2} u_j(\mathbf{x}^b) = -\mathcal{D}_{j-1,j}[u_{j-1,j-1}](\mathbf{x}^b) + \beta_{j-1} \mathcal{S}_{j-1,j}[\partial_n u_{j-1,j-1}](\mathbf{x}^b) & \\ \qquad + \mathcal{D}_{j,j}[u_{j,j}](\mathbf{x}^b) - \mathcal{S}_{j,j}[\partial_n u_{j,j}](\mathbf{x}^b) & \text{on } \Gamma_j, \\ \\ \frac{1}{2} u_N(\mathbf{x}^b) = -D_{N-1,N}[u_{N-1,N}](\mathbf{x}^b) + \beta_{N-1} S_{N-1,N}[\partial_n u_{N-1,N-1}](\mathbf{x}^b) & \text{on } \Gamma_{N-1} \end{array} $$

There exist various approaches to discretizing this system. Since we consider multilayer problems defined by smooth periodic boundaries, we choose to parameterize the boundary integrals in terms of $\mathbf{x}(\theta) = (x_1(\theta), x_2(\theta))$ for $0 \le \theta \le 2\pi$. A well-known approach for the defined boundary discretization is proposed by Kress [11]. It consists of expressing the boundary integrals in terms of trigonometric polynomials equi-spaced on $\theta_m$ such that

$$ (9) \qquad \theta_m := \frac{2\pi m}{M} \quad \text{for } m = 0, \dots, M-1. $$

A key point for BEM is the evaluation of the (singular) boundary integrals. In the case of 2D Helmholtz, layer potentials in (8) contain a kernel with a logarithmic singularity. An effective treatment method is the Kress Quadrature Rule (KQR) [11]. KQR splits each singular kernel into an analytically integrable part plus a smooth remainder, allowing the use of corrected trapezoidal weights on the latter.

From an algebraic point of view, the size of the system of equations to solve is determined by the number of boundaries $N$ and the number of equispaced points $M$ on each interface. Once we determine the unknown traces and normal traces, the solution is obtained using the boundary integral representation in the volume (7), approximated with the periodic trapezoid rule (PTR) where we evaluate on the same uniform distribution of $\theta_m$ in (9):

$$ (10) \qquad \int_0^{2\pi} f(\theta) d\theta \approx \Delta\theta \sum_{m=0}^{M} f(\theta_m) = \frac{2\pi}{M} \sum_{m=0}^{M-1} f(\theta_m) $$

For smooth, closed boundaries, the approach gives spectral convergence with respect to the number of quadrature points $M$.

A key question for this problem is how to determine the optimal number of points $M$ to use on each interface. Traditionally, a general rule based on the maximum wavenumber $k_{max}$ in the domain is used to set $M$, either uniformly across all interfaces or scaled according to the size of each boundary. When $k_{max}$ is large, the field becomes highly oscillatory, requiring a greater number of points to be accurately represented. One of our objectives is to go beyond the standard approach and determine the optimal number of points $M_i$ on each interface $\Gamma_i$, where optimal is defined as the minimal $M_i$ that achieves a desired order of accuracy in the boundary data at that interface. In particular, we aim to investigate whether a rule that accounts for the contrast at each interface, $\beta_j := \frac{k_{j+1}^2}{k_j^2}$, can be inferred. This is a largely unexplored area for wave propagation problems in multilayer media.

In the next section, we introduce an approach to identify the optimal $M_i$ for various contrasts, including comparisons to the sound-hard case, in order to examine any trends observed in multilayer problems.

### 3. BOUNDARY MESH OPTIMIZATION FOR BEM
### 3.1. A Preliminary Approach Based on the Reference Solution.
We focus on configurations with circular boundaries, for which analytic solutions can be derived (see Section A) and used as reference solutions to gain further insight into the problem. We aim to determine the minimum number of points $M_i$ on each interface $\Gamma_i$ required to achieve a target accuracy. Here, accuracy is measured as the absolute $L^2$ error between the BIE solution of (8) on the boundaries and the corresponding analytic solution $u^{ana}$ (21). Our naive preliminary approach to determine the optimal discretization consists of iteratively increasing $M_i$ until the inequality

$$ (11) \qquad \|u_j - u_j^{ana}\|_{L^2(\Gamma_i)} \le \epsilon $$

is satisfied, with $\epsilon$ the desired accuracy (between $10^{-6}$ and $10^{-8}$ in this work). We stop the iterative process when (11) is satisfied for all boundaries.

In the context of multilayer transmission problems, evaluating errors on the interfaces is essential for assessing how accurately the discretization captures the complex behavior of the solution in the surrounding volume. To evaluate the absolute $L^2$ error (11), we use the explicit parameterization of the interfaces $\Gamma_i$ for $i = j - 1$ or $j$. It follows that

$$ (12) \qquad \|u_j - u_j^{ana}\|_{L^2(\Gamma_i)} = \left( \int_{\Gamma_i} |u_j(\mathbf{x}_j(\theta)) - u_j^{ana}(R_i, \theta)|^2 d\theta \right)^{1/2} $$

for $\mathbf{x}_j(\theta) := (R_i \cos \theta, R_i \sin \theta)$.

Since the interfaces $\Gamma_i$ are periodic and parameterized in terms of $\theta$ over the interval $[0, 2\pi)$, we use PTR to evaluate the integral in (12). The boundary is discretized using $M_i$ equispaced points $\mathbf{x}_j(\theta_m)$ for $m = 0, \dots, M_i$ and $\theta_m$ defined as in (9), so that (12)

$$ (13) \qquad \|u_j - u_j^{ana}\|_{L^2(\Gamma_i)} \approx \left( \frac{2\pi}{M_i} \sum_{m=0}^{M_i-1} |u_j(\mathbf{x}_j(\theta_m)) - u_j^{ana}(R_i, \theta_m)|^2 \right)^{1/2} $$

### 3.2. Observations from Analytic Solutions.
We consider both the sound-hard boundary case imposing $\partial_n u_{0,0} = 0$ on $\Gamma_0$ and the transmission boundary case (see A and B for expressions of the associated analytic solutions). The goal of this comparison is to reveal how the properties of an additional layer affect the solution in the previous layer. Rather than focusing solely on the layer with the highest wavenumber, where oscillations are strongest, this analysis helps to understand how the solution evolves across multiple layers.

We consider a simple case with only one interface (a circular boundary at $R_0 = 4$) defining two materials. We analyze the accuracy for different contrasts, $\beta_0 = k_1^2/k_0^2$, where $k_0$ and $k_1$ are the wavenumbers in the outer and inner regions respectively as shown in Figure 2(b). We

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Schematics of the scattering domain with radius $R_0 = 4$, showing two types of boundary conditions: (a) a sound-hard obstacle and (b) a penetrable interface with transmission conditions.</figcaption></figure>
<table><caption>TABLE 1. Summary of single-boundary cases used in numerical experiments, with wavenumbers **k**, boundary radius $R_0$, and contrast $\beta_0$.</caption><thead><tr><th>Case</th><th>k = (k<sub>0</sub>, k<sub>1</sub>, ..., k<sub>N</sub>)</th><th>r = (R<sub>0</sub>, ..., R<sub>N-1</sub>)</th><th>b = (β<sub>0</sub>, ..., β<sub>N-1</sub>)</th></tr></thead><tbody><tr><td>1</td><td>(2, 6)</td><td>4</td><td>9</td></tr><tr><td>2</td><td>(6, 2)</td><td>4</td><td><sup>1</sup>&frasl;<sub>9</sub></td></tr><tr><td>3</td><td>(2, 10)</td><td>4</td><td>25</td></tr><tr><td>4</td><td>(10, 2)</td><td>4</td><td><sup>1</sup>&frasl;<sub>25</sub></td></tr><tr><td>5</td><td>(8, 2, 6)</td><td>(6, 2)</td><td>(<sup>1</sup>&frasl;<sub>16</sub>, 9)</td></tr></tbody></table>

examine all the cases in Table 1. In particular, we explore cases obtained by swapping the interior and exterior wavenumbers. The aim is to examine how the contrast influences the wave propagation and scattering behaviors. All cases are compared to the corresponding sound-hard case, using their respective $k_0$ in the exterior region (see Figure 2(a)).

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Real part of the total field analytic solutions for two contrast cases with $R_0 = 4$. The left plot on each sub-figure shows the solution for the sound-hard boundary condition (Neumann), while the right plot of each sub-figure shows the solution for the transmission boundary condition.</figcaption></figure>

In Figure 3, we plot the corresponding total field solutions for **Cases 1 and 2**. In **Case 1**, we observe that for the transmission case the solution exhibits strong oscillations within the interior region, although much of this complexity is constrained by the boundary and the lower wavenumber of the exterior region. Although this suggests that finer discretization may be necessary, the oscillatory behavior along the boundary remains comparable to the reference solution for the sound-hard boundary condition. In contrast, in **Case 2**, the transmission boundary condition produces wave fields with oscillations comparable to those of the exterior for the sound-hard case. This supports the expectation that resolution should be governed by the higher wavenumber. However, the decrease in oscillations near the boundary indicates a more subtle effect.

For the sharper contrasts in **Cases 3 and 4**, shown in Figure 4, we draw similar conclusions. However, for subplots with the transmission boundary condition, much of the solution behavior is preserved in regions where $k = 2$ despite the sharp interface contrast. For example, comparing the solution in $\Omega_0$ for **Case 1** and **Case 3** reveals little to no visual difference. This suggests that the level of discretization required to achieve the same accuracy on $\Gamma_0$ for $\beta_0 = 9$ and $\beta_0 = 25$ is more similar than one might expect given the difference in contrast. To determine whether this is indeed the case, we examine the optimal number of boundary points required to achieve a fixed target accuracy.

### 3.3. Observations from Boundary Optimization.
Now we check the number of points on the boundary required to achieve an accuracy of $10^{-6}$ using solutions approximated with
BEM. Beginning with **Case 1** with the sound-hard boundary condition, we determine the required number of discretization points on the boundary, denoted $M_0^{SH}$. $M_0^{SH}$ is obtained by incrementally refining the boundary discretization until the target $L^2$ error, in (13), falls below $10^{-6}$. Our tests indicate that a minimum of $M_0^{SH} = 46$ points is required to reach this accuracy. To estimate the number of points per wavelength $N_{ppw}$ a discretization on the boundary $\Gamma_j$ corresponds to, we use the relation

$$ (14) \qquad M_j = N_{ppw} \cdot R_j \cdot \max(k_j) \quad \text{on } \Gamma_j. $$

We need to determine $N_{ppw}$ for the target accuracy. We obtain an estimate of it for the sound-hard case using:

$$ M_0^{SH} \approx N_{ppw} \cdot R_0 \cdot k_0, $$

which gives $N_{ppw} \approx 6$. While this value is derived for a specific case, we use it as a practical reference point. The goal is not to impose a precise standard for all cases, but to guide expectations for the discretization of the boundary. Applying this to **Case 1** for transmission, when the highest wavenumber is $k_1 = 6$, we expect to need approximately $M_0 = 144$ points to achieve a similar accuracy. However, the optimal value, as determined by iterating through increasing values of $M$, is $M_0^{tar} = 84$.

For **Case 2**, the optimal interface discretization for the transmission boundary condition is $M_0^{tar} = 108$. This value is higher than that obtained in the previous case, which involved a low-to-high wavenumber contrast. This indicates that the exterior wavenumber influences the solution on both sides of the interface, so a higher exterior wavenumber leads to higher resolution demand (and the opposite is also true). However, this result is also lower than what would be predicted by the common guideline for wave propagation problems based only on the highest wavenumber. From $M_0^{SH} = 112$ we estimate $N_{ppw} \approx 5$, which yields that the approximate number of points we should expect to obtain a similar accuracy for the transmission case is $M_0 = 120$. Adopting a more conservative baseline for the same pair of wavenumbers (2, 6) would suggest using our previously obtained $M_0 = 144$. Regardless of which assumption we follow, the observed optimal resolution falls below the estimates. The general rule systematically overestimates the required boundary resolution in the multilayer setting.

To illustrate how accuracy on the boundary impacts the accuracy in the volume, we also determine the solution in the layers and compare the resulting error for different evaluation distances. To do so, we determine the $L^2$ error of the numerical solutions in the layers at various radial distances at radii $r$, both inward and outward from the boundary $\Gamma_i$ (see Figure 5). The boundary solution is evaluated at the radii $R_i$ using the BIE. Then the fields in each layer $u_j$ are evaluated on the boundary-fitted grid using (10). At each sampling radius for fixed

<figure><img src="image_5.png" /><figcaption>**Case 3:** $k = (2, 10)$, $R_0 = 4$, $\beta_0 = 25$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Case 4:** $k = (10, 2)$, $R_0 = 4$, $\beta_0 = \frac{1}{25}$</figcaption></figure>

<caption>**FIGURE 4.** Real part of the total field analytic solutions for a higher contrast case with sharper contrasts $\beta_0$ as given.</caption>
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Schematic of the samplings for $||u_j - u_j^{ana}||_{L^2(r)}$. The boundary solution is computed on the black ring at radius $R_0$, and the error between the off-boundary PTR and analytic solutions is evaluated at each sampling radii (red rings) separately.</figcaption></figure>

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

Radial errors for Case 1: $k = (2, 6)$, $R_0 = 4$, $\beta_0 = 9$

<caption>FIGURE 6. Semilog plot of the absolute $L_2$ error of solutions over radii $r$. The vertical dashed black line marks $R_0$, while the horizontal dashed gray line indicates the target order of error. Dashed segments denote skipped values due to close evaluation [16].</caption>

values of $r$, the numerical solution is compared to the analytic solution (21) on $M_i$ values of $\theta_m$. Specifically, the radial $L^2$ error is computed using (13),

$$ ||u_j - u_j^{ana}||_{L^2(r)} \approx \left( \frac{2\pi}{M_i} \sum_{m=0}^{M_i-1} |u_j(\mathbf{x}_j(\theta_m)) - u_j^{ana}(r, \theta_m)|^2 \right)^{1/2} $$

for $i = j$, or otherwise $i = N - 1$ if within the region $\Omega_N$, which is not defined by an interior boundary. Here, $\mathbf{x}_j(\theta_m) := (r \cos \theta_m, r \sin \theta_m)$.

In Figure 6, we examine the radial error behavior for **Case 1**. This figure shows that our BEM framework delivers accuracy in the volume that is comparable to, or even better than, that obtained on the boundary. This behavior comes from the smoothing effect of the boundary integral representation, combined with the spectral accuracy of the PTR for smooth and periodic interfaces, resulting in highly accurate solutions across the domain [11]. Using the estimate $M_0 = 144$, we achieve machine precision at $10^{-14}$ throughout the domain, rather than the intended target accuracy achieved with the optimal $M_0^{tar} = 84$. While obtaining machine precision is not undesirable, this result indicates that the estimate over-resolves the solution. The optimal number of points to reach machine precision is $M_0^* = 106$, which lies closer to the discretization needed for a lower order accuracy than our estimate $M_0$. This example clearly
demonstrates that only accounting for $k_{max}$ can lead to unnecessary levels of discretization for multilayered media.

The radial errors for **Case 2** are shown in Figure 7. Estimating $M_0$ from the sound-hard boundary discretization using $N_{ppw} \approx 5$ results in a solution that falls below the target error threshold. A more conservative choice of $N_{ppw} = 6$ (which is the average $N_{ppw}$ estimated for every case in Table 1) yields the optimal number of boundary points required for machine precision, $M_0^* = 144$. In either case, the discretization remains overestimated, even though the solution behavior observed in Figure 3 shows little visual difference between sound-hard and transmission boundary conditions for Case 2.

<figure><img src="image_9.png" /><figcaption>Radial errors for **Case 2**: $k = (6, 2)$, $R_0 = 4$, $\beta_0 = \frac{1}{9}$</figcaption></figure>

FIGURE 7. Semilog plot of the absolute $L_2$ error of solutions over radii $r$.

<figure><img src="image_10.png" /><figcaption>Radial errors for **Case 3**: $k = (2, 10)$, $R_0 = 4$, $\beta_0 = 25$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Radial errors for **Case 4**: $k = (10, 2)$, $R_0 = 4$, $\beta_0 = 25$</figcaption></figure>

FIGURE 8. Semilog plot of the absolute $L_2$ error of solutions for the given cases computed at different radii $r$.
To see whether these observations are more pronounced for sharper contrasts, we also plot the radial errors for **Case 3** and **Case 4**. In Figure 8, $M_0^{tar}$ in both cases deviates even further from the value obtained by the general rule, $M_0$, in (14). Even when targeting machine precision, the discretization is substantially overestimated relative to the true optimal value $M_0^*$. For example, in **Case 3**, $M_0 = 240$ is nearly twice $M_0^* = 140$. These results indicate that as the contrast increases, the general rule becomes progressively less effective at predicting the optimal boundary discretization.

Choosing an optimal boundary discretization becomes more complex as we introduce additional interfaces $\Gamma_j$ for the multilayer problem. We present the result for **Case 5** (with 2 interfaces) in Figure 9. We obtain two optimal discretizations on $\Gamma_0$ and $\Gamma_1$ respectively, $M_0^{tar} = 212$ and $M_1^{tar} = 50$, to achieve the chosen accuracy. Although we may attribute the difference in number of points to the varying boundary sizes for $R_0 > R_1$, the ratio $R_0/R_1 = 3$ does not fully account for the fact that $M_0^{tar}$ is over four times $M_1^{tar}$. With this boundary discretization, the $L^2$ error seems uniform across the domain.

These numerical experiments suggest that as long as the interfaces are discretized with a number of points large enough to reach the target accuracy, the solutions within each layer can be captured with comparable or better precision. However, this naive adaptation of the discretization has the prohibitive need for the analytic solutions. Furthermore, identifying a rule to account for multiple contrasts $\beta_j$ and their influence on the interface where they are defined is challenging. It is even more difficult to quantify their effect on other interfaces. We would like to develop an *automatic adaptive* approach, i.e., a method that is applicable to cases where we do not have access to analytic solutions.

In the following section, we introduce an adaptive mesh strategy that optimizes the discretization to minimize the interpolation error. Unlike the boundary-fitted approach it refines the mesh only where it is needed, enabling accurate and efficient computation of the field solutions.

# 4. ADAPTIVE MESH REFINEMENT

## 4.1. From Interpolation Error for Layered Domains to Optimal Mesh Design.

Deriving the optimal volume mesh. We use an automatic mesh adaptation algorithm to determine the optimal discretization for each domain (*Adaptive Method*). Our goal is to determine a mesh (with size *K*) that minimizes the interpolation error, estimated from a truncated Taylor expansion together with geometric considerations. However, solving this minimization problem

<figure><img src="image_12.png" /><figcaption>FIGURE 9. Semilog plot of the absolute $L_2$ error of BIE solutions case for $\mathbf{k} = (8, 2, 6)$, $\mathbf{r} = (6, 2)$. Only the results from the optimal $M^{tar}$ values $\mathbf{M} = (212, 50)$ are shown.</figcaption></figure>
directly at the discrete level is challenging, which motivates the use of an equivalent continuous formulation.

We follow the approach of [25, 26]. An anisotropic adapted mesh $T^*$ is generated with respect to a metric $\mathcal{T}$, so that all elements have unit size in the $\mathcal{T}$-metric. The optimal metric $\mathcal{T}$ is obtained by minimizing the continuous interpolation error for a given mesh complexity $\mathcal{N}$, thereby defining the corresponding adapted mesh $T^*$ and controlling the overall mesh density.

Given a domain $\Omega$ with smooth closed boundary $\Gamma$, the optimal metric $\mathcal{T}$, in the $L^p$ norm (for $p \ge 1$), is given by

$$ (15) \qquad \mathcal{T} = \mathcal{N}^{\frac{2}{3}} \left( \int_{\Gamma} \det\left( |H_R(u)| \right)^{\frac{p}{2p+3}} \right)^{-\frac{2}{3}} \det\left( |H_R(u)| \right)^{\frac{-1}{2p+3}} |H_R(u)| $$

where $H_R(u)$ is the recovered Hessian of the solution. In addition, with this approach, an explicit bound on the interpolation error associated with $T^*$ is available and given by

$$ (16) \qquad E_p(T^*) = \frac{\mathcal{N}^{-1}}{4} \left( \int_{\Gamma} \det\left( |H_R(u)| \right)^{\frac{p}{2(p+1)}} \right)^{\frac{p+1}{p}}. $$

Finally, the corresponding error estimate becomes

$$ \frac{1}{2} E_p(T^*) \le ||u - \Pi_h u||_{L^p(\Omega)} \le 2E_p(T^*) $$

where $\Pi_h u$ is the linear interpolant of the solution $u$ on the mesh $T^*$.

Illustration of the behavior of the Adaptive Method. For each interface $\Gamma_j$, we begin by defining a fine discretization: a dense, equispaced set of candidate points on the boundary. These points serve as part of the initial data for our mesh adaptation process. From there, we refine or coarsen the mesh based on a mesh complexity input requirement. However, a key constraint is imposed: the interface geometry must be preserved. This ensures the interface remains geometrically consistent and is not altered during mesh refinement.

Although we begin with an equispaced set of points on the interfaces, the algorithm does not necessarily preserve that spacing. Thus, the distribution of retained interface vertices may be nonuniform.

As an illustration, Figure 10 presents a representative comparison between a boundary-fitted mesh and its adapted (anisotropic) counterpart for a configuration with three interfaces and varying material contrasts. As expected, the adaptive mesh concentrates points in regions where higher resolution is required, leading to highly nonuniform discretizations along each interface compared with the uniform spacing observed in Figure 10(a).

Automatic mesh adaptation. The adaptive method follows an iterative process: we initialize the algorithm with a coarse target complexity, typically $\mathcal{N}_0 = 200$, and compute the solution in the volume. From this solution, we estimate the Hessian and construct the corresponding optimal metric $\mathcal{T}$ as well as the corresponding estimate of the interpolation error $E_p(T^*)$. A new mesh is then generated from this metric for a chosen $\mathcal{N}$, and this process is repeated.

Rather than increasing the mesh complexity at every iteration, we raise it only intermittently, e.g., by a fixed percentage every other iteration. This strategy offers the possibility to optimize the interpolation error at a given complexity before refining the mesh. At each iteration, the interpolation error associated with the current metric is measured. The relative improvement between iterations serves as an indicator of whether further adaptation at the same complexity is needed. The algorithm may be stopped once the improvement in the interpolation error becomes negligible or when a prescribed error bound estimate or maximum complexity is reached.

The adaptive method is coupled with a two-stage procedure to evaluate solutions. We first find the boundary data by solving the system of corresponding boundary integral equations (8). The solution is then propagated into the volume via boundary integral representation (7).
<figure><img src="image_13.png" /><figcaption>FIGURE 10. Comparison of meshes at initial and later iterations of mesh generation versus the solution for $\mathbf{r} = (6, 4, 2.5)$ and $\mathbf{k} = (2, 8, 4, 2)$.</figcaption></figure>

To utilize boundary element methods, the boundary geometry is preserved by the adaptive method.

A particular constraint of our implementation arises from our choice of quadrature rule used for the boundary integrals: they require a discretization that is equi-spaced, defined at $M_j$ values of $\theta_m$ (9). To satisfy this requirement, at each iteration we extract from the current mesh its closest approximation to a uniform boundary mesh. This is used to evaluate accurate boundary data for the solutions. The volume mesh, however, is adapted strictly according to the optimal metric and is not required to be uniform. It is our chosen quadrature that requires an equi-spaced distribution. In return it yields spectral accuracy in the volume, so as long as boundary solutions are sufficiently accurate, the solution remains highly accurate throughout the domain (as shown in section 3.3). This enables accurate wave field solutions with reduced computational cost for reduced, optimal discretizations.

## 4.2. Illustration of the adaptive process on circular boundaries.
We now regard the mesh adaptation process on a series of representative examples involving circular boundaries (presented in Table 2). Importantly, we want to determine whether a general rule such as (14) is appropriate.

<table><caption>TABLE 2. Summary of multiple-boundary cases used in numerical experiments with the Adaptive Method, with wavenumbers $\mathbf{k}$, boundary radii $\mathbf{r}$, and contrasts $\mathbf{b}$.</caption><thead><tr><th>Case</th><th><b>k</b> = (k<sub>0</sub>, k<sub>1</sub>, ..., k<sub>N</sub>)</th><th><b>r</b> = (R<sub>0</sub>, ..., R<sub>N-1</sub>)</th><th><b>b</b> = (β<sub>0</sub>, ..., β<sub>N-1</sub>)</th></tr></thead><tbody><tr><td>1</td><td>(2, 6)</td><td>4</td><td>9</td></tr><tr><td>5</td><td>(2, 6, 1)</td><td>(4, 2)</td><td>(9, <sup>1</sup>&frasl;<sub>36</sub>)</td></tr><tr><td>6</td><td>(2, 6, 10)</td><td>(4, 2)</td><td>(9, <sup>25</sup>&frasl;<sub>9</sub>)</td></tr><tr><td>7</td><td>(6, 2, 4, 1)</td><td>(4, 2, 1)</td><td>(<sup>1</sup>&frasl;<sub>9</sub>, 4, <sup>1</sup>&frasl;<sub>16</sub>)</td></tr></tbody></table>

We first consider **Case 1**, for which we already gained some knowledge in previous section.

Recall that the adaptive method relies on an iterative process. To obtain the results presented in Figure 11, we fix the initial complexity across the complete domain at $N_0 = 400$. The automatic mesh adaptation process undergoes 14 iterations, where we increase the complexity $N_n$ by 30% every other iteration and otherwise re-mesh on the same complexity once. As the mesh refines, the number of boundary points $M_0$ changes.
We have an adapted mesh, but we look at the accuracy of the BIE solutions by computing the corresponding absolute $L^2$ error using the known analytic solution as in (13) at each iteration. Figure 11 shows the evolution of this error as a function of the mesh complexity $\mathcal{N}$ and the number of boundary points $M_0$.

Examining Figure 11(a), we observe that the $L^2$ error is highly sensitive to the choice of $M_0$. In particular, looking at where the discretization is refined by approximately 10 additional points, increasing $M_0$ from 78 to 88 decreases the error by roughly one order of magnitude, while increasing $M_0$ to $98 \approx 100$ results in a reduction of several orders. This highlights that even slight deviations from the optimal number of points, particularly under-resolving, can have substantial impacts on the BIE solution's accuracy.

We recall that our goal is to verify the accuracy of the general rule. Figure 11(b) plots the error versus the ratio $\frac{M_0}{2\pi R_0}$, which is proportional to the points-per-wavelength metric $N_{\text{ppw}}$ from (14). The results show that this ratio is not constant; if it were, it would provide a simple guideline for selecting the optimal boundary discretization. To achieve a target accuracy of $10^{-6}$, the adaptive mesh algorithm reaches an optimal resolution with only $\frac{M_0}{2\pi R_0} \approx 3$, well below the typical rule of 5–10 points per wavelength for $k_{\text{max}} = 6$. This can be explained by the lower exterior wavenumber $k_0 = 2$, which dominates the mesh requirements and reduces the number of boundary points needed for $k_1 = 6$. Even for errors on the order of $10^{-14}$, the ratio remains below the standard 5–10 points-per-wavelength guideline.

After this introductory example, we turn to more complex cases with additional circular boundaries, as given in Table 2.

Figure 12 presents results for the 2-boundary cases from Table 2 (Case 5 and Case 6). The Adaptive Method achieves a comparable order of accuracy on both boundaries, $\Gamma_0$ and $\Gamma_1$, with different discretizations. Again, in both cases there is no clear trend when comparing to the ratios $\frac{M_i}{2\pi R_i}$ ($i = 0, 1$), reinforcing the absence of a general rule for optimal boundary discretization in multilayer transmission problems. Increasing $k_{\text{max}}$ from 6 to 10 in **Case 6** results in only slightly increased discretization demand on $\Gamma_1$ while it remains relatively unchanged on $\Gamma_0$. The discretization range that might otherwise be assumed from the general rule (14) for $5 \le N_{\text{ppw}} \le 10$, given $k_{\text{max}} = 10$, is between $M = 200$ to 400.

Figure 13 shows that similar trends are observed for **Case 7**, with three boundaries. Additionally, though a maximum wavenumber of $k_0 = 6$ is introduced in the outermost subdomain,

<figure><img src="image_14.png" /><figcaption>FIGURE 11. For **Case 1**: $k = (2, 6)$ and $R_0 = 4$: (a) Log of the absolute $L_2$ error (13) obtained with the BIE solution on $\Gamma_0$ versus the mesh complexity $\mathcal{N}$, for various boundary resolutions $M_0$ as labeled. (b) The same error plotted against $\frac{M_0}{2\pi R_0}$.</figcaption></figure>
the presence of subdomains with lower wavenumbers substantially reduces the overall discretization requirements on the interior interfaces. The transmission conditions couple the interface behavior in such a way that improved accuracy on $\Gamma_1$ and $\Gamma_2$ is obtained with little additional discretization effort, except on $\Gamma_0$. Notably, increasing the discretization on $\Gamma_0$ to improve its boundary accuracy also improves the accuracy on the other boundaries.

From these circular boundary examples examined, the Adaptive Method seems to maintain similar accuracy across all boundaries by concentrating resolution where needed, often reducing mesh requirements due to lower wavenumbers in sub-regions $\Omega_j$. Although no simple rule applies universally, it reliably identifies boundary discretizations that preserve uniform accuracy, even in multilayer configurations with high material contrasts.

Having validated the Adaptive Method on simple geometries where an analytical solution can be found, we now turn to more complex boundaries. In the next section, we assess its performance on exotic shapes, where analytic solutions are unavailable and local refinement, curvature, and material contrasts pose greater challenges.

### 4.3. Adaptive Mesh Algorithm with BIE Solutions.
In general, one does not have access to an analytic solution to investigate how the Adaptive Method is performing. Since the

<figure><img src="image_15.png" /><figcaption>Case 5: $k = (2, 6, 1)$, $r = (4, 2)$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Case 6: $k = (2, 6, 10)$, $r = (4, 2)$</figcaption></figure>

<caption>FIGURE 12. (a) Log absolute $L_2$ error of the BIE solutions on $\Gamma_0$ and $\Gamma_1$ versus total mesh complexity $N$, for various boundary resolutions $M_0$ and $M_1$ as labeled. (b) The same error plotted against the ratio $\frac{M_i}{2\pi R_i}$.</caption>
volumetric interpolation error bound (16) should consistently improve over increasing mesh complexities if the algorithm is performing well, we simply rely on computing (16) using the approximated solution, $u$. Here we make use of the BIE solutions and the boundary integral representation to obtain the recovered Hessian $H_R(u)$.

In order to validate our approach, we first test on circular boundaries so we can compare with previous results. We call **ADAPT-ANA** the error bound in terms of the recovered Hessian $H_R(u^{ana})$ which defines the mesh $T_{ana}^*$, and **ADAPT-BIE** the one where the recovered Hessian is determined from the BIE solution, defining the mesh $T^*$. We summarize the two approaches in Table 3.

Figure 14 presents the interpolation error bounds computed with both methods for mesh complexities increasing from $10^3$ to $10^5$ for **Case 1** and **Case 2**. We find linear convergence with **ADAPT-ANA**, which is in accordance with theoretical expectations [23]. We observe convergence with **ADAPT-BIE**, however with a slower rate for larger mesh complexities. This behavior is likely caused by increased refinement near the boundaries, where the evaluation of the boundary potentials becomes challenging due to nearly-singular kernels. This *close-evaluation* error is well-known and can impact the convergence behaviors of standard quadrature rules [16]. This highlights an important research direction when combining adaptive meshing with BEM solvers: the integration of close-evaluation treatment methods.

The first two error estimates give an upper bound but we would like to check if we can rely on it. Here, we also estimate the error in the volume by comparing the BIE solution with the analytic solution: this is what we call **BIE** in Figure 14. This error is computed via the $L^2$

<figure><img src="image_17.png" /><figcaption>FIGURE 13. For **Case 7**: $k = (6, 2, 4, 1)$ and $r = (4, 2, 1)$: (a) Log absolute $L^2$ error of the BIE solution on three boundaries versus total mesh complexity $N$, for various boundary resolutions $M_0, M_1$, and $M_2$ as labeled. (b) The same errors plotted against the ratio $\frac{M_i}{2\pi R_i}$.</figcaption></figure>

<table><caption>TABLE 3. Description of methods plotted in Figure 14.</caption><thead><tr><td>Method</td><td>Description</td><td>Computation</td></tr></thead><tbody><tr><td>ADAPT-ANA</td><td>Interpolation error bound for analytic solution on adapted mesh <b>T</b><sup>*</sup><sub>ana</sub></td><td>E<sub>2</sub>(<b>T</b><sup>*</sup><sub>ana</sub>)</td></tr><tr><td>ADAPT-BIE</td><td>Interpolation error bound for BIE solution on adapted mesh <b>T</b><sup>*</sup></td><td>E<sub>2</sub>(<b>T</b><sup>*</sup>)</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>FIGURE 14. Plot of the error bounds described in Table 3 for Case 1 and 2, over increasing mesh complexities $\mathcal{N}$.</figcaption></figure>

error

$$ (17) \qquad \|u - u^{ana}\|_{L^2(\tilde{\Omega})} = \sum_j \|u_j - u_j^{ana}\|_{L^2(\Omega_j \setminus \omega_j)} $$

where $\tilde{\Omega}$ is the “perforated" domain $\Omega$, where we remove interfaces vicinities (called $\omega_j$) to avoid the so-called close evaluation problem. This yields a more accurate measure of the accuracy of the solution in the volume as we do not treat the nearly-singular behaviors. Figure 14 shows spectral convergence of **BIE**, as expected theoretically.

For **Case 1** in particular, we see that the error in the volume begins to decrease exponentially at higher orders of complexity $\mathcal{N}$. This can likely be attributed to the increasing number of nodes discarded due to their proximity to the boundaries. At very high complexities, the Adaptive Method places additional nodes near the boundaries to resolve the jumps in the solution magnitude caused by the nearly-singular behavior in those regions. It is more prominent for the case when $k_0 = 2$ because the lower exterior wavenumber drives lower discretization demand to capture the overall solution behavior in the domain, so the adaptive method focuses on close evaluation behavior near the boundary. Figure 15 compares optimal meshes (according to the mesh adaptation algorithm) from ADAPT-BIE versus ADAPT-ANA. In Figure 15(a), the mesh densities are similar and relatively uniform through each region $\Omega_0$ and $\Omega_1$ in the volume. In Figure 15(b), for a higher complexity, the ADAPT-BIE mesh shows increased node concentration near the boundary, while the ADAPT-ANA mesh adds more nodes throughout the whole interior region $\Omega_1$. This near-boundary refinement difference only becomes significant at sufficiently high complexities relative to the domain size. In this example, it does not appear until $\mathcal{N} > 10^4$. Therefore, our results should remain optimal when using BIE solutions for the cases examined so far, and subsequent cases, provided the mesh complexity stays below this threshold and the relative domain size is not reduced.

Having established the interpolation error bound may tell us the effectiveness of the adaptive process where we do not have reference solutions, we now turn to more complex geometries to assess the robustness of the mesh adaptation strategy.

We use the following parameterization to model star shape objects:

$$ (18) \qquad \begin{aligned} x_i(\theta) &= (R_i + a \cos(n\theta)) \cos\theta, \\ y_i(\theta) &= (R_i + a \cos(n\theta)) \sin\theta, \end{aligned} \qquad i = 0, \dots, N-1, $$

where $n$ is the number of star points and $a \ge 0$ controls their relative size with respect to $R_i$.
We can leverage the error bound associated with the BIE solutions to determine whether additional re-meshing iterations are necessary, and how many. If the error bound continues to improve with each iteration, the current number of re-meshing iterations is likely sufficient. Conversely, if the error bound plateaus or degrades, this signals the need to adjust the number of re-meshing iterations to achieve better results.

We examine two cases for amplitudes $a = 0.05$ and $a = 0.1$ with $n = 10$, to compare the behavior of a slight versus more pronounced star point amplitude under two meshing strategies: no re-meshing versus two re-meshings, as shown in Figures 16 and 17. We observe that fewer re-meshing iterations tend to yield better convergence of the error bound.

For these cases, we also consider that the adaptive strategy may prioritize regions of high curvature and sharper geometries, allocating more points in these regions where interpolation error is likely to be higher regardless of the contrasts. In particular, cases with high exterior wavenumbers consistently demand finer resolution across all boundaries, while otherwise lower exterior wavenumbers lessen these demands. Since we cannot quantify the error on the boundaries directly in the absence of analytic solutions, we rely on visual inspection of the corresponding BIE solutions to comment on whether the meshes produced are sufficient to capture the key solution features. Typically they do, though there is improvements to be made to leverage the same optimality of results we achieve from analytic solutions. Treatment of close evaluation error near the boundaries while processing solutions, for example, will improve usage of this adaptive algorithm by improving the accuracy of the method. It will also allow

<figure><img src="image_19.png" /><figcaption>FIGURE 15. Comparison of ADAPT-BIE and ADAPT-ANA meshes for varying mesh complexities $\mathcal{N}$, with $\mathbf{k} = (2, 6)$ and $R_0 = 4$ after (a) 5 iterations versus (b) 11 iterations of adaptation.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIGURE 16. Error bounds and reference slopes of the interpolation versus complexity $\mathcal{N}$ for $\mathbf{k} = (2, 6)$, $R_0 = 4$, $a = 0.05$, $n = 10$ for different re-meshing schemes.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 17. Error bounds and reference slopes of the interpolation versus complexity $\mathcal{N}$ for $\mathbf{k} = (2, 6)$, $R_0 = 4$, $a = 0.1$, $n = 10$ for different re-meshing schemes.</figcaption></figure>

re-meshing to be used as intended, to optimally redistribute the mesh points to capture the actual solution behavior rather than try to resolve problematic solution features. With treated solutions, the adaptive method is expected to limit unnecessary over-resolution near the boundaries. Although some refinement is unavoidable due to boundary curvature in the cases shown in this section, untreated solutions likely amplify this effect, as shown in Figure 18.

Our qualitative observations reinforce the key insight from the previous section: there is no simple rule linking the optimal number of points $M_i$ on $\Gamma_i$ to the maximum wavenumber. Instead, the discretization depends on both the local geometry and material properties within the overall multilayered structure. However, as the boundary shapes become more complex, the algorithm continues to produce stable and visually reasonable discretizations. In particular, the background medium, $k_0$, strongly influences the discretization behavior. This complexity highlights the value of automatic boundary discretization approaches in multilayer transmission problems.
<figure><img src="image_22.png" /><figcaption>FIGURE 18. Real part of the PTR solutions and generated meshes for $\mathbf{k} = (2, 6)$, $R_0 = 4$, $a = 0.1$, $n = 10$. We note for the same mesh complexity $N = 2000$, after re-meshing from (a) the new mesh generated in (b) adds more refinement to the boundary.</figcaption></figure>

## 5. CONCLUSIONS AND PERSPECTIVES

This work has shown that achieving optimal interface discretizations in multilayer transmission problems requires mesh adaptivity, as no simple linear rule—whether based on the contrasts $\beta_j$ or the highest wavenumber—can reliably determine optimal discretizations. Such adaptivity proves essential to obtain the best cost-accuracy balance, particularly due to the interaction between high- and low-wavenumber regimes across material interfaces.

In particular, we have observed that the exterior wavenumber has a significant influence on the overall solution behavior. A very low exterior wavenumber can damp oscillations within the interior media, whereas a high wavenumber requires a correspondingly finer resolution [14]. Because the interfaces often differ not only in physical behavior, size, and shape, using the same number of boundary elements or scaling all interface meshes in proportion to the boundary with the highest wavenumber is not ideal. Such uniform or proportional discretizations increase computational cost and lead to over-refinement of subdomains with less complex solution behavior.

Another challenge in BEM formulations for multilayered media is the strong coupling between interfaces, which causes errors on one boundary to propagate to others. For computational efficiency, it is therefore crucial to achieve a comparable accuracy on all interfaces: over-resolving one boundary while others remain coarse does not improve the global error and only increases the computational cost. The adaptive algorithm employed in this work addresses this issue by iteratively refining the boundary meshes based on interpolation error indicators. In this way, the system automatically identifies where refinement is needed, leading to a balanced distribution of error across interfaces and an overall more efficient discretization.

Our automatic mesh adaptation algorithm for multilayer problems has been shown to be optimal for simple configurations where analytical solutions are available. However, our experiments also reveal a critical point: the need for highly accurate evaluation of nearly singular integrals. These integrals receive relatively little attention outside the BEM community as the solution is usually not required with high precision near boundaries. For multilayered media, particularly when we may want to model metamaterials where high contrasts are prevalent, it is important to address but they are often treated as post-processing options. In the context of mesh adaptation, however, inaccuracies in these integrals can dominate the local error
indicators near interfaces, misleading the refinement process and causing unnecessary mesh enrichment driven by quadrature errors rather than by the solution dynamics. Several methods for handling near-singular integrals exist in the literature, but we have not yet incorporated them into our implementation.

Though the results discussed here are in 2D, we expect that for multilayered media in 3D the reduction of computational costs would be more pronounced and the potential benefit even greater. For BEM, our focus on the interfaces is natural, but this approach and our conclusions could be useful for other mesh-based methods for wave problems as well, including FEM or FD methods. To sum up, less is more: the solution itself should dictate the discretization; it is not needed to refine everywhere for optimal accuracy.

## APPENDIX A. ANALYTIC SOLUTION FOR MULTILAYER TRANSMISSION PROBLEMS

The analytic solution for the multilayer transmission problem is available in the case of concentric, circular interfaces [27]. In these cases, the boundary $Γ_j$ for $j = 0, \dots, N-1$, is parameterized by $\mathbf{x}_j(\theta) = (R_j \cos \theta, R_j \sin \theta)$ with $R_N < R_{N-1} < \dots < R_1 < R_0$. The incident field, for the choice $\mathbf{a} = (0, 1)$, may also be expressed in terms of $r$ and $\theta$ such that $u^{in}(y) = u^{in}(r, \theta) = e^{ik_0 r \sin \theta}$.

By the Jacobi-Anger expansion [28],

$$ (19) \qquad e^{iz \sin \theta} = \sum_{m=0}^{\infty} J_m(z) e^{im\theta} $$

the expression of the incident field may be given by

$$ u^{in}(r, \theta) = e^{ik_0 r \sin \theta} = \sum_{m=0}^{\infty} J_m(k_0 r) e^{im\theta} $$

where $J_m$ is the $m$th order Bessel function. In $\Omega_0$, the scattered field must satisfy both the Helmholtz equation (1) and the Sommerfeld radiation condition (3). It takes the form of the following ansatz

$$ (20) \qquad u_0^{sc}(r, \theta) = \sum_{m=0}^{\infty} A_m H_m^{(1)}(k_0 r) e^{im\theta} $$

where $H_m^{(1)}$ are $m$th order Hankel functions of the first kind, which represent outgoing wave behavior, and $A_m \in \mathbb{C}$ are coefficients to determine. Similarly, the behavior of interior fields can be captured by expansions in terms of $J_m$ and $H_m^{(1)}$ in the regions $\Omega_j$ for $j = 1, \dots, N$ [24]. Altogether, we can construct the analytic expression for the total field solution $u_j^{ana}$ in each subdomain $\Omega_j$ as:

$$ u_0^{ana}(r, \theta) = \sum_{m=0}^{\infty} [A_m H_m^{(1)}(k_0 r) + J_m(k_0 r_0)] e^{im\theta} \quad R_0 \le r, \quad \forall \theta, $$

$$ (21) \qquad u_j^{ana}(r, \theta) = \sum_{m=0}^{\infty} [B_{j,m} H_m^{(1)}(k_j r) + C_{j,m} J_m(k_j r)] e^{im\theta} \quad R_{j-1} \le r \le R_j, \quad \forall \theta, $$

$$ u_N^{ana}(r, \theta) = \sum_{m=0}^{\infty} E_m J_m(k_N r) e^{im\theta} \qquad r \le R_N, \quad \forall \theta, $$

with $B_{j,m}$, $C_{j,m}$, $E_m \in \mathbb{C}$. By enforcing the transmission conditions (2) at each interface $Γ_j$ between regions $\Omega_j$ and $\Omega_{j+1}$, we yield a solvable linear system for the unknown coefficients $A_m$, $B_{j,m}$, $C_{j,m}$, and $E_m$. The unknown coefficients are then used in (21) to determine the solution $u_j^{ana}$ anywhere in the domain.
APPENDIX B. ANALYTIC SOLUTIONS FOR THE SOUND-HARD BOUNDARY CONDITION

We consider here the sound-hard scattering problem. For the one-boundary problem, the solutions $u_0^{ana}$ and $u_1^{ana}$ must satisfy the Helmholtz equation as defined in (1) in the regions $\Omega_0 \cup \Omega_1$. The sound-hard boundary condition on the interface $\Gamma_0$ defined at radius $R_0$ is $\frac{\partial u_0}{\partial n_0} = 0$.

As similarly stated in (21), the solutions in $\Omega_0$ and $\Omega_1$ respectively can be expressed as:

$$ \begin{cases} u_0^{ana}(r, \theta) = \sum_{m=0}^{\infty} \left[ A_m H_m^{(1)}(k_0 r) + J_m(k_0 r) \right] e^{im\theta} \\ u_1^{ana}(r, \theta) = 0 \end{cases} $$

The coefficient $A_m$ must satisfy:

$$ A_m k_0 H_m^{(1)}(k_0 R_0) + k_0 J_m'(k_0 R_0) = 0. $$

Thus the coefficient $A_m$ can be found as:

$$ A_m = \frac{-J_m'(k_0 R_0)}{H_m^{(1)}(k_0 R_0)} $$

Finally, the analytic solution for the sound-hard boundary case for a single boundary is given as:

$$ \begin{cases} u_0^{ana}(r, \theta) = \sum_{m=0}^{\infty} \left[ J_m(k_0 r) - \frac{J_m'(k_0 R_0)}{H_m^{(1)}(k_0 R_0)} H_m^{(1)}(k_0 r) \right] e^{im\theta} \\ u_1^{ana}(r, \theta) = 0 \end{cases} $$

REFERENCES

[1] J. Bertolotti, E. G. van Putten, C. Blum, A. Lagendijk, W. L. Vos, A. P. Mosk, Non-invasive imaging through opaque scattering layers, Nature 491 (7423) (2012) 232–234.

[2] J. P. Elisee, Innovative boundary integral and hybrid methods for diffuse optical imaging.

[3] J. Cai, M. Deng, A methodology for evaluating damage in thin-walled cylindrical shells using bounded ultrasonic beam scattering 204 112352.

[4] S. Caorsi, A. Massa, M. Pastorino, M. Raffetto, A. Randazzo, Microwave imaging of cylindrical inhomogeneities based on an analytical forward solver and multiple illuminations, in: 2004 IEEE International Workshop on Imaging Systems and Techniques (IST), pp. 100–105.

[5] B. Wood, Metamaterials and invisibility 10 (5) 379–390.

[6] J. B. Pendry, D. Schurig, D. R. Smith, Controlling Electromagnetic Fields 312 (5781) 1780–1782.

[7] F. Ihlenburg (Ed.), Finite Element Error Analysis and Control for Helmholtz Problems, Springer New York, New York, NY, 1998, pp. 101–188.

[8] Q. Feng, B. Han, M. Michelle, Sixth order compact finite difference method for 2d helmholtz equations with singular sources and reduced pollution effect (2023). arXiv:2112.07154.

[9] F. Nataf, E. Parolin, Coarse spaces for non-symmetric two-level preconditioners based on local extended generalized eigenproblems (2025). arXiv:2404.02758.

[10] M. Bonnet, Boundary Integral Equation Methods for Solids and Fluids 34 (4) 301–302.

[11] R. Kress, Boundary integral equations in time-harmonic acoustic scattering, Mathematical and Computer Modelling 15 (3) (1991) 229–243.

[12] J. Galkowski, E. A. Spence, Does the helmholtz boundary element method suffer from the pollution effect?, SIAM Review 65 (3) (2023) 806–828.

[13] X. Chen, Y. Liu, Thermal Stress Analysis of Multi-layer Thin Films and Coatings by an Advanced Boundary Element Method, CMES - Computer Modeling in Engineering and Sciences 2 (3) (2001) 337–349.

[14] Y. Liu, On the BEM for acoustic wave problems, Engineering Analysis with Boundary Elements 107 (2019) 53–62.

[15] A. Bendali, M. Fares, Boundary Integral Equations Methods in Acoustic Scattering, in: F. Magoulès (Ed.), Computational Science, Engineering & Technology Series, Vol. 18, Saxe-Coburg Publications, Stirlingshire, UK, 2008, pp. 1–36. doi:10.4203/csets.18.1.

[16] A. H. Barnett, Evaluation of Layer Potentials Close to the Boundary for Laplace and Helmholtz Problems on Analytic Planar Domains, SIAM J. Sci. Comput. 36 (2) (2014) A427–A451, publisher: Society for Industrial and Applied Mathematics.
[17] A. Klöckner, A. Barnett, L. Greengard, M. O'Neil, Quadrature by expansion: A new method for the evaluation of layer potentials, Journal of Computational Physics 252 (2013) 332–349.

[18] C. Carvalho, S. Khatri, A. D. Kim, Asymptotic approximations for the close evaluation of double-layer potentials, SIAM Journal on Scientific Computing 42 (1) (2020) A504–A533.

[19] E. A. Cortes, C. Carvalho, S. Chaillat, C. Tsogka, Fast and accurate boundary integral equation methods for the multi-layer transmission problem, in: WAVES 2024 - The 16th International Conference on Mathematical and Numerical Aspects of Wave Propagation, Edmond.

[20] W. A. Wartman, G. N. Ponasso, Z. Qi, J. Haueisen, B. Maess, T. R. Knösche, K. Weise, G. M. Noetscher, T. Raij, S. N. Makaroff, Fast EEG/MEG BEM-based forward problem solution for high-resolution head models 306 120998.

[21] W. A. Wartman, K. Weise, M. Rachh, L. Morales, Z.-D. Deng, A. Nummenmaa, S. N. Makaroff, An adaptive h-refinement method for the boundary element fast multipole method for quasi-static electromagnetic modeling 69 (5) 055030.

[22] J.-M. Mirebeau, Adaptive and anisotropic finite element approximation: Theory and algorithms (Mar. 2011). arXiv:1101.1555.

[23] S. Chaillat, S. P. Groth, A. Loseille, Metric-based anisotropic mesh adaptation for 3D acoustic boundary element methods, Journal of Computational Physics 372 (2018) 473–499.

[24] D. Colton, R. Kress, The Helmholtz Equation, in: D. Colton, R. Kress (Eds.), Inverse Acoustic and Electromagnetic Scattering Theory, Springer International Publishing, Cham, 2019, pp. 15–42.

[25] A. Loseille, F. Alauzet, Continuous Mesh Framework Part II: Validations and Applications.

[26] A. Loseille, Chapter 10 - Unstructured Mesh Generation and Adaptation, in: R. Abgrall, C.-W. Shu (Eds.), Handbook of Numerical Analysis, Vol. 18 of Handbook of Numerical Methods for Hyperbolic Problems, Elsevier, 2017, pp. 263–302.

[27] L. L. Marrec, Investigation ultrasonore, qualitative et quantitative, d'objets à fort contraste, thesis, Aix-Marseille 2 (Jan. 2004).

[28] NIST Digital Library of Mathematical Functions, https://dlmf.nist.gov/, Release 1.2.4 of 2025-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.

INSA LYON, CNRS, ECOLE CENTRALE DE LYON, UNIVERSITE CLAUDE BERNARD LYON 1, UNIVERSITÉ JEAN MONNET, ICJ UMR5208, 69621 VILLEURBANNE, FRANCE.

Email address: camille.carvalho@insa-lyon.fr

LABORATOIRE POEMS, CNRS-ENSTA-INRIA, INSTITUT POLYTECHNIQUE DE PARIS, 828 BD DES MARECHAUX, 91129 PALAISEAU, FRANCE.

Email address: stephnie.chaillat@ensta.fr

UNIVERSITY OF CALIFORNIA, MERCED, 5200 LAKE RD, MERCED CA 95343, USA.

Email address: ecortes7@ucmerced.edu

UNIVERSITY OF CALIFORNIA, MERCED, 5200 LAKE RD, MERCED CA 95343, USA.

Email address: ctsogka@ucmerced.edu