# Reconstruction of averaging indicators for highly heterogeneous media

Lorenzo Audibert, Houssem Haddar, Fabien Pourre

► To cite this version:

Lorenzo Audibert, Houssem Haddar, Fabien Pourre. Reconstruction of averaging indicators for highly heterogeneous media. Inverse Problems, 2024, 40 (4), pp.045028. 10.1088/1361-6420/ad2f64. hal-04223245

HAL Id: hal-04223245

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

Submitted on 29 Sep 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Reconstruction of averaging indicators for highly heterogeneous media

Lorenzo Audibert<sup>1,2</sup>, Houssem Haddar<sup>2</sup> and Fabien Pourre<sup>2</sup>

<sup>1</sup> Departement PRISME, EDF R&D, 6 quai Watier BP 49 Chatou, 78401 Cedex, France

<sup>2</sup> Inria, UMA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France

**Abstract.** We propose a new imaging algorithm capable of producing quantitative indicator functions for unknown and possibly highly oscillating media from multi-static far field measurements of scattered fields at a fixed frequency. The algorithm exploits the notion of modified transmission eigenvalues and their determination from measurements. We propose in particular the use of a new modified background obtained as the limit of a metamaterial background. It has the specificity of having a unique non trivial eigenvalue, which is particularly suited for the proposed imaging procedure. We show the efficiency of this new algorithm on some 2D experiments and emphasize its superiority with respect to some clasical approaches such as the Linear Sampling Method.

*Keywords:* Transmission Eigenvalues, Inverse Scattering, Generalized Linear Sampling Method, Steklov Eigenvalues
# 1. Introduction

In the context of inverse scattering problems at a fixed frequency, we are interested in the design of an imaging algorithm that provides interpretable indicator functions for highly fluctuating scatterers. This is motivated for instance by non destructive testing of composite materials, such as concrete, composed of unknown randomly distributed components.

We particularly target methods that do not rely on numerical modeling of the scattering problem. This is justified by the impossibility or difficulty of designing a sufficiently accurate numerical model when randomness is present in the probed media.

A first attempt to address this type of problems was done in [5]. The authors were able to design an imaging functional for highly fractured media capable of providing quantitative indicator of the crack's density. Their methodology is based on the introduction of an artificial background containing a resonator. Considering the scattering problem relative to this background, one can define some resonant frequencies that can be reconstructed from measured data. These frequencies would monotonically depend on the intersection of the artificial resonators with the set of cracks, which furnishes an indicator function on their density.

We shall follow here a similar route but improving the procedure in [5] in two important aspects. The first one is to consider measured data at a fixed frequency. This is done by encoding the spectral parameter in the design of the background. The second important aspect is the choice of the background model. We start with the metamaterial background proposed in [3] and related modified transmission eigenvalues that play the role of a spectral parameter [9]. We prove that the associated spectrum converges, as the metamaterial parameter goes to infinity, to a simplified spectrum made of only one non trivial eigenvalue. This is very appealing from the numerical point of view since it greatly simplifies the identification of the non trivial eigenvalue from measured data. In addition, the spectral parameter can also be equivalently interpreted as a weighted average of the refractive index with respect to the solution of a boundary value problem posed inside the artificial resonator. This allows for a better quantitative interpretation of the indicator function.

In this proof of concept we shall restrict ourselves to the scalar Helmholtz equation and to scatterers made of penetrable isotropic inhomogeneities (that for instance can be a collection of small inclusions randomly distributed in the probed domain). Our methodology is quite general and can be easily extended to other inverse scattering problems, such as electromagnetic or elastodynamic inverse problems. In addition to acoustic problems, the two dimensional version of the scalar problem can also be seen as a model for electromagnetic problems invariant in one direction.

The use of artificial backgrounds in imaging algorithms has been proposed in a number of recent works in the literature ([13, 12, 4, 16, 11]). Their use is different than in the imaging functional proposed in [3] and that we generalize in the present work.

The new artificial background (obtained as a limiting problem of the one introduced
first in [3]) can be seen as a special form of impedance boundary value problems. We analyse the properties of associated spectral problem (that shares similarities with Steklov eigenvalue problem [10]) and prove a monotonicity property of its unique non trivial eigenvalue with respect to the index of refraction. We then prove that this eigenvalue can be reconstructed from far field data at a fixed frequency. We employ the framework of the Generalized Linear Sampling Method [7] (GLSM) to prove this result in a constructive way. The proposed algorithm for a quantitative indicator function associated with the refractive index is then as follows: 1) choose a disc of given radius that we shall sweep over some sampling positions in the probed domain; 2) for each sampling position, identify the eigenvalue from far field data using the GLSM approach; 3) assign to the sampling position the difference between this eigenvalue and the eigenvalue associated with an index of refraction equals to one. The obtained indicator function would then be monotonically dependent with respect to the values of the true refractive index intersecting the disc. It also represent an averaged value of the contrast as we indicated above.

We explicit the details of this procedure in the numerical section and discuss the choice of the disc radius in terms of the used frequency. We then show how the obtained indicator function gives superior results to those given by geometrical approaches such as the Linear Sampling Method [14]. Moreover, the sign of the indicator function would give an indication on how the mean value of the refractive index compares with respect to 1. This is also discussed and illustrated through some numerical tests in 2D.

The article is organized as follows. In section 2, we introduce the direct problem and associated far field operator. We discuss in section 3 the metamaterial transmission eigenvalue problem and its limit as the metamaterial parameter goes to $\infty$. Then we analyse the properties of this limit problem. We introduce in section 4 the artificial background associated with the new limit eigenvalue problem. We explain in section 5 how to identify the associated modified transmission eigenvalue from far field data. Based on these materials, we propose in section 6 a new indicator function for the reconstruction some quantitative information related to the refractive index. Finally, this section is concluded with some validating numerical examples in 2D. Some complimentary results are given in the Appendix.

## 2. Settings of the direct problem and some related key results

Let $D$ be a bounded domain in $\mathbb{R}^m$, $m = 2, 3$ with piecewise smooth boundary $\partial D$ and connected complement. We denote by $\nu$ the outward normal field on $\partial D$. We shall consider the following scattering problem, with inhomogeneity supported in $D$. For $k > 0$ the wavenumber, let $u \in H_{\text{loc}}^1(\mathbb{R}^m)$ be the total field and $u^s$ be the scattered field
solutions of:

$$
\left\{
\begin{array}{l}
\Delta u + k^2 n u = 0 \text{ in } \mathbb{R}^m, \\
u = u^s + u^i \\
\lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial u^s}{\partial r} - i k u^s \right|^2 = 0,
\end{array}
\right.
\qquad (1)
$$

where the incident field $u^i = e^{ikx \cdot d}$, denoted $u^i(\cdot, d)$ for $d \in S = \{x \in \mathbb{R}^m, |x| = 1\}$ the unit sphere. We extend this notation to $u^s(\cdot, d)$ and $u(\cdot, d)$. The real refractive index $n \in L^\infty(\mathbb{R}^m)$ is such that $n > 0$ and $n = 1$ in $\mathbb{R}^m \setminus \bar{D}$. It is known that problem (1) is well posed and the scattered field $u^s$ has the following expansion for all $\hat{x} \in S$,

$$
u^s(r\hat{x}, d) = \frac{e^{ikr}}{r^{\frac{m-1}{2}}} u^\infty(\hat{x}, d) + O\left(\frac{1}{r^{\frac{m+1}{2}}}\right)
$$

with $u^\infty(\cdot, \cdot)$ being the so called far field pattern of the scattered field. We define the
far field operator $F : L^2(S) \to L^2(S)$ by

$$
(Fg)(\hat{x}) := \int_S g(d)u^{\infty}(\hat{x}, d)ds(d), \quad \hat{x} \in S. \tag{2}
$$

We remark that $u_g^\infty := Fg$ is the far field pattern of the scattered field $u_g^s$, solution of
(1) with $u^i = v_g$ the Herglotz wave function defined for $g \in L^2(S)$ by

$$
v_g(x) := \int_S e^{ikxd} g(d)ds(d). \tag{3}
$$

For a given $\psi \in L^2(D)$, we define the unique function $w^s \in H_{\text{loc}}^1(\mathbb{R}^m)$ satisfying

$$
\left\{
\begin{array}{l}
\Delta w^s + k^2 n w^s = k^2 (1-n) \psi \text{ in } \mathbb{R}^m, \\
\lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial w^s}{\partial r} - i k w^s \right|^2 = 0.
\end{array}
\right.
\qquad (4)
$$

Observe that if $\psi = v_g$, then $w^s = u_g^s$. We introduce the operator $H : L^2(S) \to L^2(D)$
defined by

$$
Hg := v_g|_{D}, \tag{5}
$$

along with its $L^2$-adjoint

$$
H^*(\psi) := \int_D \psi(y) e^{-ik\hat{x}\cdot y} dy. \qquad (6)
$$

Then, the far field operator $F$ has the following factorization [18, 7]

$$
Fg = H^*TH, \tag{7}
$$

where $T : L^2(D) \to L^2(D)$ is defined by

$$
T(\psi) = -\gamma k^2 (1-n)(\psi + w^s), \qquad (8)
$$

with $w^s$ the solution of (4), $\gamma = k^2/4\pi$ for $m=3$ and $\gamma = e^{i\pi/4}/\sqrt{8\pi k}$ for $m=2$.
For later use, we introduce the definition of transmission eigenvalues.
**Definition 1.** Values of $k > 0$ for which the problem:

$$ egin{cases} 
\Delta u + k^2 nu = 0 \text{ in } D, \\
\Delta v + k^2 v = 0 \text{ in } D, \\
v = u \text{ on } \partial D, \\
\partial_{\nu} v = \partial_{\nu} u \text{ on } \partial D,
\end{cases} 
\qquad (9) $$

has nontrivial solutions $(u, v) \in L^2(D) \times L^2(D)$, such that $w - v \in H^1(D_b)$ are called transmission eigenvalues.

## 2.1. Inverse problem and a motivating example

The inverse problem we would like to address is to retrieve information about the geometry of $D$ and the refractive index $n$ from the knowledge of the far field operator $F$. The developments in this article are motivated by the limitation of standard imaging methods to handle highly heterogeneous scatterers. Consider for instance $D$ as a union of small and close scatterers such as presented in the left part of Figure 1. The main objective is to recover the local distribution of those small inhomogeneities. The high density of the inclusions challenges classical methods such as the Linear Sampling Method (LSM) that try to construct an indicator function for the geometry of the scatterers. We display in the right part of Figure 1 the indicator result provided by the LSM for a fixed frequency. We clearly observe that not only the method fails to recover the geometry, but also gives no information about their distribution. To overcome this issue, we shall develop an imaging method based on transmission eigenvalues as introduced in [5] in order to build an indicator function that has connections with the scatterers density, which is the most relevant information one can obtain for this type of media. As opposed to [5], we shall work at a fixed frequency and replace transmission eigenvalues by a modified version based on the introduction of an artificial background.

## 3. An artificial metamaterial background

The starting point of our method is to consider so called modified transmission eigenvalues as introduced in [3] with an artificial background media being a metamaterial. We sketch in the following this idea. Let $D_b \subset \mathbb{R}^m$ be a regular bounded simply connected domain. Let $\Omega$ be an open domain such that $\bar{\Omega} = \bar{D} \cup \bar{D}_b$ and such that $\Omega$ has a connected complement. Let $u_b \in H_{\text{loc}}^1(\mathbb{R}^m)$ be the total field and $u_b^s$ be the scattered field solutions of

$$ egin{cases} 
\text{div}(a_b \nabla u_b) + k^2 n_b(\lambda) u_b = 0 \text{ in } \mathbb{R}^m, \\
u_b = u^i + u_b^s, \\
\lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial u_b^s}{\partial r} - i k u_b^s \right|^2 = 0,
\end{cases} 
\qquad (10) $$

for some incident field $u^i$. The coefficients $a_b, n_b(\lambda) \in L^\infty(\mathbb{R}^m)$ are such that

$$ a_b = \begin{cases} 1 & \text{in } \mathbb{R}^m \setminus \bar{D}_b, \\ -a & \text{in } D_b, \end{cases} \quad \text{and} \quad n_b(\lambda) = \begin{cases} 1 & \text{in } \mathbb{R}^m \setminus \bar{D}_b, \\ \lambda & \text{in } D_b, \end{cases} $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left: the domain $D$ constituted by 2250 small circles of radius 0.02 having a constant index of refraction $n = 2$. Right: indicator function provided by the LSM with 1% added noise to the far field operator $F$ and for frequency $k = 4$ (the wavelength is roughly 80 times higher than the radius of the circles.)</figcaption></figure>

for $a > 0$ and $a \neq 1$ and $\lambda \in \mathbb{R}$ that will play the role of a spectral parameter. Similarly to the forward scattering problem, we use the notation $u_b(\cdot, d)$, $u_b^s(\cdot, d)$, $u_b^\infty(\cdot, d)$ for $u^i = u^i(\cdot, d)$ and consider the far field operator $F_b : L^2(S) \to L^2(S)$

$$ (F_b g)(\hat{x}) := \int_S g(d) u_b^\infty(\hat{x}, d) ds(d), \quad \hat{x} \in S. $$

We then define the modified far field operator $\mathcal{F} : L^2(S) \to L^2(S)$

$$ \mathcal{F}g := Fg - F_b g. $$

Modified transmission eigenvalues $\lambda$ correspond to the values of $\lambda$ such that there exists a non trivial incident field in $L^2(\Omega) \cap H^1(D_b)$ such that, $u^s = u_b^s$ in $\mathbb{R}^m \setminus \Omega$. In particular, we obtain that for these values, there exists a non trivial solution $w := u|_{\Omega}$ and $v := u_b|_{\Omega}$ of the following problem

$$ \left\{ \begin{array}{l} \Delta w + k^2 n w = 0 \text{ in } \Omega, \\ \operatorname{div}(a_b \nabla v) + k^2 n_b(\lambda)v = 0 \text{ in } \Omega, \\ v = w \text{ on } \partial\Omega, \\ a_b \partial_\nu v = \partial_\nu w \text{ on } \partial\Omega. \end{array} \right. \qquad (11) $$

In the case where $u|_{\Omega} = w$ and $u_b|_{\Omega} = v$ respectively coincide with the solutions of (11) and (10) with $u^i = v_g$, we have $\mathcal{F}g = 0$. This means that for this particular Herglotz wave, the scattered wave produced by the original media and the background media are the same outside $\Omega$.

Let us denote by $\eta_0(n, D_b)$ the first Dirichlet eigenvalue of the problem:

$$ \left\{ \begin{array}{l} w \in H^1(D_b), \\ \Delta w + \eta n w = 0 \text{ in } D_b, \quad w = 0 \text{ on } \partial D_b. \end{array} \right. \qquad (12) $$
The interior transmission eigenvalue problem (11) has been studied in [3] in the case $D \subset D_b$. It is proved that if $k^2$ is not an eigenvalue of (12) then the spectrum is discrete without finite accumulation point and there exists at least one positive eigenvalue to (11). In addition, if $k^2 < \eta_0(n, D_b)$, then the eigenvalues are monotonically increasing with respect to $n$. In particular, the largest positive eigenvalue $\lambda_1$ satisfies:

$$ \lambda_1(n, D_b) = \sup_{w, v \in \mathbf{H}_*(D_b)} \frac{k^2 \int_{D_b} n|w|^2 dx - \int_{D_b} |\nabla w|^2 dx - a \int_{D_b} |\nabla v|^2 dx}{k^2 \int_{D_b} |v|^2 ds}, \quad (13) $$

where $\mathbf{H}_*(D_b) := \{w, v \in H^1(D_b), w|_{\partial D_b} = v|_{\partial D_b}, v \neq 0\}$.

The imaging algorithm relies on identifying these modified transmission eigenvalues from knowledge of the operator $F$ and numerical computations of the operator $F_b$ associated with different locations of the domain $D_b$. In practice, $D_b$ will be a circle that sweeps over the probed domain. For numerical efficiency, the parameter $a$ in (11) has to be chosen in order to have the simplest possible structure of the spectrum and the highest sensitivity with respect to $n$. These eigenvalues can be computed analytically in the case of a ball and constant $n$ as the zero of a determinant. Indeed, in dimension 2, for $D = D_b$ a centered ball at the origin circle of radius $\rho$, the eigenvalues are solutions to:

$$ \mathcal{P}_j(a, \lambda) = \det \begin{pmatrix} J_j(k\rho\sqrt{n}) & -J_j(k\rho\sqrt{\frac{-\lambda}{a}}) \\ \sqrt{n}J'_j(k\rho\sqrt{n}) & a\sqrt{\frac{-\lambda}{a}}J'_j(k\rho\sqrt{\frac{-\lambda}{a}}) \end{pmatrix} = 0, $$

where $j \in \mathbb{Z}$ and $J_j$ denotes the Bessel function of the first kind of order $j$. We display in Figure 2 the values of $\mathcal{P}_0(a, \lambda)^{-1}$. In particular, we observe that as $a$ goes to infinity, negative eigenvalues go to $-\infty$ and the largest one converges to a unique positive value. This suggests that only the positive eigenvalues would be the most stably reconstructed one. To avoid playing with the parameter $a$, it would be interesting to consider the spectral problem (11) as $a$ goes to infinity. This limit is studied in the following section.

### 3.1. The limit spectral problem as $a \to +\infty$

Let us consider the interior transmission problem (11) where $D$ does not intersect $\partial D_b$ as presented in Figure 3.

In this case, solving the spectral problem (11) is equivalent to solving the following one posed on $D_b$ by setting $w = v = 0$ in $\Omega \setminus D_b$ (assuming that $k$ is not a transmission eigenvalue for classical transmission problem in $D$ [9]): $w \in H^1(D_b), v \in H^1(D_b)$,

$$ \left\{ \begin{array}{l} \Delta w + k^2nw = 0 \text{ in } D_b, \\ -a\Delta v + k^2\lambda v = 0 \text{ in } D_b, \\ v = w \text{ on } \partial D_b, \\ -a\partial_\nu v = \partial_\nu w \text{ on } \partial D_b. \end{array} \right. \qquad (14) $$

**Theorem 1.** Assume that $k^2 < \eta_0(n, D_b)$. Let $(a_j)_j$ be an unbounded increasing sequence of positive numbers larger than 1. Let $\lambda_j$ be the largest positive eigenvalue of (14) for $a = a_j$ and $(w_j, v_j) \in H^1(D_b) \times H^1(D_b)$ be an associated eigenvector such that
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Values of $P_0(a, \lambda)^{-1}$ for constant $n|_{D_b} = 1$ and product $k\rho = 1$. The bright regions indicate where the function approaches 0.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Example of a configuration where $D \cap \partial D_b = \emptyset$</figcaption></figure>

$\|w_j\|_{H^1(D_b)}^2 + \|v_j\|_{H^1(D_b)}^2 = 1$. Then, $\lambda_j$ converges to some non negative value $\frac{\mu}{k^2|D_b|}$ and (up to a subsequence) $w_j$ converges weakly in $H^1(D_b)$ and strongly in $L^2(D_b)$ to a non trivial function $w \in H^1(D_b)$ that satisfies

$$ \left\{ \begin{array}{l} \Delta w + k^2nw = 0 \text{ in } D_b, \\ \mu w + \int_{\partial D_b} \partial_\nu w ds = 0 \text{ on } \partial D_b, \end{array} \right. \qquad (15) $$

where $|D_b| := \int_{D_b} dx$. The corresponding sequence $(v_j)_j$ strongly converges to a constant in $H^1(D_b)$.

*Proof.* According to (13), the positive sequence $(\lambda_j)$ is decreasing and therefore converges to some $\lambda \ge 0$. Moreover, since the sequence $(w_j, v_j)$ is bounded in $H^1(D_b) \times H^1(D_b)$, up to a subsequence, we can assume that $(w_j, v_j)$ converges to some $(w, v)$ weakly in $H^1(D_b) \times H^1(D_b)$ and strongly in $L^2(D_b) \times L^2(D_b)$. Problem (14) is
equivalent to the variational formulation:

$$ \int_{D_b} \nabla w_j \nabla w' dx + a_j \int_{D_b} \nabla v_j \nabla v' dx - k^2 \int_{D_b} n w_j w' dx = -k^2 \lambda_j \int_{D_b} v_j v' dx, (16) $$

for all $w', v' \in H^1(D_b)$, $w' = v'$ on $\partial D_b$. Taking $v' = v_j$ and $w' = w_j$ in (16) shows that the sequence $a_j ||\nabla v_j||_{L^2(D_b)}^2$ is bounded. Therefore $\nabla v_j$ converges strongly to 0 in $L^2(D_b)$ which proves that $v = K$ a constant in $D_b$. We then infer that $w = K$ on $\partial D_b$. Taking $v' = 0$ in $D_b$ and $w' = H_0^1(D_b)$ and passing to the limit as $j$ goes to infinity in (21) shows that

$$ \Delta w + k^2 n w = 0 \text{ in } D_b. \qquad (17) $$

Since $-a_j \Delta v_j + k^2 \lambda_j v_j = 0$ in $D_b$, we obtain, by integration over $D_b$, the Green formula and $-a_j \partial_\nu v_j = \partial_\nu w_j$ on $\partial D_b$,

$$ \int_{\partial D_b} \partial_\nu w_j ds + \lambda_j k^2 \int_{D_b} v_j dx = 0, $$

where the first integral must be understood as $< H^{-\frac{1}{2}}(\partial D_b), H^{\frac{1}{2}}(\partial D_b) >$ duality pairing. This notation convention will be constantly used in the rest of the product. Passing to the limit in the previous equation leads to

$$ \mu w + \int_{\partial D_b} \partial_\nu w ds = 0, $$

using the fact that $\partial_\nu w_j$ weakly converges to $\partial_\nu w$ in $H^{-\frac{1}{2}}(\partial D_b)$. This proves that $w$ satisfies (15).

We now prove by contradiction that $w$ is not a trivial function. Assume that $w = 0$. Then $K = 0$ and therefore $v = 0$. Taking $v' = v_j$ and $w' = w_j$ in (16) gives $\int_{D_b} |\nabla w_j|^2 dx + a_j \int_{D_b} |\nabla v_j|^2 dx \to 0$ as $j$ goes to infinity. This proves that $w_j$ converges strongly to 0 in $H^1(D_b)$ and we already have $v_j$ converges to 0 strongly in $H^1(D_b)$. This contradicts $\|w_j\|_{H^1(D_b)}^2 + \|v_j\|_{H^1(D_b)}^2 = 1$. $\square$

**Remark 1.** In dimension 2, in the case where $D_b$ is a disk of radius $\rho > 0$ and $n$ is constant, one can compute explicitly the eigenvalue $\mu$ of problem (15). We get that $\mu = \mu_{\text{ref}}(n)$ where

$$ \mu_{\text{ref}}(n) := 2\pi\rho k \sqrt{n} \frac{J_1(k\rho\sqrt{n})}{J_0(k\rho\sqrt{n})}. \qquad (18) $$

A similar expression can be derived in space dimension 3 for balls.

Theorem 1 suggests to replace the spectral problem (14) by (15) as $a$ goes to infinity. As we shall demonstrate in the sequel, the eigenvalue problem (15) has a very simple spectrum constituted by at most one non trivial eigenvalue. This is why it also constitutes a good candidate for the imaging algorithm and that will be the subject of Section 5, after introducing the associated direct scattering problem.
## 3.2. Analysis of the spectrum of (15)

This section is dedicated to the analysis of the spectral problem (15). This analysis requires to exclude values of $k^2$ that are eigenvalues of the following problem.

$$ egin{cases} w \in H_c^1(D_b), \\ \Delta w + k^2nw = 0 \text{ in } D_b, \\ \int_{\partial D_b} \partial_\nu w ds = 0 \text{ on } \partial D_b, \end{cases} \quad (19) $$

where $H_c^1(D_b) := \{w \in H^1(D_b), w|_{\partial D_b} \text{ is constant}\}$. The eigenvalue problem (19) has a similar structure as Dirichlet or Neumann eigenvalue problems. It is equivalent to the variational writing: $w \in H_c^1(D_b)$

$$ \int_{D_b} \nabla w \nabla \overline{w'} dx - k^2 \int_{D_b} nw \overline{w'} dx = 0 \text{ for all } w' \in H_c^1(D_b). \quad (20) $$

It follows from standard results on the spectrum of self-adjoint compact operators that the spectrum of (19) is discrete and formed by a real sequence that accumulates at $+\infty$.

**Proposition 2.** Assume that $k^2$ is not an eigenvalue of (19), then problem (15) has at most one eigenvalue. Moreover, this eigenvalue is non trivial.

*Proof.* Problem (15) is equivalent to the variational formulation: $w \in H_c^1(D_b)$

$$ \int_{D_b} \nabla w \nabla \overline{w'} dx - k^2 \int_{D_b} nw \overline{w'} dx = - \frac{\mu}{|\partial D_b|} \int_{\partial D_b} w \overline{w'} ds \quad (21) $$

for all $w' \in H_c^1(D_b)$. Thanks to the Riesz representation theorem, we define two selfadjoint operators $\mathcal{A}$ and $\mathcal{T}: H_c^1(D_b) \to H_c^1(D_b)$ by:

$$ \begin{aligned} (\mathcal{A}w, w')_{H^1(D_b)} &= \int_{D_b} \nabla w \nabla \overline{w'} dx - k^2 \int_{D_b} nw \overline{w'} dx, \quad \forall w' \in H_c^1(D_b), \\ (\mathcal{T}w, w')_{H^1(D_b)} &= \frac{1}{|\partial D_b|} \int_{\partial D_b} w \overline{w'} ds, \quad \forall w' \in H_c^1(D_b). \end{aligned} \quad (22) $$

The eigenvalue problem (15) is then equivalent to:

$$ \mathcal{A}w = -\mu \mathcal{T}w, \quad \text{in } H_c^1(D_b). \quad (23) $$

Notice that since $w|_{\partial D_b}$ is constant,

$$ (\mathcal{T}w, w')_{H^1(D_b)} = \frac{w|_{\partial D_b}}{|\partial D_b|} \int_{\partial D_b} \overline{w'} ds = w|_{\partial D_b} (\mathcal{T}w_1, w')_{H^1(D_b)}, $$

where $w_1 := 1$ in $D_b$. Therefore, the range of $\mathcal{T}$ is of dimension 1. The hypothesis on $k^2$ shows that $\mathcal{A}$ is injective and therefore there exists at most one non trivial eigenvalue of (23). □

**Theorem 3.** The spectrum of problem (15) can be described as follows:

1) Assume that $k^2$ is not an eigenvalue of (19).

a) If $k^2$ is not a Dirichlet eigenvalue of (12), then there exists an unique non trivial eigenvalue $\mu$ to (15) and an associated eigenvector is $w_1 \in H^1(D_b)$ the unique solution of $\Delta w_1 + k^2nw_1 = 0$ in $D_b$ and $w_1 = 1$ on $\partial D_b$. In particular, we have $\mu = -\int_{\partial D_b} \partial_\nu w_1 ds = k^2 \int_{D_b} nw_1 dx$.
b) If $k^2$ is a Dirichlet eigenvalue of (12), then there are no eigenvalues to (15).

2) Assume that $k^2$ is an eigenvalue of (19). Then 0 is an eigenvalue of (15) and the associated eigenvectors are the eigenvectors of (19). If $k^2$ is not a Dirichlet eigenvalue, 0 is the unique eigenvalue.

*Proof.* We distinguish the various cases as in the Theorem.

1) Assume that $k^2$ is not an eigenvalue of (19). Then according the Proposition 2, we have at most one eigenvalue.

a) If $k^2$ is not a Dirichlet eigenvalue of (12). Then the function $w_1$ as in item a) can be defined and one easily verifies that it is an eigenvector associated with $\mu = -\int_{\partial D_b} \partial_\nu w_1$.

b) If $k^2$ is a Dirichlet eigenvalue of (12), assume by contradiction that there exists an eigenvalue $\mu$ of (15) associated with an eigenvector $w$. Since $k^2$ is not an eigenvalue of (19), $\mu \neq 0$ and $w|_{\partial D_b} \neq 0$. Let $w_0 \in H^1(D_b)$ the eigenvector verifying:

$$ \begin{cases} \Delta w_0 + k^2 n w_0 = 0 & \text{in } D_b, \\ w_0 = 0 & \text{on } \partial D_b. \end{cases} $$

Multiplying the previous equation by $w_1$, and integrating by part gives:

$$ \int_{D_b} (\Delta w_0 + k^2 n w_0) w dx = \int_{\partial D_b} w \partial_\nu w_0 ds = w|_{\partial D_b} \int_{\partial D_b} \partial_\nu w_0 ds = 0. $$

Therefore, $w_0$ is an eigenvector of (19) which is a contradiction.

2) Assume that $k^2$ is an eigenvalue of (19). One can easily verify that 0 is an eigenvalue of (15) by inserting $\mu = 0$ in equation (15). If $k^2$ is not a Dirichlet eigenvalue of (12), assume by contradiction that there exists a non trivial eigenvalue $\mu$ of (15) and an associated eigenvector $w_1$. We denote by $w_0$ an eigenvector of (19). Notice that $w_0|_{\partial D_b} \neq 0$ and $w_1|_{\partial D_b} \neq 0$. We obtain a contradiction following a similar reasoning as in item 1-b.

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

We have proved the existence of a unique eigenvalue of (15) under certain conditions on the wavenumber. The following theorem proves the monotonicity property of the eigenvalue with respect to the refractive index.

**Theorem 4.** Assume that $0 < k^2 < \eta_0(n, D_b)$. Then $k^2$ is not an eigenvalue of (19). Let $\mu(k, n, D_b)$ be the unique eigenvalue of (15) as defined in Theorem 3 1-a. Then we also have

$$ \mu(k, n, D_b) = \sup_{w \in H_c^1(D_b), w|_{\partial D_b} = 1} (k^2 \int_{D_b} n|w|^2 dx - \int_{D_b} |\nabla w|^2 dx). \quad (24) $$

This shows in particular that $\mu(k, n, D_b)$ is positive and is monotonically increasing with respect to $n$. Moreover, we have that $\mu(k, n, D_b) \to +\infty$ as $k^2 \to \eta_0(n, D_b)$.
Proof. We start by proving that $k^2$ is not an eigenvalue of (19). First of all, we notice that the first eigenvalue of (19) is 0 associated with the constant eigenvector 1. Then the second eigenvalue is given by:

$$k_1^2 = \inf_{E_2 \subset W_2} \sup_{w \in E_2, w \neq 0} \frac{\int_{D_b} |\nabla w|^2 dx}{\int_{D_b} n|w|^2 dx}, \qquad (25)$$

where $W_2$ is the set of 2 dimensional subspaces of $H_c^1(D_b)$. For all $w \in H_0^1(D_b)$, there exists $E_2$ such that $w \in E_2$ since $H_0^1(D_b) \subset H_c^1(D_b)$. Reciprocally any $E_2 \subset W_2$ contains at least a non trivial element of $H_0^1(D_b)$, for instance $w = u_2|_{\partial D_b}u_1 - u_1|_{\partial D_b}u_2$ where $u_1, u_2$ are two linear independant vectors of $E_2$. Therefore:

$$k_1^2 \geq \inf_{w \in H_0^1(D_b)} \frac{\int_{D_b} |\nabla w|^2 dx}{\int_{D_b} n|w|^2 dx} = \eta_0(n, D_b). \qquad (26)$$

Let $\Lambda > 0$ be such that the operator $\Lambda'$ defined by

$$(\Lambda' w, w')_{H^1(D_b)} = \int_{D_b} \nabla w \overline{\nabla w'} dx - k^2 \int_{D_b} n w \overline{n w'} dx + \frac{\Lambda}{|\partial D_b|} \int_{\partial D_b} w \overline{w} ds, (27)$$

for all $w' \in H_c^1(D_b)$ is coercive. The existence of such $\Lambda$ follows from Wirtinguer's type inequalities and the proof is given in Appendix A, Theorem 12 for the reader's convenience. The eigenvalue problem (15) can be equivalently written as:

$$\Lambda' w = -(\mu - \Lambda)Tw, \qquad (28)$$

where $\Lambda'$ is a selfadjoint coercive operator and $T$ is a selfadjoint compact non negative operator. Identity (24) is then a direct consequence of the Courant-Fischer principle applied to the eigenvalue problem (28).

We now prove that $\mu(k, n, D_b) \to +\infty$ as $k^2 \to \eta_0(n, D_b)$. Let $j \in \mathbb{N}^*$ and denote by $(k_j, w_j)$ the first eigenvalue and eigenvector with $||w_j||_{H^1(D_b)} = 1$ of the following spectral problem (j is fixed):

$$\left\{ \begin{array}{l} w_j \in H_c^1(D_b), \\ \Delta w_j + k_j^2 n w_j = 0 \text{ in } D_b, \\ j w_j + \int_{\partial D_b} \partial_n w_j ds = 0 \text{ on } \partial D_b. \end{array} \right. \qquad (29)$$

We obtain by the Courant-Fischer principle:

$$k_j^2 = \inf_{w \in H_c^1(D_b), w \neq 0} \frac{\int_{\partial D_b} \frac{j}{|\partial D_b|} |w|^2 + \int_{D_b} |\nabla w|^2}{\int_{D_b} n|w|^2} \leq \inf_{w \in H_0^1(D_b), w \neq 0} \frac{\int_{D_b} |\nabla w|^2}{\int_{D_b} n|w|^2}.$$

Therefore, $(k_j^2)_j$ an increasing sequence and bounded by $\eta_0(n, D_b)$. This implies $k_j^2 \to \eta^*$ as $j$ goes to infinity for some $\eta^* \leq \eta_0(n, D_b)$. The eigenvector $w_j$ verifies:

$$\int_{D_b} |\nabla w_j|^2 dx - k_j^2 \int_{D_b} n|w_j|^2 dx = - \frac{j}{|\partial D_b|} \int_{\partial D_b} |w_j|^2 ds. \qquad (30)$$

This equality proves that the right hand side remains bounded as $j$ goes to infinity, which means that $w_j|_{\partial D_b} \to 0$. We then have that (up to a subsequence) $(w_j)_j$ converges weakly in $H^1(D_b)$ (and strongly in $L^2(D_b)$) to some $w$ in $H_c^1(D_b)$ that verifies
$\Delta w + \eta^*nw = 0$ in $D_b$ and $w = 0$ on $\partial D_b$. Let us show that $w$ is non trivial by using a contradiction argument. Suppose $w = 0$. From (30), we have:

$$\int_{D_b} |\nabla w_j|^2 dx \leq k_j^2 \int_{D_b} n|w_j|^2 dx,$$

and therefore $\nabla w_j \to 0$ in $L^2(D_b)$. Consequently $w_j$ strongly converges to 0 in $H^1(D_b)$. This contradicts $\|w_j\|_{H^1(D_b)} = 1$.

The fact that $w$ is non trivial implies that $\eta^*$ is a Dirichlet eigenvalue of (12) and therefore $\eta^* \geq \eta_0(n, D_b)$. We then have $\eta^* = \eta_0(n, D_b) = \lim_{j \to +\infty} k_j^2$. Observe that $w_j$ is not in $H_0^1(D_b)$ thanks to (25). Using (24) and (30), we then get:

$$\mu(k_j, n, D_b) \geq \frac{k_j^2 \int_{D_b} n|w_j|^2 dx - \int_{D_b} |\nabla w_j|^2 dx}{\int_{\partial D_b} |w_j|^2 ds} = \frac{j}{|\partial D_b|}.$$

$$\text{Thus, } \lim_{j \to +\infty} \mu(k_j, n, D_b) \geq \frac{j}{|\partial D_b|} \to_{j \to +\infty} +\infty \quad \square$$

### 4. An artificial background associated with (15)

The goal of this section is to introduce a scattering problem for which the spectral problem (15) would play the role of interior transmission problem similarly to the role played by (11) for (10). Let $D_b$ and $\Omega$ be as in section 3. The new artificial scattering problem can be formulated as: the total field $u_b \in H_{loc}^1(\mathbb{R}^m)$ satisfies

$$\left\{
\begin{array}{l}
\Delta u_b + k^2 u_b = 0 \text{ in } \mathbb{R}^m \setminus \overline{D_b}, \\
\mu u_b + \int_{\partial D_b} \partial_\nu u_b ds = 0 \text{ on } \partial D_b, \\
u_b = u_b^s + u^i, \\
\lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial u_b^s}{\partial r} - ik u_b^s \right|^2 = 0,
\end{array}
\right.
\quad (31)$$

for some incident field $u^i$ and $\mu \in \mathbb{R}$ is the spectral parameter. The analysis of this scattering problem can be done using a standard technique, for instance a variational approach based on the introduction of the Dirichlet to Neumann operator. Let $B_R$ (with boundary $S_R$) be a ball containing $D_b$ for $R$ large enough. We define the Dirichlet to Neumann map $\Lambda: H^{\frac{1}{2}}(S_R) \to H^{-\frac{1}{2}}(S_R)$ by $\Lambda\phi := \partial_\nu v|_{S_R}$, where $v$ is the radiating solution to the Helmholtz equation in $\mathbb{R}^d \setminus B_R$ with $v = \phi$ on $S_R$, and where $\nu$ is the outward unit normal to $S_R$. Consider the following space

$$H_c^1(B_R \setminus D_b) := \{w \in H^1(B_R \setminus D_b), w|_{\partial D_b} \text{ is constant}\}. \quad (32)$$

Then, solving the scattering problem (31) is equivalent to solving the following variational problem: $u \in H_c^1(B_R \setminus D_b)$ and for all $v \in H_c^1(B_R \setminus D_b)$,

$$\int_{B_R \setminus D_b} (\nabla u \nabla v - k^2 uv) dx - \int_{S_R} \Lambda uv ds - \frac{\mu}{|\partial D_b|} \int_{\partial D_b} uv ds = \ell(v) \quad (33)$$

where

$$\ell(v) := \int_{S_R} (\partial_\nu u^i - \Lambda u^i) v ds. \quad (34)$$
Using similar arguments as in [9], one can prove that this problem is well posed for $\mu \in \mathbb{R}$ (or more generally, for non trivial complex value $\mu$ with $\Im(\mu) \ge 0$). We summarize the wellposedness in the following theorem and leave the proof to the reader.

**Theorem 5.** The scattering problem (31) is well posed for $u^i$ such that $(\partial_\nu u^i - \Lambda u^i) \in H^{-\frac{1}{2}}(S_R)$.

Similarly to the forward scattering problem, we use the notation $u_b(\cdot, d)$, $u_b^s(\cdot, d)$, $u_b^\infty(\cdot, d)$ to respectively denote the total, scattered and far field when the incident field $u^i = u^i(\cdot, d)$, the plane wave of direction $d \in S$. Consider the far field operator $F_b^\mu : L^2(S) \to L^2(S)$ defined by

$$ (F_b^\mu g)(\hat{x}) := \int_S g(d)u_b^\infty(\hat{x}, d)ds(d), \quad \hat{x} \in S. $$

We then define the modified far field operator $\mathcal{F}^\mu : L^2(S) \to L^2(S)$ as

$$ \mathcal{F}^\mu g := Fg - F_b^\mu g. \qquad (35) $$

For further use, especially in the proofs of section 5, we notice that similarly to (4) the scattered field $u_b^s \in H_{\text{loc}}^1(\mathbb{R}^m \setminus D_b)$ satisfies:

$$ \left\{ \begin{array}{l} \Delta u_b^s + k^2 u_b^s = 0 \text{ in } \mathbb{R}^m \setminus D_b, \\ \mu u_b^s + \int_{\partial D_b} \partial_\nu u_b^s = -\mu u^i - \int_{\partial D_b} \partial_\nu u^i ds \text{ on } \partial D_b, \\ \lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial u_b^s}{\partial r} - iku_b^s \right|^2 = 0, \end{array} \right. \qquad (36) $$

for $u^i$ satisfying $\Delta u^i + k^2 u^i = 0$ in $\mathbb{R}^m$. We denote by $u_{b,g}^s$ the solution of (36) with $u^i = v_g$ defined in (3) and $u_{b,g} = u_{b,g}^s + v_g$ in $\mathbb{R}^m \setminus D_b$.

**Definition 2.** We say that $\mu \in \mathbb{R}$ is a modified transmission eigenvalue if problem (15) admits a non trivial solution.

## 5. An indicator function for the refractive index based on the eigenvalue of (15)

This section develops a method to recover some information on the refractive index $n$ from the far field operator $F$. This method is inspired by the one introduced in [5] and exploit the monotonicity property of the modified eigenvalue of the spectral problem (15) with respect to $n$ (see Theorem 4). The first step is to explain how one can identify the modified transmission eigenvalue from $F^\mu$. We shall use the GLSM framework as introduced in [7] for which we summarize the theoretical background in Apprendix, Theorem 14 for the reader's convenience. In order to avoid unnecessarily technical complications, we shall restrict ourselves to the case where $\partial D_b \cap D = \emptyset$.

### 5.1. Factorization of $\mathcal{F}^\mu$

We define $\Omega_D := D \setminus D_b$ the component of $\Omega$ outside of $D_b$. Let $\mathcal{H}^\mu : L^2(S) \to \mathbb{R} \times H^{-\frac{1}{2}}(\partial D_b) \times L^2(\Omega_D)$ defined as:

$$ \mathcal{H}^\mu g := (u_{b,g}|_{\partial D_b}, \partial_\nu u_{b,g}|_{\partial D_b}, u_{b,g}|_{\Omega_D}), $$
where $u_{b,g}$ is solution of (31) with $u^i = v_g$ and

$$V_{\text{inc}}(\Omega_D) = \{u \in L^2(\Omega_D) \text{ such that } \Delta u + k^2 u = 0 \text{ in } \Omega_D\}.$$

**Lemma 6.** We introduce the subspace $V$ defined as:

$$V = \{(\varphi, \psi, u) \in \mathbb{R} \times H^{-\frac{1}{2}}(\partial D_b) \times V_{\text{inc}}(\Omega_D) \text{ such that } \mu\varphi + \int_{\partial D_b} \psi ds = 0\}.$$

The closure of the range of $\mathcal{H}^\mu$, denoted $\mathcal{R}(\mathcal{H}^\mu)$, in $\mathbb{R} \times H^{-\frac{1}{2}}(\partial D_b) \times L^2(\Omega_D)$ satisfies:

$$\overline{\mathcal{R}(\mathcal{H}^\mu)} = V. \qquad (37)$$

*Proof.* One can easily verify that $\overline{\mathcal{R}(\mathcal{H}^\mu)} \subset V$. To prove the reverse inclusion, it is sufficient to prove that the orthogonal of $\mathcal{R}(\mathcal{H}^\mu)$ in $V$ is reduced to the trivial vector. Let $(\varphi, \psi, u) \in V$, such that $(\mathcal{H}^\mu g, (\varphi, \psi, u)) = 0$, for all $g \in L^2(S)$, which can be rewritten as follow:

$$\varphi u_{b,g} + (\psi, \partial_\nu \overline{u_{b,g}})_{H^{-\frac{1}{2}}(\partial D_b)} + (u, \overline{u_{b,g}}|_{\Omega_D})_{L^2(\Omega_D)} = 0. \qquad (38)$$

Let $\tilde{\psi} \in H^{\frac{1}{2}}(\partial D_b)$ be the unique vector such that for all $f \in H^{-\frac{1}{2}}(\partial D_b)$, we have

$$(\psi, f)_{H^{-\frac{1}{2}}(\partial D_b)} = \int_{\partial D_b} \tilde{\psi} f ds, \qquad (39)$$

where the integral on the right hand side refers to the $<H^{-\frac{1}{2}}(\partial D_b), H^{\frac{1}{2}}(\partial D_b)>$ duality product that can be replaced by the integral if $f \in L^2(\partial D_b)$ for instance. In particular, from equation (38), we get for all $g \in L^2(S)$:

$$\varphi u_{b,g} + \int_{\partial D_b} \tilde{\psi} \partial_\nu u_{b,g} ds + \int_{\Omega_D} u u_{b,g} dx = 0. \qquad (40)$$

Let $w \in H_{\text{loc}}^1(\mathbb{R}^m \setminus D_b)$ be a solution to

$$\left\{ \begin{array}{l} \Delta w + k^2 w = u|_{\Omega_D} \text{ in } \mathbb{R}^m \setminus D_b, \\ \frac{1}{\mu} \int_{\partial D_b} \partial_\nu w + w = \varphi - \tilde{\psi} \text{ on } \partial D_b, \\ \lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial w}{\partial r} - ikw \right|^2 = 0. \end{array} \right. \qquad (41)$$

Thus, we can replace the last integral in (40) by an integral on a centered ball $B_R$ for $R > 0$.

$$(u, \overline{u_{b,g}|_{\Omega_D}})_{L^2(\Omega_D)} = \int_{B_R \setminus D_b} (\Delta w + k^2 w) u_{b,g}^s dx + \int_{\Omega_D} uv_g dx. \qquad (42)$$

Applying the Green formula twice and using the relation $u_{b,g} = v_g + u_{b,g}^s$ in $\mathbb{R}^m \setminus D_b$, one obtains:

$$\int_{B_R \setminus D_b} (\Delta w + k^2 w) u_{b,g}^s dx = \int_{\partial B_R} (\partial_\nu w u_{b,g}^s - w \partial_\nu u_{b,g}^s) ds \\ - \int_{\partial D_b} (\partial_\nu w u_{b,g} - w \partial_\nu u_{b,g}) ds + \int_{\partial D_b} (\partial_\nu w v_g - w \partial_\nu v_g) ds.$$
The first term on the right hand side goes to 0 as $R \to +\infty$ thanks to the radiation condition. Therefore, inserting the previous equality in (42) gives:

$$ (u, \overline{u_{b,g}^s|_{\Omega_D}})_{L^2(\Omega_D)} = - \int_{\partial D_b} (\partial_\nu w u_{b,g} - w \partial_\nu u_{b,g}) ds + \int_{\partial D_b} (\partial_\nu w v_g - w \partial_\nu v_g) ds. \quad (43) $$

Inserting (43) in equation (40) and using the boundary condition of $w$, we obtain:

$$ \int_{\partial D_b} (\partial_\nu w v_g - \partial_\nu v_g w) ds + \int_{\Omega_D} uv_g dx = 0. \quad (44) $$

Replacing $v_g$ by its expression and inverting the integrals shows:

$$ \int_S g(d)w^\infty(-d)ds(d) = 0, \quad (45) $$

where

$$ w^\infty(d) = - \int_{\partial D_b} (\partial_{\nu(x)} w(x) e^{-ikx \cdot d} - \partial_{\nu(x)} e^{-ikx \cdot d} w) ds(x) - \int_{\Omega_D} u(x) e^{-ikx \cdot d} dx, $$

is the far field pattern of $w$. By the Rellich lemma, we then get that $w = 0$ in $\mathbb{R}^m \setminus D_b$. Hence, $w \in H_0^2(\Omega_D)$ and $\frac{\varphi}{\mu} - \tilde{\psi} = 0$ on $\partial D_b$. Since $u \in L^2(\Omega_D)$ and satisfies $\Delta \bar{u} + k^2 \bar{u} = 0$ in $\Omega_D$, we have in particular:

$$ \int_{\Omega_D} \bar{u}(\Delta w + k^2 w) dx = 0, $$

which gives $\|u\|_{L^2(\Omega_D)}^2 = 0$ therefore $u = 0$. On the other hand, taking the duality product between $\frac{\varphi}{\mu} - \tilde{\psi}$ and $f \in H^{-\frac{1}{2}}(\partial D_b)$, and using the definition of $\tilde{\psi}$, we get

$$ \frac{\varphi}{\mu} \int_{\partial D_b} \bar{f} ds - (\psi, f)_{H^{-\frac{1}{2}}(\partial D_b)} = 0, $$

for all $f \in H^{-\frac{1}{2}}(\partial D_b)$. Taking $f = \psi$ gives us $\frac{\varphi}{\mu} \int_{\partial D_b} \bar{\psi} ds - \|\psi\|_{H^{-\frac{1}{2}}(\partial D_b)}^2 = 0$. Recalling that $\mu\varphi + \int_{\partial D_b} \psi dx = 0$, we end up with $\|\psi\|_{H^{-\frac{1}{2}}(\partial D_b)}^2 + |\varphi|^2 = 0$, hence $\psi = 0$ and $\varphi = 0$. This concludes the proof. $\square$

We introduce the operator $\mathcal{G}^\mu : \overline{\mathcal{R}(\mathcal{H}^\mu)} \to L^2(S)$ such that: $\mathcal{G}^\mu(\varphi, \psi, u) := w^\infty$ where $w^\infty$ is the far field of $w^s \in H_{\text{loc}}^1(\mathbb{R}^m \setminus D_b)$ solution to:

$$ \left\{ \begin{array}{l} \Delta w + k^2 n w = 0 \text{ in } D_b, \\ \Delta w^s + k^2 w^s = k^2(1-n)(u + w^s) \text{ in } \mathbb{R}^m \setminus D_b, \\ w - w^s = \varphi \text{ on } \partial D_b, \\ \partial_\nu w - \partial_\nu w^s = \psi \text{ on } \partial D_b, \\ \lim_{r \to +\infty} \int_{|x|=r} \left| \frac{\partial w^s}{\partial r} - ikw^s \right|^2 = 0. \end{array} \right. \quad (46) $$

Observing that $w := u_g^s + v_g$ in $D_b$ and $w^s = u_g^s - u_{b,g}^s$ in $\mathbb{R}^m \setminus D_b$ satisfy (46) with $u = u_{b,g}|_{\Omega_D}$, $\varphi = u_{b,g}|_{\partial D_b}$ and $\psi = \partial_\nu u_{b,g}|_{\partial D_b}$, we then deduce the following factorization:

$$ \mathcal{F}^\mu = \mathcal{G}^\mu \mathcal{H}^\mu. \quad (47) $$
## 5.2. Properties of the range of $G^\mu$

For $z \in \mathbb{R}^m$, we set $\phi(., z) \in L_{\text{loc}}^2(\mathbb{R}^m)$, the radiating fundamental solution to the equation $\Delta\phi(., z) + k^2\phi(., z) = -\delta_z$ in $\mathbb{R}^m$ and denote by $\phi_z^\infty$ its far field pattern. We recall that:

$$ \phi(x, z) := \begin{cases} \frac{e^{ik|x-z|}}{4\pi|x-z|} & \text{in } \mathbb{R}^3, \\ \frac{i}{4} H_0^{(1)}(k|x-z|) & \text{in } \mathbb{R}^2, \end{cases} \qquad (48) $$

where $H_0^{(1)}$ denotes the Hankel function of the first kind of order zero.

**Lemma 7.** Assume $\mu$ is not an eigenvalue of problem (15). Then $\phi_z^\infty \in \mathcal{R}(G^\mu)$, the range of $G^\mu$, for all $z \in D_b$.

*Proof.* Let $z \in D_b$. Let $\chi_z$ be a $C^\infty(\mathbb{R}^m)$ function such that $\chi_z = 1$ in $\mathbb{R}^m \setminus D_b$ and $\chi_z|_V = 0$ where $V$ is a neighbourhood of $z$. Consider $\Theta_z := \chi_z\phi(., z) \in H_{\text{loc}}^2(\mathbb{R}^m)$. Let $W_z \in H_c^1(D_b)$ be the unique solution of

$$ \left\{ \begin{array}{l} \Delta W_z + k^2 n W_z = -(\Delta + k^2 n) \Theta_z \text{ in } D_b, \\ \mu W_z + \int_{\partial D_b} \partial_\nu W_z ds = 0 \text{ on } \partial D_b. \end{array} \right. \qquad (49) $$

We define $w_z = W_z + \Theta_z \in H^1(D_b)$. We notice that $w_z$ is solution of

$$ \left\{ \begin{array}{l} \Delta w_z + k^2 n w_z = 0 \text{ in } D_b, \\ \mu w_z + \int_{\partial D_b} \partial_\nu w_z ds = \mu \phi(., z) + \int_{\partial D_b} \partial_\nu \phi(., z) ds \text{ on } \partial D_b. \end{array} \right. \qquad (50) $$

Let us set $w := w_z$ in $D_b$ and $w^s := \phi(., z)$ in $\mathbb{R}^m \setminus D_b$. We shall prove that $w$ and $w^s$ satisfy (46) with $\varphi = (u_{b,z}^s + v_z)|_{\partial D_b}$, $\psi = (\partial_\nu u_{b,z}^s + \partial_\nu v_z)|_{\partial D_b}$ and $u = -\phi(., z)|_{\partial D_b}$ where

$$ v_z(x) := \int_{\partial D_b} (\partial_{\nu(y)} w_z(y) \phi(x, y) - \partial_{\nu(y)} \phi(x, y) w_z(y)) ds(y), $$

and $u_{b,z}^s$ the solution of problem (36) with $u^i = v_z$. Indeed, the first and second equations in (46) are clearly verified together with the radiation condition. It only remains to prove that the $w - w^s = \varphi$ and $\partial_\nu w - \partial_\nu w^s = \psi$ on $\partial D_b$.

From representation formula of solutions to the Helmholtz equation, we get the decomposition $w_z = w_z^s + v_z$ where

$$ w_z^s(x) := -k^2 \int_{D_b} (1-n)w_z(y)\phi(x,y)dy, \quad x \in D_b. \qquad (51) $$

We extend $w_z^s(x)$ for $x \in \mathbb{R}^m \setminus D_b$ using the same formula. Set $u_{b,z}^s := w_z^s - \phi(., z)$ in $\mathbb{R}^m \setminus D_b$. We notice that $u_{b,z}^s$ is solution to (36) with $u^i = v_z$. Thus we get on the boundary $\partial D_b$:

$$ \left\{ \begin{array}{l} w_z = w_z^s + v_z = u_{b,z}^s + \phi(., z) + v_z, \\ \partial_\nu w_z = \partial_\nu w_z^s + \partial_\nu v_z = \partial_\nu u_{b,z}^s + \partial_\nu \phi(., z) + \partial_\nu v_z, \end{array} \right. $$

which correspond to the desired transmission conditions since $w_s = \phi(., z)$. $\square$
**Lemma 8.** Assume $\mu$ is an eigenvalue of problem (15). Then the set of points $z$ for which $\phi_z^\infty \in \mathcal{R}(\mathcal{G}^\mu)$ is nowhere dense in $D_b$.

*Proof.* Let $w_0 \in H_c^1(D_b)$ be the eigenvector of (15) such that $\|w_0\|_{H^1(D_b)} = 1$. We will proceed by contradiction. Assume there exists $(\varphi, \psi, u) \in V$, such that,

$$ \mathcal{G}^{\mu}(\varphi, \psi, u) = \phi_z^{\infty}, $$

for all $z \in \mathcal{Z}$, a dense set in $D_b$. From the definition of $\mathcal{G}^\mu$ and the Rellich lemma, we obtain that $w_z := w$, where $w$ and $w^s$ verify the scattering problem (46) is a solution of (50). Let $\Theta_z \in H^2(D_b)$ defined as in the previous lemma and consider $W_z := w_z - \Theta_z \in H_c^1(D_b)$ that satisfies (49).

Multiplying the first equation in (49) by $\overline{w_0}$ and using twice the Green formula, yields the compatibility condition

$$ \int_{D_b} (\Delta \Theta_z + k^2 n \Theta_z) \overline{w_0} dx = 0, \qquad (52) $$

which gives, by applying again twice the Green formula, that

$$ v_{\mu}(z) := \int_{\partial D_b} (\partial_{\nu(x)}\phi(x, z)\overline{w_0}(x) - \partial_{\nu(x)}\overline{w_0}(x)\phi(x, z))ds(x) = 0, \qquad (53) $$

for all $z \in \mathcal{Z}$. Since $v_\mu$ satisfies the Helmholtz equation in $D_b$ and vanishes on a dense subset of $D_b$, we get $v_\mu = 0$ in $D_b$.

From representation formula of $w_0$, we have the following decomposition $\overline{w_0} = v_\mu + w_0^s$ in $D_b$, where

$$ w_0^s(x) := -k^2 \int_{D_b} (1-n) \overline{w_0(y)} \phi(x, y) dy, \quad x \in \mathbb{R}^m. \qquad (54) $$

Hence, $w_0^s$ defines a radiating solution of (4) with $D = D_b$ and $\psi = v_\mu$ in $D_b$. Since $v_\mu = 0$, we then conclude by uniqueness of solutions to (4) that $w_0^s = 0$ in $\mathbb{R}^m$ and then $\overline{w_0} = 0$, which gives a contradiction. $\square$

### 5.3. The GLSM Theorem

The previous lemmas allow us to apply Theorem (14) with $X = X^* = L^2(S)$ and $Y = \mathbb{R} \times H^{-\frac{1}{2}}(\partial D_b) \times L^2(\Omega_D)$ to the modified far field operator $\mathcal{F}^\mu$.

However, we need to construct an operator $B$ that satisfies Assumption 13, with $H := \mathcal{H}^\mu$ defined in (5.1). Keeping an implementable inverse method in mind, the operator $B$ must use computable quantities, such as $Fg$, $v_g|_{D_b}$ or $F_b^\mu g$.

**Lemma 9.** If $k$ is not a transmission eigenvalue and there exists $\beta > 0$ such that $n - 1 \ge \beta$ (or $1 - n \ge \beta$) in a neighborhood of $\partial D$, then

$$ B(g) = |(Fg, g)_{L^2(S)}| + \|v_g\|_{H^1(D_b)}^2, \qquad (55) $$

satisfies Assumption 13.
Proof. We start our proof by mentioning that under the assumption on $k^2$, the operator $T$ given by (8) is coercive on $V_{\text{inc}}(D)$ (see [2]).

Let us start with the first implication. Assume that the quantity $|(Fg_n, g_n)_{L^2(S)}| + \|v_{g_n}\|_{H^1(D_b)}^2$ is bounded for a sequence $\{g_n\}$. The latter implies in particular (thanks to classical trace theorem) that $\|v_{g_n}\|_{H^{\frac{1}{2}}(\partial D_b)}$ and $\|\partial_\nu v_{g_n}\|_{H^{-\frac{1}{2}}(\partial D_b)}$ are bounded sequences. The well-posedness of (36) yields that the sequence $u_{b,g_n}^s$ is bounded in $H^1(K\setminus D_b)$ for any compact set $K$ containing $D_b$. Then, since $u_{b,g_n} = u_{b,g_n}^s + v_{g_n}$ in $\mathbb{R}^m\setminus D_b$, we also get that $\|u_{b,g_n}\|_{H^{\frac{1}{2}}(\partial D_b)}$ and $\|\partial_\nu u_{b,g_n}\|_{H^{-\frac{1}{2}}(\partial D_b)}$ are bounded sequences. From the coercivity of $T$ and the factorization $F = H^*TH$, we get:

$$ \|v_g\|_{L^2(\Omega_D)}^2 \le C |(Fg, g)_{L^2(S)}|, \qquad (56) $$

where $C$ is a constant independent from $g$. This implies that $\|v_{g_n}\|_{L^2(\Omega_D)}$ is bounded and again from the identity $u_{b,g_n} = u_{b,g_n}^s + v_{g_n}$ in $\mathbb{R}^m\setminus D_b$ we also deduce that $\|u_{b,g_n}\|_{L^2(\Omega_D)}$. This proves that $\|\mathcal{H}^\mu g_n\|_Y^2$ is bounded.

Now, for the reverse implication, assume $\|\mathcal{H}^\mu g_n\|_Y^2$ bounded. From the representation formula on the scattered field and from the fact that $v_g$ satisfies Helmholtz equation, we have for $x \in \mathbb{R}^m\setminus D_b$

$$ u_{b,g}(x) = \int_{\partial D_b} (u_b(y)\partial_{\nu(y)}\phi(x, y) - \partial_\nu u_b(y)\phi(x, y))ds(y). $$

Hence, from continuity properties of single and double layers potentials, ([19]) we infer that the scattered field $u_{b,g_n}^s$ is bounded in $H^1(K\setminus D_b)$, for any compact set $K$ containing $D_b$, and particularly $(u_{b,g_n}^s, \partial_\nu u_{b,g_n}^s)$ are bounded in $H^{\frac{1}{2}}(\partial D_b) \times H^{-\frac{1}{2}}(\partial D_b)$. Therefore, $(v_{g_n}, \partial_\nu v_{g_n})|_{\partial D_b}$ is a bounded sequence in $H^{\frac{1}{2}}(\partial D_b) \times H^{-\frac{1}{2}}(\partial D_b)$. The representation formula for solution to the Helmholtz equation in $D_b$ applied on $v_{g_n}$ gives the boundedness of $v_{g_n}$ in $H^1(D_b)$.

In addition, since $u_{b,g_n}$ is bounded in $L^2(\Omega_D)$, we also get $\|v_{g_n}\|_{L^2(\Omega_D)}$ bounded. We then deduce that $\|v_{g_n}\|_{L^2(D)}$ is bounded. From the continuity of the operator $T$ and the factorization (7), we finally conclude that $|(Fg_n, g_n)_{L^2(S)}|$ is bounded. $\square$

We are ready to apply Theorem 14 in Appendix Appendix B to the operator $\mathcal{F}^\mu$. We observe that if $\mu$ is not an eigenvalue of (15), then $\mathcal{F}^\mu$ is injective and by reciprocity relations for far fields, we also obtain that it has dense range.

Combining Lemma 8, 7 and Theorem 14, we conclude with the main theorem of this section:

**Theorem 10.** Assume that the hypothesis of Lemma 9 hold, with $B$ defined by (55). Consider the functional

$$ J_{\alpha}(\phi_{z}^{\infty}, g) := \alpha B(g) + \|\mathcal{F}^{\mu}g - \phi_{z}^{\infty}\|_{L^2(S)}^2 \quad \text{and set} \quad j_{\alpha}(\phi_{z}^{\infty}) := \inf_{g \in L^2(S)} J_{\alpha}(\phi_{z}^{\infty}, g). $$

Let $g_\alpha^z$ be a minimizing sequence defined by

$$ J_{\alpha}(\phi_{z}^{\infty}, g_{\alpha}^{z}) \le j_{\alpha}(\phi_{z}^{\infty}) + C\alpha, $$
where $C > 0$ is a constant independent of $\alpha$. Then a real number $\mu > 0$ is not an eigenvalue of (15) if and only if the set of points $z$ such that $B(g_z^\alpha)$ is bounded and $\|\mathcal{F}^\mu g_z^\alpha - \phi_z^\infty\|_{L^2(S)} \to 0$ as $\alpha \to 0$ is nowhere dense in $D_b$.

In the case of noisy data, yielding a noisy operator $F^\delta$ satisfying $\|F^\delta - F\| \le \delta$, the theorem above also applies under the following modifications (insertcitation): the functional $J_\alpha$ should be replaced by

$$J_\alpha^\delta(\phi_z^\infty, g) = \alpha(B^\delta(g) + \delta\|g\|^2) + \|(F^\delta - F_b^\mu)g - \phi_z^\infty\|_{L^2(S)}^2,$$
with $B^\delta(g) := |(F^\delta g, g)_{L^2(S)}|^2 + \|v_g\|_{H^1(D_b)}^2$.

**Theorem 11.** Assume that $F^\delta$ is compact with dense range and the hypothesis of Lemma 9 hold. The functional $J_\alpha^\delta(\phi_z^\infty, \cdot)$ admits a minimizer $g_\alpha^{z,\delta}$. Moreover a real number $\mu > 0$ is not an eigenvalue of (15) if only if the set of points $z$ such that

$$\lim_{\alpha \to 0} \lim_{\delta \to 0} (B^\delta(g_\alpha^{z,\delta}) + \delta\|g_\alpha^{z,\delta}\|^2) < \infty \quad \text{and} \quad \lim_{\alpha \to 0} \lim_{\delta \to 0} \|\mathcal{F}^\mu g_\alpha^{z,\delta} - \phi_z^\infty\|_{L^2(S)} = 0$$
is nowhere dense in $D_b$.

# 6. Numerical inversion algorithm and validation

## 6.1. The inversion algorithm

Motivated by the analysis above and inspired by the inversion algorithm suggested in [5], we propose the following scheme to build an indicator function for the refractive index $n$ of the media from the knowledge of the far field operator $F$ at a fixed frequency $k$. It can be formally synthesised by the following three steps:

(i) Let $\rho > 0$ be a given parameter. Choose $D_b$ to be the ball $B_\rho^y$ of center $y \in \mathbb{R}^m$ and radius $\rho$ satisfying the condition $k^2 < \eta_0(n, D_b)$.

(ii) Evaluate the eigenvalue $\mu(y, n)$ of problem (15) from the measurements $F$ and the analytically computed $F_b^\mu$ using the GLSM method described by Theorem 14.

(iii) Plot the function $\mathcal{I}: y \to \mu(y, n) - \mu_{\text{ref}}(1)$ for $y$ sampling the probed domain containing the inclusions.

We now describe how to numerically implement the steps 2 and 3 for 2D problems. We will discuss the choice of $\rho$ and the influence of $k$ in the numerical validation subsection.

### 6.1.1. Synthetic generation of the data
Our numerical validating examples (in $\mathbb{R}^2$) use synthetic data $F$ that is numerically generated by solving the scattering problem (4) using a finite element method implemented using the Freefem++ [17] package. We use the Perfectly Matched Layer ([8]) technique to bound the computational domain and model the Sommerfeld radiation condition.

The outcome of our numerical solver is the matrix $F$ with entries

$$F_{pq} = u^\infty(\hat{x}_p, d_q), \quad 1 \le p, q \le N, \qquad (57)$$

where $\hat{x}_p = d_p = (\cos(\theta_p), \sin(\theta_p))$ with $\theta_p = \frac{p}{N}2\pi$ and where $u^\infty$ is the numerically computed far field.
### 6.1.2. Evaluation of the background operator
In order to evaluate the far field pattern $u_{b,y}^\infty$ associated with the background domain $D_b = B_\rho^y$, we shall rely on the following transformation of the far field pattern under translation [15, eq (5.3)]:

$$u_{b,y}^{\infty}(\hat{x}, d) = e^{iky \cdot (d-\hat{x})} u_{b}^{\infty}(\hat{x}, d) \qquad (58)$$

where $u_b^\infty$ is the far field associated with $D_b = B_\rho^0$ (i.e the ball of radius $\rho$ centered at the origin). The numerical evaluation of $u_b^\infty$ can be done using the numerical solver as described above or also analytically using separation of variables. We chose to use the latter option. Indeed, for $D_b = B_\rho^0$, the solution of equation (31) can be analytically expressed as:

$$u_b(r\hat{x}, d) = \sum_{j \in \mathbb{Z}} (\alpha_j H_j^{(1)}(kr) + i^j J_j(kr)) e^{ij(\theta_{\hat{x}} - \theta_d)} \text{ for } r \ge \rho,$$

with

$$\left\{ \begin{array}{l} \alpha_j = -i^j \frac{J_j(k\rho)}{H_j^{(1)}(k\rho)} \quad j \neq 0, \\ \alpha_0 = - \frac{\mu J_0(k\rho) + 2\pi\rho k J_0'(k\rho)}{\mu H_0^{(1)}(k\rho) + 2\pi\rho k H_0^{(1)'}(k\rho)}, \end{array} \right.$$

where $H_j^{(1)}$ denotes the Hankel function of the first kind of order $j$. Then the far field pattern assumes the following analytic expression:

$$u_{b,0}^{\infty}(\hat{x}, d) = \sqrt{\frac{2}{\pi k}} e^{i(-j\frac{\pi}{2} - \frac{\pi}{4})} \sum_{j \in \mathbb{Z}} \alpha_j e^{ij(\theta_{\hat{x}} - \theta_d)}, \qquad (59)$$

with $\hat{x} = (\cos(\theta_{\hat{x}}), \sin(\theta_{\hat{x}}))$ and $d = (\cos(\theta_d), \sin(\theta_d))$.

The numerical outcome of (59) and (58) is the background far field matrix $\mathbb{F}_{b,y}^\mu$ which entries are:

$$(\mathbb{F}_{b,y}^\mu)_{pq} = u_{b,y}^\infty(\hat{x}_p, d_q), \quad 1 \le p, q \le N, \qquad (60)$$

where $u_{b,y}^\infty$ is evaluated using (58) and where $u_{b,0}^\infty$ is approximated by a truncating of the sum in (59), keeping the indices $j \in \mathbb{Z}$ such that $|j| < M$.

### 6.1.3. Implementation of the GLSM method to identify $\mu(y, n)$
For the implementation of the GLSM method, one needs to introduce a discretization of $\|v_g\|_{H^1(D_b)}^2$ where $D_b = B_\rho^y$ for some given $y$ in $\mathbb{R}^2$. This norm can be computed using separation of variables. We first expand $g \in L^2(S)$ in Fourier series. More precisely, setting

$$g(\theta) := \sum_{m \in \mathbb{Z}} \hat{g}_m^y e^{im\theta} e^{-ik\hat{x}(\theta)y},$$

with $\hat{x}(\theta) = (\cos(\theta), \sin(\theta))$ and

$$\hat{g}_m^y := \frac{1}{2\pi} \int_0^{2\pi} g(\theta) e^{ik\hat{x}(\theta)y} e^{-im\theta} d\theta,$$

we have:

$$v_g(x + y) = 2\pi \sum_{m \in \mathbb{Z}} i^m \hat{g}_m^y J_m(kr) e^{im\theta}, \text{ for } x = r(\cos(\theta), \sin(\theta)). \qquad (61)$$
Then we deduce:

$$ \|v_g\|_{H^1(D_b)}^2 = 8\pi^3 \sum_{m \in \mathbb{Z}} |\hat{g}_m^y|^2 \left( (1+m^2) \int_0^R J_m(kr)^2 r dr + k^2 \int_0^R J'_m(kr)^2 r dr \right) $$

In order to simulate noise in the data, we change the values of the synthetic data $\mathbb{F}$ by adding random noise of level $\delta$ to construct the noisy far field matrix $\mathbb{F}^\delta = \mathbb{F} \cdot (1 + \delta(A + iB))$ where the entries of the matrices $A$ and $B$ are uniformly distributed real values in $[-1, 1]$ and where $\cdot$ denotes the element-wise product of matrices.

Since in our experiments, the refractive index $n$ is real, when $k^2$ is not a classical transmission eigenvalue, it is proved in [1] that we can substitute the term $|(Fg, g)_{L^2(S)}|$ in $B(g)$ by $((F^*F)^{\frac{1}{2}}g, g)_{L^2(S)} = \|(F^*F)^{\frac{1}{4}}g\|_{L^2(S)}^2$. For $z \in D_b$, considering $\phi_z^\infty := (\phi_z^\infty(\theta_1), \dots, \phi_z^\infty(\theta_N)) \in \mathbb{C}^N$, the cost function $J_\alpha$ can be approximated as

$$ J_\alpha(\phi_z^\infty, g) := \alpha(\|\mathbb{H}_y g\|_{\mathbb{C}^{2M+1}}^2 + \|((\mathbb{F}^\delta)^*\mathbb{F}^\delta)^{\frac{1}{4}}g\|_{\mathbb{C}^N}^2 + \delta\|g\|_{\mathbb{C}^N}^2) + \|(\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)g - \phi_z^\infty\|_{\mathbb{C}^N}^2 $$

for $g \in \mathbb{C}^N$ where $\|\cdot\|_{\mathbb{C}^N}$ denotes the euclidian norm in $\mathbb{C}^N$ and the matrix $\mathbb{H}_y$ in $\mathbb{C}^{(2M+1)\times N}$ defined by:

$$ (\mathbb{H}_y g)_m = \hat{g}_m^y \sqrt{8\pi^3((1+(m-M)^2) \int_0^R J_{m-M}(kr)^2 r dr + k^2 \int_0^R J'_{m-M}(kr)^2 r dr)}, $$

with:

$$ \hat{g}_m^y := \frac{1}{2\pi N} \sum_{\ell=1}^{N} g_\ell e^{ik\hat{x}(\theta_\ell) \cdot y} e^{-i(m-M)\theta_\ell}, \quad (62) $$

for $m = 0, \dots, 2M$.

For $y$ fixed, the minimizing sequence of Theorem 11 is actually computed as the minimizer $g_y^*(\alpha, z, \mu)$ of the function $J_\alpha(\phi_z^\infty, \cdot)$ which is also solution to the normal equation:

$$ (\alpha(D_y^*\mathbb{H}_y + ((\mathbb{F}^\delta)^*\mathbb{F}^\delta)^{\frac{1}{2}} + \delta I) + (\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)^*(\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu))g_y^*(\alpha, z, \mu) = (\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)^*\phi_z^\infty. $$

The theoretical result of Theorem 11 does not indicate how to choose $\alpha$ in practice. Inspired by the classical Linear Sampling Method (insertref), we employ the Morozov principle (insertref) by seeking $\alpha$ that ensures the equality:

$$ \|(\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)g_\alpha - \phi_z^\infty\|^2 = \delta\|g_\alpha\|^2, $$

for $g_\alpha$ being the solution to the normal equation:

$$ (\alpha I + (\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)^*(\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)) g_\alpha = (\mathbb{F}^\delta - \mathbb{F}_{b,y}^\mu)^* \phi_z^\infty. \quad (63) $$

### 6.1.4. Implementation of the inversion algorithm
Let $\mathcal{Y}$ be an uniform grid in the probed domain and $\Lambda$ a uniform discretization of an interval of positive values $\mu$ that contains $\mu_{\text{ref}}(1)$ (see (18)), the reference eigenvalue (ie, the eigenvalue of problem (15) with $n|_{D_b} = 1$). Beforehand, we precompute a collection of far field matrices $\mathbb{F}_{b,y}^\mu$ for each $(y, \mu) \in \mathcal{Y} \times \Lambda$ (recall that, for fixed $\mu$, $\mathbb{F}_{b,y}^\mu$ can be computed from $\mathbb{F}_{b,0}^\mu$ using the relation (58)).
For each $y \in \mathcal{Y}$, we compute the minimizer $g_y^*(\alpha, z)$ for $z \in Z$ a finite set of points in inside $D_b$ of cardinal $|Z|$. We then set as an approximation of the eigenvalue $\mu(y, n)$ the value of $\mu$ where the function

$$P_y(\mu) := \frac{1}{|Z|} \sum_{z \in Z} \|g_y^*(\alpha, z, \mu)\|_{\mathbb{C}^N}^2, \quad (64)$$

attains its maximum. To give an illustration of the behavior the function $\mu \to P_y(\mu)$, we display in Figure 4 this function in the case where $D = D_b = B_\rho^0$ with $\rho = 0.5$, $n|_{D_b} = 1$ and $k = 2$. We observe that the function has a sharp peak at the expected eigenvalue $\mu_{\text{ref}}(1)$ for $k = 2$ and $\rho = 0.5$.

After obtaining the numerical approximation of $\mu(y, n)$, we plot the function $\mathcal{I}(y) =$

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The indicator function $P_y(\mu)$ given by (64) for $D = D_b = B_\rho^0$ with $\rho = 0.5$, $n|_{D_b} = 1$ and $k = 2$.</figcaption></figure>

$$\mu(y, n) - \mu_{\text{ref}}(1).$$

We recall that thanks to the monotonicity property of $\mu(y, n)$ with respect to $n$ as stated in Theorem 4, $\mathcal{I}(y) > 0$ if $n|_{B_\rho^y} > 1$ and $\mathcal{I}(y) < 0$ if $n|_{B_\rho^y} < 1$. Furthermore, the larger is the value of $n|_{B_\rho^y}$, the larger is the value of $\mathcal{I}(y)$.

## 6.2. Numerical validation

### 6.2.1. The case of extended scatterers

Before treating the motivating example given in the introduction, we first show that the inverse method can also work in the case of extended scatterers. We consider in this example four circular inclusions of radii 0.3 and having four different constant refractive indices as presented in Figure 5 (right). The indicator function $\mathcal{I}$ is computed on an uniform grid $\mathcal{Y}$ of $50 \times 50$ sampling points of $[-2, 2] \times [-2, 2]$ and we take $\rho = 0.3$ for the artificial background (see Figure 5 (left)). The interval $\Lambda$ for identifying $\mu(y, n)$ is set to be equal to $[0, \mu_{\text{ref}}(2)]$. Figure 6 (left) presents the values of the function $\mathcal{I}: \mathcal{Y} \to \mathbb{R}^+$ for $k = 4$. We indeed distinguish 4 scattering circular objects. We observe from the sign of the indicator function, that the objects on the left have indeed a refractive index smaller than 1 and greater than
1 on the right. Thanks to the monotonicity property of the eigenvalue (Theorem 4), it is then possible to ordinate the 4 different refractive indices from the values of the indicator function $\mathcal{I}$.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left: Sampling grid $\mathcal{Y}$ formed by 50 × 50 points (blue points) and $D_b = B_\rho^0$ with $\rho = 0.3$ (green). Right: Four diffracting circles of radius 0.3 associated with four different values of the refractive index $n = 0.25$ (bottom left), $n = 0.5$ (top left), $n = 1.5$ (bottom right) and $n = 2$ (upper right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Indicator function $\mathcal{I}$ on the grid $\mathcal{Y}$ indicated in Figure 5 (left) for the scatterers shown in Figure 5 (right). We chose the parameter $\rho = 0.3$ and the noise level to be $\delta = 1\%$. Left: Reconstruction obtained for wavenumber $k = 4$. Right: Reconstruction obtained for wavenumber $k = 5.5$.</figcaption></figure>

*Discussion on the choice of the parameters k and ρ:* In the particular case of $D = D_b$ and constant $n$, the analytical expression of the eigenvalue (18) shows that its sensitivity
with respect to the refractive index is correlated with the value of the product $k\rho$. The behavior of the eigenvalue (18) is directly linked to the function $x \rightarrow x \frac{J_1(x)}{J_0(x)}$, which derivative is an increasing function that goes to infinity as $x \rightarrow j_0$, the first zero of $J_0$. Therefore the product $k\rho$ must be chosen with 3 constraints in mind:

1) The condition $0 < k^2 < \eta_0(n, D_b)$ which can be rewritten as $k\rho\sqrt{n} < j_0$ must be satisfied.

2) If the product $k\rho\sqrt{n}$ is close to $j_0$, the sensitivity of $\mu(n, y)$ with respect to $n$ would be very high, and only the regions with the largest values of $n$ will be visible. Furthermore, a small error in the data can lead to a significant error in the computed eigenvalue.

3) If the product $k\rho\sqrt{n}$ is far from $j_0$, the sensitivity would be smaller. Different values of $n$ may become more challenging to differentiate.

To illustrate this sensitivity, we compare the result of Figure 6 (left) with the indicator function $\mathcal{I}$ for $k = 5.5$ keeping the same radius $\rho = 0.3$ (i.e., a larger product $k\rho$) shown in Figure 6 (right). In this case, $k^2$ is close to the Dirichlet eigenvalue (12) for $D_b = B_\rho^0$ with constant $n = 2$. Although the values of $\mathcal{I}$ are coherent with the theory, it becomes challenging to visually discern the obstacles with $n < 1$ as some values of $\mathcal{I}$ get large inside the inclusion with refractive index $n = 2$.

### 6.2.2. The case of densely distributed small scatterers
We consider in this example the configuration given in the introductory section shown in Figure 1 (left) where we have a large number of small scatterers densely distributed on the probed domain with a constant refractive index $n = 2$. For the inverse algorithm, we keep the sets $\mathcal{Y}$ and $\mathcal{A}$ the same as in the previous example. Although the eigenvalue $\eta_0(n, D_b)$ is unknown, since it is monotonically decreasing with respect to $n$, using that $n|_{D_b} \le 2$, one obtains that $\eta_0(2, D_b) = \frac{j_0^2}{2\rho^2}$ is a lower bound of $\eta_0(n, D_b)$. Selecting $\rho$ such that $0 < k^2 < \eta_0(2, D_b)$ will ensure the hypothesis on the wavenumber needed in Theorem 4. We chose the radius $\rho$ of the background $D_b$ in order to have a high sensitivity of the eigenvalue to visually highlight the areas of higher concentration. Several numerical experimentations suggested that the choice $\rho = \sqrt{\frac{3\eta_0(2, D_b)}{k}}$ would give satisfactory results and this is the value of $\rho$ that is used in the following results. Figure 7 presents the function $\mathcal{I}$ for two wavenumbers: $k = 3$ (left) and $k = 4$ (right). The indicator function manages in both cases to distinguish 5 main regions that correspond to the 5 areas of higher concentration. The obtained result have indeed better qualitative aspect in relation to the local density of the small scatterers than the one given by the Linear Sampling Method in Figure 1 (right).

A second example of densely distributed obstacles is shown in Figure 8 where we have a fewer number of scatterers. We observe that the indicator function works even better in this case.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Indicator function $\mathcal{I}$ on the grid $\mathcal{Y}$ indicated in Figure 5 (left) for the scatterers shown in Figure 1 (left). We chose the parameter $\rho = \sqrt{\frac{3\eta_0(2, D_b)}{k}}$ and the noise level to be $\delta = 1\%$. Left: Reconstruction obtained for wavenumber $k = 3$ and $\rho = 0.57$. Right: Reconstruction obtained for wavenumber $k = 4$ and $\rho = 0.37$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: the domain $D$ constituted by 1257 small circles of radius 0.02 having a constant index of refraction $n = 2$. Right: Indicator function $\mathcal{I}$ on the grid $\mathcal{Y}$ indicated in Figure 5 (left). The noise level is $\delta = 1\%$.</figcaption></figure>

## Appendix A. Complementary technical results

We first give a proof of a technical result used in Theorem 4.

**Lemma 12.** If $k^2 < \eta_0(n, D_b)$, then there exists $\Lambda > 0$

$$ \int_{D_b} |\nabla w|^2 dx - k^2 \int_{D_b} n|w|^2 dx + \Lambda \int_{\partial D_b} |w|^2 \geq c \|w\|_{H^1(D_b)}^2, $$

for all $w \in H_c^1(D_b)$ and for some constant $c > 0$.
Proof. Assume by contradiction there is no $\Lambda > 0$ such that the (12) holds for $k^2 < \eta_0(n, D_b)$.

Let $(\Lambda_j)_j$ an increasing unbounded sequence of positive constant. The contradiction argument ensures that for each $\Lambda_j$, there exists $w_j \in H^1(D_b)$ with $\|w_j\|_{H^1(D_b)} = 1$ such that

$$ \int_{D_b} |\nabla w_j|^2 dx - k^2 \int_{D_b} n|w_j|^2 dx + \Lambda_j \int_{\partial D_b} |w_j|^2 ds \le \frac{1}{j}. \quad (A.1) $$

This gives for instance

$$ \int_{D_b} |\nabla w_j|^2 dx + \Lambda_j \int_{\partial D_b} |w_j|^2 ds \le k^2 \int_{D_b} n|w_j|^2 dx + \frac{1}{j}. $$

One can assume that, up to a subsequence, $w_j$ converges to some $w$, weakly in $H^1(D_b)$ and strongly in $L^2(D_b)$.

The term $\Lambda_j \int_{\partial D_b} |w_j|^2 ds$ must remain bounded, leading to $w = 0$ in $\partial D_b$. Hence $w \in H_0^1(D_b)$. Since the norm is weakly lower semi continuous, we have

$$ \|\nabla w\|_{L^2(D_b)}^2 \le \lim_{j \to +\infty} \inf \int_{D_b} |\nabla w_j|^2 dx \le \lim_{j \to +\infty} k^2 \int_{D_b} n|w_j|^2 dx = k^2 \int_{D_b} n|w|^2 dx, $$

which contradicts

$$ k^2 < \inf_{w \in H_0^1(D_b), w \ne 0} \frac{\|\nabla w\|_{L^2(D_b)}^2}{\int_{D_b} n|w|^2 dx} = \eta_0(n, D_b). $$

□

## Appendix B. Analytical Framework for GLSM

Let $X$ and $Y$ be two Hilbert spaces. We consider a bounded linear operator $F: X \to X$ that has dense range and can be factorized as $F = GH$ where $H: X \to Y$ and $G: \overline{R(H)} \subset Y \to X$ are bounded linear operators with $\overline{R(H)}$ being the closure of the range of $H$ in $Y$. In addition let $B: X \to \mathbb{R}^+$ be a continuous functional that satisfies the following assumption.

**Assumption 13.** Given a sequence $\{g_n\} \in X$, the sequence $\{B(g_n)\}$ is bounded if and only if the sequence $\{\|Hg_n\|_Y\}$ is bounded.

For a given parameter $\alpha > 0$ and $\phi \in X$, we consider the following cost functional

$$ J_{\alpha}(g, \phi) := \alpha B(g) + \|Fg - \phi\|_{X}^{2} $$

This cost functional has no minimizer in general, however its positivity implies that we can define $j_{\alpha}(\phi) := \inf_{g \in X} J_{\alpha}(g, \phi)$.

The central theorem of the GLSM is the following characterization of the range of $G$ in terms of $F$ and $B$.
**Theorem 14.** In addition to Assumption 13 we assume that $F$ has dense range. Let $C > 0$ be a given constant independent of $\alpha$ and consider a minimizing sequence $\{g_\alpha\}$ of $J_\alpha$ such that:

$$J_{\alpha}(\phi, g_{\alpha}) \leq j_{\alpha}(\phi) + C\alpha$$

Then $\phi \in \mathcal{R}(G)$ if and only if the sequence $B(g_\alpha)$ is bounded as $\alpha \to 0$.

A proof can be found in [6].

## References

[1] L. AUDIBERT, Qualitative methods for heterogeneous media, theses, Ecole polytechnique X, Sept. 2015.

[2] ——, The generalized linear sampling and factorization methods only depends on the sign of contrast on the boundary., Inverse Problems & Imaging, 11 (2017).

[3] L. AUDIBERT, F. CAKONI, AND H. HADDAR, New sets of eigenvalues in inverse scattering for inhomogeneous media and their determination from scattering data, Inverse Problems, 33 (2017), pp. 1–30.

[4] L. AUDIBERT, L. CHESNEL, AND H. HADDAR, Transmission eigenvalues with artificial background for explicit material index identification, Comptes Rendus Mathematique, 356 (2018), pp. 626–631.

[5] L. AUDIBERT, L. CHESNEL, H. HADDAR, AND K. NAPAL, Qualitative indicator functions for imaging crack networks using acoustic waves, 2020.

[6] L. AUDIBERT AND H. HADDAR, A generalized formulation of the Linear Sampling Method with exact characterization of targets in terms of farfield measurements, Inverse Problems, 30 (2014).

[7] ——, A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements, Inverse Problems, 30 (2014), p. 035011.

[8] E. BÉCACHE, S. FAUQUEUX, AND P. JOLY, Stability of perfectly matched layers, group velocities and anisotropic waves., Journal of Computational Physics, 188 (2003), pp. 399–433.

[9] F. CAKONI, D. COLTON, AND H. HADDAR, Inverse Scattering Theory and Transmission Eigenvalues, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016.

[10] F. CAKONI, D. COLTON, S. MENG, AND P. MONK, Stekloff eigenvalues in inverse scattering, SIAM journal on applied mathematics, 76 (2016), pp. 1737–1763.

[11] S. COGAR, A modified transmission eigenvalue problem for scattering by a partially coated crack, Inverse Problems, 34 (2018), p. 115003.

[12] S. COGAR, D. COLTON, S. MENG, AND P. MONK, Modified transmission eigenvalues in inverse scattering theory, Inverse Problems, 33 (2017), p. 125002.

[13] D. COLTON, F. CAKONI, AND P. MONK, Nondestructive testing and target identification, (2016).

[14] D. COLTON AND A. KIRSCH, A simple method for solving inverse scattering problems in the resonance region, Inverse Problems, 12 (1996), pp. 383–393.

[15] D. COLTON AND R. KRESS, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences, Springer Berlin Heidelberg, 2013.

[16] R. GRIESMAIER AND B. HARRACH, Monotonicity in inverse medium scattering on unbounded domains, SIAM Journal on Applied Mathematics, 78 (2018), pp. 2533–2557.

[17] F. HECHT, New development in freefem++, Journal of Numerical Mathematics, 20 (2012).

[18] A. KIRSCH AND N. GRINBERG, The Factorization Method for Inverse Problems, Oxford University Press, 12 2007.

[19] W. C. MCLEAN, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge University Press, 2000.