# Scattering in a partially open waveguide: the inverse problem

Laurent Bourgeois, Jean-François Fritsch, Arnaud Recoquillay

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

## To cite this version:

Laurent Bourgeois, Jean-François Fritsch, Arnaud Recoquillay. Scattering in a partially open waveguide: the inverse problem. Inverse Problems and Imaging, 2023, 17 (2), pp.463-469. 10.3934/ipi.2022052. hal-03619606

HAL Id: hal-03619606

https://inria.hal.science/hal-03619606

Submitted on 25 Mar 2022

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
# Scattering in a partially open waveguide: the inverse problem

March 25, 2022

## LAURENT BOURGEOIS

Laboratoire Poems, CNRS, INRIA, Ensta Paris, Institut Polytechnique de Paris,
828 Boulevard des Maréchaux, 91762 Palaiseau, France

## JEAN-FRANÇOIS FRITSCH

Laboratoire Poems, CNRS, INRIA, Ensta Paris, Institut Polytechnique de Paris,
828 Boulevard des Maréchaux, 91762 Palaiseau, France
and
Université Paris-Saclay, CEA, List
F-91120, Palaiseau, France

## ARNAUD RECOQUILLAY

Université Paris-Saclay, CEA, List
F-91120, Palaiseau, France

## Abstract

In this paper we consider an inverse scattering problem which consists in retrieving obstacles in a partially embedded waveguide in the acoustic case, the measurements being located on the accessible part of the structure. Such accessible part can be considered as a closed waveguide (with a finite cross section), while the embedded part can be considered as an open waveguide (with an infinite cross section). We propose an approximate model of the open waveguide by using Perfectly Matched Layers in order to simplify the resolution of the inverse problem, which is based on a modal formulation of the Linear Sampling Method. Some numerical results show the efficiency of our approach. This paper can be viewed as a continuation of the article [11], which was focused on the forward problem.

# 1 Introduction

The context of the present paper is Non Destructive Testing (NDT) of an elongated structure which is partially embedded in a surrounding medium. In the field of civil engineering, this elongated structure could be a cable made of steel which is partially free and partially surrounded by concrete. In the oil and gas industry, it could be a metallic tube which is partially free and partially immersed in water. Note that in these two examples, the celerity of waves is larger in the elongated structure (the core) than in the surrounding medium (the sheath), which is an
important feature in what follows. The case when the celerity is smaller in the core
than in the sheath would be interesting for other types of applications like optic
fibers, but we restrict to the first case in order to simplify the presentation. In the
NDT procedure we consider in this paper, the objective is to retrieve some defects
lying within the embedded part of the cable/tube from measurements located on
the accessible part of such cable/tube, this is the free one. More precisely, we wish
to send some incident waves from the accessible part of the structure and measure
the corresponding scattered waves which come back to this accessible part, then try
to identify the defects from those multistatic data. This inverse problem is delicate
for at least two reasons. The first one is it corresponds to a back-scattering problem:
the emitters and receivers are located on a single side of the area to probe. The
second one is that the waves coming from the free part of the cable/tube partially
leak in the surrounding medium, which from the point of view of the inverse problem
implies a loss of information. Indeed, since the ratio of the wave celerities between
the core and the sheath is larger than one, there are no propagating modes in the
core. On the contrary, in the absence of the surrounding medium these propagating
modes would transport the information at long distance without attenuation (see
[15]). Such delicate inverse problem, which is addressed in [34] and in [25], has
received very little attention as far as we can judge. It however deserves some
further investigation since it has important applications.

From a mathematical point of view, the configuration described above can be
viewed as a junction between a closed waveguide (with a *bounded* cross-section) and
an open waveguide (with an *unbounded* cross-section), the open waveguide having
itself a core and a sheath. In our paper, the model we consider is simplified in three
different ways. Firstly, the model is two-dimensional and corresponds to isotropic
antiplane elasticity, that is the sole SH waves are taken into account. Secondly,
the defect is a Dirichlet-type impenetrable obstacle denoted *O* and located in the
core or in the sheath of the open waveguide. Thirdly, the problem is addressed
in the time-harmonic regime at a given frequency *ω*. The configuration that we
consider is illustrated in Figure 1. The closed part of the waveguide is the domain
(-∞, 0) × (-h, h), with h > 0, while the open part of the waveguide is the domain
(0, +∞) × ℝ. The shear modulus and the density are denoted *μ* and *ρ*, respectively,
from which we define the speed *c* and the wave number *k* by the formulas *c* := *μ*/*ρ*
and *k* := *ω*/*c*. In the closed part of the waveguide, the shear modulus and the density
are constant and given by (*μ*, *ρ*) = (*μ*₀, *ρ*₀). In the open part of the waveguide, the
shear modulus and the density are piecewise positive constants and given by (*μ*, *ρ*) =
(*μ*₀, *ρ*₀) in the core, that is the domain (0, +∞) × (-h, h) and (*μ*, *ρ*) = (*μ*<sub>∞</sub>, *ρ*<sub>∞</sub>) in
the sheath, that is the domain (0, +∞) × ((-∞, -h) ∪ (h, +∞)). The corresponding
speed and wave number are denoted *c*₀ and *k*₀ in the core, while they are denoted
*c*<sub>∞</sub> and *k*<sub>∞</sub> in the sheath. Instead of considering the original configuration 1, which
is difficult to study as such due to the open part of the waveguide, a classical and
convenient idea consists following [7] in using some Perfectly Matched Layers in the
transverse direction: it has the effect to transform the open waveguide into a closed
one by introducing external layers which are made of an artificial material. Such
layers do not perturb the waves inside the physical medium but absorb the waves
before they reach the external boundaries of the layers. They are perfectly matched
in the sense that they do not introduce reflections at the interface with the physical
part of the medium. In the presence of those PMLs, the physical medium reduces
now to the domain (0, +∞) × (-h<sub>in</sub>, h<sub>in</sub>), with h<sub>in</sub> > h, and the PMLs consist of
<figure><img src="image_3.png" /><figcaption>Figure 1: Original configuration</figcaption></figure>

the domain $(0, +\infty) \times ((-h_{\text{out}}, -h_{\text{in}}) \cup (h_{\text{in}}, h_{\text{out}}))$, with $h_{\text{out}} > h_{\text{in}}$. The PMLs are characterized by a complex-valued function $\alpha$ such that $\alpha$ is equal to 1 in the physical domain and $\alpha$ is well chosen in the PMLs. The configuration with PMLs is represented in Figure 2, which will be the configuration of interest in the remainder of the paper.

In [11], for the configuration (2), the forward scattering problem which is gener-
ated by an incident mode which is sent from the left half-waveguide is extensively
analyzed: well-posedness of such problem is proved in an appropriate functional
space, the asymptotic behaviour at infinity of the solution in the longitudinal direc-
tion is studied and some artificial boundary conditions in such longitudinal direction
are proposed and analyzed in order to compute an approximate solution numeri-
cally. All these results are obtained by using the Kondratiev theory [31]. Concerning
the inverse problem, we propose to use the Linear Sampling Method to identify the
defects. The LSM appeared for the first time in [24], while a more sophisticated
variant of this method called the Factorization Method has emerged a couple of
years later [29]. First introduced in acoustics for the free homogeneous space, sam-
pling methods have been adapted to various situations ever since, both from the
point of view of the physics which governs the problem (electromagnetism [23, 30],
elasticity [1, 21], fluid mechanics [32], fluid-solid interaction [33]), the nature of the
defects (obstacles/inclusions, cracks [18]), the properties of the background medium
(anisotropic [19], periodic [3], heterogeneous and unknown [4]) and its geometry
(free space, half-space [22], shallow ocean [2], waveguides [20, 15, 33, 9, 8, 35]). By
shallow ocean, we mean a 3D domain which is bounded in only one direction and
unbounded in the other directions. By waveguide, we mean a 2D or 3D domain
which is unbounded in one direction and bounded in the other ones. The Linear
Sampling Method (LSM) consists in defining an integral operator, the kernel of
which is exactly formed by the multistatic scattering data. For each point Z of a
sampling grid, we ask whether a test function depending on Z belongs or not to
<figure><img src="image_4.png" /><figcaption>Figure 2: Configuration with PMLs</figcaption></figure>

the range of this integral operator. That test range enables us to plot the indicator function of the defects. It is important to note that, because the integral operator is compact and the data are noisy in practice, some regularization is required to apply the range test. In [15], the authors introduced a modal formulation of the LSM which consists in putting the incident and scattered fields in the form of an infinite linear combination of the modes of the waveguide, and to restrict those series to a finite number of terms which coincides with the number of propagating modes. This way of regularizing the problem is based on a physical argument which amounts to keep the sole information that propagates at long distance and to get rid of the evasentesent one. This is all the more justified as the support of data is far away from the defects. In a long series of papers [12, 16, 17, 19, 20, 22], this technique was generalized to more and more complicated situations. In particular, in [5] it was validated with experimental data in a realistic NDT context. An important point in the Linear Sampling Method is that the test function is nothing but the fundamental solution of the background domain at point Z. For simple geometries, a closed-form expression of the fundamental solution is known. But for junctions of half-waveguides such as in Figure 2, a closed-form expression of the fundamental solution is not available. This problem was addressed in [12] for a junction of several closed half-waveguides, that is with real parameters, where the fundamental solution was expressed in terms of the so-called reference fields, which were themselves computed with the help of a Finite Element Method. Importantly, those FEM computations do not depend on the sampling point Z and are performed once and for all. We reapply this technique in the present paper.

Our article is organized as follows. In section 2 we recall some important results concerning the forward problem and which are established in [11]. Section 3 is dedicated to the fundamental solution of the structure shown on Figure 2 in the absence of the obstacle. An explicit expression of such fundamental solution is required to solve the inverse problem by using the Linear Sampling Method, which is the subject of section 4. Some numerical experiments showing the feasibility of our method are presented in section 5. Lastly, an appendix details the derivation of the fundamental solution in a uniform PML-waveguide, which is a first step to derive the fundamental solution of the overall structure.
# 2 An overview of the forward problem

## 2.1 Description of the modes

In this paragraph we recall the expression of the modes, on the one hand in a straight homogeneous closed waveguide having the properties of the left half-waveguide of Figure 2, on the other hand in a straight stratified closed waveguide with PMLs having the properties of the right half-waveguide of Figure 2.

Let us first denote $\tilde{\Omega} := \mathbb{R} \times I$, $I := (-h, h)$. Such closed waveguide is characterized by the constant material properties $(\mu_0, \rho_0)$, from which we deduce the celerity $c_0 := \sqrt{\mu_0/\rho_0}$ and the wave number $k_0 := \omega/c_0$. We introduce the coordinates $(x, y)$, where $x$ is the coordinate along the unbounded direction of the waveguide, while $y$ is the coordinate along the transverse direction. The modes are the solutions $u$ of the form $u(x, y) = e^{\lambda x} \varphi(y)$ for some $\lambda \in \mathbb{C}$ to the problem

$$ \begin{cases} -\Delta u - k_0^2 u = 0 & \text{in } \tilde{\Omega} \\ \partial_{\nu} u = 0 & \text{on } \partial \tilde{\Omega}, \end{cases} $$

where $\nu$ is the outward unit normal vector to $\partial \tilde{\Omega}$. The modes are given by

$$ \tilde{w}_n^{\pm}(x, y) := e^{\pm \tilde{\lambda}_n x} \tilde{\varphi}_n(y), \quad (1) $$

with

$$ \tilde{\lambda}_n := i \sqrt{k_0^2 - \frac{n^2 \pi^2}{4h^2}}, \quad n \in \mathbb{N}, \quad (2) $$

where the branch cut in the square root is chosen such that $\Im(\sqrt{z}) \ge 0$, and with

$$ \tilde{\varphi}_0(y) := \frac{1}{\sqrt{2h}}, \quad \tilde{\varphi}_n(y) := \frac{1}{\sqrt{h}} \cos \left( \frac{n\pi}{2h} (y+h) \right), \quad n \in \mathbb{N}^*, \quad (3) $$

the family of functions $\tilde{\varphi}_n$, $n \in \mathbb{N}$, forming a complete orthonormal basis of $L^2(I)$.

Let us introduce the following assumption on $k_0$, which implies that none of the $\tilde{\lambda}_n$ vanish.

**Assumption 2.1.** For all $n \in \mathbb{N}$, $n\pi/2h \ne k_0$.

This implies in particular that there exists some $\tilde{N} \in \mathbb{N}^*$ such that for $n = 0, \dots, \tilde{N}-1$, $\tilde{\lambda}_n \in i(0, +\infty)$, while for $n = \tilde{N}, \dots, +\infty$, $\tilde{\lambda}_n \in (-\infty, 0)$. Hence for $n = 0, \dots, \tilde{N}-1$, the $\tilde{w}_n^{\pm}$ are called propagating modes, the $\tilde{w}_n^{+}$ propagating from the left to the right, the $\tilde{w}_n^{-}$ propagating from the right to the left. For $n = \tilde{N}, \dots, +\infty$, the $\tilde{w}_n^{\pm}$ are called the evanescent modes, the $\tilde{w}_n^{+}$ exponentially decaying from the left to the right, the $\tilde{w}_n^{-}$ exponentially decaying from the right to the left.

We now describe the modes of the closed waveguide with PMLs, that is $\Omega_{\text{out}} := \mathbb{R} \times I_{\text{out}}$, with $I_{\text{out}} := (-h_{\text{out}}, h_{\text{out}})$. The shear modulus and the density are given in the truncated domain $\Omega_{\text{out}}$ by

$$ (\mu, \rho)(y) := \begin{cases} (\mu_0, \rho_0) & \text{if } |y| < h, \\ (\mu_\infty, \rho_\infty) & \text{if } h < |y| < h_{\text{out}}. \end{cases} \quad (4) $$

From $\mu_\infty$ and $\rho_\infty$ we deduce the celerity $c_\infty := \sqrt{\mu_\infty/\rho_\infty}$ and the wave number $k_\infty := \omega/c_\infty$. We assume in what follows that the celerity in the core is larger than the celerity in the sheath, that is:
**Assumption 2.2.** We have that $c_0 > c_\infty$ or in other words $k_0 < k_\infty$.

The PML complex function $\alpha$ is defined in $\Omega_{\text{out}}$ by

$$ \alpha(y) := \begin{cases} 1 & \text{if } |y| < h_{\text{in}}, \\ \alpha_{\infty}(y) & \text{if } h_{\text{in}} < |y| < h_{\text{out}}, \end{cases} \qquad (5) $$

where $\alpha_\infty$ is a function which satisfies

$$ -\frac{\pi}{2} < \arg(\alpha_{\infty}(y)) < 0 \quad \text{for} \quad h_{\text{in}} < |y| < h_{\text{out}}. \qquad (6) $$

The properties (5) (6) of the function $\alpha$ are essential for the PMLs to play their role (see [11]). The modes are the solutions $u$ of the form $u(x, y) = e^{\lambda x} \varphi(y)$ to the problem

$$ \left\{ \begin{array}{ll} -\alpha \partial_y (\alpha \mu \partial_y u) - \mu \partial_{xx} u - \omega^2 \rho u = 0 & \text{in } \Omega_{\text{out}} \\ \partial_\nu u = 0 & \text{on } \partial \Omega_{\text{out}}, \end{array} \right. \qquad (7) $$

where $\nu$ is the outward unit normal vector to $\partial\Omega_{\text{out}}$. There exists a countable set of complex numbers $\lambda_n$ and a countable family of functions $\varphi_n$ on $I_{\text{out}}$ such that the modes are given for $n \in \mathbb{N}$ by

$$ w_n^{\pm}(x, y) = e^{\pm \lambda_n x} \varphi_n(y). \qquad (8) $$

More precisely, if we consider the operator $\mathcal{L}(\lambda) : H^1(I_{\text{out}}) \to H^1(I_{\text{out}})^*$, where $H^1(I_{\text{out}})^*$ is the topological dual of $H^1(I_{\text{out}})$, which is defined by

$$ \langle \mathcal{L}(\lambda)\varphi, \psi \rangle_{I_{\text{out}}} := \int_{I_{\text{out}}} \left( \alpha \mu d_y \varphi d_y \psi - \frac{\mu}{\alpha}(\lambda^2 + k^2)\varphi\psi \right) dy, \quad \forall \varphi, \psi \in H^1(I_{\text{out}}), \quad (9) $$

and the bracket $\langle \cdot, \cdot \rangle_{I_{\text{out}}}$ means duality between $H^1(I_{\text{out}})$ and $H^1(I_{\text{out}})^*$, the complex numbers $\pm \lambda_n$ and the corresponding functions $\varphi_n$ coincide with the pairs $(\lambda, \varphi) \in \mathbb{C} \times H^1(I_{\text{out}}) \setminus \{0\}$ such that

$$ \mathcal{L}(\lambda)\varphi = 0. \qquad (10) $$

The set of $\lambda \in \mathbb{C}$ such that there exists some non trivial $\varphi \in H^1(I_{\text{out}})$ satisfying (10) is denoted $\Lambda$. The eigenvalues $\lambda_n, n \in \mathbb{N}$, satisfy $\Re(\lambda_n) < 0$ and $\Im(\lambda_n) > 0$ and are numbered such that

$$ \dots \le \Re(\lambda_{n+1}) \le \Re(\lambda_n) \le \dots \le \Re(\lambda_0) < 0. \qquad (11) $$

Contrary to the previous homogeneous case (see the expression of the numbers $\tilde{\lambda}_n$ given by (2)), the numbers $\lambda_n$ do not have a simple expression: they are solutions to a dispersion relationship which is given in [11] (see Remark 3.1). The modes given by (8) are divided into the so-called leaky modes, which are localized in the physical part of the waveguide, and into the PML modes, which are localized in the PMLs. The physical meaning of those two kinds of modes are described in [11]. In particular, the leaky modes inherit from the fact that in the true unbounded configuration of Figure 1, some modes radiate from the core to the sheath in the transverse direction. It should be noted that the modes $w_n^+$ are exponentially decaying from the left to the right while the modes $w_n^-$ are exponentially decaying from the right to the left. In other words, and contrary to the homogeneous waveguide, all the modes $w_n^\pm$ are evanescent, as can be seen on Figure 5 of [11]. That none of the modes is
propagating is a consequence of Assumption 2.2. Depending on the parameters of the sheath and the core, some leaky modes may have a very small real part (it will be the case in our numerical section where the core is made of steel and the sheath is made of concrete), which implies that they are almost propagating. That will be useful for the resolution of the inverse problem. Let us introduce the following fundamental assumption on the eigenfunctions $\varphi$.

**Assumption 2.3.** The eigenvectors $\varphi$ are such that

$$ \int_{I_{\text{out}}} \frac{\mu}{\alpha} \varphi^2 dy \neq 0. $$

It is proved in [11] (see lemma 2 and 3) that under the Assumption 2.3, the geometrical and algebraical multiplicities of all eigenvalues $\lambda$ are both equal to 1. It is unknown to us whether the $\varphi_n$, $n \in \mathbb{N}$, form a complete basis of $L^2(I_{\text{out}})$ or not. However, again thanks to Assumption 2.3 (see proposition 4 in [11]), the eigenvectors $\varphi_n$ satisfy the biorthogonality relationship: for $n, m \in \mathbb{N}$, we have

$$ \int_{I_{\text{out}}} \frac{\mu}{\alpha} \varphi_n \varphi_m dy = J_n \delta_{mn}, \qquad (12) $$

where $\delta_{mn} = 1$ if $m = n$ and $\delta_{mn} = 0$ otherwise, which implicitly defines $J_n$.

## 2.2 The forward problem

We consider the problem illustrated by Figure 2, which models the diffraction by an obstacle of an incident wave coming from the left in a junction between a half-closed waveguide and a half-open waveguide closed by finite PMLs. We introduce the following notations: $\Omega^- := (-\infty, 0) \times (-h, h)$, $\Omega^+ := (0, +\infty) \times (-h_{\text{out}}, h_{\text{out}})$, $\Sigma_0 := \{0\} \times (-h, h)$ and $\Omega := \Omega^- \cup \Sigma_0 \cup \Omega^+$. The obstacle $O$ is a smooth bounded domain which lies either within the core or within the sheath of the open half-waveguide with PMLs, and we denote $D := \Omega \setminus \bar{O}$ and $D^+ := \Omega^+ \setminus \bar{O}$. Let us define the incident wave $u^i = \tilde{w}_{n,0}^+$ in $\Omega$ as

$$ \tilde{w}_{n,0}^+ := \begin{cases} \tilde{w}_n^+ & \text{in } \Omega^-, \\ 0 & \text{in } \Omega^+, \end{cases} \qquad (13) $$

where $\tilde{w}_n^+$ for $n \in \mathbb{N}$ is a mode coming from the left closed waveguide and defined by (1). Let us introduce the spaces $\tilde{H}_{\text{loc}}^1(D)$ (resp. $\tilde{L}_{\text{loc}}^2(D)$) as the set of distributions $v$ in $D$ such that $\chi v \in H^1(D)$ (resp. $L^2(D)$), for all $\chi \in \mathcal{C}^\infty(\mathbb{R}^2)$ vanishing as $x \to -\infty$. The scattering problem we consider is: find $u \in \tilde{H}_{\text{loc}}^1(D)$ such that

$$ \left\{ \begin{array}{ll} -\Delta u - k_0^2 u = 0 & \text{in } \Omega^-, \\ -\partial_y(\alpha \mu \partial_y u) - \frac{\mu}{\alpha} \partial_{xx} u - \frac{\mu}{\alpha} k_0^2 u = 0 & \text{in } D^+, \\ [u] = 0 & \text{on } \Sigma_0, \\ [\partial_x u] = 0 & \text{on } \Sigma_0, \\ \partial_\nu u = 0 & \text{on } \partial\Omega, \\ u = 0 & \text{on } \partial O, \\ u - u^i & \text{is outgoing,} \end{array} \right. \qquad (14) $$

where $u^i = \tilde{w}_{n,0}^+$ and $[\cdot]$ denotes the jump of the solution at the interface $\Sigma_0$. The obstacle $O$ is hence given by a Dirichlet-type boundary condition. As shown in [11],
the radiation condition on the right side of $D$ consists in the choice of a $H^1$-type functional space in view of the absence of propagating modes. On the left side of $D$, such radiation condition amounts to forcing the solution to be outgoing, that is to assuming the existence of a sequence of complex numbers $a_n^-$ such that in $\Omega^-$:

$$u - u^i = \sum_{n \in \mathbb{N}} a_n^- \tilde{w}_n^-.$$

It is proved in [11] (see Theorem 5.1) that if the scattering problem (14) satisfies the uniqueness property, that is $u = 0$ if $u^i = 0$, then it is well-posed. In the remainder of the paper, we assume that such uniqueness property is satisfied.

### 3 The fundamental solution

We first introduce the notion of reference field, which in the context of the inverse problem will play a crucial role to compute the fundamental solution in $\Omega$. Let us define the space $\tilde{H}_{\text{loc}}^1(\Omega)$ (resp. $\tilde{L}_{\text{loc}}^2(\Omega)$) as the set of distributions $v$ in $\Omega$ such that $\chi v \in H^1(\Omega)$ (resp. $L^2(\Omega)$), for all $\chi \in \mathcal{C}^\infty(\mathbb{R}^2)$ vanishing as $x \to -\infty$. For $n \in \mathbb{N}$, the reference fields $\tilde{r}_n$ and $r_n$ satisfy the following transmission problems: find $u \in \tilde{H}_{\text{loc}}^1(\Omega)$ such that

$$ \left\{ \begin{array}{ll} -\Delta u - k_0^2 u = 0 & \text{in } \Omega^-, \\ -\partial_y(\alpha \mu \partial_y u) - \frac{\mu}{\alpha} \partial_{xx} u - \frac{\mu}{\alpha} k^2 u = 0 & \text{in } \Omega^+, \\ \llbracket u \rrbracket = 0 & \text{on } \Sigma_0, \\ \llbracket \partial_x u \rrbracket = 0 & \text{on } \Sigma_0, \\ \partial_\nu u = 0 & \text{on } \partial\Omega, \\ u - u^i & \text{is outgoing,} \end{array} \right. \quad (15) $$

where the incident field $u^i$ is set either to the field $\tilde{w}_{n,0}^+$ in $\Omega$ given by (13) or the field $w_{n,0}^-$ in $\Omega$ defined by

$$ w_{n,0}^- := \begin{cases} 0 & \text{in } \Omega^-, \\ w_n^- & \text{in } \Omega^+, \end{cases} \quad (16) $$

and $w_n^-$ is given by (8). The fields $\tilde{r}_n$ and $r_n$ are hence the total fields associated with the incident fields $\tilde{w}_{n,0}^+$ and $w_{n,0}^-$ which come from the left and from the right, respectively, the scattered fields being generated by the sole junction of the two half-waveguides (that is in the absence of the defects). We have the following result, which states that the fields $\tilde{r}_n$ and $r_n$ are well-defined if the uniqueness property is satisfied. The proof is omitted since it is almost the same as the one of Theorem 5.1 in [11].

**Theorem 3.1.** Assume that the problem (15) for $u^i = 0$ only admits the trivial solution. Then the problems (15) associated with $u^i = \tilde{w}_{n,0}^+$ and $u^i = w_{n,0}^-$, $n \in \mathbb{N}$, have a unique solution.

Secondly, we introduce the fundamental solutions of the straight waveguides $\tilde{\Omega}$ and $\Omega_{\text{out}}$, the first one being the homogeneous closed waveguide, the second one being the stratified waveguide closed by PMLs. The fundamental solution of the
closed waveguide $\tilde{\Omega}$ is, for $M' = (x', y') \in \tilde{\Omega}$, the solution $\tilde{G}(\cdot; M')$ to the following problem

$$ \left\{ \begin{array}{l} -\mu_0(\Delta\tilde{G} + k_0^2\tilde{G}) = \delta_{M'} \quad \text{in } \tilde{\Omega}, \\ \partial_\nu\tilde{G} = 0 \quad \text{on } \partial\tilde{\Omega}, \\ \tilde{G} \text{ is outgoing.} \end{array} \right. \qquad (17) $$

It is well-known that the function $\tilde{G}(\cdot; M')$ is well-defined in $L^2_{\text{loc}}(\tilde{\Omega})$ and is given, for $M = (x, y)$, by

$$ \tilde{G}(M; M') = - \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n \mu_0} e^{\tilde{\lambda}_n |x-x'|} \tilde{\varphi}_n(y) \tilde{\varphi}_n(y'), \qquad (18) $$

where the $\tilde{\lambda}_n$ and the $\tilde{\varphi}_n$ are given by (2) and (3), respectively. The fundamental solution in the stratified waveguide $\Omega_{\text{out}}$ closed by PMLs is, for $M' \in \Omega_{\text{out}}$, the solution $G(\cdot; M')$ to the following problem

$$ \left\{ \begin{array}{l} -d_y(\alpha \mu d_y G) - \frac{\mu}{\alpha} d_{xx} G - \frac{\mu}{\alpha} k^2 G = \delta_{M'} \quad \text{in } \Omega_{\text{out}}, \\ \partial_\nu G = 0 \quad \text{on } \partial\Omega_{\text{out}}. \end{array} \right. \qquad (19) $$

In order to prove well-posedness of problem (19), we introduce a more general one: for a source term $f$ defined in $\Omega_{\text{out}}$, find the solution $u$ defined in $\Omega_{\text{out}}$ such that

$$ \left\{ \begin{array}{l} -d_y(\alpha \mu d_y u) - \frac{\mu}{\alpha} d_{xx} u - \frac{\mu}{\alpha} k^2 u = f \quad \text{in } \Omega_{\text{out}}, \\ \partial_\nu u = 0 \quad \text{on } \partial\Omega_{\text{out}}. \end{array} \right. \qquad (20) $$

Let us define the operator $A : H^1(\Omega_{\text{out}}) \to H^1(\Omega_{\text{out}})^*$ such that for all $(u, v) \in H^1(\Omega_{\text{out}}) \times H^1(\Omega_{\text{out}})$,

$$ \langle Au, v \rangle_{\Omega_{\text{out}}} := \int_{\Omega_{\text{out}}} \left( \alpha \mu \partial_y u \partial_y v + \frac{\mu}{\alpha} \partial_x u \partial_x v - \frac{\mu}{\alpha} k^2 u v \right) dxdy, \qquad (21) $$

where $\langle \cdot, \cdot \rangle_{\Omega_{\text{out}}}$ is the duality bracket between $H^1(\Omega_{\text{out}})^*$ and $H^1(\Omega_{\text{out}})$. It was proved in [11] that $A$ is an isomorphism, in other words that for $f \in H^1(\Omega_{\text{out}})^*$, the problem (20) has a unique solution $u$ in $H^1(\Omega_{\text{out}})$. We first complement such result by a regularity result. In this view we need to introduce the Fourier-Laplace transform $\mathcal{F}_{x \to \lambda}$ which is given, for $v \in S'(\mathbb{R})$ (here, $S'(\mathbb{R})$ denotes the classical space of tempered distributions) and $\lambda \in i\mathbb{R}$, by

$$ \hat{v}(\lambda) = (\mathcal{F}_{x \to \lambda} v)(\lambda) := \int_{-\infty}^{+\infty} e^{-\lambda x} v(x) dx. $$

The inverse $\mathcal{F}_{x \to \lambda}^{-1}$ is given by

$$ \mathcal{F}_{x \to \lambda}^{-1} \hat{v}(x) = \frac{1}{2\pi i} \int_{i\mathbb{R}} e^{\lambda x} \hat{v}(\lambda) d\lambda. $$

A well known property of $\mathcal{F}$ is that for all $v \in S'(\mathbb{R})$ and $\lambda \in i\mathbb{R}$ we have

$$ \hat{v}' = \lambda \hat{v}. $$
Using this property, one can prove that for $m \in \mathbb{N}$ and two reals $a < b$, the standard norm of $H^m(\mathbb{R} \times (a, b))$ is equivalent to the norm

$$ \|v\| := \left( \frac{1}{i} \int_{i\mathbb{R}} \|\hat{v}(\lambda, \cdot)\|_{H^m((a,b),|\lambda|)}^2 d\lambda \right)^{1/2}, $$

where

$$ \|\varphi\|_{H^m((a,b),|\lambda|)}^2 := \|\varphi\|_{H^m((a,b))}^2 + |\lambda|^{2m} \|\varphi\|_{L^2((a,b))}^2. \quad (22) $$

Let us denote $H_{\#}^2(\Omega_{\text{out}})$ the subset of functions $u \in H^1(\Omega_{\text{out}})$ such that for each $j \in \{-2, -1, 0, 1, 2\}$, the restrictions $u|_{\Omega_j}$ belong to $H^2(\Omega_j)$, where $\Omega_{-2} := \mathbb{R} \times (-h_{\text{out}}, -h_{\text{in}})$, $\Omega_{-1} := \mathbb{R} \times (-h_{\text{in}}, -h)$, $\Omega_0 := \mathbb{R} \times (-h, h) = \tilde{\Omega}$, $\Omega_1 := \mathbb{R} \times (h, h_{\text{in}})$ and $\Omega_2 := \mathbb{R} \times (h_{\text{in}}, h_{\text{out}})$, and such that the function $\alpha\mu\partial_y u$ is continuous across all the horizontal lines which separate these domains $\Omega_j$, $j \in \{-2, -1, 0, 1, 2\}$.

**Proposition 1.** For $f \in L^2(\Omega_{\text{out}})$, the problem (20) has a unique solution in $H_{\#}^2(\Omega_{\text{out}})$.

*Proof.* By applying $\mathcal{F}_{x \to \lambda}$ to the equation $Au = f$, one obtains that $\mathcal{L}(\lambda)\hat{u}(\lambda, \cdot) = \hat{f}(\lambda, \cdot)$ for $\lambda \in i\mathbb{R}$, so that one is naturally led to study the symbol $\mathcal{L}(\lambda)$ defined by (9). By using Lemma 4.1 in [11], there is $\tau_0 > 0$ such that for $\lambda = i\tau$, $\tau \in \mathbb{R}$ with $|\tau| \ge \tau_0$, $\mathcal{L}(\lambda) : H^1(I_{\text{out}}) \to H^1(I_{\text{out}})^*$ is an isomorphism. In particular, denoting

$$ \|g\|_{H^1(I_{\text{out}}, |\lambda|)^*} := \sup_{\psi \in H^1(I_{\text{out}}) \setminus \{0\}} \frac{|\langle g, \bar{\psi} \rangle_{I_{\text{out}}}|}{\|\psi\|_{H^1(I_{\text{out}}, |\lambda|)}}, \quad \forall g \in H^1(I_{\text{out}})^*, \quad (23) $$

where $\langle \cdot, \cdot \rangle_{I_{\text{out}}}$ denotes the duality pairing between $H^1(I_{\text{out}})^*$ and $H^1(I_{\text{out}})$, $\varphi := \mathcal{L}(\lambda)^{-1}g$ satisfies

$$ \|\varphi\|_{H^1(I_{\text{out}}, |\lambda|)} \le C \|g\|_{H^1(I_{\text{out}}, |\lambda|)^*}, \quad (24) $$

where $C > 0$ is independent of $g$ and $\lambda$. Since $i\mathbb{R} \notin \Lambda$, the operator $\mathcal{L}(\lambda)$ is injective for $\lambda \in i\mathbb{R}$, and the Fredholm analytic theorem guarantees that the operator $\mathcal{L}(\lambda)^{-1}$ is well defined and continuous for $\lambda \in i\mathbb{R}$. This implies that the estimate (24) is also valid for $\lambda$ in the vertical segment $i[-\tau_0, \tau_0]$, hence for all $\lambda \in i\mathbb{R}$, with a new constant $C > 0$ which depends neither on $g$ nor on $\lambda$. Integrating such estimate on $i\mathbb{R}$ implies that if $f \in L^2(\Omega_{\text{out}})$, then the solution to the problem $Au = f$, given by

$$ u(x, \cdot) := \frac{1}{2\pi i} \int_{i\mathbb{R}} e^{\lambda x} \mathcal{L}(\lambda) \hat{f}(\lambda, \cdot) d\lambda, $$

belongs to $H^1(\Omega_{\text{out}})$. Let us prove that $u \in H^2(\Omega_0)$. For $g \in L^2(I_{\text{out}})$, we deduce from (23) that for all $\lambda \in i\mathbb{R}$,

$$ (1 + |\lambda|^2) \|g\|_{H^1(I_{\text{out}}, |\lambda|)^*}^2 \le \|g\|_{L^2(I_{\text{out}})}^2. \quad (25) $$

The estimates (24) and (25) imply that for all $g \in L^2(I_{\text{out}})$, $\varphi$ satisfies for all $\lambda \in i\mathbb{R}$,

$$ \|\varphi\|_{H^1(-h, h)}^2 + |\lambda|^4 \|\varphi\|_{L^2(-h, h)}^2 \le C \|g\|_{L^2(I_{\text{out}})}^2. \quad (26) $$

In $(-h, h)$, the function $\varphi$ satisfies the equation

$$ -d_{yy}^2 \varphi - (\lambda^2 + k_0^2) \varphi = \frac{g}{\mu_0}, $$
which implies that there exists some constant $C > 0$ such that

$$ \|d_{yy}^2 \varphi\|_{L^2(-h,h)}^2 \le C(1 + |\lambda|^2)^2 \|\varphi\|_{L^2(-h,h)}^2 + C \|g\|_{L^2(-h,h)}^2. $$

We conclude from the above inequality and (26) that

$$ \|\varphi\|_{H^2((-h,h),|\lambda|)}^2 = \|\varphi\|_{H^2(-h,h)}^2 + |\lambda|^4 \|\varphi\|_{L^2(-h,h)}^2 \le C \|g\|_{L^2(I_{\text{out}})}^2, $$

and by integration on $i\mathbb{R}$, that $u \in H^2(\Omega_0)$. The same reasoning applies to all domains $\Omega_j, j \in \{-2, -1, 0, 1, 2\}$, which completes the proof. $\square$

**Remark 1.** Due to the discontinuity of the function $\partial_y u$ across the horizontal lines which separate the domains $\Omega_j, j \in \{-2, -1, 0, 1, 2\}$, the function $u$ does not belong to $H^2(\Omega_{\text{out}})$.

**Remark 2.** Lemma 4.1 in [11] relies on the fact that, due to (5) and (6), we have $-\pi/2 < \arg(\alpha_\alpha) \le 0$, which implies that $\Re(e(\alpha)) \ge C$ and $\Re(1/\alpha) \ge C$ for $C > 0$. Since such property is also satisfied by $\bar{\alpha}$, Proposition 1 holds if we replace coefficient $\alpha$ by $\bar{\alpha}$ in problem (20).

We are now in a position to prove well-posedness of problem (19).

**Theorem 3.2.** For all $M' \in \Omega_{\text{out}}$, the problem (19) has a unique solution in $L^2(\Omega_{\text{out}})$.

*Proof.* From Proposition 1 we observe that $u \in H^{1+\beta}(\Omega_{\text{out}})$, for all $\beta \in [0, 1/2)$. Indeed, since $u \in H_{\#}^2(\Omega_{\text{out}})$, we have in particular $u|_{\Omega_j} \in H^{1+\beta}(\Omega_j)$ for all $j \in \{-2, -1, 0, 1, 2\}$, and thus $\partial_x u, \partial_y u \in H^\beta(\Omega_j)$ for all $j \in \{-2, -1, 0, 1, 2\}$. We know from [26] that for $\beta \in [0, 1/2)$ the extension by 0 from $\Omega_j$ to $\Omega_{\text{out}}$ of any function in $H^\beta(\Omega_j)$ belongs to $H^\beta(\Omega_{\text{out}})$. As a consequence $\partial_x u, \partial_y u \in H^\beta(\Omega_{\text{out}})$, for all $\beta \in [0, 1/2)$. From a classical Sobolev embedding theorem in two dimensions we have that for $\beta \in (0, 1/2)$, the functions in $H^{1+\beta}(\Omega_{\text{out}})$ are continuous in $\Omega_{\text{out}}$. Hence for all $M' \in \Omega_{\text{out}}$, $\delta_{M'} \in H_{\#}^2(\Omega_{\text{out}})^*$. Proposition 1 amounts to say that the operator $A : H_{\#}^2(\Omega_{\text{out}}) \to L^2(\Omega_{\text{out}})$ is an isomorphism. In view of Remark 2, we have that $\bar{A} : H_{\#}^2(\Omega_{\text{out}}) \to L^2(\Omega_{\text{out}})$ is also an isomorphism. As a consequence, its adjoint operator $\bar{A}^* : L^2(\Omega_{\text{out}}) \to H_{\#}^2(\Omega_{\text{out}})^*$ is an isomorphism and since $A$ is a symmetric operator, namely $\bar{A}^* = A$, we conclude that $A$ is an isomorphism from $L^2(\Omega_{\text{out}})$ to $H_{\#}^2(\Omega_{\text{out}})^*$. That $\delta_{M'} \in H_{\#}^2(\Omega_{\text{out}})^*$ completes the proof. $\square$

**Remark 3.** It should be noted that Theorem 3.2 is in particular valid for a source point $M'$ located at the interface of two different domains $\Omega_j, j \in \{-2, -1, 0, 1, 2\}$.

As shown in the Appendix, the function $G(\cdot; M')$ is given by

$$ G(M; M') = - \sum_{n \in \mathbb{N}} \frac{1}{2\lambda_n J_n} e^{\lambda_n |x-x'|} \varphi_n(y) \varphi_n(y'), \quad (27) $$

where the pairs $(\lambda_n, \varphi_n)$ are the solutions to problem (10), the $\varphi_n$ satisfying (12) and the $J_n$ being defined in Assumption 2.3. The fundamental solution of the junction
of the two half-waveguides in the absence of defects is, for $M' \in \Omega$, the solution $\mathcal{G}(\cdot; M')$ to the following problem

$$ egin{cases} -\mu_0(\Delta\mathcal{G} + k_0^2\mathcal{G}) = f & \text{in } \Omega^{-}, \\ -\partial_y(\alpha\mu\partial_y\mathcal{G}) - \frac{\mu}{\alpha}\partial_{xx}u - \frac{\mu}{\alpha}k^2\mathcal{G} = g & \text{in } \Omega^{+}, \\ \llbracket \mathcal{G} \rrbracket = 0 & \text{on } \Sigma_0, \\ \llbracket \partial_x\mathcal{G} \rrbracket = 0 & \text{on } \Sigma_0, \\ \partial_\nu\mathcal{G} = 0 & \text{on } \partial\Omega, \\ \mathcal{G} \text{ is outgoing,} \end{cases} \quad (28) $$

where $(f,g) = (\delta_{M'}, 0)$ if $M' \in \Omega^-$ and $(f,g) = (0, \delta_{M'})$ if $M' \in \Omega^+$. We have the following result.

**Proposition 2.** If the problem (28) for $(f,g) = (0,0)$ only has the trivial solution, the problem (28) has a unique solution in $L^2_{\text{loc}}(\Omega)$ which is given by the following formulas:

* for $M' = (x', y') \in \Omega^-$ and $M = (x, y) \in \Omega$,

$$ \mathcal{G}(M; M') = \begin{cases} \tilde{G}(M; M') - \frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{w}_n^-(M')(\tilde{r}_n(M) - \tilde{w}_n^+(M)) & \text{for } x < x' \\ -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{w}_n^-(M')\tilde{r}_n(M) & \text{for } x > x', \end{cases} \quad (29) $$

* for $M' = (x', y') \in \Omega^+$ and $M = (x, y) \in \Omega$,

$$ \mathcal{G}(M; M') = \begin{cases} -\sum_{n \in \mathbb{N}} \frac{1}{2\lambda_n J_n} w_n^+(M') r_n(M) & \text{for } x < x' \\ G(M; M') - \sum_{n \in \mathbb{N}} \frac{1}{2\lambda_n J_n} w_n^+(M')(r_n(M) - w_n^-(M)) & \text{for } x > x', \end{cases} \quad (30) $$

where $\tilde{r}_n$ and $r_n$ are the reference fields given by the system (15) for $u^i = \tilde{w}_{n,0}^+$ and $u^i = w_{n,0}^-$, respectively. In addition, we have the symmetry relationship:

$$ \mathcal{G}(M; M') = \mathcal{G}(M'; M), \quad \forall M, M' \in \Omega. \quad (31) $$

*Proof.* Let us consider the case when $M' \in \Omega^-$ and let us define $\tilde{G}^i(\cdot, M')$ as

$$ \tilde{G}^i(M; M') := \begin{cases} \tilde{G}(M; M') & \text{if } M \in \Omega^-, \\ 0 & \text{if } M \in \Omega^+. \end{cases} \quad (32) $$

We use the decomposition $\mathcal{G}(\cdot; M') = \tilde{G}^i(\cdot; M') + \tilde{G}^s(\cdot; M')$, where $\tilde{G}^i(\cdot; M')$ plays the role of an incident field and $\tilde{G}^s(\cdot; M')$ plays the role of a scattered field by the
junction. It is readily seen that the scattered field $\tilde{G}^s(\cdot; M')$ satisfies the transmission problem

$$ egin{cases} 
\begin{aligned}
-\Delta \tilde{G}^s - k_0^2 \tilde{G}^s &= 0 & \text{in } \Omega^-, \\
-\partial_y(\alpha \mu \partial_y \tilde{G}^s) - \frac{\mu}{\alpha} \partial_{xx} \tilde{G}^s - \frac{\mu}{\alpha} k^2 \tilde{G}^s &= 0 & \text{in } \Omega^+, \\
\llbracket \tilde{G}^s \rrbracket &= \tilde{G} & \text{on } \Sigma_0, \\
\llbracket \partial_x \tilde{G}^s \rrbracket &= \partial_x \tilde{G} & \text{on } \Sigma_0, \\
\partial_\nu \tilde{G}^s &= 0 & \text{on } \partial\Omega, \\
\tilde{G}^s & \text{is outgoing,}
\end{aligned}
\end{cases} \tag{33} $$

where $\llbracket \cdot \rrbracket$ is the jump of the solution at the interface $\Sigma_0$ from the left to the right. When the problem (28) for $(f,g) = (0,0)$ only has the trivial solution, by proceeding as in [11], we obtain that the problem (33) is well-posed in $\tilde{H}_{\text{loc}}^1(\Omega)$.

Now, using the same decomposition of the total field $\tilde{r}_n = \tilde{w}_{n,0}^+ + \tilde{r}_n^s$ into the incident field $\tilde{w}_{n,0}^+$ and the scattered field $\tilde{r}_n^s$, such scattered field is solution to the transmission problem

$$ egin{cases} 
\begin{aligned}
-\Delta \tilde{r}_n^s - k_0^2 \tilde{r}_n^s &= 0 & \text{in } \Omega^-, \\
-\partial_y(\alpha \mu \partial_y \tilde{r}_n^s) - \frac{\mu}{\alpha} \partial_{xx} \tilde{r}_n^s - \frac{\mu}{\alpha} k^2 \tilde{r}_n^s &= 0 & \text{in } \Omega^+, \\
\llbracket \tilde{r}_n^s \rrbracket &= \tilde{w}_n^+ & \text{on } \Sigma_0, \\
\llbracket \partial_x \tilde{r}_n^s \rrbracket &= \partial_x \tilde{w}_n^+ & \text{on } \Sigma_0, \\
\partial_\nu \tilde{r}_n^s &= 0 & \text{on } \partial\Omega, \\
\tilde{r}_n^s & \text{is outgoing.}
\end{aligned}
\end{cases} \tag{34} $$

Since the section $\Sigma_0$ is located on the right of the source point $M'$ of the fundamental solution $\tilde{G}$, the formula (18) reduces for $x > x'$ to

$$ \tilde{G}(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} e^{\tilde{\lambda}_n x} e^{-\tilde{\lambda}_n x'} \tilde{\varphi}_n(y) \tilde{\varphi}_n(y'), $$

that is

$$ \tilde{G}(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{w}_n^+(M) \tilde{w}_n^-(M'). \tag{35} $$

Comparing the systems (33) and (34), by linearity we obtain that for all $M \in \Omega$,

$$ \tilde{G}^s(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{r}_n^s(M) \tilde{w}_n^-(M'), $$

that is

$$ \mathcal{G}(M; M') - \tilde{G}(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} (\tilde{r}_n(M) - \tilde{w}_n^+(M)) \tilde{w}_n^-(M') \tag{36} $$

for $M \in \Omega^-$ and

$$ \mathcal{G}(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{r}_n(M) \tilde{w}_n^-(M') \tag{37} $$

for $M \in \Omega^+$. It remains to remark that for $M \in \Omega^-$ and $x > x'$, due to the simplification of the fundamental solution $\tilde{G}$ given by (35), the expression (36) simply becomes (37).
To address the case when $M' \in \Omega^+$, we have to define $G^i(\cdot, M')$ as

$$ G^i(M; M') := \begin{cases} 0 & \text{if } M \in \Omega^-, \\ G(M; M') & \text{if } M \in \Omega^+, \end{cases} \qquad (38) $$

and to use the decomposition $\mathcal{G}(\cdot; M') = G^i(\cdot; M') + \mathcal{G}^s(\cdot; M')$, where $G^i(\cdot; M')$ plays the role of an incident field and $\mathcal{G}^s(\cdot; M')$ plays the role of a scattered field. The derivation of (30) for $M' \in \Omega^+$ follows the same lines as in the case when $M' \in \Omega^-$. In order to obtain the symmetry relationship, it suffices to rewrite the problem (28) as

$$ \left\{ \begin{array}{lll} -\operatorname{div}(A\nabla\mathcal{G}) - \frac{\mu}{\alpha}k^2\mathcal{G} & = & \delta_{M'} \quad \text{in } \Omega, \\ \partial_{\nu}\mathcal{G} & = & 0 \quad \text{on } \partial\Omega, \\ \mathcal{G} & \text{is} & \text{outgoing,} \end{array} \right. \quad A := \begin{pmatrix} \mu/\alpha & 0 \\ 0 & \alpha\mu \end{pmatrix}, $$

with $(\alpha, \mu, k) = (1, \mu_0, k_0)$ in $\mathbb{R} \times (-h, h)$, $(\alpha, \mu, k) = (1, \mu_\infty, k_\infty)$ in $(0, +\infty) \times ((-h_{\text{in}}, -h) \cup (h, h_{\text{in}}))$ and $(\alpha, \mu, k) = (\alpha_\infty, \mu_\infty, k_\infty)$ in $(0, +\infty) \times ((-h_{\text{out}}, -h_{\text{in}}) \cup (h_{\text{in}}, h_{\text{out}}))$. Since the differential operator $P := -\operatorname{div}(A\nabla\cdot) - \frac{\mu}{\alpha}k^2$ is symmetric (in particular the matrix $A$ is symmetric), the fundamental solution $\mathcal{G}$ satisfies (31), following a general and classical result on fundamental solutions. $\square$

# 4 The Linear Sampling Method

Let us come back to the situation presented in the introduction and illustrated by Figure 2, that is the junction of the left closed half-waveguide and the right half-waveguide closed by PMLs, in the presence of the Dirichlet obstacle $O$ in the physical part of the right half-waveguide. More precisely, such obstacle is either embedded in the core, that is $\bar{O} \subset (0, +\infty) \times (-h, h)$ or embedded in the sheath, that is $\bar{O} \subset (0, +\infty) \times ((-h_{\text{in}}, -h) \cup (h, h_{\text{in}}))$. Our method would enable us to address the case when the obstacle $O$ lies in the left half-waveguide, but we do not consider this case which is less interesting for the applications.

For $M' \in \Omega^-$, we consider the following forward problem: find $u(\cdot; M') \in \tilde{L}_{\text{loc}}^2(D)$ such that

$$ \left\{ \begin{array}{lll} -\mu_0(\Delta u + k_0^2)u & = & \delta_{M'} \quad \text{in } \Omega^-, \\ -\partial_y(\alpha\mu\partial_y u) - \frac{\mu}{\alpha}\partial_{xx}u - \frac{\mu}{\alpha}k^2 u & = & 0 \quad \text{in } D^+, \\ [[u]] & = & 0 \quad \text{on } \Sigma_0, \\ [[\partial_x u]] & = & 0 \quad \text{on } \Sigma_0, \\ \partial_{\nu}u & = & 0 \quad \text{on } \partial\Omega, \\ u & = & 0 \quad \text{on } \partial O, \\ u & \text{is} & \text{outgoing.} \end{array} \right. \qquad (39) $$

The field $u(\cdot; M')$ can be seen as a total field having the decomposition $u(\cdot; M') = \mathcal{G}(\cdot; M') + u^s(\cdot; M')$, where the fundamental solution $\mathcal{G}(\cdot; M')$ of the junction $\Omega$, which is the solution to problem (28) with $(f, g) = (\delta_{M'}, 0)$, is the incident field and the field $u^s(\cdot; M')$ is the resulting scattered field by the defect. Such field is solution
to the following problem: find $u^s(\cdot; M') \in \tilde{H}_{\text{loc}}^1(D)$ such that

$$ egin{cases} 
-\Delta u^s - k_0^2 u^s = 0 & \text{in } \Omega^{-}, \\
-\partial_y(\alpha \mu \partial_y u^s) - \frac{\mu}{\alpha} \partial_{xx} u^s - \frac{\mu}{\alpha} k^2 u^s = 0 & \text{in } D^{+}, \\
\llbracket u^s \rrbracket = 0 & \text{on } \Sigma_0, \\
\llbracket \partial_x u^s \rrbracket = 0 & \text{on } \Sigma_0, \\
\partial_\nu u^s = 0 & \text{on } \partial\Omega, \\
u^s = F & \text{on } \partial O, \\
u^s \text{ is outgoing,} & 
\end{cases} 
\quad (40) $$

for $F = -\mathcal{G}(\cdot; M')|_{\partial O}$. For all $n \in \mathbb{N}$, we will also need to consider the solutions $\tilde{u}_n^s$ to the same problem (40) for $F = -\tilde{r}_n|_{\partial O}$, where $\tilde{r}_n$ is a reference field. It is readily seen that $\tilde{u}_n^s$ satisfies such problem if and only if $(\tilde{u}_n^s + \tilde{r}_n)$ satisfies problem (14). From Theorem 5.1 in [11], if we assume that the problem (14) when replacing the incident wave by 0 only has the trivial solution, the problem (14) is well-posed, which implies that the problem (40) for $F = -\tilde{r}_n|_{\partial O}$ is also well-posed. Proceeding similarly, we can prove that the problem (39) and the problem (40) for $F = -\mathcal{G}(\cdot; M')|_{\partial O}$ are well-posed.

Let us now introduce the inverse problem. The support of the measurements is the transverse section $\Sigma_{-R} := \{-R\} \times (-h, h)$ with $R > 0$, which is a section of the left closed half-waveguide. We assume that for all $M' \in \Sigma_{-R}$, we measure the corresponding scattered field $u^s(M; M')$, that is the solution to the problem (40), for all $M \in \Sigma_{-R}$. The goal of the inverse problem is to identify the defect $O$ from those multistatic data. Since the support of both the sources and the receivers are located on a single side of the defect, our data can be viewed as back-scattering data. Let us try to identify the obstacle $O$ from the data $u^s(M; M')$ for $(M, M') \in \Sigma_{-R} \times \Sigma_{-R}$ with the help of the Linear Sampling Method, which is based on the near-field operator:

$$ egin{cases} 
\mathcal{N} : L^2(\Sigma_{-R}) \to L^2(\Sigma_{-R}) \\
h \mapsto \mathcal{N}h, \quad (\mathcal{N}h)(M) = \int_{\Sigma_{-R}} u^s(M; M')h(M') ds(M'), \quad M \in \Sigma_{-R}, 
\end{cases} 
\quad (41) $$

where $u^s(\cdot; M')$ is the solution to problem (40) for $F = -\mathcal{G}(\cdot; M')|_{\partial O}$. We will also need the operator

$$ egin{cases} 
\mathcal{H} : L^2(\Sigma_{-R}) \to H^{1/2}(\partial O) \\
h \mapsto (\mathcal{H}h)(M) = \int_{\Sigma_{-R}} \mathcal{G}(M; M')h(M') ds(M'), \quad M \in \partial O. 
\end{cases} 
\quad (42) $$

The Linear Sampling Method is justified by the following theorem, the proof of which mimics the one proved in [14]. It clarifies the reason why for all point $Z$ of a sampling grid, we will approximately solve the near-field equation

$$ \mathcal{N}h = \mathcal{G}(\cdot; Z)|_{\Sigma_{-R}}. \quad (43) $$

**Theorem 4.1.** We assume that the exterior problems (40) are well-posed and that the interior problem: find $v \in H^1(O)$ such that

$$ egin{cases} 
-\Delta v - k^2 v = 0 & \text{in } O \\
v = 0 & \text{on } \partial O 
\end{cases} $$
only has the trivial solution. Let $\mathcal{N}$ and $\mathcal{H}$ be the operators defined by (41) and (42), respectively.

* If $Z \in O$, then for all $\varepsilon > 0$ there exists a solution in $L^2(\Sigma_{-R})$ parametrized by $Z$ and denoted $h_{\varepsilon}(\cdot; Z)$, of the inequality

$$ \|\mathcal{N}h_{\varepsilon}(\cdot; Z) - \mathcal{G}(\cdot; Z)\|_{L^2(\Sigma_{-R})} \le \varepsilon $$

such that the function $\mathcal{H}h_{\varepsilon}(\cdot; Z)$ converges in $H^{1/2}(\partial O)$ as $\varepsilon \to 0$. Furthermore, for a given fixed $\varepsilon$, the function $h_{\varepsilon}(\cdot; Z)$ satisfies

$$ \lim_{Z \to \partial O} \|h_{\varepsilon}(\cdot; Z)\|_{L^2(\Sigma_{-R})} = +\infty \quad \text{and} \quad \lim_{Z \to \partial O} \|\mathcal{H}h_{\varepsilon}(\cdot; Z)\|_{H^{1/2}(\partial O)} = +\infty. $$

* If $Z \in \Omega \setminus \overline{O}$, then every solution $h_{\varepsilon}(\cdot; Z) \in L^2(\Sigma_{-R})$ of the inequality

$$ \|\mathcal{N}h_{\varepsilon}(\cdot; Z) - \mathcal{G}(\cdot; Z)\|_{L^2(\Sigma_{-R})} \le \varepsilon $$

satisfies

$$ \lim_{\varepsilon \to 0} \|h_{\varepsilon}(\cdot; Z)\|_{L^2(\Sigma_{-R})} = +\infty \quad \text{and} \quad \lim_{\varepsilon \to 0} \|\mathcal{H}h_{\varepsilon}(\cdot; Z)\|_{H^{1/2}(\partial O)} = +\infty. $$

The Linear Sampling Method consists then, for all $Z$ of a sampling grid of $\Omega$, in solving a regularized version of the near-field equation (43). Such regularization is required by the fact that the operator $\mathcal{N}$ is compact. Following [15], we introduce a modal formulation of the Linear Sampling Method: the principle is to project such near-field equation on the complete basis $(\tilde{\varphi}_n)_{n \in \mathbb{N}}$ of the transverse section $\Sigma_{-R}$. This enables us to propose a “physical regularization” which consists in replacing the series which result from these projections by the sum of their first $\tilde{N}$ terms. This amounts to keep, among the information contained in the incident and scattered waves, their propagating parts only, in other words to neglect their evanescent parts. In this view we will use, for $n \in \mathbb{N}$, the solution $\tilde{u}_n^s \in \dot{H}_{\text{loc}}^1(D)$ to the problem (40) for $F = -\tilde{r}_n|_{\partial O}$. We have the following lemma.

**Lemma 4.2.** Denoting again $M = (x, y)$ and $M' = (x', y')$, for $x' = -R$ and $x > -R$,

$$ u^s(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{u}_n^s(M) \tilde{w}_n^{-}(M'). $$

*Proof.* From Proposition 2, since $M' \in \Omega^-$ and $x > x'$, we have

$$ \mathcal{G}(M; M') = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{r}_n(M) \tilde{w}_n^{-}(M'). $$

Since $u^s(\cdot; M')$ is the solution to the problem (40) for $F = -\mathcal{G}(\cdot; M')|_{\partial O}$ and the $\tilde{u}_n^s$ are the solutions to the same problem (40) for $F = -\tilde{r}_n|_{\partial O}$, the result follows by linearity. $\square$

Now, let us project the near-field equation (43) in the complete basis $(\tilde{\varphi}_n)$, $n \in \mathbb{N}$, of the transverse section $\Sigma_{-R}$. We have the following decomposition.
**Proposition 3.** For all $h \in L^2(\Sigma_{-R})$ and $Z \in \Omega^+$, we have

$$ \mathcal{N}h = -\frac{1}{\mu_0} \sum_{m,n \in \mathbb{N}} \frac{e^{\tilde{\lambda}_n R}}{2\tilde{\lambda}_n} U_{mn} h_n \tilde{\varphi}_m \quad (44) $$

and

$$ \mathcal{G}(\cdot; Z)|_{\Sigma_{-R}} = - \sum_{m,n \in \mathbb{N}} \frac{e^{\lambda_n x_Z}}{2\lambda_n J_n} \varphi_n(y_Z)(r_n, \tilde{\varphi}_m)_{L^2(\Sigma_{-R})} \tilde{\varphi}_m, \quad (45) $$

where we use the decompositions

$$ \tilde{u}_n^s|_{\Sigma_{-R}} = \sum_{m \in \mathbb{N}} U_{mn} \tilde{\varphi}_m, \quad \forall n \in \mathbb{N}, \quad \text{and} \quad h = \sum_{m \in \mathbb{N}} h_m \tilde{\varphi}_m. $$

An alternative formula of $\mathcal{G}(\cdot; Z)|_{\Sigma_{-R}}$ is

$$ \mathcal{G}(\cdot; Z)|_{\Sigma_{-R}} = -\frac{1}{\mu_0} \sum_{m \in \mathbb{N}} \frac{e^{\tilde{\lambda}_m R}}{2\tilde{\lambda}_m} \tilde{r}_m(Z) \tilde{\varphi}_m. \quad (46) $$

*Proof.* Since for all $M \in \Sigma_{-R}$,

$$ (\mathcal{N}h)(M) = \int_{\Sigma_{-R}} u^s(M; M')h(M')ds(M'), $$

by using Lemma 4.2, we obtain

$$ (\mathcal{N}h)(M) = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_n} \tilde{u}_n^s(M) \int_{\Sigma_{-R}} \tilde{w}_n^-(M')h(M')ds(M'). $$

That $\tilde{w}_n^-(x', y') = e^{-\tilde{\lambda}_n x'} \tilde{\varphi}_n(y')$, $h = \sum_{m \in \mathbb{N}} h_m \tilde{\varphi}_m$ and $(\tilde{\varphi}_m)_m$ is a complete basis of $L^2(\Sigma_{-R})$, imply

$$ (\mathcal{N}h)(M) = -\frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{e^{\tilde{\lambda}_n R}}{2\tilde{\lambda}_n} h_n \tilde{u}_n^s(M). $$

The expression (44) is obtained by decomposing the trace on $\Sigma_{-R}$ of the fields $\tilde{u}_n^s$ in the basis $(\tilde{\varphi}_n)_n$ of $L^2(\Sigma_{-R})$. In order to obtain (45), we start from Proposition 2: in the case when $Z \in \Omega^+$ and $M \in \Omega^-$, we have

$$ \mathcal{G}(M; Z) = - \sum_{n \in \mathbb{N}} \frac{1}{2\lambda_n J_n} w_n^+(Z) r_n(M). $$

Using both $w_n^+(Z) = e^{\lambda_n x_Z} \varphi_n(y_Z)$ and the decomposition of the field $r_n$ in the basis $(\tilde{\varphi}_n)_n$, yields (45). To obtain (46), we use the symmetry relationship $\mathcal{G}(M; Z) = \mathcal{G}(Z; M)$ and Proposition 2 again: in the case when $Z \in \Omega^+$ and $M \in \Omega^-$, we have

$$ \mathcal{G}(Z; M) = -\frac{1}{\mu_0} \sum_{m \in \mathbb{N}} \frac{1}{2\tilde{\lambda}_m} \tilde{w}_m^-(M) \tilde{r}_m(Z) $$

and (46) follows from the fact that $\tilde{w}_m^-(x, y) = e^{-\tilde{\lambda}_m x} \tilde{\varphi}_m(y)$. $\square$
**Remark 4.** It is important to note from (44) that it is equivalent to measure $u^s(M; M')$ for all $(M, M') \in \Sigma_{-R} \times \Sigma_{-R}$ and to measure the projections of the fields $\tilde{u}_n^s$ on the transverse functions $\tilde{\varphi}_m$ of $\Sigma_{-R}$, for all $m, n \in \mathbb{N}$. This is why in the numerical section, the data will consist of those projections.

From Proposition 3, the near-field equation (43) is equivalent to the infinite systems: find $h = \sum_{n \in \mathbb{N}} h_n \tilde{\varphi}_n \in L^2(\Sigma_{-R})$ such that for all $m \in \mathbb{N}$,

$$ \frac{1}{\mu_0} \sum_{n \in \mathbb{N}} \frac{e^{\tilde{\lambda}_n R}}{\tilde{\lambda}_n} U_{mn} h_n = \sum_{p \in \mathbb{N}} \frac{e^{\lambda_p x_Z}}{\lambda_p J_p} \varphi_p(y_Z)(r_p, \tilde{\varphi}_m)_{L^2(\Sigma_{-R})} \quad (47) $$

or

$$ \sum_{n \in \mathbb{N}} \frac{e^{\tilde{\lambda}_n R}}{\tilde{\lambda}_n} U_{mn} h_n = \frac{e^{\tilde{\lambda}_m R}}{\tilde{\lambda}_m} \tilde{r}_m(Z). \quad (48) $$

**Remark 5.** The main difference between the near-field equations (47) and (48) is that the reference fields $r_n$ are required in the first case, that is the response of the structure without any obstacle to a mode coming from the right, while the reference fields $\tilde{r}_n$ are required in the second case, that is the response of the structure without any obstacle to a mode coming from the left. It is important to note that, for both the near-field equations (47) and (48), the reference fields $r_n$ and $\tilde{r}_n$ are computed once and for all during an offline step: those computations do not depend on the sampling point $Z$. The system (48) is particularly interesting because its right-hand side does not require the truncation of a series. Such system is very similar to the one derived for the back-scattering imaging in a straight closed waveguide, for which the reference field $\tilde{r}_m$ reduces to a right-going mode.

# 5 Numerical experiments

In this numerical section, we compute some artificial data on the section $\Sigma_{-R}$ by choosing some obstacle $O$ and by solving the corresponding forward scattering problems (14) with the help of a finite element method. In order to reduce the unbounded domain to a computational bounded domain without changing the solution, we make use of transparent boundary conditions in the longitudinal direction on each side of the obstacle $O$. Those artificial boundary conditions are based on Dirichlet-to-Neumann operators. More precisely, we use a “thin” DtN operator on a transverse section located on the left waveguide and a “thick” DtN with an overlap between two transverse sections located on the right waveguide, as described in the section 5 of [11]. The reference fields, which are the solution to the problem (15) for $u^i = \tilde{w}_{n,0}^+$, are computed using the same finite element method and the same artificial boundary conditions. To be more precise, concerning the finite elements themselves, we use standard P2 triangular elements, the mesh size being such that a wavelength contains about 10 discretization points. We also have to compute the eigenpairs $(\lambda_n, \varphi_n)$ which satisfy (10), which is done by using one-dimensional P1 finite elements based on a very refined mesh size of $10^{-5} m$. The DtN operators are approximated with the help of a series which is truncated to a finite number of terms, as described in the section 6 of [11]. In the following numerical computations, we have chosen $R = 17 \times 10^{-2} m$. The parameters which characterize the overall structure are the same as in [11], that is $h = 5 \times 10^{-2} m$, $h_{in} = 7.5 \times 10^{-2} m$ and
$h_{\text{out}} = 12.5 \times 10^{-2} m$. The core is made of steel, with $\mu_0 = 84.298 \times 10^9 Pa$ and $\rho_0 = 7932 kg.m^{-3}$, while the sheath is made of concrete, with $\mu_\infty = 15.908 \times 10^9 Pa$ and $\rho_\infty = 2300 kg.m^{-3}$. Those values ensure that $c_0 > c_\infty$, which is consistent with Assumption 2.2. All the numerical experiments were performed with the help of the FEM-BEM library XLiFE++ (see [28]). In order to make our numerical experiment closer to a real NDT experiment, we produce noisy data by combining two ideas. First of all, the forward scattering problem (14) is solved by using PMLs which are different from those which are used to compute the reference fields involved in the inverse problem. In the first case, the function $\alpha_\infty$ involved in (5) is given by the parabolic profile

$$ \alpha_{\infty}(y) = \frac{1}{1 + 3(3 + 4i)(|y| - h_{\text{in}})^2/h_{\text{out}}^2}, $$

while in the second case, $\alpha_\infty$ is the constant $e^{-i\pi/3}$. Note that a comparison between those two profile functions is proposed in [11]. Since the PMLs are part of the inverse procedure, in some sense this idea enables us to avoid committing an inverse crime. Secondly, our data are artificially corrupted following the method described in [15]. We consider the trace on the measurement section $\Sigma_{-R}$ of a scattered field $\tilde{u}_n^s$ obtained by solving (14) for $F = -\tilde{r}_n|_{\partial O}$. For simplicity, let us denote $\tilde{u}_n^s|_{\Sigma_{-R}} = U$. Decomposing $\Sigma_{-R}$ into a finite number of intervals, we compute a pointwise Gaussian noise $B$ and form a noisy data $U_\delta$ on $\Sigma_{-R}$ by setting

$$ U_{\delta} := U + \alpha B, $$

where the number $\alpha > 0$ is calibrated such that

$$ \|U_{\delta} - U\|_{L^2(\Sigma_{-R})} = \sigma \|U\|_{L^2(\Sigma_{-R})}, $$

the number $\sigma$ being a given relative amplitude of noise. Now let us explain how we discretize the near-field equations (47) and (48) in view of the numerical implementation. We recall that the infinite systems (47) and (48) are ill-posed since the underlying operator $\mathcal{N}$ we “invert” is compact. In practice it is natural to approximate the infinite systems (47) and (48) by restricting $m$ and $n$ to $0, \dots, \tilde{N} - 1$, where we recall that $\tilde{N}$ is the number of propagating modes in the left closed half-waveguide, and by restricting $p$ to the first $N$ terms, where no obvious value of $N$ is available. In order to specify $N$, we can view the modes in the waveguide $\Omega_{\text{out}}$ as the modes in the waveguide $\tilde{\Omega}$ perturbed by the surrounding medium. This is why we choose $N = \tilde{N}$, that is in $\Omega^+$ we restrict to the leaky modes which correspond to the propagating modes in $\Omega^-$. In addition to the truncation of the series, we regularize the finite systems with the help of the Tikhonov regularization combined with the Morozov principle to determine the regularization parameter as a function of the amplitude of noise. In fact we exactly follow the procedure described in [15]. This additional regularization is all the more necessary as the number of propagating modes in the left closed half-waveguide, that is $\tilde{N}$, is large. Calling $h(\cdot; Z)$ the regularized solution to either problem (47) or problem (48) (note that such solution implicitly depends on the sampling point $Z$), we hereafter plot the function

$$ \psi(Z) := \log_{10} \left( \frac{1}{\|h(\cdot; Z)\|_{L^2(\Sigma_{-R})}} \right), \qquad (49) $$

which in view of Theorem 4.1 is bounded inside the obstacle $O$ and almost $-\infty$ outside. In the numerical experiments hereafter, we analyze the effect of several
parameters on the quality of the identification. In all the following pictures, the contour of the obstacle is represented by a black solid line.

## 5.1 About the choice of the near-field equation

We first consider an obstacle in the core of the structure. Here, the frequency is $\omega = 12.56 \text{ rad.s}^{-1}$, which corresponds to $N = 13$ propagating modes, while the relative amplitude of noise is $\sigma = 10\%$. On Figure 3, we show the reconstructed obstacles with the help of the function (49), by using both near-field equations (47) and (48). Those reconstructions have to be compared with the one obtained with a wrong near-field equation, that is by doing as if the surrounding medium were neglected in the right part of the structure, in other words as if the global structure did consist of the uniform closed waveguide $\tilde{\Omega}$. We observe that both near-field equations (47) and (48) approximately produce the same images, which unsurprisingly are better than the image obtained with the wrong near-field equation. In all further identification results, the near-field equation (48) is chosen.

<figure><img src="image_5.png" /><figcaption>Figure 3: Choice of the near-field equation. Top left: near-field equation (47). Top right: near-field equation (48). Bottom: wrong near-field equation associated with the uniform closed waveguide</figcaption></figure>

## 5.2 About the distance between the defect and the interface

Let us recall that the scattering fields $u_n^s$ satisfy the problem (14). In view of Theorem 14 in [11] and the properties of the eigenvalues $\lambda_n$ given by (11), those fields are exponentially decaying when $x \to +\infty$. We hence expect that the quality of the identification result decreases when the distance between the interface between the two half-waveguides and the obstacle increases. In order to illustrate this property, we show on Figure 4 how the quality of the identification is affected by such distance, either for a single obstacle or for two obstacles (one of them is fixed). Here the frequency is $\omega = 21.99 \text{ rad.s}^{-1}$, that is the number of propagating modes is $N = 22$, and the relative amplitude of noise is $\sigma = 10\%$. We notice that the left part of the obstacle is better reconstructed than its right part, which is of course due to our back-scattering situation (measurements are located only on the left transverse section $\Sigma_{-R}$). When one obstacle is close and the other is far from the interface, it
seems that the latter is invisible while the former is visible: this is due to the high
contrast between the values of the indicator function $\psi$.

<figure><img src="image_6.png" /><figcaption>Figure 4: Influence of the distance between the defects and the interface. Left: one obstacle. Right: two obstacles</figcaption></figure>

## 5.3 About the amplitude of noise

Let us illustrate the impact of the amplitude of noise on the reconstruction. In the case when $\omega = 21.99 \text{ rad.s}^{-1}$, which corresponds to $N = 22$, we test different values of the amplitude of noise $\sigma$, that is $\sigma = 1\%$, $\sigma = 10\%$ and $\sigma = 20\%$. The results are shown on Figure 5. We observe that, as always for the type of Gaussian pointwise noise we introduce, the reconstructions with the Linear Sampling Method are very robust with respect to the amplitude of noise.

## 5.4 About the frequency

Then we analyze the impact of the frequency $\omega$ on the reconstruction, which amounts to study the impact of the number $N$ of propagating modes. The amplitude of noise being $\sigma = 10\%$, we compare in Figure 6 the results obtained for $\omega = 6.28 \text{ rad.s}^{-1}$, $\omega = 12.56 \text{ rad.s}^{-1}$ and $\omega = 25.13 \text{ rad.s}^{-1}$, which correspond to $N = 7, N = 13$ and $N = 26$, respectively. Clearly, the higher is the frequency, the better is the reconstruction, provided the size of the mesh used to solve the forward problems is chosen accordingly.
<figure><img src="image_7.png" /><figcaption>Figure 5: Impact of the amplitude of noise. Top left: 1% noise. Top right: 10% noise. Bottom: 20% noise.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Impact of the frequency. Top left: $N = 7$. Top right: $N = 13$. Bottom: $N = 26$.</figcaption></figure>

## 5.5 About obstacles in the sheath

Up to now we have considered examples of obstacles embedded in the core of the right half-waveguide. For $\omega = 21.99 \text{ rad.s}^{-1}$ ($N = 22$) and $\sigma = 10\%$, we show in Figure 7 the identification results obtained with one or two obstacles, one of them at least being embedded in the sheath of the right half-waveguide. The obstacles located in the sheath are not as well retrieved as the ones located in the core. The reason is that when an obstacle is located in the core, its reconstruction benefits from the reflections of waves between the core and the sheath, while it does not when it lies in the sheath (waves leak at infinity).
<figure><img src="image_9.png" /><figcaption>Figure 7: One obstacle at least is located in the sheath</figcaption></figure>

# 6 Appendix

In order to derive an explicit expression of the solution $G(\cdot; M')$ to problem (19), let us consider the operator $L: D(L) \subset L^2(I_{\text{out}}) \to L^2(I_{\text{out}})$ defined by

$$ \left\{ \begin{array}{l} L\varphi := -\frac{\alpha}{\mu}d_y(\alpha\mu d_y\varphi) - k^2\varphi \\ D(L) := \{\varphi \in H^1(I_{\text{out}}), \alpha\mu d_y\varphi \in H^1(I_{\text{out}}), d_y\varphi(-h_{\text{out}}) = d_y\varphi(h_{\text{out}}) = 0\}. \end{array} \right. \qquad (50) $$

Finding the pairs $(\lambda, \varphi) \in \mathbb{C} \times H^1(I_{\text{out}}) \setminus \{0\}$ satisfying (10) is equivalent to finding the eigenvalues $\gamma = \lambda^2$ and the corresponding eigenfunctions $\varphi$ of operator $L$. It is readily seen that $\mathcal{L}(\lambda) = (\mu/\alpha)(L - \lambda^2)$, hence the inverse of $\mathcal{L}(\lambda)$ is related to the resolvent of $L$ by the relationship

$$ \mathcal{L}(\lambda)^{-1} = (L - \lambda^2)^{-1} \left( \frac{\alpha}{\mu} \right). \qquad (51) $$

Using [27] and the fact that the geometric and algebraic multiplicities of all eigenvalues $\gamma_n = \lambda_n^2$ of $L$ are equal to 1, we have that in a vicinity of $\gamma_n$,

$$ (L - \gamma)^{-1} = -\frac{P_n}{\gamma - \gamma_n} + U_n(\gamma), \qquad (52) $$

where $P_n$ is the projector on the eigenfunction $\varphi_n$ and $U_n$ is a holomorphic function. From (51) and (52), we get that in a vicinity of $\lambda_n$,

$$ \mathcal{L}(\lambda)^{-1} = -\frac{P_n\left(\frac{\alpha}{\mu}\right)}{\lambda^2 - \lambda_n^2} + U_n(\lambda), \qquad (53) $$
where the function $U_n$ is holomorphic.

In [11] it is proved that the operator $\mathcal{L}(\lambda)$ is invertible for all $\lambda \in i\mathbb{R}$, hence for a source term $f$ given in $\Omega_{\text{out}}$, the unique solution to the problem (20) is given by

$$u(x, y) = \frac{1}{2\pi i} \int_{i\mathbb{R}} e^{\lambda x} \mathcal{L}(\lambda)^{-1} \hat{f}(\lambda, y) d\lambda. \quad (54)$$

We apply such formula for $f = \delta_{M'}$. The idea is now to evaluate the solution $u$ by applying the residue theorem to a well-chosen contour in the complex plane. Without loss of generality, we shall assume that $x' = 0$, since the general case is easily obtained by a translation along the $x$-axis. In the case when $x > 0$, let us consider the oriented contour $C_R = i[-R, R] \cup \Gamma_R$, where $\Gamma_R = \{Re^{i\theta}, \theta \in [\pi/2, 3\pi/2]\}$ (see Figure 8).

<figure><img src="image_10.png" /><figcaption>Figure 8: Contour $C_R$ for the residue theorem</figcaption></figure>

Since the $\lambda$ which belong to the bounded domain $D_R$ delimited by $C_R$ satisfy $\Re(\lambda) < 0$, we obtain

$$u(x, y) = \lim_{R \to +\infty} \frac{1}{2\pi i} \int_{C_R} e^{\lambda x} \mathcal{L}(\lambda)^{-1} \hat{f}(\lambda, y) d\lambda.$$

We observe that $F_{x,y}(\lambda) := e^{\lambda x} \mathcal{L}(\lambda)^{-1} \hat{f}(\lambda, y)$ is a meromorphic function, the poles of which are given by the eigenvalues of the symbol $\mathcal{L}$ in the half-plane $\Re(\lambda) < 0$, that is the $\lambda_n$. By the residue theorem, we hence have

$$\frac{1}{2\pi i} \int_{C_R} e^{\lambda x} \mathcal{L}(\lambda)^{-1} \hat{f}(\lambda, y) d\lambda = \sum_{\lambda_n \in D_R} \text{Res}(F_{x,y}(\lambda), \lambda_n).$$

For $f = \delta_{0,y'}$, we have that $\hat{f}(\lambda, y) = \delta_{y'}(y)$ and by using the biorthogonality relationship (12), the projector $P_n$ is given by

$$P_n\varphi = \frac{1}{J_n} \left( \int_{I_{\text{out}}} \frac{\mu}{\alpha} \varphi_n \varphi dy \right) \varphi_n.$$

We conclude in view of (53) that

$$\text{Res}(F_{x,y}(\lambda), \lambda_n) = -\frac{1}{2\lambda_n} e^{\lambda_n x} P_n \left( \frac{\alpha}{\mu} \hat{f}(\lambda_n, y) \right) = -\frac{1}{2\lambda_n J_n} e^{\lambda_n x} \varphi_n(y') \varphi_n(y),$$
so that

$$G(x, y; 0, y') = - \sum_{n \in \mathbb{N}} \frac{1}{2\lambda_n J_n} e^{\lambda_n x} \varphi_n(y) \varphi_n(y').$$

If $x' > x$, it suffices to replace $x$ in the above right-hand side by $(x-x')$. For $x' < x$, the same method can be applied, which will involve the eigenvalues $-\lambda_n$ of $\mathcal{L}(\lambda)$ in the half-plane $\Re(\lambda) > 0$. We finally obtain the symmetric function

$$G(x, y; x', y') = - \sum_{n \in \mathbb{N}} \frac{1}{2\lambda_n J_n} e^{\lambda_n |x-x'|} \varphi_n(y) \varphi_n(y').$$

## References

[1] T. Arens. Linear sampling methods for 2D inverse elastic wave scattering. *Inverse Problems*, 17(5):1445–1464, 2001.

[2] T. Arens, D. Gintides, and A. Lechleiter. Direct and inverse medium scattering in a three-dimensional homogeneous planar waveguide. *SIAM J. Appl. Math.*, 71(3):753–772, 2011.

[3] T. Arens and A. Kirsch. The factorization method in inverse scattering from periodic structures. *Inverse Problems*, 19(5):1195–1211, 2003.

[4] L. Audibert, A. Girard, and H. Haddar. Identifying defects in an unknown background using differential measurements. *Inverse Probl. Imaging*, 9(3):625–643, 2015.

[5] V. Baronian, L. Bourgeois, B. Chapuis, and A. Recoquillay. Linear sampling method applied to non destructive testing of an elastic waveguide: theory, numerics and experiments. *Inverse Problems*, 34(7):075006, 34, 2018.

[6] V. Baronian, L. Bourgeois, and A. Recoquillay. Imaging an acoustic waveguide from surface data in the time domain. *Wave Motion*, 66:68–87, 2016.

[7] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. *J. of Comput. Phys.*, 114(2):185–200, 1994.

[8] Liliana Borcea, Fioralba Cakoni, and Shixu Meng. A direct approach to imaging in a waveguide with perturbed geometry. *Journal of Computational Physics*, 392:556–577, 2019.

[9] Liliana Borcea and Shixu Meng. Factorization method versus migration imaging in a waveguide. *Inverse Problems*, 35(12):124006, 2019.

[10] L. Bourgeois and S. Fliss. On the identification of defects in a periodic waveguide from far field data. *Inverse Problems*, 30(9):095004, 31, 2014.

[11] L. Bourgeois, S. Fliss, J.-F. Fritsch, C. Hazard, and A. Recoquillay. Scattering in a partially open waveguide: the forward problem. preprint, October 2021.

[12] L. Bourgeois, J.-F. Fritsch, and A. Recoquillay. Imaging junctions of waveguides. *Inverse Problems and Imaging*, 15(2):285–314, 2021.
[13] L. Bourgeois, F. Le Louër, and E. Lunéville. On the use of Lamb modes in the linear sampling method for elastic waveguides. *Inverse Problems*, 27(5):055001, 27, 2011.

[14] L Bourgeois and E Lunéville. The linear sampling method in a waveguide: A formulation based on modes. *Journal of Physics: Conference Series*, 135:012023, nov 2008.

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

[16] L. Bourgeois and É. Lunéville. On the use of sampling methods to identify cracks in acoustic waveguides. *Inverse Problems*, 28(10):105011, 18, 2012.

[17] L. Bourgeois and E. Lunéville. On the use of the linear sampling method to identify cracks in elastic waveguides. *Inverse Problems*, 29(2):025017, 19, 2013.

[18] F. Cakoni and D. Colton. The linear sampling method for cracks. *Inverse Problems*, 19(2):279-295, 2003.

[19] F. Cakoni, D. Colton, and H. Haddar. The linear sampling method for anisotropic media. *J. Comput. Appl. Math.*, 146(2):285-299, 2002.

[20] A. Charalambopoulos, D. Gintides, K. Kiriaki, and A. Kirsch. The factorization method for an acoustic wave guide. In *Mathematical methods in scattering theory and biomedical engineering*, pages 120-127. World Sci. Publ., Hackensack, NJ, 2006.

[21] A. Charalambopoulos, A. Kirsch, K. A. Anagnostopoulos, D. Gintides, and K. Kiriaki. The factorization method in inverse elastic scattering from penetrable bodies. *Inverse Problems*, 23(1):27-51, 2007.

[22] D. Colton, J. Coyle, and P. Monk. Recent developments in inverse acoustic scattering theory. *SIAM Rev.*, 42(3):369-414, 2000.

[23] D. Colton, H. Haddar, and P. Monk. The linear sampling method for solving the electromagnetic inverse scattering problem. *SIAM J. Sci. Comput.*, 24(3):719-731, 2002.

[24] D. Colton and A. Kirsch. A simple method for solving inverse scattering problems in the resonance region. *Inverse Problems*, 12(4):383-393, 1996.

[25] M. Gallezot. *Numerical modelling of non-destructive testing of buried waveguides*. Theses, École centrale de Nantes, November 2018.

[26] P. Grisvard. *Singularities in boundary value problems*, volume 22 of *Recherches en Mathématiques Appliquées* [*Research in Applied Mathematics*]. Masson, Paris; Springer-Verlag, Berlin, 1992.

[27] T. Kato. *Perturbation Theory For Linear Operators*. 1980.

[28] N. Kielbasiewicz and E. Lunéville. User documentation for xlife++. https://uma.ensta-paris.fr/soft/XLiFE++/var/files/docs/usr/user documentation.pdf, 2019.
[29] A. Kirsch. Characterization of the shape of a scattering obstacle using the spectral data of the far field operator. *Inverse Problems*, 14(6):1489–1512, 1998.

[30] A. Kirsch. The factorization method for Maxwell's equations. *Inverse Problems*, 20(6):S117–S134, 2004.

[31] V.A. Kondratiev. Boundary-value problems for elliptic equations in domains with conical or angular points. *Trans. Moscow Math. Soc.*, 16:227–313, 1967.

[32] A. Lechleiter and T. Rienmüller. Factorization method for the inverse Stokes problem. *Inverse Probl. Imaging*, 7(4):1271–1293, 2013.

[33] P. Monk and V. Selgas. An inverse fluid-solid interaction problem. *Inverse Problems & Imaging*, 3(2):173, 2009.

[34] B. Pavlakovic, M. Lowe, and P. Cawley. Prediction of reflection coefficients from defects in embedded bars. In *Review of Progress in Quantitative Nonde-structive Evaluation*, pages 207–214. Springer, 1999.

[35] C. Tsogka, D. A. Mitsoudis, and S. Papadimitropoulos. Imaging extended reflectors in a terminating waveguide. *SIAM Journal on Imaging Sciences*, 11(2):1680–1716, 2018.