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

HAL
open science

# Reconstruction of smooth shape defects in waveguides using locally resonant frequencies

Angèle Niclas, Laurent Seppecher

► To cite this version:

Angèle Niclas, Laurent Seppecher. Reconstruction of smooth shape defects in waveguides using locally resonant frequencies. Inverse Problems, 2023, 39 (5), pp.055006. 10.1088/1361-6420/acc7c0. hal-03694665v2

HAL Id: hal-03694665

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

Submitted on 11 Apr 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
# Reconstruction of smooth shape defects in waveguides using locally resonant frequencies

Angèle Niclas¹,\* and Laurent Seppecher²

¹CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

²Institut Camille Jordan, École Centrale Lyon, France

\*Corresponding author: angele.niclas@ec-lyon.fr

## Abstract

This article aims to present a new method to reconstruct slowly varying width defects in 2D waveguides using locally resonant frequencies. At these frequencies, locally resonant modes propagate in the waveguide under the form of Airy functions depending on a parameter called the locally resonant point. In this particular point, the local width of the waveguide is known and its location can be recovered from boundary measurements of the wavefield. Using the same process for different frequencies, we produce a good approximation of the width in all the waveguide. Given multi-frequency measurements taken at the surface of the waveguide, we provide a $L^\infty$-stable explicit method to reconstruct the width of the waveguide. We finally validate our method on numerical data, and we discuss its applications and limits.

## 1 Introduction

This article presents a new method to reconstruct width variations of a slowly varying waveguide from multi-frequency one side boundary measurements in dimension 2. The considered varying waveguide is described by

$$ \Omega := \{(x, y) \in \mathbb{R}^2 | 0 < y < h(x)\}, \quad (1) $$

where $h \in C^2(\mathbb{R}) \cap W^{2,\infty}(\mathbb{R})$ is a positive profile function defining the top boundary. The bottom boundary is assumed to be flat (see an illustration in Figure 1) but a similar analysis could be done when both boundaries are varying. In the time-harmonic regime, the wavefield $u$ satisfies the Helmholtz equation with Neumann boundary conditions

$$ \left\{ \begin{array}{ll} \Delta u + k^2 u = -f & \text{in } \Omega, \\ \partial_{\nu} u = b & \text{on } \partial\Omega, \end{array} \right. \quad (2) $$

where $k \in (0, +\infty)$ is the frequency, $f$ is an interior source term, and $b$ is a boundary source term. In this work, a waveguide is said to be slowly varying when there exists a small parameter $\eta > 0$ such that $\|h'\|_{L^\infty(\mathbb{R})} \le \eta$ and $\|h''\|_{L^\infty(\mathbb{R})} \le \eta^2$. Such waveguides are good models of ducts, corroded pipes, or metal plates (see [16, 19]).
We focus in this work on the recovery of the function $h$ modeling the waveguide shape from the knowledge of the wavefield $d^{\text{ex}}(x) := u(x, 0)$ on one surface of the waveguide and for multiple frequencies $k$. This model and inverse problem is inspired from non destructive monitoring of plates done in [5, 18, 13]. Hence we assume the knowledge of measurements of $u(x, y)$ for all frequencies in $K \subset (0, +\infty)$ and for $(x, y) \in I \times \{0\}$ where $I$ is an interval of $\mathbb{R}$ as shown in Figure 1.

If $k$ is chosen such that $k = n\pi/h(x_k^*)$ with $n \in \mathbb{N}$ and $x_k^* \in \mathbb{R}$, the Helmholtz problem is not well posed in general. Nevertheless, we prove in [10] that there exists a unique solution to this problem as long as the waveguide is slowly varying. In the same work, we also give a suitable explicit approximation of the wavefield that explicitly depends on $x_k^*$. The aim in this article to recover the position of $x_k^*$ for different frequencies, and then to recover the shape function $h$. Using these frequencies, the proposed inverse problem is highly non linear but a unique and stable recovering of $h$ is possible up to a controllable approximation error.

<figure><img src="image_3.png" /><figcaption>Figure 1: Parametrization of a slowly variable waveguide of width $h$. A wavefield $u$ is generated by an internal source $f$ and/or a boundary source $b$. Squares represent measurements of $u$ taken on the surface $y = 0$.</figcaption></figure>

## 1.1 Scientific context

The detection and reconstruction of shape defects in a waveguide are mentioned in different works. In articles [20, 2, 1], the authors use a conformal mapping to map the geometry of the perturbed waveguide to that of a regular waveguide. This method is beneficial to understand the propagation of waves in irregular waveguides. Still, it is not easily adaptable to the inverse problem and the reconstruction of defects since the transformation to a regular waveguide is not explicit and proves numerically expensive. Another approach based on the scattering field treatment is developed in [22]. Finally, articles [23, 15] study the forward problem and give leads on how to use one side boundary measurements to reconstruct the width of a slowly varying waveguide.

Our work concerns a different approach, also used in [7, 6], where we assume the data to be available for a whole interval of frequencies. This provides additional information that should help to localize and reconstruct the shape of the defect. Moreover, the use of multi-frequency data often provides uniqueness of the reconstruction (see [4]) and a better stability (see [8, 17, 24]).

In a previous work [11], we already presented a method to reconstruct small width variations using an inverse scattering procedure. However, we avoided all the cut-off frequencies of the waveguide, which are frequencies $k > 0$ such that $k = n\pi/h(x_k^*)$ for a mode $n \in \mathbb{N}$ and a longitudinal position $x_k^*$. Since experimental works presented in [5, 13] suggest that these resonant frequencies are helpful to reconstruct width variations, we choose in this article to work only with these frequencies. Using the study of the forward problem already done in [10], we know that if $k$ is a locally resonant frequency, the wavefield $u$ strongly depends on $x_k^*$. The main idea of our reconstruction method is to use measurements of $u$ to find back $x_k^*$. Since $h(x_k^*) = n\pi/k$, it then gives up the information
about the width of the waveguide in one point. By taking different locally resonant frequencies $k$, we show that one can obtain a complete approximation of the width $h$ of the waveguide.

## 1.2 Outline of the paper

The key result of this paper is Theorem 2, which proves that $u$ is close to a three parameters Airy function, one parameter being $x_k^*$. As explained in Proposition 3, it enables us to find the value of $x_k^*$ for every locally resonant frequency, and to prove that our reconstruction method is $L^\infty$-stable. The paper is organized as follows. In section 2, we briefly recall results on the modal decomposition and the study of the forward problem. In section 3, we study the inverse problem with measurements taken at the surface of the waveguide and we provide a stability result for the reconstruction of the width of the waveguide. Finally, in section 4, we provide numerical reconstruction of different width defects.

## 1.3 Notations

The varying waveguide is denoted by $\Omega$, its boundary by $\partial\Omega$ and the subscript “top” (resp. “bot”) indicates the upper boundary of the waveguide (resp. lower). We denote $\nu$ the outer normal unit vector. For every $\ell > 0$, we set $\Omega_\ell = \{(x, y) \in \Omega \mid |x| < \ell\}$. Spaces $H^1$, $H^2$, $W^{1,1}$, $H^{1/2}$ over $\Omega$ or $\mathbb{R}$ are classic Sobolev spaces. The Airy function of the first kind (resp. second kind) is denoted by $\mathcal{A}$ (resp. $\mathcal{B}$). They are linear independent solutions of the Airy equation $y'' - xy = 0$ (see [3] for more results about Airy functions). See in Figure 2 the graph of these two functions. The term $\delta_{x=s}$ denotes the Dirac distribution at the point $s \in \mathbb{R}$ and the function $1_E$ is the indicator function of the set $E$. Finally, the notation $\{a : b : \ell\}$ designates the uniform discretization of the interval $[a, b]$ with $\ell$ points.

<figure><img src="image_4.png" /><figcaption>Figure 2: Representation of the Airy functions $\mathcal{A}$ and $\mathcal{B}$.</figcaption></figure>

# 2 Brief study of the forward problem

Before studying the inverse problem associated with the reconstruction of the width in a varying waveguide, we need to study the forward problem in order to find an approximation of the available data. In this section, we briefly recall all the main results on the study of the forward problem. These results and their proofs can be found in [11, 12].
A useful tool when working in waveguides is the modal decomposition. The following definition provides a modal decomposition in varying waveguides:

**Definition 1.** We define the sequence of functions $(\varphi_n)_{n \in \mathbb{N}}$ by

$$ \forall (x, y) \in \Omega, \quad \varphi_n(x, y) := \begin{cases} 1/\sqrt{h(x)} & \text{if } n = 0, \\ \frac{\sqrt{2}}{\sqrt{h(x)}} \cos\left(\frac{n\pi y}{h(x)}\right) & \text{if } n \ge 1, \end{cases} \quad (3) $$

which for any fixed $x \in \mathbb{R}$ defines an orthonormal basis of $L^2(0, h(x))$ (see [25, th 5.11])

Hence, a solution $u \in H_{\text{loc}}^2(\Omega)$ of (2) admits a unique modal decomposition

$$ u(x, y) = \sum_{n \in \mathbb{N}} u_n(x)\varphi_n(x, y) \quad \text{where} \quad u_n(x) := \int_0^{h(x)} u(x, y)\varphi_n(x, y)dy. \quad (4) $$

Note that $u_n$ does not satisfy in general any nice equation. However, when $h$ is constant (outside of supp($h'$)), it satisfies an equation of the form $u_n'' + k_n^2 u_n = -g_n$ where $k_n^2 = k^2 - n^2\pi^2/h^2$ is the wavenumber. More precisely, we have the following expression for $g_n$ and $u_n$ (see [11] for more details about the proof):

$$ g_n(x) = b_{\text{top}}(x)\phi_n(1) + b_{\text{bot}}(x)\phi_n(0) + f_n(x), \quad u_n(x) = \frac{i}{2k_n} \int_{\mathbb{R}} g_n(z)e^{ik_n|x-z|}dz. \quad (5) $$

When $h$ is variable, the decomposition (4) motivates the following definition:

**Definition 2.** The local wavenumber function of the mode $n \in \mathbb{N}$ is the complex function $k_n : \mathbb{R} \to \mathbb{C}$ defined by

$$ k_n^2(x) := k^2 - \frac{n^2\pi^2}{h(x)^2}, \quad (6) $$

with $\text{Re}(k_n)$, $\text{Im}(k_n) \ge 0$.

In this work, as $h(x)$ is non constant, $k_n(x)$ may vanish for some $x \in \mathbb{R}$ and change from a positive real number to a purely imaginary one. We then distinguish three different situations:

**Definition 3.** A mode $n \in \mathbb{N}$ falls in one of these three situations:

1. If $n < kh(x)/\pi$ for all $x \in \mathbb{R}$ then $k_n(x) \in (0, +\infty)$ for all $x \in \mathbb{R}$ and the mode $n$ is called propagative.

2. If $n > kh(x)/\pi$ for all $x \in \mathbb{R}$ then $k_n(x) \in i(0, +\infty)$ for all $x \in \mathbb{R}$ and the mode $n$ is called evanescent.

3. If there exists $x_k^* \in \mathbb{R}$ such that $n = kh(x_k^*)/\pi$ the mode $n$ is called locally resonant. Such points $x_k^*$ are called resonant points, and they are simple if $h'(x_k^*) \ne 0$, and multiple otherwise.

A frequency $k > 0$ for which there exists at least a locally resonant mode is called a locally resonant frequency.
Using the wavenumber function, one can adapt the classic Sommerfeld (or outgoing) condition, defined in [11] for regular waveguides, to general varying waveguides $\Omega$. This condition is used to guarantee uniqueness for the source problem given in equation (2).

**Definition 4.** A wavefield $u \in H_{\text{loc}}^2(\Omega)$ is said to be outgoing if it satisfies

$$ \left| u'_{n}(x) \frac{x}{|x|} - ik_{n}(x)u_{n}(x) \right| \underset{|x| \to +\infty}{\longrightarrow} 0 \quad \forall n \in \mathbb{N}, \qquad (7) $$

where $u_n$ is given in (4).

In all this work, we make the following assumptions:

**Assumption 1.** We assume that $h \in C^2(\mathbb{R}) \cap W^{2,\infty}(\mathbb{R})$ with $h'$ compactly supported and we then define two constants $0 < h_{\min} \le h_{\max} < +\infty$ such that

$$ \forall x \in \mathbb{R} \quad h_{\min} \le h(x) \le h_{\max}, $$

and $h(x) = h_{\min}$ or $h(x) = h_{\max}$ if $x \notin \text{supp}(h')$. For such a function, we define a parameter $\eta > 0$ that satisfies

$$ \|h'\|_{L^{\infty}(\mathbb{R})} < \eta \quad \text{and} \quad \|h''\|_{L^{\infty}(\mathbb{R})} < \eta^2. $$

Such a waveguide is represented in Figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: Representation of an increasing slowly and compactly varying waveguide.</figcaption></figure>

The forward source problem is defined for every frequency $k > 0$ by

$$ (\mathcal{H}_k) : \quad \left\{ \begin{array}{ll} \Delta u + k^2 u = -f & \text{in } \Omega, \\ \partial_{\nu} u = b_{\text{top}} & \text{in } \partial \Omega_{\text{top}}, \\ \partial_{\nu} u = b_{\text{bot}} & \text{in } \partial \Omega_{\text{bot}}, \\ u \text{ is outgoing,} \end{array} \right. \qquad (8) $$

As explained in [12], this problem is not well-posed when $\{x | k_n(x) = 0\}$ is a non-trivial interval of $\mathbb{R}$. This especially happens when $k = n\pi/h_{\min}$ or $k = n\pi/h_{\max}$. We then avoid these two forbidden situations and we set

$$ \delta(k) := \min_{n \in \mathbb{N}} \left( \sqrt{|k^2 - \frac{n^2 \pi^2}{h_{\min}^2}|}, \sqrt{|k^2 - \frac{n^2 \pi^2}{h_{\max}^2}|} \right). \qquad (9) $$

From now on, we define $(f_n)_{n \in \mathbb{N}}$ the modal decomposition of $f$, and

$$ g_n(x) := \frac{f_n(x)}{\sqrt{h(x)}} + \varphi_n(1)b_{\text{top}}(x) \frac{\sqrt{1 + (h'(x))^2}}{\sqrt{h}} + \varphi_n(0)b_{\text{bot}}(x) \frac{1}{\sqrt{h}}. \qquad (10) $$

Using the work done in [10], we are able to provide an approximation of the solution of (8). If $h$ is increasing, we can state the following result using Theorem 1 and Remark 3 in [10].
**Theorem 1.** Let $h$ be an increasing function defining a varying waveguide $\Omega$ that satisfies Assumption 1 with a variation parameter $\eta > 0$. Consider sources $f \in L_c^\infty(\Omega)$, $b := (b_{\text{bot}}, b_{\text{top}}) \in (\text{H}^{1/2}(\mathbb{R}))^2 \cap (L_c^\infty(\mathbb{R})^2)$. Assume that there is a unique locally resonant mode $N \in \mathbb{N}$, associated with a simple resonant point $x_k^* \in \mathbb{R}$. Let $I \subset \mathbb{R}$ be an interval of length $R > 0$, and $\Omega_I := \{(x, y) \in \Omega | x \in I\}$.

There exists $\eta_0 > 0$ depending only on $h_{\min}$, $h_{\max}$, $\delta(k)$ and of the length of supp($h'$) such that if $\eta \le \eta_0$, then the problem $(\mathcal{H}_k)$ admits a unique solution $u \in H_{\text{loc}}^2(\Omega)$. Moreover, this solution is approached by $u^{\text{app}}$ defined for almost every $(x, y) \in \Omega$ by

$$u^{\text{app}}(x, y) := \sum_{n \in \mathbb{N}} \left( \int_{\mathbb{R}} G_n^{\text{app}}(x, s) g_n(s) \mathrm{d}s \right) \varphi_n(y), \quad (11)$$

where $(g_n)_{n \in \mathbb{N}}$ are defined in (10), $\varphi_n$ is defined in (3) and $G_n^{\text{app}}$ is given by

$$G_n^{\text{app}}(x, s) := \begin{cases} \frac{i}{2\sqrt{k_n(s)k_n(x)}} \exp\left(i\left|\int_s^x k_n\right|\right), & \text{if } n < N, \\ \frac{1}{2\sqrt{|k_n(s)|k_n(x)}} \exp\left(-\left|\int_s^x k_n\right|\right), & \text{if } n > N, \end{cases} \quad (12) \\ \begin{cases} \frac{\pi(\xi(s)\xi(x))^{1/4}}{\sqrt{k_n(s)k_n(x)}} (i\mathcal{A} + \mathcal{B}) \circ \xi(s)\mathcal{A} \circ \xi(x) & \text{if } x < s, \\ \frac{\pi(\xi(s)\xi(x))^{1/4}}{\sqrt{k_n(s)k_n(x)}} (i\mathcal{A} + \mathcal{B}) \circ \xi(x)\mathcal{A} \circ \xi(s) & \text{if } x > s, \end{cases} \quad \text{if } n = N.$$

Function $k_n$ is the wavenumber function defined in Definition 2 and the function $\xi$ is given by

$$\xi(x) := \begin{cases} \left(-\frac{3}{2}i \int_x^{x_k^*} k_N(t) \mathrm{d}t\right)^{2/3} & \text{if } x < x_k^*, \\ -\left(\frac{3}{2} \int_{x_k^*}^x k_N(t) \mathrm{d}t\right)^{2/3} & \text{if } x > x_k^*. \end{cases} \quad (13)$$

Precisely, there exist a constant $C_1 > 0$ depending only on $h_{\min}$, $h_{\max}$ and $N$ such that

$$\|u - u^{\text{app}}\|_{\text{H}^1(\Omega_I)} \le C_1 \eta R^2 \delta(k)^{-8} \left( \|f\|_{\text{L}^2(\Omega)} + \|b\|_{(\text{H}^{1/2}(\mathbb{R}))^2} \right). \quad (14)$$

This result provides an approximation of the measurements of the wavefield for every frequency, and a control of the approximation error. We represent in Figure 4 the wavefield for different frequencies. We also point out in Figure 5 that the source should be located in an area where $h > h(x^*)$ in order to generate a significant locally resonant mode.

We notice that at locally resonant frequencies, the wavefield strongly depends on the position of $x_k^*$, which justify the idea of using it to develop an inverse method to reconstruct the width $h$. An illustration is provided in Figure 6 with a representation of the wavefield $u$ and the one side boundary measurement $u(x, 0)$ when $k$ is a locally resonant frequency.

### 3 Shape Inversion using a monochromatic source

In this section, we describe the method to recover the width $h$ from one side boundary measurements of the wavefield at locally resonant frequencies. We focus here on the reconstruction of the shape of
<figure><img src="image_6.png" /><figcaption>Figure 4: Representation of the wavefield amplitude $|u|$ solution of (8) with a monochromatic source $f(x, y) := \delta_{x=s}\varphi_1(y)$ for different frequencies $k$. Data are generated using the finite element method described in Section 5.1. When $k = 30.9$, the mode $n = 1$ is evanescent, and the wavefield decreases very fast around the source. When $k = 31.1$ and $k = 31.4$, the mode $n = 1$ is locally resonant and the wavefield propagates in the waveguide as an Airy function until it reaches the point $x_k^*$. When $k = 32$, the mode $n = 2$ is propagative and the wavefield propagates in all the waveguide. In all these representations, the width $h$ is defined by the function (52).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Representation of the wavefield amplitude $|u|$ solution of (8) with a monochromatic source $f(x, y) := \delta_{x=s}\varphi_1(y)$ for different positions $s$. On top, $h$ is a dilation defined in (55) and on bottom, $h$ is a shrinkage defined in (56). Depending on the location of the source, we observe different behaviors: while top left, bottom left and bottom right show locally resonant modes, the picture on top right show an evanescent mode. In order to generate a significant locally resonant mode, the source $s$ should be placed at a width where $h(s) > h(x_k^*)$ and the mode only propagates until it reaches the point $x_k^*$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Illustration of the amplitude of the wavefield measurements in a varying waveguide. The wavefield $u$ is solution of (8) where the profile $h$ is given in (52) and $k = 31.5$ is a locally resonant frequency. Data are generated using the finite element method described in Section 5.1. Top: amplitude of $|u|$ in the whole waveguide $\Omega$ (non scaled). Bottom: amplitude of the measurements of $|u|$ on the surface $y = 0$.</figcaption></figure>

$h$ on supp($h'$), assuming the a priori knowledge of the constants $h_{min}$ and $h_{max}$ and of an interval containing supp($h'$). We detail in Appendix A how these constants can be estimated. To recover the shape function $h$ on supp($h'$), we start by studying the simpler case of a source term $f$ generating only a single locally resonant mode in the waveguide. Hence we assume that $f$ takes the form

$$ \forall (x, y) \in \Omega \quad f(x, y) = f_N(x)\varphi_N(y), \quad f_N \in L^2(\mathbb{R}), \qquad (15) $$

where $f_N$ is compactly supported. We also assume the absence of boundary source term, meaning that $b = 0$. This simplified situation is useful to understand the method of reconstruction. It will be generalized to any kind on internal and boundary sources in section 4.

In the simplest case of a single internal source (15), we know from the study of the forward problem in Theorem 1 that the measured data without noise $d^{\text{ex}}$ satisfies

$$ d^{\text{ex}}(x) := u(x, 0) \approx u_N^{\text{app}}(x)\varphi_N(0) \quad \text{where} \quad u_N^{\text{app}}(x) = \int_{\mathbb{R}} G_N^{\text{app}}(x, s)f_N(s)ds. \qquad (16) $$

Our reconstruction method is based on the recovery of the resonant point $x_k^*$ for every locally resonant frequencies $k$. It can be seen in Theorem 1 that the approached Green function $G_N^{\text{app}}$ depends on $x_k^*$ through the function $\xi$. However, this dependence is intricate and hardly usable to find a direct link between $d^{\text{ex}}$ and $x_k^*$. Thus, we need to find a simpler approximation of the measurements. In a first part, we provide an approximation of the measured data $d^{\text{ex}}$ using a three parameters model function

$$ d_{z,\alpha,x_k^*}^{\text{app}} : x \mapsto z\mathcal{A}(\alpha(x_k^* - x)), \qquad (17) $$

where $z \in \mathbb{C}^*$ is a complex amplitude, $\alpha > 0$ is a scaling parameter and $x_k^*$ is the resonant point playing the role of a longitudinal shift. In a second part, we first control the approximation error
between this function and the exact measurements and then, we develop a stable way to reconstruct $x_k^*$ from this approximated data.

## 3.1 Wavefield approximation and measurements approximation

In order to find a reconstruction of $x_k^*$, we need to find an exploitable link between $d^{\text{ex}} := u(x, 0)$ and $x_k^*$. To do so, we make a Taylor expansion of $G_N^{\text{app}}$ around the point $x_k^*$. For every frequency $k > 0$ and $R > 0$, we denote

$$ \Omega_R(x_k^*) := \{(x, y) \in \Omega \mid |x - x_k^*| < R\}, \quad \Gamma_R(x_k^*) := (x_k^* - R, x_k^* + R), \quad (18) $$

and we consider the Taylor expansion on the interval $\Gamma_R$. Moreover, we assume that the source is located at the right of the interval $\Gamma_R$, and we thus we define

$$ \Omega_R^+(x_k^*) := \{(x, y) \in \Omega \mid x - x_k^* > R\}, \quad \Gamma_R^+(x_k^*) = (x_k^* + R, +\infty). \quad (19) $$

Both these sets are represented in Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7: Representation of $\Gamma_R(x_k^*)$, $\Omega_R(x_k^*)$, $\Gamma_R^+(x_k^*)$ and $\Omega_R^+(x_k^*)$. The source $f$ is assumed to be compactly supported in $\Gamma_R(x_k^*)$.</figcaption></figure>

The following Proposition shows that $u_N^{\text{app}}$ can be approached by a three parameters function of type $d_{z,\alpha,x_k^*}^{\text{app}}$ if the shape function $h$ is steep enough at $x_k^*$.

**Proposition 1.** Assume that $\Omega$ satisfies assumptions 1, let $R > 0$ and $k > 0$ be a locally resonant frequency associated to the mode $N \in \mathbb{N}$ and the locally resonant point $x_k^*$ that satisfies $h'(x_k^*) \ge \theta\eta$ for some $\theta > 0$. There exist $z \in \mathbb{C}^*$, $\alpha > 0$ and a constant $C_2 > 0$ depending only on $h_{\min}$, $h_{\max}$, $N$ and $\theta$ such that

$$ \|u_N^{\text{app}} - d_{z,\alpha,x_k^*}^{\text{app}}\|_{L^2(\Gamma_R(x_k^*))} \le C_2 (R^{3/2}\eta^{5/6} + R^{5/2}\eta^{7/6}), \quad (20) $$

where $d_{z,\alpha,x_k^*}^{\text{app}}$ is defined in (17).

*Proof.* Using the information about the support of the source term, we know that

$$ u_N^{\text{app}}(x) = \int_{x_k^*+R}^{+\infty} G_N^{\text{app}}(x, s)f_N(s)ds. $$

Using the definition of $G_N^{\text{app}}$ given in (12), there exists a function $q_k$ such that for every $x \in \Gamma_R(x_k^*)$,

$$ G_N^{\text{app}}(x, s) = q_k(s) \frac{(-\xi(x))^{1/4}}{\sqrt{k_n(x)}} A(\xi(x)). $$
It follows that

$$u_N^{\text{app}}(x) = \frac{(-\xi(x))^{1/4}}{\sqrt{k_n(x)}} \mathcal{A}(\xi(x)) \int_{x_k^*+R}^{+\infty} q_k(s) f_N(s) ds.$$

In the following, we denote $\mathcal{O}(\cdot)$ bounds depending only on $h_{\max}$, $h_{\min}$ and $N$. We see that

$$k_N(x)^2 = \frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3}(x - x_k^*) + \mathcal{O}(\eta^2(x - x_k^*)^2).$$

From now on, we assume that $x > x_k^*$, which leads to

$$k_N(x) = \sqrt{\frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3}}(x - x_k^*)^{1/2} + \mathcal{O}\left(\frac{\eta^2(x - x_k^*)^{3/2}}{\sqrt{h'(x_k^*)}}\right), \quad (21)$$

and so

$$\xi(x) = \left( \frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3} \right)^{1/3} (x_k^* - x) + \mathcal{O} \left( \frac{\eta^2(x - x_k^*)^2}{h'(x_k^*)^{2/3}} \right).$$

Then,

$$\mathcal{A}(\xi(x)) = \mathcal{A}\left(\left(\frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3}\right)^{1/3}(x_k^* - x)\right) + \mathcal{O}\left(\frac{\eta^2(x - x_k^*)^2}{h'(x_k^*)^{2/3}}\right),$$

and

$$\frac{(-\xi(x))^{1/4}}{\sqrt{k_N(x)}} = \left( \frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3} \right)^{-1/6} + \mathcal{O} \left( \frac{\eta^2(x - x_k^*)}{h'(x_k^*)^{7/6}} \right).$$

We set

$$z = \left( \frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3} \right)^{-1/6} \int_{x_k^*+R}^{+\infty} q_k(s) f_N(s) ds, \quad \alpha = \left( \frac{2N^2\pi^2h'(x_k^*)}{h(x_k^*)^3} \right)^{1/3}, \quad (22)$$

and it follows that

$$u_N^{\text{app}}(x) = z\mathcal{A}(\alpha(x_k^* - x)) + \mathcal{O}\left(\frac{\eta^2(x - x_k^*)}{h'(x_k^*)^{7/6}}\right) + \mathcal{O}\left(\frac{\eta^2(x - x_k^*)^2}{h'(x_k^*)^{5/6}}\right).$$

The exact same study can be done in the case $x < x_k^*$. Since $|x - x_k^*| \le R$, we conclude that

$$\|u_N^{\text{app}} - d_{z,\alpha,x_k^*}^{\text{app}}\|_{L^2(\Gamma_R(x_k^*))} = \eta^2 \mathcal{O} \left( R^{3/2} (h'(x_k^*))^{-7/6} + R^{5/2} (h'(x_k^*))^{-5/6} \right).$$

To conclude, we use the fact that $h'(x_k^*) \ge \theta\eta$. $\square$

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

**Corollary 1.** Assume that $\Omega$ satisfies assumptions 1, let $R > 0$ and $k > 0$ be a locally resonant frequency associated to the mode $N \in \mathbb{N}$ and the locally resonant point $x_k^*$ that satisfies $h'(x_k^*) \ge \theta\eta$ for some $\theta > 0$. There exist $\eta > 0$ such that if $\eta \le \eta_0$ then there exist $z \in \mathbb{C}^*$, $\alpha > 0$ and a constant $C_3 > 0$ depending only on $h_{\min}$, $h_{\max}$, $N$ and $\theta$ such that

$$\|d^{\text{ex}} - d_{z,\alpha,x_k^*}^{\text{app}}\|_{L^2(\Gamma_R(x_k^*))} \le C_3 (\delta(k)^{-8} R^2 \eta + R^{3/2} \eta^{5/6} + R^{5/2} \eta^{7/6}). \quad (23)$$
*Proof.* We apply Proposition 1, Theorem 1 and trace results which prove that there exists a constant $\gamma > 0$ depending only on $R$, $h_{\min}$ and $h_{\max}$ such that

$$ \|u(\cdot, 0) - u^{\text{app}}(\cdot, 0)\|_{L^2(\Gamma_R(x_k^*))} \le \|u - u^{\text{app}}\|_{H^{1/2}(\Gamma_R(x_k^*))} \le \gamma \|u - u^{\text{app}}\|_{H^1(\Omega_R(x_k^*))}. \quad \square $$

This result indicates a strategy to determine $x_k^*$ from the exact data $d^{\text{ex}}$. Assuming that $\eta$ is small enough to see $d_{z,\alpha,x_k^*}^{\text{app}}$ as a good approximation of $d^{\text{ex}}$, one may fit the three parameters $(\alpha, z, x_k^*)$ that minimize the misfit $d^{\text{ex}} - d_{z,\alpha,x_k^*}^{\text{app}}$ for some norm. We also see from this result that the error in the Taylor expansion strongly depends on the values of $\theta$, $R$ and $\delta(k)$. Since we plan on using different locally resonant frequencies, we need to get a uniform control over $k$ on the Taylor expansion. Moreover, we need to control the length $R$ of the measurement interval. If $R$ is too large, the quality of the Taylor expansion diminishes and if $R$ is too small, the data on the interval $\Gamma_R(x_k^*)$ may not contain enough information to fit the three parameters $(z, \alpha, x_k^*)$ with stability.

To prevent this, $R$ needs to be scaled depending on the parameter $\alpha$. Indeed, the change of variable $x \mapsto \alpha(x_k^* - x)$ must cover a large enough interval to perform the desired fitting. Using the expression of $\alpha$ given in (22), we say that $\alpha(x_k^* - x)$ covers an interval of fixed radius $\sigma > 0$ if

$$ R \ge \frac{\sigma}{\alpha} \ge \frac{\sigma h_{\max}}{2\beta^{1/3}N^2\pi^2}\eta^{-1/3}. $$

This means that $R$ needs to be scaled as $\eta^{-1/3}$ and we assume now that

$$ R := r\eta^{-1/3}. \quad (24) $$

where $r > 0$ is a constant. We can now give a uniform approximation result between $d^{\text{ex}}$ and $d_{z,\alpha,x_k^*}^{\text{app}}$ on the interval $\Gamma_R(x_k^*)$.

**Theorem 2.** Assume the same hypotheses than in Theorem 1, and fix $N \in \mathbb{N}^*$ and $\theta > 0$. There exists $\eta_0 > 0$ such that for any $\eta < \eta_0$ and for any $k > 0$ associated to the mode $N$ and the locally resonant point $x_k^*$ satisfying $h'(x_k^*) \ge \theta\eta$, there exist $\alpha > 0$ and $z \in \mathbb{C}$ such that the following approximation holds on the interval $\Gamma_R(x_k^*)$:

$$ d^{\text{ex}} = d_{z,\alpha,x_k^*}^{\text{app}} + \mathcal{O}(\eta^{1/3}). $$

More precisely, there exists a constant $C_4 > 0$ depending only on $h_{\min}$, $h_{\max}$, $\theta$ and $r$ such that

$$ \|d^{\text{ex}} - d_{z,\alpha,x_k^*}^{\text{app}}\|_{L^2(\Gamma_R(x_k^*))} \le C_4 \eta^{1/3}. \quad (25) $$

*Proof.* To prove this result, we need to provide a uniform lower bound for $\delta(k)$ under the hypothesis $h'(x_k^*) \ge \theta\eta$. From the definition of $\delta(k)$, we assume without loss of generality that $\delta(k)$ is reached by the mode number $N \in \mathbb{N}$ and that $\delta(k)^2 = |k^2 - N^2\pi^2/h_{\max}^2|$. Hence, as $k$ is locally resonant, we know that $k = N\pi/h(x_k^*)$ and that

$$ \delta(k)^2 = N^2\pi^2 \left| \frac{1}{h(x_k^*)^2} - \frac{1}{h_{\max}^2} \right| \le \frac{2N^2\pi^2}{h_{\min}^2h_{\max}}(h_{\max} - h(x_k^*)). $$
Let us call $t \ge 0$ such that $h(x_k^*+t) = h_{\max}$, we necessarily have $h'(x_k^*+t) = 0$. Using the hypothesis $h'' \ge -\eta^2$, we know that $h'(x_k^*+s) \ge \eta\theta - \eta^2s$ for all $s \in (0,t)$. This implies that $t \ge \theta/\eta$. Then,

$$h_{\max} - h(x_k^*) = \int_0^t h'(x_k^* + s)ds \geq \eta\theta t - \eta^2 \frac{t^2}{2} \geq \frac{\theta^2}{2}.$$

Then,

$$\delta(k) \geq \frac{N\pi}{h_{\min}\sqrt{h_{\max}}}\theta.$$

The proof of Theorem 2 is now straightforward, we simply replace $R$ by $r\eta^{-1/3}$ in Corollary 1 and use the lower bound on $\delta(k)$. $\square$

**Remark 1.** We underline the fact that the constant in this estimation depends on the value of $\theta$, and tends toward infinity if $\theta$ tends to zero. This result is illustrated in Figure 8 where we clearly notice that the error between $d^{\text{ex}}$ and $d_{z,\alpha,x_k^*}^{\text{app}}$ deteriorates when $\theta$ become too small.

<figure><img src="image_11.png" /><figcaption>Figure 8: Error of approximation $\|d^{\text{ex}} - d_{z,\alpha,x_k^*}^{\text{app}}\|_{L^2(\Gamma_R(x_k^*))}$ for different values of $\theta$ with a fixed value of $\eta$. The width $h$ is represented on the left of the picture and exact data $d^{\text{ex}}$ are generated as explained in section 5.1 with a locally resonant mode $N = 1$ and a source $f(x, y) = \delta_{x=0.6}\varphi_N(y)$. Then, $d^{\text{ex}}$ is compared with $d_{z,\alpha,x_k^*}^{\text{app}}$ where $\alpha$ and $z$ are defined using the expression (22). Here, $r = 0.2$ and $\eta = 8.10^{-4}$.</figcaption></figure>

## 3.2 Stable reconstruction of $x_k^*$

We proved so far that the one side boundary measurements are close to the three parameters function

$$d_{z,\alpha,x_k^*}^{\text{app}} : x \mapsto z\mathcal{A}(\alpha(x_k^* - x)), \qquad (26)$$

on the interval $\Gamma_{r\eta^{-1/3}}(x_k^*) = (x_k^* - r\eta^{-1/3}, x_k^* + r\eta^{-1/3})$. The question is now to understand if one can get a stable evaluation of $x_k^*$ from this approached data. We define then the following forward operator

$$F: \begin{cases} \mathbb{C}^* \times (0, +\infty) \times \mathbb{R} & \longrightarrow \mathcal{C}^0(\Gamma_R(\beta/\alpha)) \\ (z, \alpha, \beta) & \longmapsto (x \mapsto z\mathcal{A}(\beta - \alpha x)) \end{cases}. \qquad (27)$$
If the three parameters are uniquely defined, we can deduce an approximation of the position $x_k^*$ from the solution of $F(z, \alpha, \beta) = d_{z,\alpha,x_k^*}^{\text{app}}$ through the identification $x_k^* = \beta/\alpha$. The following result guarantee the uniqueness of the solution to this problem.

**Proposition 2.** Let $d_0 := F(z_0, \alpha_0, \beta_0)$, there exists $r_0 > 0$ such that if $r > r_0$ then the problem $F(z, \alpha, \beta) = d_0$ has a unique solution $(z_0, \alpha_0, \beta_0)$.

*Proof.* If $r$ is high enough, the interval $\Gamma_{r\eta^{-1/3}}$ is large enough to contain the maximum amplitude of the function $d_0$ at position $x_{\max}$, and the two first zeros of the function $d_0$ called $x_1$ and $x_2$. Hence $(z_0, \alpha_0, \beta_0)$ are uniquely determined by

$$z_0 = \frac{d_0(x_{\max})}{\|\mathcal{A}\|_{\infty}}, \quad \alpha_0 = \frac{y_1 - y_2}{x_2 - x_1}, \quad \beta_0 = \frac{y_1 x_2 - y_2 x_1}{x_2 - x_1}, \qquad (28)$$

where $y_1$ and $y_2$ are the two first zeros of the Airy function $\mathcal{A}$.

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

This result of uniqueness and the corresponding inversion formulas are not directly applicable. Indeed, these formulas are not robust to noise or data uncertainties and we remind the reader that an approximation is already made between $d^{\text{ex}}$ and $d_{z,\alpha,x_k^*}^{\text{app}}$.

A first strategy would be to use the explicit formula (28) and apply some data regularization with a low-pass filter before doing the inversion (see [21] for more details). This method works perfectly with exact data, and will be called from now the “direct parameters reconstruction method”.

However, few issues can be raised with this method. Firstly, it is not sufficiently robust if data are very noisy, as illustrated in Figure 10 where we use this method to reconstruct parameters with a random additive noise of increasing amplitude. Secondly, as mentioned in the previous subsection, we will only make measurements of the wavefield on the interval $\Gamma_R$, and we cannot be sure that neither the first two zeros nor the maximum of $\mathcal{A}$ will occur in this interval. Finally, it is more realistic to assume that we only have access to $d_{z,\alpha,x_k^*}^{\text{app}}$ on a finite number of receivers positions. For all this reasons, a least squares approach is introduced. On can keep the previous method to find an initial guess of the parameters $\mathbf{p} := (z, \alpha, \beta)$.

We now assume that we have access to the data

$$d_i := F(p)(t_i), \quad i = 1, \dots, n, \qquad (29)$$

where $\mathbf{p} := (z, \alpha, \beta)$ and $\mathbf{t} := (t_i)_{i=1,...,n} \in \Gamma_R$ is an increasing subdivision of $\Gamma_R$. We denote by $\mathbf{d} := (d_1, \dots, d_n) \in \mathbb{C}^n$ the discrete data on the subdivision $\mathbf{t}$. We aim to minimize the least squares energy functional

$$J_d(\mathbf{p}) := \frac{1}{2} \|F(\mathbf{p})(\mathbf{t}) - \mathbf{d}\|_{\ell^2}^2. \qquad (30)$$

In this expression, the norm $\ell^2$ is the normalized euclidean norm defined by

$$\|\mathbf{d}\|_{\ell^2}^2 = \frac{1}{n} \sum_{i=1}^{n} d_i^2. \qquad (31)$$

Denoting by $\tau$ the step of the subdivision $\mathbf{t}$, we know using quadrature results (see [14] for more details) that for all functions $f \in \mathbb{H}^2(\Gamma_R)$,

$$\left| \|f\|_{\text{L}^2(\Gamma_R)} - \|f(\mathbf{t})\|_{\ell^2} \right| \le 2\tau \|f\|_{\text{H}^1(\Gamma_R)}. \qquad (32)$$
It shows that the $\ell^2$ norm is a good approximation of the $L^2(\Gamma_R)$ norm if the step of discretization $\tau$ is small enough.

We assume the knowledge of an open set $U \subset \mathbb{C}^* \times (0, +\infty) \times \mathbb{R}$ of the form

$$ U := B_{\mathbb{C}}(0, z_{\max}) \times (\alpha_{\min}, \alpha_{\max}) \times (\beta_{\min}, \beta_{\max}), \quad z_{\max}, \alpha_{\min}, \alpha_{\max} \in \mathbb{R}_+^*, \quad \beta_{\min}, \beta_{\max} \in \mathbb{R}, \quad (33) $$

containing the solution $p$. The following proposition shows that the least squares problem is locally well-posed if the sampling size $n$ is large enough, and it quantifies the error on the recovered parameters.

**Proposition 3.** There exists $n_0 \in \mathbb{N}^*$ such that if $n \ge n_0$, then for every $p_0 \in U$ and $d_0 := F(p_0)(x)$, there exist $\varepsilon > 0$ and $U' \subset U$ such that for every $d \in \mathbb{R}^n$ satisfying $\|d - d_0\|_{\ell^2} < \varepsilon$ then the functional $J_d$ is strictly convex on $U'$ and admits a unique minimizer denoted $p_{LS} = (z_{LS}, \alpha_{LS}, \beta_{LS})$. Moreover, there exists a constant $C_5 > 0$ depending only on $U$, $p_0$ and $n$ such that

$$ \|p_{LS} - p_0\|_2 \le C_5 \|d - d_0\|_{\ell^2}. \quad (34) $$

Finally, if we denote

$$ \Lambda: \left\{ \begin{array}{rcl} B_2(d_0, \varepsilon) & \rightarrow & \mathbb{R} \\ d & \mapsto & \frac{\beta_{LS}}{\alpha_{LS}} \end{array} \right. , \quad (35) $$

the operator which approach the value of $x_k^*$, there exists a constant $C_6 > 0$ depending only on $U$, $p_0$ and $n$ such that

$$ |\Lambda(d) - \Lambda(d_0)| \le C_6 \|d - d_0\|_{\ell^2}. \quad (36) $$

The proof of this result is given in Appendix B. In this proof we can see that the choice of the data discretization points $t$ is important in order to improve the accuracy of the reconstruction. We illustrate that by choosing $p_0 = (2 + i, 1.4, -2.8)$ and comparing the direct reconstruction of the parameters (28) with the least squares method for three different sets $t$:

* $t_1$ is a discretization of $[-6, -1]$ with 100 points, where the tail of the Airy function is located.

* $t_2$ is a discretization of $[-2, 6]$ with 100 points, where the Airy function varies.

* $t_3$ is a discretization of $[-6, 6]$ with 200 points.

In Figure 9, we represent the four reconstructions where $d = d_0 + \mathcal{N}$ and $\mathcal{N}$ is a normal random additive noise of amplitude 0.3. We see that reconstructions with $t_2$ and $t_3$ seems more accurate than the one with $t_1$ and the direct reconstruction. We detail this in Figure 10, where we compare the reconstruction errors with different noise amplitudes. We clearly see that using a least squares method improves the reconstruction if $t$ is well-chosen.

We now have a stable method to reconstruct an approximation $x_k^{*,\text{app}}$ of $x_k^*$ from given measurements of $d(x) := u(x, 0)$. Since $h(x_k^*) = N\pi/k$, we can approximate the width at positions $x_k^{*,\text{app}}$ with the formula

$$ h^{\text{app}}(x_k^{*,\text{app}}) = \frac{N\pi}{k}. \quad (37) $$

This provides an approximation of $h$ in one point. Then, we change the frequency $k > 0$ to get approximations of $h$ all along the support of $h'$. From now on, we denote $u_k$ the wavefield
<figure><img src="image_13.png" /><figcaption>Figure 9: Comparison of four different reconstructions of $F(p_0)$ where $d-d_0$ is a normal random additive noise of amplitude 0.3. Blue: direct reconstruction. Green: least square reconstruction with $t = t_1$. Red: least square reconstruction with $t = t_2$. Purple: least square reconstruction with $t = t_3$.</figcaption></figure>

propagating at frequency $k$. As mentioned before, we assume that we already have an approximation of $\text{supp}(h')$, $h_{\min}$ and $h_{\max}$ (see Appendix A). We denote

$$k_{\max} := \frac{N\pi}{h_{\min}}, \quad k_{\min} := \frac{N\pi}{h_{\max}}. \qquad (38)$$

and we take a finite set of frequencies $K \subset (k_{\min}, k_{\max})$. For every frequency $k \in K$, we introduce $t_k$ a discretization of $\Gamma_R$ and we set

$$\forall k \in K \quad d_k := u_k(t_k, 0) + \zeta_k, \qquad (39)$$

the measurement of $u_k$ corrupted by an error term $\zeta_k$. For every $k \in K$, using Proposition 3, we have an $\varepsilon_k$ and a constant $C_6^k$ associated with $d_{k,0} = d_{z,\alpha,x_k^*}^{\text{app}}$. We define

$$\varepsilon := \min_{k \in K}(\varepsilon_k), \quad C_6 := \max_{k \in K}(C_6^k), \qquad (40)$$

and

$$X^{*,\text{app}} := \{\Lambda(d_k), k \in K\}. \qquad (41)$$

Using the known approximation of $\text{supp}(h')$ provided by Appendix A, we set the coordinates $x_0$ and $x_{n+1}$ such that $\text{supp}(h') \subset (x_0, x_{n+1})$ and $X^{*,\text{app}} \subset (x_0, x_{n+1})$. We then define the function $h^{\text{app}}$ as the piecewise linear function such that

$$h^{\text{app}}(x_k^{*,\text{app}}) = \frac{N\pi}{k} \quad \forall k \in K, \quad h^{\text{app}}(x_0) = \frac{N\pi}{k_{\max}}, \quad h^{\text{app}}(x_{n+1}) = \frac{N\pi}{k_{\min}}. \qquad (42)$$

Using all the previous results, we are able to quantify the error of reconstruction between $h^{\text{app}}$ and $h$:
<figure><img src="image_14.png" /><figcaption>Figure 10: Error of reconstruction of $p_0$ and $\Lambda(d_0)$ with respect to the noise on the data for different method of reconstruction. On the left, $\|p_{LS} - p_0\|_2$ with respect to $\|d - d_0\|_{\ell^2}/\|d\|_{\ell^2}$. On the right, $|\Lambda(d_0) - \Lambda(d)|$ with respect to $\|d - d_0\|_{\ell^2}/\|d\|_{\ell^2}$.</figcaption></figure>

**Theorem 3.** Let $K$ be a finite subset of $(k_{\min}, k_{\max})$. Assume the same hypotheses than in Theorem 1, and fix $N \in \mathbb{N}^*$ and $\theta > 0$ such that for every $k \in K$, $h'(x_k^*) \ge \theta\eta$. There exist $\eta_1 > 0$ and $\zeta_0 > 0$ such that if $\eta < \eta_1$ and $\max_{k \in K} \|\zeta_k\|_{\ell^2} \le \zeta_0$, then there exists a constant $C_7 > 0$ depending only on $h_{\min}$, $h_{\max}$, $N$, $r$, $\mu$ and $\theta$ such that

$$ \|h^{\text{app}}(X^{\star, \text{app}}) - h(X^{\star, \text{app}})\|_{\infty} \le \eta C_7 \left( \eta^{1/3} + \max_{k \in K} \|\zeta_k\|_{\ell^2} \right). \quad (43) $$

*Proof.* For every $k \in K$, we notice that

$$ |h^{\text{app}}(x_k^{\star, \text{app}}) - h(x_k^{\star, \text{app}})| = |h(x_k^{\star}) - h(x_k^{\star, \text{app}})| \le \eta|x_k^{\star} - x_k^{\star, \text{app}}| = \eta|\Lambda(d_0) - \Lambda(d_{k,0})|. $$

Using Theorem 2 combined with the control (32), we know that

$$ \|d_k - d_{0,k}\|_{\ell^2} \le 2\eta^{1/3}\tau C_4 \left( \|d^{\text{ex}}\|_{\text{H}^1(\Gamma_R)} + \|d_{z,\alpha,x_k^*}^{\text{app}}\|_{\text{H}^1(\Gamma_R)} \right). $$
Using Theorem 1, we can control both these norm by a constant $c_2 > 0$ which does not depend on $k$. Then, if $\eta$ and $\max_{k \in K} ||\zeta_k||_{\ell^2}$ are small enough then,

$$2\tau C_4 c_2 \eta^{1/3} + \max_{k \in K} ||\zeta_k||_{\ell^2} \le \varepsilon,$$

and using Proposition 3,

$$|h^{\text{app}}(x_k^{\star, \text{app}}) - h(x_k^{\star, \text{app}})| \le \eta C_6 \left( 2\tau C_4 c_2 \eta^{1/3} + \max_{k \in K} ||\zeta_k||_{\ell^2} \right).$$

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

The first error term in this theorem is a consequence of the approximation of the data $d^{\text{ex}}$ by the Airy function $d^{\text{app}}$, while the second error term is caused by the eventual presence of measurements noises. We illustrate this reconstruction in Figure 11 where we choose a subset $K$ with 20 points, and we compare $h^{\text{app}}(x_k^{\star, \text{app}})$ and the exact values $h(x_k^{\star})$.

<figure><img src="image_16.png" /><figcaption>Figure 11: Representation of $x_k^{\star}$ and $x_k^{\star, \text{app}}$ when $K = \{30.9 : 31.95 : 20\}$. Here, $h$ is defined in (50), $N = 1$, and data are generated as explained in section 5.1 with a source $f(x, y) = \delta_{x=6\pi/5}(y)$.</figcaption></figure>

To conclude, we have proved in this section that we are able to reconstruct in a stable way a set of resonant positions $x_k^{\star}$, which leads to a stable reconstruction of the function $h$. We present in the next section the general idea needed to generalize this reconstruction to more realistic non-monochromatic sources terms.

# 4 Inversion using a general source term

We now consider the general case without any a priori assumptions on the source terms $f$ and $b$. We use the same arguments as before to reconstruct the width using the locally resonant mode $N$. However, other modes may also be present in the wavefield, and the previous approximation of $d^{\text{ex}}$ by $u_N^{\text{app}}$ (see (16)) may not be valid. We provide here two methods to exploit the framework developed in the previous section.
The first idea is to treat all resonant and evanescent modes as an added noise to the locally resonant mode. Indeed, we know using Theorem 1 that

$$d^{\text{ex}} := u(x, 0) \approx \sum_{n \in \mathbb{N}} \int_{\mathbb{R}} G_n^{\text{app}}(x, s) g_n(s) \mathrm{d}s \varphi_n(y). \quad (44)$$

Defining the noise

$$\zeta_k(x) := \sum_{n \neq N} \int_{\mathbb{R}} G_n^{\text{app}}(x, s) g_n(s) \mathrm{d}s \varphi_n(y), \quad (45)$$

we see that

$$d^{\text{ex}} \approx u_N^{\text{app}}(x)\varphi_N(0) + \zeta_k(x). \quad (46)$$

To make a successful use of Theorem 3, the noise $\zeta_k$ needs to be smaller than the parameter $\zeta_0$. Using the expression of $G_n^{\text{app}}$ given in (12), we notice that

$$\forall n \in \mathbb{N}, \quad |G_n^{\text{app}}(x)| = \mathcal{O}\left(\frac{1}{\min(|k_n|)}\right). \quad (47)$$

Since

$$\frac{|k_N|}{|k_n|} \le \frac{|k_N|^2 h_{\max}^2}{\sqrt{2N - 1}\pi}, \quad (48)$$

the amplitude of $\zeta_k$ seems to be neglectable compared to the amplitude of $u_N^{\text{app}}$ if the width $h_{\max}$ is small enough and if $\|f_n\|_{L^2(\mathbb{R})} = \mathcal{O}(\|f_N\|_{L^2(\mathbb{R})})$ for every $n \neq N$. We illustrate this in Figure 12 where we compare the amplitude of $u(x, 0)$ with the amplitude of $u_N^{\text{app}}(x)\varphi_N(0)$ for source terms $f(x, y) = y^2\delta_{x=6}$ and $b_{\text{top}}(x) = \delta_{x=6}$. We see that it is possible to fit directly the Airy function on $u(x, 0)$, since the noise is very small.

<figure><img src="image_17.png" /><figcaption>Figure 12: Representation of $u(x, 0)$ and $u_N^{\text{app}}(x)\varphi_N(0)$ generated by a source $f(x, y) = y^2\delta_{x=6}$ and $b_{\text{top}}(x) = \delta_{x=6}$. Here $N = 1$, $k = 31.7$ and $h$ is defined as in (50).</figcaption></figure>

However, in the general case, the noise of the other modes can not always be neglected. In this case, we can filter the measurements to keep only the locally resonant part. Evanescent modes vanish away from the source and so their contribution is negligible in $u(x, 0)$. As for propagative
modes, we notice that $k_n(x)$ is almost constant when $n < N$. Using (12), it means that every propagative mode is oscillating with a frequency almost equal to $k_n$. We illustrate it in Figure 13 with the Fourier transform of $u(x, 0)$ and the contribution of each mode. Using a filter cutting all frequencies around $k_n(x)$ for $n < N$, we can clean the signal and get a good approximation of $u_N(x)\varphi_N(0)$ (see [21] for more details). We illustrate this in Figure 14 where we represent $u(x, 0)$, $u_N(x)\varphi_n(0)$ and the approximation obtained using this Fourier filtering.

<figure><img src="image_18.png" /><figcaption>Figure 13: Influence of each mode on the Fourier transform of the measurements real part. Fourier transform of Real($u_n\varphi_n(0)$) are plotted for every non evanescent mode ($n = 0, 1, 2$). For comparison purpose, $k_n(0)$ is represented for every propagative modes ($n = 1, 2$). Here, $k = 63.4$, $N = 2$, $f(x, y) = \delta_{x=6}(3\varphi_0(y) + 2\varphi_1(y) + \varphi_2(y))$ and $b_{top} = \delta_{x=6}$.</figcaption></figure>

We can now fit the three parameters Airy function on the measurements, and find an approximation of $x_k^*$ as before. It proves that our method can be used in the case of general sources, providing a filtering of propagative mode, and does not need any *a priori* information on sources as long as their locally resonant part does not vanish. Using the extension of Theorem 1 to non monotone waveguides provided in Section 4 of [10], all the results presented in this section remain true in any kind of slowly varying waveguide.

# 5 Numerical computations

In this section, we show some numerical applications of our reconstruction method on slowly varying waveguides. We simulate one side boundary measurements using numerically generated data, and we provide reconstructions of increasing an non monotone waveguides with different shape profiles.

## 5.1 Generation of data

In the following, numerical solutions of (8) are generated using the software Matlab to solve numerically the equation in the waveguide $\Omega$. In every numerical simulation, we assume that $h'$ is supported between $x = -7$ and $x = 7$. To generate the solution $u$ of (8) on $\Omega_7$, we use a self-coded finite element method and a perfectly matched layer [9] on the left side of the waveguide between $x = -15$ and $x = -8$ and on the right side between $x = 8$ and $x = 15$. The coefficient of absorption
<figure><img src="image_19.png" /><figcaption>Figure 14: Fourier filtering of Real($u_k(x, 0)$) and comparison with Real($u_{k,N}^{\text{app}}(x)\varphi_N(0)$). Here, $k = 63.4$, $N = 2$, $f(x, y) = \delta_{x=6}(3\varphi_0(y) + 2\varphi_1(y) + \varphi_2(y))$ and $b_{\text{top}} = \delta_{x=6}$.</figcaption></figure>

for the perfectly matched layer is defined as $\alpha = -k((x-8)\mathbf{1}_{x \ge 8} - (x+8)\mathbf{1}_{x \le -8})$ and $k^2$ is replaced in the Helmholtz equation with $k^2 + i\alpha$. The structured mesh is built with a stepsize of $10^{-3}$.

## 5.2 Method of reconstruction

Using all the previous results, we summarize here all the steps to reconstruct the approximation $h^{\text{app}}$ of the width $h$.

1. Find an approximation of supp($h'$), $k_{\min}$ and $k_{\max}$ using the method described at the beginning of Appendix A.

2. Choose a set of frequencies $K \subset (k_{\min}, k_{\max})$ and sources $f, b_{\text{top}}, b_{\text{bot}}$. Then, for every frequency, measure the wavefield $d(x) := u(x, 0)$ solution of (8).

3. Filter the data by eliminating in the Fourier transform responses around $k_n$ for every a prop- agative mode $n$.

4. Find an approximation of the coordinate $x_{\max}$ where $|u(x, 0)|$ is maximal. Then, choose a length $R > 0$, and discretize the interval $[x_{\max} - R, x_{\max} + R]$ into the set $t_k$. The available data $t_k$ are then the measurements of $u(t_k, 0)$.

5. For every frequency, minimize the quantity $\|d_k - F(p)(t_k)\|_{\ell^2}$ using a gradient descent to find the approximation $x_k^{\star, \text{app}}$ of $x_k^*$. The direct reconstruction method can be used to initialize the gradient descent method.

6. Compute $h^{\text{app}}$ using expression (42).

In step 3, we should normally discretize the interval $\Gamma_R(x_k^*)$. However, since we do not know yet the value of $x_k^*$, we use the fact that $x_k^*$ is located near $x_{\max}$, and that the interval $[x_{\max}-R, x_{\max}+R]$ can be included into a bigger interval $\Gamma_{R'}(x_k^*)$ where $R' > R$.
## 5.3 Numerical results

We now apply this method to reconstruct different profiles of slowly varying waveguides. Firstly, we present in Figure 15 the reconstruction $h^{\text{app}}$ obtained for different increasing functions $h$. We choose four different waveguide profiles

$$ \Omega_1: \quad h_1(x) = 0.1 + \gamma_1 \left( \frac{x^5}{5} - 32 \frac{x^3}{3} + 256x \right) \mathbf{1}_{-4 \le x \le 4} - \gamma_2 \mathbf{1}_{x<-4} + \gamma_2 \mathbf{1}_{x>4}, \quad (49) $$

$$ \Omega_2: \quad h_2(x) = 0.1 + \gamma_3 \left( \frac{x^5}{5} - 2x^4 + 16 \frac{x^3}{3} \right) (\mathbf{1}_{0 \le x \le 4} - \mathbf{1}_{-4 \le x < 0}) + \gamma_4 (\mathbf{1}_{x>4} - \mathbf{1}_{x<-4}), \quad (50) $$

$$ \Omega_3: \quad h_3(x) = 0.1 + \gamma_5 x \mathbf{1}_{-4 \le x \le 4} + 4\gamma_5 \mathbf{1}_{x>4} - 4\gamma_4 \mathbf{1}_{x<-4}, \quad (51) $$

$$ \Omega_4: \quad h_4(x) = 0.1 - 4\gamma_5 + 4\gamma_5 \frac{\sqrt{x+4}}{\sqrt{2}} \mathbf{1}_{-4 \le x \le 4} + 8\gamma_5 \mathbf{1}_{x>4}. \quad (52) $$

where $\gamma_1 = 3.10^{-6}$, $\gamma_2 = 8192/5.10^{-6}$, $\gamma_3 = 5.10^{-5}$, $\gamma_4 = 53/3.10^{-5}$, $\gamma_5 = 0.01/30$. All these profiles are represented in black in Figure 15. We impose that sources of (8) need to be located at the right of the waveguide in order to generate significant locally resonant mode (see Figure 5), and we define

$$ f(x, y) = y\delta_{x=6}, \quad b_{\text{top}}(x) = \delta_{x=6}, \quad b_{\text{bot}} = 0. \quad (53) $$

We choose to work with the following sets of frequencies

$$ K_1 = \{30.92 : 31.93 : 20\}, \quad K_2 = \{30.9 : 31.95 : 20\}, \quad K_3 = K_4 = \{31.01 : 31.83 : 20\}. \quad (54) $$

The profile $h_1$ and the set $K_1$ satisfy all the assumptions of Theorem 2, while the derivative of $h_2 \in C^2(\mathbb{R})$ vanishes once. The last two profiles $h_3$ and $h_4$ are not in $C^2(\mathbb{R})$ and show corners where the derivative of the profile is not continuous. The function $h'_3$ is piecewise continuous and bounded, contrary to $h'_4$ which explodes at $x = -4$.

We plot in red in Figure 15 the reconstructions $h^{\text{app}}$, slightly shifted, obtained using our method of reconstruction. We also compute the $L^\infty$ norm of $h - h^{\text{app}}$ in each situation. We clearly see that the reconstruction deteriorates when $h'(x_k^*)$ is too small or when the function $h$ is not sufficiently smooth.

Secondly, we present in Figure 16 the reconstruction $h^{\text{app}}$ obtained for different non monotoneous widths. We choose three different profiles defined by

$$ \Omega_5: \quad h_5(x) = 0.1 + \gamma_6 \sin \left( \frac{\pi}{10} (x+5) \right) \mathbf{1}_{-5 \le x \le 5}, \quad (55) $$

$$ \Omega_6: \quad h_6(x) = 0.1 - \gamma_7(x+5)\mathbf{1}_{-5 \le x \le 0} + \frac{\gamma_6}{4}(x-4)\mathbf{1}_{0 < x \le 4}, \quad (56) $$

$$ \Omega_7: \quad h_7(x) = 0.1 - \gamma_8 \sqrt{3} + 2\gamma_8 \sin \left( \frac{4\pi\sqrt{x+5}}{3} \right) \mathbf{1}_{-3.5 \le x \le 4} + 2\gamma_8 \sin \left( \frac{4\pi\sqrt{1.5}}{3} \right) \mathbf{1}_{x<-3.5}, \quad (57) $$

where $\gamma_6 = 25.10^{-4}$, $\gamma_7 = 5.10^{-4}$, $\gamma_8 = 4.10^{-4}$. All these profiles are represented in black in Figure 16. The profile $h_5$ represent a dilation of the waveguide, while $h_6$ represent a compression of the waveguide. The profile $h_7$ is the more general one with both compressions and dilations. Again,
<figure><img src="image_20.png" /><figcaption>Figure 15: Reconstruction of four different increasing profiles. In black, the initial shape of $\Omega_5$ (not scaled), and in red the reconstruction, slightly shifted for comparison purposes. In each case, $K = K_i$ is defined in (54), $h = h_i$ is defined in (49), (50), (51), (52), and the sources of (8) are defined in (53). We also compute the relative error $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}}$ between the exact reconstruction and the approximation. Top left: $i = 1$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 0.49\%$. Top right: $i = 2$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 0.94\%$. Bottom left: $i = 3$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 0.40\%$. Bottom right: $i = 4$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 1.6\%$.</figcaption></figure>

sources in (8) are located in every large area of the waveguide to generate significant locally resonant modes and are defined by

$$f^5(x, y) = \delta_{x=0}y, \quad b_{\text{top}}^5(x) = \delta_{x=0}(x), \quad b_{\text{bot}}^5 = 0, \qquad (58)$$

$$f^6(x, y) = (\delta_{x=6} + \delta_{x=-6})y, \quad b_{\text{top}}^6(x) = (\delta_{x=6} + \delta_{x=-6}), \quad b_{\text{bot}}^6 = 0, \qquad (59)$$

$$f^7(x, y) = (\delta_{x=-1.5} + \delta_{x=6})y, \quad b_{\text{top}}^7(x) = (\delta_{x=-1.5} + \delta_{x=6}), \quad b_{\text{bot}}^7 = 0. \qquad (60)$$

We also define the frequency sets

$$K_5 = \{30.65 : 31.4 : 20\}, \quad K_6 = \{31.42 : 32.21 : 20\}, \quad K_7 = \{30.97 : 31.43 : 20\}. \qquad (61)$$

We plot in red in Figure 16 the reconstructions $h^{\text{app}}$, slightly shifted, obtained using our method of reconstruction. We also compute the $L^{\infty}$ norm of $h - h^{\text{app}}$ in each situation.

To conclude, we have provided in this section a method to reconstruct slowly varying width defects given one side boundary measurements of the wavefield at local resonant frequencies. This method works for every type of source and does not require any a priori information on the source except the fact that it is located away from the defect. This reconstruction method is very sensitive to small defects and works numerically to reconstruct dilatations or compressions of the waveguide.

# 6 Conclusion

In this article, we have used the study of the forward problem in slowly varying waveguide presented in [10] and the approximation of the solutions as combination of Airy functions to develop a new inverse method to reconstruct the width of slowly varying waveguides. Given wavefield measurements at the surface of the waveguide for different locally resonant frequencies, we reconstruct the associated locally resonant points which lead to a good approximation of the width of the waveguide.
<figure><img src="image_21.png" /><figcaption>Figure 16: Reconstruction of three different general profiles. In black, the initial shape of $\Omega_6$, and in red, the reconstruction, slightly shifted for comparison purposes. In each case, $K = K_i$ is defined in (61), $h = h_i$ is defined in (55), (56), (57), and sources of (8) are defined in (58), (59), (60). We also compute the relative error $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}}$ between the exact reconstruction and the approximation. Top left: $i = 5$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 0.57\%$. Top right: $i = 6$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 0.81\%$. Bottom: $i = 7$, $\|h - h^{\text{app}}\|_{\infty}/h_{\text{max}} = 0.97\%$.</figcaption></figure>

One main advantage of this new method is that it does not require any *a priori* information on the sources, and we believe that it could be applied to develop new non destructive passive monitoring methods. Moreover, using locally resonant modes, this method can detect small variations of the width with a high sensibility. More importantly, when we compare this new method with the usual back scattering method developed for instance in [11], we notice that this new method seems a lot more precise: while the best relative reconstruction errors are of the order of 8% in [11], we reach in this papers relative errors of the order of less than 1%. Even if measurements are not taken in the same area, this improvement of the reconstruction precision must be underlined.

We believe that this work could be extended to elastic waveguides in two dimensions, using the modal decomposition in Lamb modes presented in [23]. After generalizing it to elastic waveguide, we plan in future works to test the method developed in this article on experimental data to see if the good numerical results obtained using data generated by finite element methods can be reproduced. Indeed, different physical experiments have already been conducted to recover width defects using locally resonant frequencies (see [5, 13]) and we hope that the present reconstruction method could both justify and improve these width reconstructions.

## Appendix A: Identification of $\text{supp}(h)$, $k_{\text{min}}$ and $k_{\text{max}}$

Giving a compactly perturbed waveguide $\Omega$, we describe here how one side boundary measurements enable to approximate very precisely the quantities $\text{supp}(h)$, $k_{\text{min}}$ and $k_{\text{max}}$. The article [12] mentions that the problem (8) is not well-defined when $k_n(x) = 0$ in a non-trivial interval, which especially happens when $k = n\pi/h_{\text{min}}$ or $k = n\pi/h_{\text{max}}$. Numerically, this results in an explosion of the solution when $k$ tends to $n\pi/h_{\text{min}}$ (resp. $n\pi/h_{\text{max}}$) with a source term located in the area where $h(x) = h_{\text{min}}$ (resp. $h(x) = h_{\text{max}}$). This behavior is illustrated in Figure 17.

Measuring the surface wavefield while $k$ varies and detecting its explosions provides a good approximation of the width at the left and the right of the waveguide. Then, we choose a frequency $k = n\pi/h_{\text{max}}$ or $k = n\pi/h_{\text{min}}$ and we move the sources while measuring the amplitude of the
<figure><img src="image_22.png" /><figcaption>Figure 17: L²-norm of $u(x, 0)$ on the interval $(-8, 8)$ with respect to $k$ for a source $b_{\text{top}} = \delta_{x=-5}$ at the left of $\text{supp}(h)$, and a source $b_{\text{top}} = \delta_{x=5}$ at the right of $\text{supp}(h)$. For comparison purposes, $n\pi/h_{\text{max}}$ and $n\pi/h_{\text{min}}$ are plotted for $n = 1$ and $n = 2$.</figcaption></figure>

wavefield. It the source is located outside of the support of $h$, the wavefield is supposed to explode, which provides a good approximation of the support of $h$. This behavior is illustrated in Figure 18.

<figure><img src="image_23.png" /><figcaption>Figure 18: L²-norm of $u(x, 0)$ on the interval $(-8, 8)$ with respect to the position $s$ of the source $b_{\text{top}} = \delta_{x=s}$, for a frequency $k = pi/h_{\text{max}}$ and $k = \pi/h_{\text{min}}$. For comparison purposes the support of $h$ is plotted.</figcaption></figure>

## Appendix B: Proof of Proposition 3

*Proof.* For every $\boldsymbol{p} \in U$, we compute

$$ \nabla J_d(\boldsymbol{p}) = \left( \begin{array}{c} \displaystyle \sum_{i=1}^{n} \mathcal{A}(\beta - \alpha t_i)(z\mathcal{A}(\beta - \alpha t_i) - d_i) \\ \displaystyle \sum_{i=1}^{n} -zt_i\mathcal{A}'(\beta - \alpha t_i)(z\mathcal{A}(\beta - \alpha t_i) - d_i) \\ \displaystyle \sum_{i=1}^{n} z\mathcal{A}'(\beta - \alpha t_i)(z\mathcal{A}(\beta - \alpha t_i) - d_i) \end{array} \right), \quad (62) $$
$$ \nabla^2 J_d(\mathbf{p}) = \sum_{i=1}^{n} M_i(\beta - \alpha t_i), \qquad (63) $$

with

$$ M_i = \begin{pmatrix} \mathcal{A}^2 & -t_i \mathcal{A}' \times (2z\mathcal{A} - d_i) & \mathcal{A}'(2z\mathcal{A} - d_i) \\ -t_i \times \mathcal{A}' \times (2z\mathcal{A} - d_i) & z t_i^2 \mathcal{A}'' \times (z\mathcal{A} - d_i) & -zt_i \times \mathcal{A}''(z\mathcal{A} - d_i) \\ \mathcal{A}' \times (2z\mathcal{A} - d_i) & -zt_i \times \mathcal{A}''(z\mathcal{A} - d_i) & z \mathcal{A}'' \times (z\mathcal{A} - d_i) \\ & -t_i z^2 \times (\mathcal{A}')^2 & z^2(\mathcal{A}')^2 \end{pmatrix}. \quad (64) $$

For every $h \in \mathbb{C} \times \mathbb{R}^2$, $h \neq 0$,

$$ (\nabla^2 J_d(\mathbf{p})h|h) = A(\mathbf{p}, h) + B(\mathbf{p}, h), \qquad (65) $$

where

$$ A(\mathbf{p}, h) = \sum_{i=1}^{n} (\mathcal{A}(\beta - \alpha t_i)h_1 + z\mathcal{A}'(\beta - \alpha t_i)(h_3 - h_2t_i))^2, \qquad (66) $$

$$ B(\mathbf{p}, h) = \sum_{i=1}^{n} (z\mathcal{A}(\beta - \alpha t_i) - d_i)(2\mathcal{A}'(\beta - \alpha t_i)h_1(h_3 - h_2t_i) + z\mathcal{A}''(\beta - \alpha t_i)(h_3 - h_2t_i)^2). \quad (67) $$

We want to prove that $A(\mathbf{p}, h) > 0$. To do so, we use the following Lemma.

**Lemma 1.** For every $(h_1, h_2, h_3) \in \mathbb{C} \times \mathbb{R} \times \mathbb{R}$ and $(z, \alpha, \beta) \in U$, the set of zeros of

$$ x \mapsto h_1 \mathcal{A}(\beta - \alpha x) + z \mathcal{A}'(\beta - \alpha x) (h_3 - h_2 x), \qquad (68) $$

is finite, and at most equal to $3\ell + 3$ where $\ell$ is the number of zeros of $\mathcal{A}'$ on

$$ I := [\beta_{\min} - t, \beta_{\max} + t], \quad \text{where} \quad t = \alpha_{\max} \max(|t_1|, |t_n|). \qquad (69) $$

*Proof.* We do a change of variable $u = \beta - \alpha x$, and we see that for every $x \in [x_1, x_n]$, $u \in I$. We now look for $u \in I$ satisfying

$$ g_1(u) := h_1 \mathcal{A}(u) + z \mathcal{A}'(u) \left( h_3 - \frac{h_2}{\alpha} (\beta - u) \right) = 0. \qquad (70) $$

We notice that if $h_1 = 0$ then $g_1(u) = 0$ if and only if $\mathcal{A}'(u) = 0$ or $x = h_3/h_2$, which gives $\ell + 1$ zeros of $g_1$. Otherwise, we notice that if $\mathcal{A}'(u) = 0$ then $g_1(u) \neq 0$ since every zero of the Airy function is simple (see [3]). It means that there exist $\alpha \in \mathbb{C}$ and $\beta \in \mathbb{R}$ such that

$$ g_1(u) = 0 \quad \Leftrightarrow \quad \frac{\mathcal{A}}{\mathcal{A}'}(u) = \alpha u + \beta. \qquad (71) $$

We define

$$ g_2(u) = \frac{\mathcal{A}}{\mathcal{A}'}(u) - \alpha u - \beta, \qquad (72) $$

and then

$$ g_2'(u) = \left(1 - u \left(\frac{\mathcal{A}(u)}{\mathcal{A}'(u)}\right)^2\right), \quad g_2''(u) = \frac{\mathcal{A}(u)}{\mathcal{A}'(u)} \underbrace{\left(2 - 2u \left(\frac{\mathcal{A}(u)}{\mathcal{A}'(u)}\right)^2 - \frac{\mathcal{A}(u)}{\mathcal{A}'(u)}\right)}_{>0}. \qquad (73) $$

Between two zeros of $\mathcal{A}', \mathcal{A}$ vanishes only once, meaning that depending of the value of $\alpha, g_2'$ vanishes at most twice, and so depending of the value of $\beta, g_2$ vanishes at most three times. $\square$
Back to the proof of Proposition 3, we now set $n_0 = 3\ell + 3$, and if $n > n_0$ then it shows that $A(\mathbf{p}, h) > 0$. We denote $\lambda_1(\mathbf{p}) = \min_{h \in \mathbb{R}^3, h \neq 0} a(\mathbf{p}, h) / \|h\|_2^2$. This function is continuous on a subset $U_1 \subset U$, and we denote by $m$ the minimum of $\lambda_1$ on $U_1$. We also notice that

$$|b(\mathbf{p}, h)| \le \|z\mathcal{A}(\beta - \alpha X) - \mathbf{d}\|_1 (\|\mathcal{A}'\|_{\infty} + 2z\|\mathcal{A}''\|_{\infty})(1 + \|\mathbf{t}\|_{\infty})\|h\|_2^2. \quad (74)$$

We see that

$$\|z\mathcal{A}(\beta - \alpha X) - \mathbf{d}\|_1 = \|F(\mathbf{p}) - \mathbf{d}\|_1 \le \|F(\mathbf{p}) - F(\mathbf{p}_0)\|_1 + \|\mathbf{d}_0 - \mathbf{d}\|_1. \quad (75)$$

There exists a constant $M > 0$ depending on $U$ such that for every $\mathbf{p} \in U$,

$$\|F(\mathbf{p}) - F(\mathbf{p}_0)\|_1 \le M\|\mathbf{p} - \mathbf{p}_0\|_1. \quad (76)$$

We define $U_2 = U_1 \cap B_1(p_0, \varepsilon_1/M)$ where

$$\varepsilon_1 := \frac{m}{4(\|\mathcal{A}'\|_{\infty} + 2z_{\max}\|\mathcal{A}''\|_{\infty})(1 + \max(|x_1|, |x_n|))}. \quad (77)$$

It follows that for every $\mathbf{p} \in U_2$ and $\mathbf{d} \in B_1(\mathbf{d}_0, \varepsilon_1)$,

$$\|F(\mathbf{p}) - F(\mathbf{p}_0)\|_1 \le \varepsilon_1, \quad |b(\mathbf{p}, h)| \le \frac{m}{2}\|h\|_2^2. \quad (78)$$

The operator $J_d$ is then strictly convex on $U_2$ since

$$\|\nabla^2 J_d(\mathbf{p})\|_2 \ge \min_{h \in \mathbb{C} \times \mathbb{R}^2, h \ne 0} \frac{|A(\mathbf{p}, h)| - |B(\mathbf{p}, h)|}{\|h\|_2^2} \ge \frac{m}{2} > 0 \quad (79)$$

We now need to prove that the minimum of $J_d$ is located inside of $U_2$ and not on its boundary. To do so, we look for a point $\mathbf{p} \in U_2$ such that $\nabla J_d(\mathbf{p}) = 0$. We already know that $\nabla J_{d_0}(\mathbf{p}_0) = 0$. Using the implicit function theorem, there exists an open set $V \subset \mathbb{R}^n$ containing $\mathbf{d}_0$ such that there exists a unique continuously differentiable function $G_1: V \to \mathbb{C} \times \mathbb{R} \times \mathbb{R}$ such that $\nabla F_d(G_1(\mathbf{d})) = 0$. We define $U' = G_1(V) \cap U_2$ and there exists $\varepsilon > 0$ such that $B_2(\mathbf{d}_0, \varepsilon) \subset G_1^{-1}(U') \cap B_1(\mathbf{d}_0, \varepsilon_1)$. It shows that the application

$$G: \left\{ \begin{array}{l} B_2(\mathbf{d}_0, \varepsilon) \rightarrow U' \\ \mathbf{d} \mapsto \underset{p \in U'}{\operatorname{argmin}} J_d(\mathbf{p}) \end{array} \right. , \quad (80)$$

is well-defined and continuously differentiable. Moreover, we also know that

$$\partial_{d_j} G(\mathbf{d}) = -[\nabla^2 J_d(G(\mathbf{d}))]^{-1}[\partial_{d_j} \nabla J_d(G(\mathbf{d}))]. \quad (81)$$

We denote $\mathbf{p}_{LS} = G(\mathbf{d})$, and

$$\|\partial_{d_j} \nabla J_d(\mathbf{p}_{LS})\|_2 = \left\| \begin{pmatrix} -\mathcal{A}(\beta_{LS} - \alpha_{LS}t_j) \\ z_{LS}t_j \mathcal{A}'(\beta_{LS} - \alpha_{LS}t_j) \\ -z_{LS} \mathcal{A}'(\beta_{LS} - \alpha_{LS}t_j) \end{pmatrix} \right\|_2 \le \sqrt{\|\mathcal{A}\|_{\infty}^2 + \alpha_{\max} \|\mathcal{A}'\|_{\infty}^2 (1 + \|\mathbf{t}\|_{\infty})} := c_1. \quad (82)$$

If follows that

$$\|\mathbf{p}_{LS} - \mathbf{p}_0\|_2 = \|G(\mathbf{d}) - G(\mathbf{d}_0)\|_2 \le \frac{2\sqrt{n}c_1}{m} \|\mathbf{d} - \mathbf{d}_0\|_2. \quad (83)$$

Finally,

$$|\Lambda(\mathbf{d}) - \Lambda(\mathbf{d}_0)| \le \frac{1}{\alpha_{\min}} |\beta_{LS} - \beta_0| + \frac{\beta_{\max}}{(\alpha_{\min})^2} |\alpha_{LS} - \beta_{LS}| \le \left( \frac{1}{\alpha_{\min}} + \frac{\beta_{\max}}{(\alpha_{\min})^2} \right) \frac{2\sqrt{n}c_1}{m} \|\mathbf{d} - \mathbf{d}_0\|_2. \quad (84)$$

□
References

[1] L. Abrahamsson. Orthogonal grid generation for two-dimensional ducts. *Journal of Computational and Applied Mathematics*, 34(3):305-314, 1991.

[2] L. Abrahamsson and H. O. Kreiss. Numerical solution of the coupled mode equations in duct acoustics. *Journal of Computational Physics*, 111(1):1-14, 1994.

[3] M. Abramowitz and I. A. Stegun. *Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables*. Applied mathematics series. Dover Publications, 1965.

[4] S. Acosta, S. Chow, J. Taylor, and V. Villamizar. On the multi-frequency inverse source problem in heterogeneous media. *Inverse Problems*, 28(7):075013, 2012.

[5] O. Balogun, T. W. Murray, and C. Prada. Simulation and measurement of the optical excitation of the s1 zero group velocity lamb wave resonance in plates. *Journal of Applied Physics*, 102(6):064914, 2007.

[6] G. Bao and P. Li. Inverse medium scattering problems for electromagnetic waves. *SIAM Journal on Applied Mathematics*, 65(6):2049-2066, 2005.

[7] G. Bao and F. Triki. Reconstruction of a defect in an open waveguide. *Science China Mathematics*, 56(12):2539-2548, 2013.

[8] G. Bao and F. Triki. Stability for the multifrequency inverse medium problem. *Journal of Differential Equations*, 269(9):7106-7128, 2020.

[9] J. P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. *Journal of Computational Physics*, 114(2):185-200, 1994.

[10] E. Bonnetier, A. Niclas, L. Seppecher, and G. Vial. The helmholtz problem in slowly varying waveguides at locally resonant frequencies. *submitted in Wave Motion*, 2022.

[11] E. Bonnetier, A. Niclas, L. Seppecher, and G. Vial. Small defects reconstruction in waveguide from multifrequency one-side scattering data. *Inverse Problems and Imaging*, 16(2):417-450, 2022.

[12] L. Bourgeois and E. Lunéville. The linear sampling method in a waveguide: A modal formulation. *Inverse Problems*, 24(1), 2008.

[13] M. Cès. *Etude théorique et expérimentale des résonances mécaniques locales de modes guidés par des structures complexes*. Theses, Université Paris-Diderot - Paris VII, 2012.

[14] P. J. Davis and P. Rabinowitz. *Methods of Numerical Integration*. Academic Press, second edition, 1984.

[15] A. Folguera and J. G. Harris. Coupled raleigh surface waves in slowly varying elastic waveguide. *Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences*, 455(1983):917-931, 1999.
[16] F. Honarvar, F. Salehi, V. Safavi, A. Mokhtari, and A. N. Sinclair. Ultrasonic monitoring of erosion/corrosion thinning rates in industrial piping systems. *Ultrasonics*, 53(7):1251-1258, 2013.

[17] V. Isakov and S. Lu. Increasing stability in the inverse source problem with attenuation and many frequencies. *SIAM J. Appl. Math.*, 78:1-18, 2018.

[18] F. Legrand. *Ondes de Lamb et Réfraction Négative*. PhD thesis, Sorbonne University, UPMC, 2020.

[19] F. Legrand, B. Gérardin, J. Laurent, C. Prada, and A. Aubry. Negative refraction of lamb modes: A theoretical study. *Physical Review B*, 98(21), 2018.

[20] Y. Y. Lu. Exact one-way methods for acoustic waveguides. *Mathematics and Computers in Simulation*, 50(5):377-391, 1999.

[21] S. Mallat. *A Wavelet Tour of Signal Processing, Chapter 1: Sparse Representations*. Academic Press, Boston, third edition, 2009.

[22] M. Norgren, M. Dalarsson, and A. Motevasselian. Reconstruction of boundary perturbations in a waveguide. *2013 International Symposium on Electromagnetic Theory*, pages 934-937, 2013.

[23] V. Pagneux and A. Maurel. Lamb wave propagation in elastic waveguides with variable thickness. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 462(2068):1315-1339, 2006.

[24] M. Sini and N. T. Thanh. Inverse acoustic obstacle scattering problems using multifrequency measurements. *Inverse Problems and Imaging*, 6(4):749-773, 2012.

[25] G. Teschl. *Ordinary differential equations and dynamical systems*, volume 140. American Mathematical Soc., 2012.