# Fast Imaging of Local Perturbations in a Unknown Bi-Periodic Layered Medium

Fioralba Cakoni, Houssem Haddar, Thi-Phong Nguyen

► To cite this version:

Fioralba Cakoni, Houssem Haddar, Thi-Phong Nguyen. Fast Imaging of Local Perturbations in a Unknown Bi-Periodic Layered Medium. Journal of Computational Physics, 2024, 501, pp.112773. 10.1016/j.jcp.2024.112773. hal-04223269

HAL Id: hal-04223269

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

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
Fast Imaging of Local Perturbations in a Unknown
Bi-Periodic Layered Medium

Fioralba Cakoni* Houssem Haddar † and Thi-Phong Nguyen ‡

Abstract

We discuss a novel approach for imaging local faults inside an infinite bi-periodic layered medium in $\mathbb{R}^3$ using acoustic measurements of scattered fields at the bottom or the top of the layer. The faulted area is represented by compactly supported perturbations with erroneous material properties. Our method reconstructs the support of perturbations without knowing or reconstructing the constitutive material parameters of healthy or faulty bi-period layer; only the size of the period is needed. This approach falls under the class of non-iterative imaging methods, known as the generalized linear sampling method with differential measurements, first introduced in [2] and adapted to periodic layers in [25]. The advantage of applying differential measurements to our inverse problem is that instead of comparing the measured data against measurements due to healthy structures, one makes use of periodicity of the layer where the data operator restricted to single Floquet-Bloch modes plays the role of the one corresponding to healthy material. This leads to a computationally efficient and mathematically rigorous reconstruction algorithm. We present numerical experiments that confirm the viability of the approach for various configurations of defects.

# 1 Formulation of the Problem

We consider nondestructive testing of an infinite bi-periodic penetrable layer in $\mathbb{R}^3$ by means of acoustic waves. This is an important problem with growing interest since periodic structures are part of many fascinating modern technological designs with applications in (bio)engineering and material sciences. In many sophisticated devises the periodic structure is complicated or difficult to model mathematically, hence evaluating its Green's function, which is the fundamental tool of many imaging methods, is computationally expensive or even impossible. On the other hand, when looking for faults in such complex media, the option of reconstructing everything, i.e. both periodic structure and the defects, may not

<footnote>*Department of Mathematics, Rutgers University, Piscataway, NJ 08854-8019, USA (fc292@math.rutgers.edu)</footnote>

<footnote>†INRIA, CMAP, Ecole polytechnique, Université Paris Saclay, Route de Saclay, 91128 Palaiseau, France (houssem.haddar@polytechnique.edu)</footnote>

<footnote>‡Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102, USA (tphong.nguyen.k4@gmail.com)</footnote>
be viable. Here we propose an approach that reconstructs the support of local anomalies without knowing explicitly or reconstructing the constitutive material properties of the periodic layer, except for the size of the period. The support of local perturbations is visualized by means of the indicator function computable from scattering data, which leads to computationally efficient imaging method. Such a connection is made through a rigorous mathematical analysis of the scattering problem by the bi-periodic layer without and with local perturbations present.

To be more specific, let $\boldsymbol{x} := (x_1, x_2, x_3) \in \mathbb{R}^3$, and assume that the scattering media is a penetrable infinite layer that is periodic in the $x_1$ and $x_2$ variable with period $L_1$ and $L_2$, respectively. Given a 2-dimensional vector $\boldsymbol{L} = (L_1, L_2)$, we call a function $w$ defined in $\mathbb{R}^3$ $\boldsymbol{L}$-periodic if $w$ is periodic in $x_1$ and $x_2$ with periods $L_1$ and $L_2$, respectively. With this notation, we assume that the refractive index of the periodic layer $n_p \in L^\infty(\mathbb{R}^3)$ is $\boldsymbol{L}$-periodic with $\mathrm{Re}(n_p) > 0$ and $\mathrm{Im}(n_p) \ge 0$. Furthermore, we assume that there exists an $h > 0$ such that $n_p = 1$ for $|x_3| > h$, hence the support of $(n_p - 1)$ represents the bi-periodic layer of width $2h$. The scattering of a time harmonic incident field $u^i$ (to become precise later on) is governed by

$$ \left\{ \begin{array}{l} \Delta u + k^2 n_p u = 0 \quad \text{in } \mathbb{R}^3, \\ u \text{ is } \boldsymbol{L}\text{-periodic} \end{array} \right. \qquad (1) $$

where $u := u^i + u^s$ is the total field, $u^s$ is the scattered field and $k > 0$ is the wave number proportional to interrogating frequency. Scattering of time-harmonic acoustic or electromagnetic waves by periodic gratings, inhomogeneous layers or waveguides is a research topic that has received enormous amount of attention. Some monographs and influential papers among the vast literature in the topic are [1, 3, 4, 5, 12, 14, 15, 19, 20, 21, 24, 25, 27]. In these type of scattering problems, both the incident field and scatterer field are assumed to be periodic in the horizontal directions, but the periods need not be the same. In such situations, one may multiply the fields by a quasi-periodicity-factor to restore overall periodicity. This allows to pose the problem on the unit cell of periodicity, where an application of well-known techniques from analysis such as analytic Fredholm theory [14] and variational methods provide existence and uniqueness of solution. For the purpose of this study we will always assume that the scattering problem (1) is well-posed, and refer the reader to aforementioned references for more details.

However, our main interest is in the case when local perturbations are present inside the bi-periodic layer. We denote by $\omega$ the support of local perturbations, and suppose that $\omega$ is not necessarily connected, but has connected complement. Without loss of generality (up to a possible rearrangement of the cell), we may assume that $\omega$ is located in one period (see Figure 1). The refractive index of the bi-periodic layer together with perturbations, which is no longer periodic function, is denoted by $n \in L^\infty(\mathbb{R}^3)$, and it satisfies $\mathrm{Re}(n) \ge n_0 > 0$, $\mathrm{Im}(n) \ge 0$. Note that $\omega$ is the support of $n-n_p$. The well-posedness of the scattering problem for the locally perturbed bi-periodic layer is handled by considering it as rough layer due to loss of periodicity. For the analysis and numerical implementations of the scattering of waves by rough penetrable layers or gratings we referee the reader to [11, 9, 10, 22, 25, 26]. Our goal is to determine the support of the damaged region $\omega$ by using the measured scattered field outside the layer due to appropriate incident fields. The challenging task however is to resolve $\omega$ without an explicit knowledge of $n_p$ (which in practice can have a compli-
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the geometry for the L-periodic problem. The healthy bi-periodic layer consists of a homogenous layer occupying $-h \le x_3 \le h$ with periodically distributed inhomogeneities indicated by red balls. Blue balls indicate the compactly supported perturbations $\omega$ located in one period denoted by $\Omega_0$.</figcaption></figure>

Fig:ML-periodic

cated form) nor reconstructing it, but just using the fact that $n_p$ is L-periodic. In this context, the response of the periodic background does not need to be measured. It is replaced by the extraction of measurements associated with a single Floquet-Bloch mode to encode some differential behavior for an appropriately designed indicator functions. This extraction requires information only on the period size of the background. To achieve this, for technical reasons, we must replace the infinite bi-periodic layer with M actual periods truncation (containing the defective period) for $M = (M_1, M_2) \in \mathbb{N}^2$ sufficiently large and then extend it as ML-periodic layer. As it is shown in [18], this truncation is equivalent to approximating the problem in the Floquet-Bloch domain using uniform discretization of the Floquet-Bloch variable and a trapezoidal rule to approximate the discretized solution. However, it is important to notice that this technical process of truncation and ML-periodic extension is needed only for the analysis of derivation of the indicator function of $\omega$, but is not involved in the computation of this indicator function. We also remark that many available inversion methods (see e.g. [6] and [23]) do not require this technical step in the analysis, but all of them assume an explicit knowledge of the Green's function of the bi-periodic layered background, which is not the case for our method. Although it is desirable to remove this technical step in the analysis since it is not used in the algorithm, it could be seen as trade-off for not using any a priori information on $n_p$, except for the fact that it L-periodic. Hence, from now on the scattering by perturbed bi-periodic layer is replaced by the following problem for the total field $u := u^i + u^s$

$$ \left\{ \begin{array}{l} \Delta u + k^2 nu = 0 \quad \text{in } \mathbb{R}^3, \\ u \text{ is } ML\text{-periodic} \end{array} \right. \qquad (2) $$

ML-periodic

where $u^i$ is the probing incident field and $u^s$ is the scattered field. Here $M = (M_1, M_2) \in \mathbb{N}^2$
with the natural numbers $M_1$ and $M_2$ sufficiently large, refers to the number of periods we consider in the $x_1$ and $x_2$ directions respectively. Thanks to **ML**-periodicity, solving (2) is equivalent to solving it in the period

$$ \Theta := \bigcup_{\mathbf{m}=(m_1, m_2) \in \mathbb{Z}_M^2} \left[ -\frac{L_1}{2} + m_1 L_1, \frac{L_1}{2} + m_1 L_1 \right] \times \left[ -\frac{L_2}{2} + m_2 L_2, \frac{L_2}{2} + m_2 L_2 \right] \times \mathbb{R}. $$

Adapting the notation $[\mathbf{a}, \mathbf{b}] := [a_1, b_1] \times [a_2, b_2]$ for any two generic vectors $\mathbf{a}$ and $\mathbf{b}$ in $\mathbb{R}^2$, we can rewrite $\Theta$ equivalently as

$$ \Theta = [\mathbf{M}_L^{-}, \mathbf{M}_L^{+}] \times \mathbb{R}, $$

where $\mathbf{M}_L^{-} := [\lfloor -\frac{M_1}{2} \rfloor L_1 + \frac{L_1}{2}, \lfloor -\frac{M_2}{2} \rfloor L_2 + \frac{L_2}{2}]$, $\mathbf{M}_L^{+} := [\lfloor \frac{M_1}{2} \rfloor L_1 + \frac{L_1}{2}, \lfloor \frac{M_2}{2} \rfloor L_2 + \frac{L_2}{2}]$, and $\mathbb{Z}_M^2 := \{\mathbf{m} = (m_1, m_2) \in \mathbb{Z}^2, \lfloor -\frac{M_\ell}{2} \rfloor + 1 \le m_\ell \le \lfloor \frac{M_\ell}{2} \rfloor, \ell = 1, 2\}$, with $\lfloor \cdot \rfloor$ denoting the floor function and $\mathbb{Z}$ the set of integers.

Note that in (2) we still call $n$ the **ML**-periodic extension of $n|_{\Theta}$ and without loss of generality assume that the defective period is

$$ \Omega_0 := \left[-\frac{\mathbf{L}}{2}, \frac{\mathbf{L}}{2}\right] \times \mathbb{R} = \left[-\frac{L_1}{2}, \frac{L_1}{2}\right] \times \left[-\frac{L_2}{2}, \frac{L_2}{2}\right] \times \mathbb{R}. \quad (3) \quad \boxed{\text{Oz}} $$

We now specify the incident wave $u^i$ we will use in our algorithm. To this end, we consider down-to-up or up-to-down incident plane waves of the form

$$ u^{i, \pm}(\mathbf{x}, \mathbf{j}) = \frac{-i}{2\bar{\beta}_{\#}(\mathbf{j})} e^{i\alpha_{\#}(\mathbf{j})\cdot\bar{\mathbf{x}} \pm i\bar{\beta}_{\#}(\mathbf{j})x_3}, \quad \text{with } \mathbf{x} = (\bar{\mathbf{x}}, x_3) \in \mathbb{R}^2 \times \mathbb{R} \quad (4) \quad \boxed{\text{inc}} $$

where for each mode $\mathbf{j} = (j_1, j_2) \in \mathbb{Z}^2$

$$ \alpha_{\#}(\mathbf{j}) := \left( \frac{2\pi}{M_1 L_1} j_1, \frac{2\pi}{M_2 L_2} j_2 \right) \in \mathbb{R}^2 \quad \text{and} \quad \beta_{\#}(\mathbf{j}) := \sqrt{k^2 - |\alpha_{\#}(\mathbf{j})|^2}, \quad \text{Im}(\beta_{\#}(\mathbf{j})) \ge 0. $$

Note that, given expression (17) below, considering these plane waves is formally equivalent to illuminating the media with periodic point sources (see also [16]). In addition the scattered field $u^s$ is outgoing which is expressed by imposing a radiation condition in the form of Rayleigh expansions:

$$ \left\{ \begin{array}{l} u^s(\bar{\mathbf{x}}, x_3) = \sum_{\ell \in \mathbb{Z}^2} \hat{u}^s(\ell) e^{i(\alpha_{\#}(\ell) \cdot \bar{\mathbf{x}} + \beta_{\#}(\ell)(x_3-h))}, \quad \forall x_3 > h, \\ u^s(\bar{\mathbf{x}}, x_3) = \sum_{\ell \in \mathbb{Z}^2} \hat{u}^s(-\ell) e^{i(\alpha_{\#}(\ell) \cdot \bar{\mathbf{x}} - \beta_{\#}(\ell)(x_3+h))}, \quad \forall x_3 < -h, \end{array} \right. \quad (5) \quad \boxed{\text{C4:RDC}} $$

where the Rayleigh coefficients $\hat{u}^s(\ell)$ are given by

$$ \begin{align} \hat{u}^{s+}(\ell) &:= \frac{1}{M_1 L_1 M_2 L_2} \int_{[M_L^-, M_L^+]} u^s(\bar{\mathbf{x}}, h) e^{-i\alpha_{\#}(\ell)\cdot\bar{\mathbf{x}}} \mathrm{d}\bar{\mathbf{x}}, \tag{6} \\ \hat{u}^{s-}(\ell) &:= \frac{1}{M_1 L_1 M_2 L_2} \int_{[M_L^-, M_L^+]} u^s(\bar{\mathbf{x}}, -h) e^{-i\alpha_{\#}(\ell)\cdot\bar{\mathbf{x}}} \mathrm{d}\bar{\mathbf{x}}. \nonumber \end{align} $$
Recall that the $[M_L^-, M_L^+]$ is a rectangle on $x_1x_2$-plane which is restricted by $M_1$ periods along $x_1$ and $M_2$ periods along $x_2$ directions. Evidently, the area of the rectangle $[M_L^-, M_L^+]$ is $M_1L_1M_2L_2$. We shall use the notation

$$ \Theta^h := [M_L^-, M_L^+] \times ] - h, h[ \quad \text{and} $$

$$ \Gamma_M^h := [M_L^-, M_L^+] \times \{h\}, \quad \Gamma_M^{-h} := [M_L^-, M_L^+] \times \{-h\}. $$

We denote by $H_{\#}^1(\Theta^h)$ the restrictions to $\Theta^h$ of functions that are in the Sobolev space $H_{\text{loc}}^1(|x_3| \le h)$ and are **ML**-periodic. The space $H_{\#}^{1/2}(\Gamma_M^h)$ is then defined as the space of traces on $\Gamma_M^h$ of functions in $H_{\#}^1(\Theta^h)$ and the space $H_{\#}^{-1/2}(\Gamma_M^h)$ is defined as the dual of $H_{\#}^{1/2}(\Gamma_M^h)$ with similar definitions for $H_{\#}^{\pm 1/2}(\Gamma_M^{-h})$.

More generally we consider the following direct problem: given $f \in L^2(\Theta^h)$ find $w \in H_{\#}^1(\Theta^h)$ satisfying

$$ \Delta w + k^2nw = k^2(1-n)f \quad (7) \quad \text{eq:w} $$

together with the Rayleigh radiation condition (5). We remark that the solution $w \in H_{\#}^1(\Theta^h)$ of (7) can be extended to a function in $\Theta$ satisfying $\Delta w + k^2nw = k^2(1-n)f$, using the Rayleigh expansion (5), and hence by **ML**-periodicity to a solution in the entire $\mathbb{R}^3$. Note that the scattering problem for **L**-periodic layer (1) is equivalent to (7) where $w := u^s$, $f := u^i|_{\Theta}^h$ and $n := n_p$, whereas the scattering problem for **ML**-periodic layer (2) is equivalent to (7) where $w := u^s$ and $f := u^i|_{\Theta}^h$. Throughout the paper we make the following assumption:

Ass:nk

**Assumption 1.** The refractive index $n$ and $k > 0$ are such that (7) as well as (7) with $n$ replaced by $n_p$ are both well-posed for all $f \in L^2(\Theta^h)$.

For sufficient conditions that guarantee Assumption 1 we refer the reader to [26], [20], [25] and the references therein. If $\Phi(n_p; \cdot)$ is the fundamental solution to

$$ \left\{ \begin{array}{l} \Delta\Phi(n_p; \cdot) + k^2 n_p \Phi(n_p; \cdot) = -\delta_0(\cdot), \\ \Phi(n_p; \cdot) \text{ is } ML \text{-periodic,} \\ \text{and the Rayleigh radiation condition (5).} \end{array} \right. \quad (8) \quad \text{phi} $$

then the solution $w$ of (7) has the representation as

$$ w(x) = - \int_D \left( k^2(n_p - n)w + k^2(1-n)f \right) (y) \Phi(n_p; x - y) dy. \quad (9) \quad \text{eq:2forw} $$

Finally as it will become clear latter on, our inversion algorithm is well suited to the case when a period of healthy **L**-periodic layer contains several compactly supported inhomogeneities sitting in homogeneous structure see Figure 1, which is the case in many applications in periodic gratings. For sake of presentation we assume that the homogenous base structure of the **L**-period media has refractive index one. In this case, we introduced the following notations which will be used throughout the paper (we refer the reader to Figure 1 for visualization of some of these regions). From now on, we denotes by **ab** the element wise multiplication of two generic vectors $\mathbf{a} = (a_1, a_2)$ and $\mathbf{b} = (b_1, b_2)$, that is

$$ \mathbf{ab} = (a_1 b_1, a_2 b_2). $$
**not** Notation 1. Recalling that $\Omega_0$ is the period containing $\omega$, we make the notations:

$$D_p := \text{Supp}(n_p - 1) \quad D := \text{Supp}(n - 1) \quad \text{and} \quad \omega := \text{Supp}(n - n_p)$$

and assume that the exterior of each of $D_p$, $D$ and $\omega$ is connected. We call $O$ union of components of $D_p \cap \Omega_0$ that intersect $\omega$, and $O^c$ its complement in $D_p \cap \Omega_0$.

$$\Lambda := O \cup \omega, \quad \hat{D} := O^c \cup \Lambda.$$

Let $\nu_m := (mL, 0) \in \mathbb{R}^3$ be the translate vector $\Omega_0 \mapsto \Omega_m$ ($m$-th period) for $m \in \mathbb{Z}^2$, then

$$O_p := \bigcup_{m \in \mathbb{Z}_M^2} O + \nu_m, \quad O_p^c := \bigcup_{m \in \mathbb{Z}_M^2} O^c + \nu_m, \quad \Lambda_p := \bigcup_{m \in \mathbb{Z}_M^2} \Lambda + \nu_m, \quad \hat{D}_p := \bigcup_{m \in \mathbb{Z}_M^2} \hat{D} + \nu_m$$

which are $L$-periodic copies of the respective aforementioned regions. Finally we define $\omega^{\text{mis}} := \omega \setminus \bar{D}$ as the missing (possibly part of) components of $D_p$ and denote by $\omega_p^{\text{mis}}$ its $L$-periodic copies. Note that $D_p$ is periodic and $\hat{D}_p \supseteq D \cup D_p$.

We refer the reader to Figures 2-3 for an illustration of the type of defect $\omega$ considered in this paper and for an illustration of the notation introduced above.

<figure><img src="image_3.png" /><figcaption>Figure 2: Two dimensional illustration for $D$, $D_p$ and $\omega$</figcaption></figure>

F1

**Remark 1.** In the special case when the defect $\omega$ only consists of missing components of $D_p$ in $\Omega_0$, where $n = 1$, we have $\hat{D}_p = D_p$ and $\Lambda = O = \omega^{\text{mis}} \cap \Omega_0$.

**Remark 2.** We remark that the results presented can readily be extended to the case when the homogeneous base structure of the $L$-period media has refractive index given by a constant different from one. In this case, the free space with refractive index one is replaced by a flat layer with refractive index one in $|x_3| > h$ and constant different from one in $|x_3| \le h$.
<figure><img src="image_4.png" /><figcaption>Figure 3: Two dimensional illustration for the domains $O_p^c$, $\Lambda_p$ and $\hat{D}_p$ associated with the configuration in Figure 2.</figcaption></figure>

The imaging method discussed in this paper falls into the class of recently developed qualitative approaches to inverse scattering or otherwise referred to as non-iterative methods [7]. The underlining idea is to design an indicator function of the support of perturbation solely from the scattering data without involving any partial differential model. Hence the support of perturbation is reconstructed without knowing the physical properties of the perturbation. In their standard presentation these methods however require an expression for the Green's function. However, in many applications an accurate modeling of the background is difficult to compute, and one way to avoid this is to use differential measurements. This idea was first introduced in [2] where two sets of scattering data, one for the healthy structure and the other in a latter time, where mathematically analyzed to arrive at an indicator function for the support of any perturbation to the healthy based structure. This idea was adapted to periodic layers in [25] (see also [17]). The advantage of applying differential measurements to the inverse problem for periodic layer described above is that instead of comparing the measured data against measurements corresponding to healthy structure, one makes use of the periodic structure of the background. More precisely, the measurement operator restricted to single Floquet-Bloch modes plays the role of data operator corresponding to healthy background. As a result information on the support $\omega$ of the perturbation is obtained without knowing or reconstructing the geometry and physical properties of the periodic background. As it is seen later in this paper, the analysis of a non-standard boundary value problem for two elliptic partial differential equations, referred to as the interior transmission problem, is at the core of this comparison. The analysis of the interior transmission problem for the problem introduced above was completed in [8]. The same problem for anisotropic media is discussed in [26].

In this paper we provide a comprehensive presentation of this new imaging method for a locally perturbed bi-periodic layer in $\mathbb{R}^3$ and present first three-dimensional numerical reconstructions. In Section 2 we discuss the analytical aspects of the inversion method for the scattering problem we introduced above. In Section 3 we remark on the analytical main steps of the direct and inverse scattering problem for an anisotropic bi-periodic layer. Section
4 is dedicated to a comprehensive numerical study of reconstructions for various configuration of local perturbations in both cases of isotropic and anisotropic bi-periodic media.

## 2 The Inverse Problem

sec2

We start by defining precisely the scattering data. As described above we have two choices of interrogating waves. If we use down-to-up (scaled) incident plane waves $u^{i,+}(x; j)$ defined by (4), then the (measured) scattering data is given by the Rayleigh coefficients

$$ \hat{u}^{s^{+}}(\ell; j), \quad (j, \ell) \in \mathbb{Z}^2 \times \mathbb{Z}^2, $$

(that means the transmitters are under the layer whereas the receivers are above the layer), whereas if we use up-to-down (scaled) incident plane waves $u^{i,-}(x; j)$ defined by (4) then the (measured) scattering data is given by the Rayleigh coefficients

$$ \hat{u}^{s^{-}}(\ell; j), \quad (j, \ell) \in \mathbb{Z}^2 \times \mathbb{Z}^2, $$

(that means the transmitters are above the layer whereas the receivers are under the layer). The inverse problem reads: from a knowledge of Rayleigh sequences $\{\hat{u}^{s^{+}}(\ell; j)\}_{\ell \in \mathbb{Z}^2}$ due to all incident waves $u^{i,+}(x; j)$ for $j \in \mathbb{Z}^2$ (or Rayleigh sequences $\{\hat{u}^{s^{-}}(\ell; j)\}_{\ell \in \mathbb{Z}^2}$ due to all incident waves $u^{i,-}(x; j)$ for $j \in \mathbb{Z}^2$) determine $\omega = \text{Supp}(n - n_p)$ without knowing $n_p$ but only the fact that $n_p$ is bi-periodic with period $L = (L_1, L_2)$ and that the layer is situated between $x_3 = -h$ and $x_3 = h$, for some $h > 0$.

To fix the idea we consider the Rayleigh sequences $\{\hat{u}^{s^{+}}(\ell; j)\}_{\ell, j \in \mathbb{Z}^2 \times \mathbb{Z}^2}$ due to incident waves $u^{i,+}(x; j)$, and to simplify the notation we let $\hat{u}^{s}(\ell; j) := \hat{u}^{s^{+}}(\ell; j)$ and $u^{i}(x; j) := u^{i,+}(x; j)$ for all $\ell \in \mathbb{Z}^2$ and $j \in \mathbb{Z}^2$. This scattering data defines the so-called near field (or data) operator $N: \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$ (recall that $\ell^2(\mathbb{Z}^2)$ is the Hilbert space of square summable sequences in $\mathbb{Z}^2$) by

$$ \{N(\boldsymbol{a})\}_{\ell \in \mathbb{Z}^2} = \sum_{j \in \mathbb{Z}^2} \boldsymbol{a}(j) \hat{u}^s(\ell; j), \quad \text{for } \boldsymbol{a} = \{\boldsymbol{a}(j)\}_{j \in \mathbb{Z}^2} \in \ell^2(\mathbb{Z}^2). \qquad (10) \qquad \text{nfo} $$

This operator is a main ingredient of our imaging method, whose properties are used to derive an indicator function for $D$. To understand its connection to the bi-periodic layer, we define the (Herglotz) operator $\mathcal{H}: \ell^2(\mathbb{Z}^2) \to L^2(D)$ by

$$ \mathcal{H}\boldsymbol{a} := \sum_{j \in \mathbb{Z}^2} \boldsymbol{a}(j)u^i(\cdot; j)|_D. \qquad (11) \qquad \text{2defH} $$

Its adjoint $\mathcal{H}^*: L^2(D) \to \ell^2(\mathbb{Z}^2)$ is given by

$$ \mathcal{H}^*\varphi := \{\hat{\varphi}(j)\}_{j \in \mathbb{Z}^2}, \quad \text{where} \quad \hat{\varphi}_j := \int_D \varphi(x) \overline{u^i(\cdot; j)}(x) \,dx. \qquad (12) \qquad \text{2adjointH} $$
Let us denote by $H_{\text{inc}}(D)$ the closure of the range of $\mathcal{H}$ in $L^2(D)$. We then consider the (compact) operator $G: H_{\text{inc}}(D) \to \ell^2(\mathbb{Z}^2)$ defined by

$$ G(f) := \{\hat{w}(\ell)\}_{\ell \in \mathbb{Z}^2}, \qquad (13) \qquad \text{defG} $$

where $\{\hat{w}(\ell)\}_{\ell \in \mathbb{Z}^2}$ is the Rayleigh sequence of $w \in H_{\#}^1(\Theta^h)$ the solution of (7). Then by linearity $N: \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$ we have

$$ N(a) = G \mathcal{H}(a). \qquad (14) $$

The following properties of $G$ and $\mathcal{H}$ are crucial to our inversion method. To state them, we must recall the standard interior transmission problem: find $(u, v) \in L^2(D) \times L^2(D)$ such that $u - v \in H^2(D)$ and satisfying

$$ \left\{ \begin{array}{ll} \Delta u + k^2 n u = 0 & \text{in } D, \\ \Delta v + k^2 v = 0 & \text{in } D, \\ u - v = \varphi & \text{on } \partial D, \\ \partial(u - v)/\partial v = \psi & \text{on } \partial D, \end{array} \right. \qquad (15) \qquad \text{oitp} $$

for given $(\varphi, \psi) \in H^{3/2}(\partial D) \times H^{1/2}(\partial D)$ where $\nu$ denotes the outward normal on $\partial D$. The wave number $k$ is called a transmission eigenvalue if the homogeneous problem (15), i.e. with $\varphi = 0$ and $\psi = 0$, has non-trivial solutions. Up-to-date results on this problem can be found in [7, Chapter 3] where in particular one finds sufficient solvability conditions. Without loss of generality we may assume that $\partial D \cap \partial \Omega_0 = \emptyset$ where $\Omega_0$ is given by (3). If the boundary of $D$ intersects with the vertical sides of the boundary $\Omega_0$, then the previous interior transmission problem should be augmented with periodicity conditions on $\partial D \cap \partial \Theta$ (intersection of $D$ with horizontal sides of the boundary $\Omega_0$ causes no problems). Since this condition does not affect the assumptions on the solvability of the interior transmission problem (in $H^2(D)$ with periodic conditions on $\partial D \cap \partial \Theta$), we make the choice of simplifying this technicality. Thus in the sequel we make the following assumption.

**HypoLSM** **Assumption 2.** $\partial D \cap \partial \Omega_0 = \emptyset$, and the refractive index $n$ and the wave number $k > 0$ are such that (15) has a unique solution.

In particular, if $\text{Re}(n - 1) > 0$ or $-1 < \text{Re}(n - 1) < 0$ uniformly in a neighborhood of $\partial D$ inside $D$ the interior transmission problem (15) satisfies the Fredholm alternative, and the set of real standard transmission eigenvalues is discrete (possibly empty). Thus Assumption 2 holds as long as $k > 0$ is not a transmission eigenvalue.

The next three proposition proven in [17] summarize the main properties of the operators $\mathcal{H}$ and $G$ as they relate to the operator $N$ which is available from the measurements.

**lemHerg** **Proposition 1** (Lemma 3.3 in [17]). The operator $\mathcal{H}: \ell^2(\mathbb{Z}^2) \to L^2(D)$ is compact and injective. The closure of the range of $\mathcal{H}$ in $L^2(D)$, denoted by $H_{\text{inc}}(D)$, is given as

$$ H_{\text{inc}}(D) := \{v \in L^2(D) : \Delta v + k^2 v = 0 \text{ in } D\}. \qquad (16) \qquad \text{eq:forHinc} $$
Let us denote by $\hat{\Phi}(\cdot; z) := \{\hat{\Phi}(\ell; z)\}_{\ell \in \mathbb{Z}^2}$ for $z \in \Theta^h$ the Rayleigh sequences of $\Phi(1, z)$ given by (8) in the special case when $n_p \equiv 1$ (i.e. corresponding to the free space), whose entries are given by

$$ \hat{\Phi}(\ell; z) := \frac{i}{2M_1L_1M_2L_2\beta_{\#}(\ell)} e^{-i(\alpha_{\#}(\ell)\cdot\bar{z} - \beta_{\#}(\ell)|z_3-h|)}, \quad z = (\bar{z}, z_3). \quad (17) \quad \text{lunch} $$

**Theorem 2** (Theorem 3.5 in [17]). Under Assumptions 1 and 2 the operator $G: H_{\text{inc}}(D) \to \ell^2(\mathbb{Z}^2)$ defined by (13) is injective with dense range. Moreover $\hat{\Phi}(\cdot; z)$ belongs to $\mathcal{R}(G)$ if and only if $z \in D$.

Finally we introduce the solution operator $T: L^2(D) \to L^2(D)$ given by

$$ T f := k^2(n-1)(f + w|_D) \quad (18) \quad \text{defT} $$

with $w$ being the solution of (7). By construction we have the following relation

$$ G f = \mathcal{H}^* T f $$

which leads to the following symmetric factorization of $N$

$$ N = \mathcal{H}^* T \mathcal{H}, \quad (19) \quad \text{symfac} $$

This symmetric factorization applied to an appropriate operator given in terms of $N$ allows us to characterize $D$ in terms of the scattering data. To this end, let us define

$$ N_{\#} := |\text{Re}(N)| + |\text{Im}(N)| \quad (20) \quad \text{sharp} $$

where $\text{Re}(N) := \frac{1}{2}(N + N^*)$, $\text{Im}(N) := \frac{1}{2i}(N - N^*)$, and $N^* : \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$ denotes the adjoint of $N : \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$. Similarly, letting $T_{\#} := |\text{Re}(T)| + |\text{Im}(T)|$ we have the following result:

**Proposition 3** (Theorem 4.2 in [17]). Under Assumptions 1 and 2 we have that

$$ N_{\#} = \mathcal{H}^* T_{\#} \mathcal{H}, \quad (21) \quad \text{titifactN} $$

where $T_{\#}: L^2(D) \to L^2(D)$ is self-adjoint and coercive on $H_{\text{inc}}(D)$. Moreover, $z \in D$ if and only if $\hat{\Phi}(\cdot; z) \in \mathcal{R}((N_{\#})^{1/2})$.

Provided that Assumption 1 holds true, Proposition 3 provides a mathematically rigorous criteria to determine $D$ which includes the support of periodic inhomogeneities in the healthy bi-periodic layer together with the defective region from the data operator $N^\pm$. In particular the Picard's series expressing the range of $N_{\#}$ converges if and only if the sampling point $z$ is in $D$. For complicated periodic structure and relatively small defects or defects with peculiar location such as $\omega \subset D_p$ this is unsatisfactory. Our aim is to find an indication function for $\omega$ without knowing or recovering $D_p$ in the same spirit as above. To achieve this, in [25] introduced the idea of using a data (or near field) operator corresponding to one Floquet-Bloch mode and perform the same type of analysis to it as above. This operator plays the role of the data due to the scattering by the background used in the sampling
method with differential measurements in [2]. The main difference here is that this operator is computationally extracted from the measurement operator N without extra measurements. Indeed for a fixed $q \in \mathbb{Z}_M^2$ we consider only $q + \ell M$ Rayleigh coefficients of scattered waves which are generated by only incident waves $u^i(\cdot; q + jM)$, i.e., $\hat{u}^s(q + \ell M; q + jM)$.

Thus, single Floquet-Bloch mode data operator: $N_q: \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$ is defined by

$$ \{N_q(\mathbf{a})\}_{\ell \in \mathbb{Z}^2} = \sum_{j \in \mathbb{Z}^2} \mathbf{a}(j) \hat{u}^s(\mathbf{q} + \ell \mathbf{M}; \mathbf{q} + j \mathbf{M}), \quad \text{for } \mathbf{a} = \{\mathbf{a}(j)\}_{j \in \mathbb{Z}^2}. \qquad (22) \qquad \text{nfo} $$

The operators $N_q$ and $N$ are related through the projection operator $I_q: \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$ given by

$$ (I_q \mathbf{a})(\ell) = \begin{cases} \mathbf{a}(j), & \text{if } \ell = \mathbf{q} + j\mathbf{M} \\ 0, & \text{else.} \end{cases} \qquad (23) \qquad \text{iop} $$

with its adjoint $I_q^*: \ell^2(\mathbb{Z}^2) \to \ell^2(\mathbb{Z}^2)$ given by $(I_q^* b)(j) = b(q + jM)$.

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

Hence we have

$$ N_q \mathbf{a} = I_q^* N I_q \mathbf{a}. \qquad (24) \qquad \text{2defNq} $$

(In the figure on the left the Rayleigh coefficients that define $N_q$ are indicated by red squares).

From physical point of view, $N_q$ is associated with $\alpha_q$-quasi-periodic fields with period $L$, where $\alpha_q = \frac{2\pi}{ML}q$, since the sequence $N_q \mathbf{a}$ corresponds to the Fourier coefficients of the $\alpha_q$-quasi-periodic component of the scattered field. Recall that a function $u$ is called $\alpha_q$-quasi-periodic fields with period $L = (L_1, L_2)$ if

$$ u(\bar{\mathbf{x}} + jL, x_3) = e^{i\alpha_q \cdot (jL)} u(\bar{\mathbf{x}}, x_3), \quad \text{for all } j \in \mathbb{Z}^2. $$

To understand the operator $N_q$ we need the $\alpha_q$-quasi-periodic fundamental solution $\Phi_q(\cdot)$

$$ \Delta \Phi_q(\cdot) + k^2 \Phi_q(\cdot) = -\delta_0 \quad \text{in } \Omega_0. \qquad (25) \qquad \text{phiq} $$

Then $\Phi_q(\cdot - z)$ solves

$$ \Delta \Phi_q(\cdot - z) + k^2 \Phi_q(\cdot - z) = -\delta_z \quad \text{in } \Omega_0 $$

for all $z \in \mathbb{R}^d$, and the Rayleigh coefficients $\hat{\Phi}_q(\cdot; z)$ of $\Phi_q(\cdot - z)$ are given by

$$ \hat{\Phi}_q(j; z) = \begin{cases} \frac{i}{2L_1 L_2 \beta_{\#}(q+M\ell) \bar{z} - \beta_{\#}(q+M\ell)|z_3 \mp h|} e^{-i(\alpha_{\#}(q+M\ell)\bar{z} - \beta_{\#}(q+M\ell)|z_3 \mp h|)\ell} & \text{if } j = q + M\ell, \ell \in \mathbb{Z}^2, \\ 0 & \text{if } j \neq q + M\ell, \ell \in \mathbb{Z}^2. \end{cases} \qquad (26) \qquad \text{hhh} $$

Similarly to the Herglotz operator $\mathcal{H}$, the single Floquet-Bloch mode Herglotz operator $\mathcal{H}_q: \ell^2(\mathbb{Z}^2) \to L^2(D)$ is defined by

$$ \mathcal{H}_q \mathbf{a} := \mathcal{H} I_q \mathbf{a} = \sum_j \mathbf{a}(j) u^i(\cdot; q + jM)|_D, \qquad (27) \qquad \text{2defHq} $$
where we note that $\mathcal{H}_q a$ restricted to $D_p$ is also $\alpha_q$-quasi-periodic function with period $L$ (See Remark 3). Then the factorization (19) along with (24) immediately implies

$$N_q = \mathcal{H}_q^* T \mathcal{H}_q, \qquad (28) \qquad \textcolor{red}{2\text{FactfNq}}$$

where $T$ is defined by (18). The role of $G$ given by (13) with respect to $N_q$ is now played by $G_q : \overline{\mathcal{R}(\mathcal{H}_q)} \to \ell^2(\mathbb{Z}^2)$ by

$$G_q = \mathcal{H}_q^* T |_{\overline{\mathcal{R}(\mathcal{H}_q)}}. \qquad (29) \qquad \textcolor{red}{2\text{defGq}}$$

Note that $\overline{\mathcal{R}(\mathcal{H}_q)}$ is characterized in Lemma 4. Observing that

$$\varphi(\mathbf{j}; \bar{\mathbf{x}}) := e^{i\alpha_{\#}(\mathbf{j})\cdot\bar{\mathbf{x}}} = e^{\frac{2\pi}{ML}\mathbf{j}\cdot\bar{\mathbf{x}}}, \quad \mathbf{j} \in \mathbb{Z}$$

is a Fourier basic of $ML$ periodic function in $L^2(\Omega_M)$, we have that any $w \in L^2(\Omega_M)$ which is $ML$ periodic, has the expansion

$$w(\mathbf{x}) = \sum_{\mathbf{j} \in \mathbb{Z}} \hat{w}(\mathbf{j}, x_3)\varphi(\mathbf{j}; \bar{\mathbf{x}}), \quad \text{where} \quad \hat{w}(\mathbf{j}, x_3) := \frac{1}{M_1 L_1 M_2 L_2} \int_{\Omega_M} w(\mathbf{x}) \overline{\varphi(\mathbf{j}; \bar{\mathbf{x}})} d\bar{\mathbf{x}}, \qquad (30)$$

(here the line over $\varphi$ denotes the conjugation whereas $\bar{\mathbf{x}} = (x_1, x_2)$). Splitting $\mathbf{j}$ by module $M$ component by component (i.e. splitting $j_\ell$ by module $M_\ell$, $\ell = 1, 2$) we can arrange the expansion of $w$ as

$$w(\mathbf{x}) = \sum_{\mathbf{q} \in \mathbb{Z}_M^2} \sum_{\ell \in \mathbb{Z}^2} \hat{w}(\mathbf{q} + M\ell, x_3)\varphi(\mathbf{q} + M\ell; \bar{\mathbf{x}}), \qquad (31)$$

where $\varphi(\mathbf{q} + M\ell; \bar{\mathbf{x}})$ is $\alpha_q$-quasi-periodic with period $L$. We also have that

$$w_q := \sum_{\ell \in \mathbb{Z}^2} \hat{w}(\mathbf{q} + M\ell, x_3)\varphi(\mathbf{q} + M\ell; \bar{\mathbf{x}})$$

is $\alpha_q$-quasi-periodic with period $L$. Thus any $ML$-periodic function $w \in L^2(\Omega_M)$ can be decomposed

$$w = \sum_{\mathbf{q} \in \mathbb{Z}_M^2} w_q \qquad (32) \qquad \textcolor{red}{\text{for:decomp}}$$

where $w_q$ is $\alpha_q$-quasi-periodic with period $L$. Moreover, by the orthogonality of the Fourier basis $\{\varphi(\mathbf{j}; \cdot)\}_{\mathbf{j} \in \mathbb{Z}}$, we have that

$$\hat{w}_q(\mathbf{j}) = 0 \quad \text{if } \mathbf{j} \neq \mathbf{q} + M\ell, \quad \ell \in \mathbb{Z}^2 \quad \text{and} \quad \hat{w}(\mathbf{q} + M\ell) = \hat{w}_q(\mathbf{q} + M\ell) \qquad (33)$$

where $\hat{w}_q(\mathbf{j})$ the Rayleigh sequence of $w_q$ defined in (6). From the definition of $G_q$, we see that $G_q(f)$ is the Rayleigh subsequence of $\hat{w}(\mathbf{j})$ corresponding to the indices $\mathbf{j} = \mathbf{q} + M\ell, \ell \in \mathbb{Z}^2$, where $w$ is solution of (7).

The above discussion is helpful to proving the following properties for $\mathcal{H}_q$ and $G_q$ that are the counterpart results to Proposition 1 and Proposition 2 needed to analyze the range of the operator $N_q$.
**Proposition 4.** With reference to Notation 1, the operator $\mathcal{H}_q : \ell^2(\mathbb{Z}^2) \to L^2(D)$ is injective. Furthermore the closure of its range is given by $\overline{\mathcal{R}(\mathcal{H}_q)} = H_{\text{inc}}^q(D)$ where

$$H_{\text{inc}}^q(D) := \{v \in L^2(D), \Delta v + k^2 v = 0 \text{ in } D \text{ and } v \text{ is } \alpha_q\text{-quasi-periodic in } D_p\}. \quad (34) \quad \text{hq1.2}$$

**Remark 3.** We remark that if $\omega$ contains missing (part of) components of $D_p$ in the above we mean that $v$ has an extension as $\alpha_q$-quasi-periodic in $D_p$. More specifically this extension is defined as

$$H_{\text{inc}}^q(D) := \left\{ \begin{array}{ll} v \in L^2(D), \quad \Delta v + k^2 v = 0 \text{ in } D \quad \text{such that } \tilde{v}(\boldsymbol{x}) \in L^2(D \cup D_p) \text{ defined by:} \\ \tilde{v}(\boldsymbol{x}) = \left\{ \begin{array}{ll} v(\boldsymbol{x}), & \forall \boldsymbol{x} \in D \\ v(\boldsymbol{x} + \boldsymbol{L})e^{-i\boldsymbol{\alpha}_q \cdot \boldsymbol{L}}, & \forall \boldsymbol{x} \in D_p \setminus D \end{array} \right. \quad \text{is } \alpha_q\text{-quasi-periodic in } D_p \end{array} \right\}. \quad (35) \quad \text{hq1.2.3}$$

Note also that, in the case $D_p \subseteq D$ (that is studied in [8]), i.e. there are no missing components, the definition in (34) becomes

$$H_{\text{inc}}^q(D) := \{v \in L^2(D), \Delta v + k^2 v = 0 \text{ in } D \text{ and } v|_{D_p} \text{ is } \alpha_q\text{-quasi-periodic}\}. \quad (36) \quad \text{hq1}$$

*Proof.* Injectivity of $\mathcal{H}_q$ follows from injectivity of the operators $\mathcal{H}$ and $I_q$. Hence it suffices to show that $\mathcal{H}_q^*$ is injective in $H_{\text{inc}}^q(D)$. The case when $\omega \subset D$ corresponds exactly to Lemma 3.1 in [8]. We sketch here this proof to confirm that it works also in the case when $\omega \not\subset D$. In particular, $\omega$ is considered to be the whole component $\mathcal{O}$ or a part of the component $\mathcal{O}$. To this end, let $\varphi \in H_{\text{inc}}^q(D)$ and assume $\mathcal{H}_q^*(\varphi) = 0$. Then we define

$$u(\boldsymbol{x}) := \frac{1}{M_1 M_2} \int_D \Phi_q(\boldsymbol{x} - \boldsymbol{y})\varphi(\boldsymbol{y}) \, d\boldsymbol{y}.$$

where $\Phi_q$ has the Rayleigh coefficients given by (26). By definition of $u$ and using (26) we see that the Rayleigh coefficients of $u$ are given by $\hat{u}(\boldsymbol{j}) = 0$ for all $\boldsymbol{j} \neq \boldsymbol{q} + \boldsymbol{M}\boldsymbol{l}$ and $\hat{u}(\boldsymbol{q} + \boldsymbol{M}\boldsymbol{l}) = (\mathcal{H}^*(\varphi))(\boldsymbol{q} + \boldsymbol{M}\boldsymbol{l}) = (\mathcal{H}_q^*(\varphi))(\boldsymbol{l}) = 0$. Therefore $u$ has all zero Rayleigh coefficients, which implies that $u = 0$ for $x_3 > h$ and $x_3 < h$. We now observe that for all $\boldsymbol{y} \in D$, $\Delta\Phi_q(\cdot; \boldsymbol{y}) + k^2\Phi_q(\cdot; \boldsymbol{y}) = 0$ in the complement of $\hat{D}_p$. This implies that

$$\Delta u + k^2 u = 0 \quad \text{in } \mathbb{R}^3 \setminus \hat{D}_p.$$

Using a unique continuation argument we infer that $u = 0$ in $\Omega_M \setminus \hat{D}_p$. Therefore, $u \in H_0^2(\hat{D}_p)$ by the regularity of volume potentials. Remark that $\omega$ can contain more than one components for all possible locations corresponding to $\mathcal{O}$ (See Figure 2 for an illustration). For instance, $\omega$ may include components that are part of $\mathcal{O}$ but are missing in $D$. For configurations where this does not happen, we can proceed exactly in the same way as [8, pages 11-12 in the proof of Lemma 3.1]. This is why we concentrate here on the case where the defect occupies a region in $\mathcal{O}$ that is missing (or partially missing) in $D$. Specifically in this case, $\omega = D_p \setminus D$, and $n = 1$ in $\omega$. Therefore $\hat{D}_p = D_p$ and $D \cap \omega = \emptyset$. We then obtain $u \in H_0^2(D_p)$, and it can be rewritten $u$ as

$$u(\boldsymbol{x}) = \frac{1}{M_1 M_2} \int_{D_p} \Phi_q(\boldsymbol{x} - \boldsymbol{y})\varphi(\boldsymbol{y}) \, d\boldsymbol{y} - \frac{1}{M_1 M_2} \int_{\omega} \Phi_q(\boldsymbol{x} - \boldsymbol{y})\varphi(\boldsymbol{y}) \, d\boldsymbol{y},$$
where $\varphi$ is $\alpha_q$-quasi-periodic in $D_p$. We first see from the definition of $u(x)$ that $u(x)$ is $\alpha_q$-quasi-periodic. In fact,

$$ 
\begin{aligned}
u(x + mL) &:= \frac{1}{M_1 M_2} \int_D \Phi_q(x + mL - y)\varphi(y) \, dy \\
&= e^{i\alpha_q \cdot mL} \frac{1}{M_1 M_2} \int_D \Phi_q(x - y)\varphi(y) \, dy = e^{i\alpha_q \cdot mL} u(x),
\end{aligned}
$$

for all $m \in \mathbb{Z}_M^2$. We now consider $x \in D_p \cap \Omega_0$. By the $\alpha_q$-quasi-periodic of $\Phi_q(x)$ and $\varphi$ with period $L$, we have

$$ \frac{1}{M_1 M_2} \int_{D_p} \Phi_q(x - y)\varphi(y) \, dy = \int_{D_p \cap \Omega_0} \Phi_q(x - y)\varphi(y) \, dy, $$

and hence

$$ 
\begin{aligned}
u(x) &= \int_{D_p \cap \Omega_0} \Phi_q(x - y)\varphi(y) \, dy - \frac{1}{M_1 M_2} \int_{\omega} \Phi_q(x - y)\varphi(y) \, dy \\
&= \int_{\hat{D}_p \cap \Omega_0} \Phi_q(x - y)\varphi(y) \, dy + \frac{M_1 M_2 - 1}{M_1 M_2} \int_{\omega} \Phi_q(x - y)\varphi(y) \, dy, \quad \forall x \in D_p \cap \Omega_0
\end{aligned}
$$

Define $\tilde{\varphi}(x)$ for all $x \in D_p \cap \Omega_0$ as

$$ \tilde{\varphi}(x) = \begin{cases} \varphi(x), & \forall x \in \hat{D}_p \cap \Omega_0 \\ \frac{M_1 M_2 - 1}{M_1 M_2} \varphi(x), & \forall x \in \omega \end{cases} $$

we have

$$ \Delta u + k^2 u = -\tilde{\varphi}(x), \quad \forall x \in D_p \cap \Omega_0. $$

Evidently, extending $\tilde{\varphi}$ by $\alpha_q$-quasi-periodicity to $D_p$ we have

$$ \Delta u + k^2 u = -\tilde{\varphi}(x), \quad \forall x \in D_p. $$

Furthermore, the linearity of the Laplace operator and the fact that $\varphi(x)$ satisfies equation $\Delta\varphi + k^2\varphi = 0$ in $D_p$ implies that

$$ \Delta \tilde{\varphi} + k^2 \tilde{\varphi} = 0, \quad \forall x \in D_p \qquad (37) $$

Since $u \in H_0^2(D_p)$ we then have

$$ \|\tilde{\varphi}\|_{L^2(D_p)}^2 = - \int_{D_p} (\Delta u + k^2 u) \bar{\tilde{\varphi}} = 0, $$

which implies $\tilde{\varphi} = 0$ in $D_p$. In particular, $\tilde{\varphi} = 0$ in $D_p \cap \Omega_0$. If $M_1 M_2 = 1$, then the definition of $\tilde{\varphi}$ implies $\varphi = 0$ in $\hat{D}_p \cap \Omega = D$, which ends the proof. Otherwise, when $M$ is such that $M_1 M_2 \ge 2$ we have $\varphi = 0$ in $D_p \cap \Omega_0$, and therefore by the quasi-periodicity $\varphi = 0$ in $D_p \supset D$. This proves the injectivity of $(H^\pm)^*$ on $H_{\text{inc}}^q(D)$ and hence proves the Lemma. $\square$
In the case of missing component, we assume that $f \in H_{\text{inc}}^q(D)$ then $\tilde{f}$ (defined by the same way as $\tilde{v}$ in (34)) is $\alpha_q$-quasi-periodic in $D_p$. Furthermore, since $n=1$ in $\omega$, solution $w$ to equation (7) associated with data $f$ equals

$$ \Delta w + k^2 n w = k^2(1-n)\tilde{f} = k^2(1-n_p)\tilde{f} + k^2(n_p - n)\tilde{f} \quad (38) $$

Using decomposition (32) for $w$ the solution of (7) along with the facts that $n_p$ is periodic, $\tilde{f}$ is $\alpha_q$-quasi-periodic and that $n-n_p$ is compactly supported in one period $\omega \subset \Omega_0$, then (7) in terms of the coefficients $w_q$ in (32) for this $w$ takes the form

$$ \Delta w_q + k^2 n_p w_q = k^2(n_p - n)w + k^2(1 - n)\tilde{f} \quad \text{in } \Omega_0. \quad (39) $$

eq:wq0

However, since $f = \tilde{f}$ in the support of $n-1$, equation (39) is equivalent to or equivalent to

$$ \Delta w_q + k^2 n_p w_q = k^2(n_p - n)w + k^2(1 - n)f \quad \text{in } \Omega_0 \quad (40) $$

eq:wq

Therefore, operator $G_q: \overline{\mathcal{R}(\mathcal{H}_q)} = H_{\text{inc}}^q(D) \to \ell^2(\mathbb{Z}^2)$ can be equivalently defined as

$$ G_q(f) := I_q^* \{ \widehat{w_q}(\ell) \}_{\ell \in \mathbb{Z}^2}, \quad (41) $$

defGqequiv

where $w_q$ solves (40), $w$ is solution to (7), and $\{\widehat{w_q}(\ell)\}_{\ell \in \mathbb{Z}^2}$ is the Rayleigh sequence of $w_q$.

We need one more interior transmission problem to introduce, which is related to the characterization of the range of $G_q$.

nitp

**Definition 1.** Seek $(u, f) \in L^2(\Lambda) \times L^2(\Lambda)$ such that $u - f \in H^2(\Lambda)$ and

$$ \left\{ \begin{array}{ll} \Delta u + k^2 n u = k^2(n_p - n)\tilde{S}_k(f) & \text{in } \Lambda, \\ \Delta f + k^2 f = 0 & \text{in } \Lambda, \\ u - f = \varphi & \text{on } \partial\Lambda, \\ \partial(u - f)/\partial\nu = \psi & \text{on } \partial\Lambda, \end{array} \right. \quad (42) \quad \text{eqnew} $$

for given $(\varphi, \psi) \in H^{3/2}(\partial\Lambda) \times H^{1/2}(\partial\Lambda)$ where

$$ \tilde{S}_k : L^2(\Lambda) \to L^2(\Lambda) : $$

$$ f \mapsto - \int_{\Lambda} k^2(1-n_p)f(y) \left( \sum_{0 \neq m \in \mathbb{Z}_M^2} e^{i\alpha_q \cdot mL} \Phi(n_p; x - mL - y) \right) dy \quad (43) \quad \text{def:2Stild} $$

and $\Phi(n_p; \cdot)$ is the **ML**-periodic fundamental solution given by (8) and where $\nu$ denotes the unit normal on $\partial\Lambda$ outward to $\Lambda$.

This problem is first introduced in [17] and we refer the reader to this paper for the analysis of its solvability. We here make the following assumption on its solvability.

as-nitp

**Assumption 3.** The refractive index $n$ and $k > 0$ are such that the interior transmission problem in Definition 1 has a unique solution.
**Theorem 1.** Suppose that Assumptions 1, 2, 3 hold, and Assumption 2 holds when $(n, D)$ is replaced by $(n_p, D_p)$. Then the operator $G_q: H_{\text{inc}}^q(D) \to \ell^2(\mathbb{Z}^2)$ is injective with dense range.

*Proof.* As in the previous proof, we only consider the case where the defect is only constituted by components in $D_p$ that are missing in $D$, i.e. $\omega \subseteq \mathcal{O}$ and $\omega \cap D = \emptyset$. The other cases can be treated as in [8, Theorem 3.2]. Assume that $f \in H_{\text{inc}}^q(D)$ such that $G_q(f) = 0$ and $\tilde{f} \in L^2(D_p)$ is an extension of $f$ given in (34). Let $w$ be solution of (7) with data $f$. We have that

$$G_q(f) := I_q^* \{\widehat{w_q}(\ell)\}_{\ell \in \mathbb{Z}^2},$$

where and $w_q$ is solution to

$$\Delta w_q + k^2 n_p w_q = k^2(n_p - n)w + k^2(1 - n)f \quad \text{in } \Omega_0 \quad (44) \qquad \text{eq:2wq1}$$

Remark that $\Delta w_q + k^2 w_q = 0$ in $\Omega_M \setminus D_p$. Using a similar unique continuation argument as at the beginning of the proof of Lemma 4 we deduce that

$$w_q = 0 \quad \text{in } \Omega_M \setminus D_p.$$

In other words, $w_q = 0$ outside $D_p$. From now to the end of the proof, we consider only period $\Omega_0$. Recall that $\text{Supp}(n - n_p) \cap \text{Supp}(1 - n) = \omega \cap D = \emptyset$, and $n = n_p$ in $D$, equation (44) can be split into two equations

$$\Delta w_q + k^2 n_p w_q = k^2(n_p - 1)w \quad \text{in } \omega, \quad (45) \qquad \text{eq:2wq1.1}$$

where $w \in H_{\text{loc}}^1(\Theta)$ is solution to (7) and

$$\Delta w_q + k^2 n_p w_q = k^2(1 - n_p)f \quad \text{in } \Omega_0 \cap D. \quad (46) \qquad \text{eq:2wq1.2}$$

Observe that $w$ satisfies $\Delta w + k^2 w = 0$ in $\omega$. Therefore, from (45) $(w_q, -w) \in H_0^2(\omega) \times L^2(\omega)$ and verifies

$$\left\{ \begin{array}{ll} \Delta w_q + k^2 n_p w_q = k^2(n_p - 1)w & \text{in } \omega, \\ \Delta w + k^2 w = 0 & \text{in } \omega. \end{array} \right. \quad (47) \qquad \text{eq:2itpwq1}$$

Assumption 2 holds when $(n, D)$ is replaced by $(n_p, D_p)$ implies that equation (47) has a trivial solution and therefore

$$w_q = -w = 0 \quad \text{in } \omega.$$

We again see that, in the domain $\Omega_0 \cap D$, $(w_p, f) \in H_0^2(\Omega_0 \cap D) \times L^2(\Omega_0 \cap D)$ and verifies

$$\left\{ \begin{array}{ll} \Delta w_q + k^2 n_p w_q = k^2(1 - n_p)f & \text{in } \Omega_0 \cap D, \\ \Delta f + k^2 f = 0 & \text{in } \Omega_0 \cap D. \end{array} \right. \quad (48) \qquad \text{eq:2itpwq1}$$

Assumption 2 holds when $(n, D)$ is replaced by $(n_p, D_p)$ implies that equation (48) has a trivial solution, and therefore

$$w_q = f = 0 \quad \text{in } \Omega_0 \cap D.$$

This proves the lemma.

<figure><img src="image_6.png" /></figure>
Finally we can prove exactly in the same way as in the proof of [8, Theorem 3.5] the following range test related to the operator $G_q$ which will play an important role in the design of our imaging method.

**Theorem 2.** Suppose that Assumptions 1, 2 and 3 hold. Then, $I_q^*\hat{\Phi}_q(\cdot; z) \in \mathcal{R}(G_q)$ if and only if $z \in \hat{D}_p$.

We now have all the ingredients to arrive at our imaging method in the next section.

# 3 A Differential Imaging Algorithm

We now apply all the results of Section 2 to design an algorithm that provides us an indicator function of the support of the defect $\omega$ without reconstructing $D_p$ or computing the Green's function of the periodic media. Such an indicator function is based on the analysis of the range of operators $N$ and $N_q$ in relation to interior transmission problems discussed above. Roughly speaking we can construct three appropriate sequences $a^{\alpha,z}$, $a_q^{\alpha,z}$ and $\tilde{a}_q^{\alpha,z}$ (to become precise latter), as nearby solutions to

$$ \|N a^{\alpha,z} - \hat{\Phi}(\cdot; z)\|_{\ell^2} \le \alpha, \quad \|N a^{\alpha,z} - \hat{\Phi}_q(\cdot; z)\|_{\ell^2} \le \alpha, \quad \|N_q a^{\alpha,z} - I_q^* \hat{\Phi}_q(\cdot; z)\|_{\ell^2} \le \alpha $$

as $\alpha \to 0$, where $\hat{\Phi}(\cdot; z)$ are the Rayleigh coefficients of $\Phi(1, z)$ (i.e. $\Phi(n_p; z)$ defined by (8) with $n_p = 1$) given by (17) and $\hat{\Phi}_q(\cdot; z)$ are the Rayleigh coefficients of $\Phi_q(\cdot - z)$ given by (26). Then we show that these nearby solutions satisfy:

* $z \in D$ if and only if $\lim_{\alpha \to 0} \langle N_{\#} a^{\alpha,z}, a^{\alpha,z} \rangle < \infty$.

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

* $z \in D_p \setminus \overline{\omega_p^{mis}}$ if and only if $\lim_{\alpha \to 0} \langle N_{\#} a_q^{\alpha,z}, a_q^{\alpha,z} \rangle < \infty$.

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

* $z \in \hat{D}_p$ if and only if $\lim_{\alpha \to 0} \langle N_{q,\#} \tilde{a}_q^{\alpha,z}, \tilde{a}_q^{\alpha,z} \rangle < \infty$.

An appropriate combination of these three indicator functions yield a visualization of the support of local perturbation. In the following we introduce and mathematically justify the algorithm.
The domain $\hat{D}_p := \bigcup_{m \in \mathbb{Z}_M^2} \hat{D} + \nu_m$

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

Next we rigorously show how to obtain the nearby solutions with the behavior as in the above. To this end, let $N_\sharp$ be defined by (20) and $N_{q,\sharp} := I_q^* N_\sharp I_q$, then for given $\phi$ and $a$ in $\ell^2(\mathbb{Z}^2)$ we define the functionals

$$J_{\alpha}(\phi, a) := \alpha \langle N_{\sharp}a, a \rangle + \|N_{\alpha} - \phi\|^2, \quad (49)$$

$$J_{\alpha,q}(\phi, a) := \alpha \langle N_{q,\sharp}a, a \rangle + \|N_q a - \phi\|^2,$$

here $\langle \cdot, \cdot \rangle$ denotes the inner product in $\ell^2(\mathbb{Z}^2)$.

Let $a^{\alpha,z}$, $a_q^{\alpha,z}$ and $\tilde{a}_q^{\alpha,z}$ in $\ell(\mathbb{Z}^2)$ verify (i.e. are minimizing sequences)

$$\begin{align}
J_{\alpha}(\hat{\Phi}(\cdot; z), a^{\alpha,z}) &\le \inf_{a \in \ell^2(\mathbb{Z}^2)} J_{\alpha}(\hat{\Phi}(\cdot; z), a) + c(\alpha) \nonumber \\
J_{\alpha}(\hat{\Phi}_q(\cdot; z), a_q^{\alpha,z}) &\le \inf_{a \in \ell^2(\mathbb{Z}^2)} J_{\alpha}(\hat{\Phi}_q(\cdot; z), a) + c(\alpha) \tag{50} \\
J_{\alpha,q}(I_q^* \hat{\Phi}_q(\cdot; z), a_q^{\alpha,z}) &\le \inf_{a \in \ell^2(\mathbb{Z}^2)} J_{\alpha,q}(I_q^* \hat{\Phi}_q(\cdot; z), a) + c(\alpha) \nonumber
\end{align}$$

with $\frac{c(\alpha)}{\alpha} \to 0$ as $\alpha \to 0$. The standard analysis of the generalized linear sampling method (see e.g. [7, Section 2.2]) making use of the factorization of $N_\sharp$ in Theorem 3 along with all the properties of the involved operators developed in Section 2 can be used to prove the lemma below.

glsm

**Lemma 1.** Under assumptions 1, 2, 3 the following results hold true:

(i) $z \in D$ if and only if $\lim_{\alpha \to 0} \langle N_\sharp a^{\alpha,z}, a^{\alpha,z} \rangle < \infty$. Moreover, if $z \in D$ then $H a^{\alpha,z} \to v_z$ in $L^2(D)$ where $(u_z, v_z)$ is the solution of problem (15) with $\varphi = \Phi(1; z)$ and $\psi = \partial\Phi(1; z)/\partial\nu$ on $\partial D$.

(ii) $z \in D_p \setminus \omega_p^{\text{mis}}$ if and only if $\lim_{\alpha \to 0} \langle N_\sharp a_q^{\alpha,z}, a_q^{\alpha,z} \rangle < \infty$. Moreover, if $z \in D_p \setminus \omega_p^{\text{mis}}$ then $H a_q^{\alpha,z} \to v_z$ in $L^2(D)$ where $(u_z, v_z)$ is the solution of problem (15) with $\varphi = \Phi_q(\cdot - z)$ and $\psi = \partial\Phi_q(\cdot - z)/\partial\nu$ on $\partial D$.

(iii) $z \in \hat{D}_p$ if and only if $\lim_{\alpha \to 0} \langle N_{q,\sharp} \tilde{a}_q^{\alpha,z}, \tilde{a}_q^{\alpha,z} \rangle < \infty$. Moreover, if $z \in \hat{D}_p$ then $H_q \tilde{a}_q^{\alpha,z} \to h_z$ in $L^2(D)$ where $h_z$ is defined by

$$\begin{align}
h_z &= \begin{cases} -\Phi_q(\cdot - z) & \text{in } \Lambda_p \\ v_z & \text{in } \mathcal{O}_p^c \end{cases} & \text{if } z \in \mathcal{O}_p^c \nonumber \\
h_z &= \begin{cases} \hat{v}_z & \text{in } \Lambda_p \\ -\Phi_q(\cdot - z) & \text{in } \mathcal{O}_p^c \end{cases} & \text{if } z \in \Lambda_p \tag{51} \end{align}$$

def:fz
where $(u_z, v_z)$ is the solution of problem (15) with $\varphi = \Phi_q(\cdot - z)$ and $\psi = \partial\Phi_q(\cdot - z)/\partial\nu$ on $\partial D$ and $(\hat{u}_z, \hat{v}_z)$ is $\alpha_q$-quasi-periodic extension of the solution $(u, f)$ of the interior transmission problem in Definition 1 with $\varphi = \Phi_q(\cdot - z)$ and $\psi = \partial\Phi_q(\cdot - z)/\partial\nu$ on $\partial\Lambda$.

Here $\mathcal{H}: \ell^2(\mathbb{Z}^2) \to L^2(D)$ is defined by (11) and $\mathcal{H}_q: \ell^2(\mathbb{Z}^2) \to L^2(D)$ is defined by (27).

*Proof.* For the proof of the items (i) and (ii), we refer to Theorem 3.5 and Lemma 4.7 in [17]. The proof of items (iii) is a direct application of Theorem A.4 in [17] in combination with Theorem 2. □

We then consider the following imaging functional, that characterizes the defects and the defective components of the periodic background,

$$ \mathcal{I}_{\alpha}(z) = \left( \langle N_{\sharp} \mathbf{a}^{\alpha,z}, \mathbf{a}^{\alpha,z} \rangle \left( 1 + \frac{\langle N_{\sharp} \mathbf{a}^{\alpha,z}, \mathbf{a}^{\alpha,z} \rangle}{D(\mathbf{a}_q^{\alpha,z}, \tilde{\mathbf{a}}_q^{\alpha,z})} \right) \right)^{-1} \quad (52) \quad \text{IndFuncDiff} $$

where for $\mathbf{a}$ and $\mathbf{b}$ in $\ell^2(\mathbb{Z}^2)$,

$$ D(\mathbf{a}, \mathbf{b}) := \langle N_{\sharp}(\mathbf{a} - I_q \mathbf{b}), (\mathbf{a} - I_q \mathbf{b}) \rangle. $$

**Theorem 3.** Under Assumptions 1, 2 and 3, we have that

fImagFunc

$z \in D \setminus \mathcal{O}_p^c$ if and only if $\lim_{\alpha \to 0} \mathcal{I}_\alpha(z) > 0$.

The domain $D \setminus \mathcal{O}_p^c$

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

*Proof.* This theorem is an improvement of Theorem 5.2 in [8] where only the configuration where $D \setminus \mathcal{O}_p^c = \Lambda_p$ had been investigated. The same arguments as in the proof of Theorem 5.2 in [8] show that $\lim_{\alpha \to 0} \mathcal{I}_\alpha(z) = 0$ for all $z \notin D$ and $z \in \mathcal{O}_p^c$. Therefore we only need to show that

$$ \lim_{\alpha \to 0} \mathcal{I}_\alpha(z) > 0 \quad \text{for all} \quad z \in D \setminus \mathcal{O}_p^c $$

(note that $D \setminus \mathcal{O}_p^c \subset \Lambda_p$ - see the figure above and Figure 3 for an illustration). By Lemma 1, the factorization of $N_{q,\sharp}$ and the identity $N_{q,\sharp} = I_q^* N_{\sharp} I_q$ we have

$$ (N_{\sharp} I_q \tilde{\mathbf{a}}_q^{\alpha,z}, I_q \tilde{\mathbf{a}}_q^{\alpha,z}) \to (T_{\sharp} h_z, h_z) < +\infty, $$

where $h_z \in L^2(\hat{D}_p)$ is defined by (52). Let us split domain $D \setminus \mathcal{O}_p^c$ into $(D \setminus \mathcal{O}_p^c) \cap D_p$ and $(D \setminus \mathcal{O}_p^c) \setminus D_p$ and treat each domain separately.

Case 1: We assume that $z \in (D \setminus \mathcal{O}_p^c) \setminus D_p$ (the defects that do not intersect the periodic domain). By Lemma 1(ii) $(N_{\sharp} \mathbf{a}_q^{\alpha,z}, \mathbf{a}_q^{\alpha,z}) \to +\infty$ as $\alpha \to 0$. This implies,

$$ D(\mathbf{a}_q^{\alpha,z}, \tilde{\mathbf{a}}_q^{\alpha,z}) \geq |(N_{\sharp} \mathbf{a}_q^{\alpha,z}, \mathbf{a}_q^{\alpha,z}) - (N_{\sharp} I_q \tilde{\mathbf{a}}_q^{\alpha,z}, I_q \tilde{\mathbf{a}}_q^{\alpha,z})| \to +\infty. $$
We then conclude that $\lim_{\alpha \to 0} \mathcal{I}_\alpha(\boldsymbol{z}) > 0$. Case 2: We assume that $\boldsymbol{z} \in (D \setminus \mathcal{O}_p^c) \cap D_p$. The case $(D \setminus \mathcal{O}_p^c) \cap D_p = \mathcal{O}_p$ is treated in [8]. The case $(D \setminus \mathcal{O}_p^c) \cap D_p \neq \mathcal{O}_p$ corresponds to the case where $\omega^{\text{mis}} \neq \emptyset$. In this domain $n = 1$ and $n_p \neq 1$. If $\boldsymbol{z}$ is not in $\omega_p^{\text{mis}}$ (i.e. one of the periodic copies of $\omega^{\text{mis}}$) then the same arguments as in [8] apply and show that $D(\boldsymbol{a}_q^{\alpha,z}, \tilde{\boldsymbol{a}}_q^{\alpha,z})$ remains bounded from 0 as $\alpha \to 0$. This is due to the fact that $h_z \neq v_z$ where $v_z$ and $h_z$ are respectively defined in Lemma 1 (ii) and (iii). This implies $\lim_{\alpha \to 0} \mathcal{I}_\alpha(\boldsymbol{z}) > 0$. The case where $\boldsymbol{z} \in \omega_p^{\text{mis}} \setminus \omega^{\text{mis}}$ can treated similarly to Case 1. By Lemma 1 (ii) and (iii), $\langle N_\sharp \boldsymbol{a}_q^{\alpha,z}, \boldsymbol{a}_q^{\alpha,z} \rangle \to \infty$ while $\langle N_\sharp \boldsymbol{a}_q^{\alpha,z}, \boldsymbol{a}_q^{\alpha,z} \rangle$ remains bounded as $\alpha \to 0$. Consequently $D(\boldsymbol{a}_q^{\alpha,z}, \tilde{\boldsymbol{a}}_q^{\alpha,z}) \to \infty$ as $\alpha \to 0$. Lemma 1 (i) indicates that $\langle N_\sharp \boldsymbol{a}^{\alpha,z}, \boldsymbol{a}^{\alpha,z} \rangle$ remains bounded as $\alpha \to 0$. The last two items show that $\lim_{\alpha \to 0} \mathcal{I}_\alpha(\boldsymbol{z}) > 0$, which ends the proof of the theorem. $\square$

## 3.1 Numerical Studies

This part will discuss some numerical examples using simulated data. The probed domain for the search is $\Omega_M$ that contains $M_1 \times M_2$ periods containing the defective cell. The data are generated by the set of incident plane waves described in (4) with the index $\boldsymbol{j} = (j_1, j_2)$ forms a set

$$ \mathbb{Z}_{inc}^2 := \{\boldsymbol{j} = (j_1, j_2) | \boldsymbol{j} = \boldsymbol{q} + M\boldsymbol{l} \text{ with } \boldsymbol{l} \in [-N_{min}, N_{max}]\}. $$

Here, $\boldsymbol{q} \in \mathbb{Z}_M^2$ is fixed, $N_{min}, N_{max} \in \mathbb{Z}_+^2$ are given. The scattered wave associated with each individual incident wave is solved numerically by implementing the Floquet-Bloch transform and volume integral method (see [18]). The collection of Rayleigh coefficients at index $\boldsymbol{l} \in \mathbb{Z}_{inc}^2$ of the scattered wave associated with each individual incident wave form the data for the inverse problem. Thus, if $N_{inc}$ denotes the number of incident waves, then the measurements $\{\hat{u}^s(\boldsymbol{l}; \boldsymbol{j})\}_{\boldsymbol{l},j \in \mathbb{Z}_{inc}^2}$ form a $N_{inc} \times N_{inc}$ matrix

$$ N := [\hat{u}^s(\boldsymbol{l}; \boldsymbol{j})]_{\boldsymbol{l},j \in \mathbb{Z}_{inc}^2}. \qquad (53) $$

Note that matrix $N$ in (54) is the matrix of measured data and is also known as the near-field operator. In practice, the measured data are always perturbed by some noise. To allow the presence of noise in data with the noisy level is denoted by $\delta > 0$ such that

$$ \|N^\delta - N\| \leq \delta \|N\|, $$

we introduce the noisy data, represented by noisy near-field operator $N^\delta$ by adding some random noise, that is

$$ N^{\delta}(\boldsymbol{j}, \boldsymbol{l}) := N(\boldsymbol{j}, \boldsymbol{l})(1 + \delta A(\boldsymbol{j}, \boldsymbol{l})), \quad \forall (\boldsymbol{j}, \boldsymbol{l}) \in \mathbb{Z}_{inc}^2 \times \mathbb{Z}_{inc}^2 \qquad (54) $$

where $A = [A(\boldsymbol{j}, \boldsymbol{l})]_{N_{inc} \times N_{inc}}$ is a matrix of uniform complex random variables with real and imaginary parts in $[-1, 1]^2$. Associated with the noisy data, one defines functionals $J_\alpha^\delta$ and $J_{\alpha,q}^\delta$ as

$$ J_{\alpha}^{\delta}(\boldsymbol{\phi}, \boldsymbol{a}) := \alpha \left( \langle N_{\sharp}^{\delta} \boldsymbol{a}, \boldsymbol{a} \rangle + \delta \|N_{\sharp}^{\delta}\| \|\boldsymbol{a}\|^2 \right) + \|N^{\delta} \boldsymbol{a} - \boldsymbol{\phi}\|^2, \qquad (55) $$

$$ J_{\alpha,q}^{\delta}(\boldsymbol{\phi}, \boldsymbol{a}) := \alpha \left( \langle N_{\sharp}^{\delta} I_q \boldsymbol{a}, I_q \boldsymbol{a} \rangle + \delta \|N_{\sharp}^{\delta}\| \|\boldsymbol{a}\|^2 \right) + \|N_{q}^{\delta} \boldsymbol{a} - \boldsymbol{\phi}\|^2. $$
We then consider $a_{\delta}^{\alpha,z}$, $a_{q,\delta}^{\alpha,z}$ and $\tilde{a}_{q,\delta}^{\alpha,z}$ in $\ell^2(\mathbb{Z}^2)$ as the minimizing sequences of, respectively,

$$J_{\alpha}^{\delta}(\hat{\Phi}(\cdot; z), \mathbf{a}), \quad J_{\alpha}^{\delta}(\hat{\Phi}_{q}(\cdot; z), \mathbf{a}), \quad \text{and} \quad J_{\alpha,q}^{\delta}(\hat{\Phi}_{q}(\cdot; z), \mathbf{a}).$$

The noisy indicator function takes the form

$$\mathcal{I}_{\alpha}^{\delta}(z) = \left( \mathcal{G}^{\delta}(a_{\delta}^{\alpha,z}) \left( 1 + \frac{\mathcal{G}^{\delta}(a_{\delta}^{\alpha,z})}{D^{\delta}(a_{q,\delta}^{\alpha,z}, a_{q,\delta}^{\tilde{\alpha},z})} \right) \right)^{-1} \quad (56)$$

where for $\mathbf{a}$ and $\mathbf{b}$ in $\ell^2(\mathbb{Z}^2)$,

$$D^{\delta}(\mathbf{a}, \mathbf{b}) := \langle N_{\#}^{\delta}(\mathbf{a} - I_q \mathbf{b}), (\mathbf{a} - I_q \mathbf{b}) \rangle$$

and

$$\mathcal{G}^{\delta}(\mathbf{a}) := \langle N_{\#}^{\delta}\mathbf{a}, \mathbf{a} \rangle + \delta \|N_{\#}^{\delta}\| \|\mathbf{a}\|^2.$$

To determine the minimizing sequences above we approximate the inner product $\langle N_{\#}^{\delta}\mathbf{a}, \mathbf{a} \rangle$, which is positive defined, by the norm $\|N_{\#}^{\delta}\|\|\mathbf{a}\|^2$ and then seek the minimiser of the functional where $\langle N_{\#}^{\delta}\mathbf{a}, \mathbf{a} \rangle$ is replaced by $\|N_{\#}^{\delta}\|\|\mathbf{a}\|^2$. For instance, to find $\mathbf{a}_{\delta}^{\alpha,z}$, we find the minimiser of the functional

$$J_{\alpha,mod}^{\delta}(\phi, \mathbf{a}) := \alpha(1 + \delta)\|N_{\#}^{\delta}\|\|\mathbf{a}\|^2 + \|N^{\delta}\mathbf{a} - \phi\|^2, \quad (57)$$

which can be viewed as a Tikhonov regularization with the regularization parameter $\alpha_{app} := \alpha(1+\delta)\|N_{\#}^{\delta}\|$. We then employ the singular value decomposition technique with the Morozov discrepancy parameter choice as described in [13] to solve $\mathbf{a}_{\delta}^{\alpha,z}$ and parameter $\alpha_{app}$. Finally, parameter $\alpha$ is chosen as $\alpha = \frac{\alpha_{app}}{(1+\delta)\|N_{\#}^{\delta}\|}$. The processes are similar for solving $\mathbf{a}_{q,\delta}^{\alpha,z}$ and $\tilde{\mathbf{a}}_{q,\delta}^{\alpha,z}$.

In the following examples, we normalize the the imaging functional $\mathcal{I}_{\alpha}^{\delta}(z)$ and consider all the domain $z$ in the domain $\{z : \mathcal{I}_{\alpha}^{\delta}(z) \ge \kappa\}$ for 3D isosurface plot. In most examples, $\kappa = 0.45$. Furthermore, we mainly focus on taking into account different type of defects rather than varying the refractive index or other parameters. Thus, for all examples, values of the following parameters are fixed and set to be

$$n_p = 2, n = 4, k = \pi, L_x = L_y = 2\pi, M = (3,2), q = (1,1), \delta = 5\% \quad (58)$$

and the width of the layer $h = 1.5\lambda$, where $\lambda := 2\pi/k$ denotes the wavelength. The information of the periodic structure is provided only in the cell $\Omega_0$ unless there is other remarks.

### Example 1:

In the first example, we consider the periodic domain made of two cubes in each period. In particular, in the period $\Omega_0$, such two cubes are of the edge sizes $0.6\lambda$ and $0.5\lambda$ and centered at $(0.6\lambda, -0.6\lambda, -0.3\lambda)$ and $(-0.5\lambda, 0, 0)$ respectively. The local perturbation is a ball of the radius $r = 0.3\lambda$ and non-overlapping periodic components (see Figure 4(a) - 4(b) for 3D plot and a slice through the center of the ball and perpendicular with $Oy$ respectively). The reconstruction using the indicator function $\mathcal{I}_{\alpha}^{\delta}(z)$ is displayed in Figure 4(c) - Figure 4(d).
<figure><img src="image_11.png" /><figcaption>Figure 4: Example 1: The defect (a ball) does not overlap any components of the periodic background. (a) - the true background in 3D; (b) - the trace of the true background cut by a plane passing through the center of the ball and perpendicular with *Oy*. (c) - the reconstruction of the local perturbation using the indicator function $I_\alpha(z)$ plotted in 3D; (d) - the trace of the reconstruction cut by a plane passing through the center of the ball and perpendicular with *Oy*.</figcaption></figure>

Example 2:

The periodic background in this example is as same as the one in Example 1. The perturbation $\omega$ is a ball that overlaps the smaller cube (see Figures 5(a) for 3D plot and 5(b) for the trace of the true background cut by a plane passing through the center of the ball and perpendicular with $Oy$). Here $n = 4$ inside the overlapped domain and $n = 3$ in the rest of the ball (see Figure 5(b)). The reconstruction using indicator function $I_\alpha^\delta(z)$ is displayed in Figures 5(c) - 5(d).
<figure><img src="image_12.png" /><figcaption>Figure 5: Example 2: The defect (a ball) overlaps a component of the periodic background. (a) - the true background in 3D; (b) - the trace of the true background cut by a plane passing through the center of the ball and perpendicular with *Oy*. (c) - the reconstruction of the local perturbation using the indicator function *I*<sub>α</sub>(*z*) plotted in 3D; (d) - the trace of the reconstruction cut by a plane passing through the center of the ball and perpendicular with *Oy*.</figcaption></figure>

intersection

Theorem 3 shows that values of indicator function I_α^δ(z) is positive for all point z belongs to the small component in all periods as well as the local perturbation. However, we do not have quantitative conclusion about value of I_α(z) at each point in the probed domain. Figure 5 - right can be interpreted that the value of I(z) at such point z belongs to the cube that overlaps the defect, but not in the overlapped domain, is much smaller than ones belongs to the overlapped domain. Moreover, comparing the reconstruction of defects illustrated in Figure 4 with Figure 5, we can point out the case of overlapping and non-overlapping thanks to the reconstruction of (partial) components of periodic background.

**Example 3:**

In this example, we consider the case that each period of the periodic background includes
two cubes and one ball (see Figure 6(a)) and the local perturbation is included in the ball
of period Ω₀ (see Figure 6(b)). Even the theory justifies that the indicator function should
provide a periodic domain that includes only the ball component, the numerical result allows
identifying the period containing the defect (See Figures 6(c) - 6(d)).
<figure><img src="image_13.png" /><figcaption>Figure 6: Example 3: The defect is included in the ball component of the periodic domain. (a) - the true background in 3D; (b) - the trace of the true background cut by a plane passing through the center of the ball and perpendicular with $Oy$. (c) - the reconstruction of the local perturbation using the indicator function $I_\alpha(z)$ plotted in 3D; (d) - the trace of the reconstruction cut by a plane passing through the center of the ball and perpendicular with $Oy$.</figcaption></figure>

## Example 4.

We will end with showing three more numerical tests for the case that there is a complete component missing (see Figure 7), a partial component missing (see Figure 8, and a the partial component missing but the rest of it also overlaps with another defect (see Figure 9 and we call this case by “missing and overlapping”).

For the case of complete missing component, the numerical results justify exactly what the theory stated.

For the case of partial missing component only, the reconstruction using $I_\alpha(z)$, showed in see Figure 8, seems to be like the case of complete missing component. This can be explained, because the value of $I_\alpha(z)$ in the rest of the partial missing component and its periodic copies in other periods may be much smaller than the value of $I_\alpha(z)$ in the periodic copies of the missing part in other periods.

We can summary from both cases that the method allow to identify the cases of missing component (in the condition that the missing component does not overlap with other defects – see the last case for more explanation), but cannot determine whether the missing part is a whole or a partial component.
<figure><img src="image_14.png" /><figcaption>Figure 7: Example 4: The complete missing component. (a) - the true background in 3D. (b) - the reconstruction of the local perturbation using the indicator function $I_\alpha(z)$ plotted in 3D.</figcaption></figure>

missingcomp

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

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

fig:DLSMmissi

fig:BGMisP2D

fig:DLSMMisP2D

<caption>Figure 8: Example 4: The defect (a ball) overlaps a component of the periodic background. (a) - the true background in 3D; (b) - the trace of the true background cut by a plane passing through the center of the ball and perpendicular with $Oy$. (c) - the reconstruction of the local perturbation using the indicator function $I_\alpha(z)$ plotted in 3D; (d) - the trace of the reconstruction cut by a plane passing through the center of the ball and perpendicular with $Oy$.</caption>

Iso-MisP

For the case of “missing and overlapping”, it shows in the period $\Omega_0$ the perturbation that overlaps the partial missing component. In other periods, it shows the partial component
corresponding to the missing part in the period $Ω_0$. This is suitable the theoretical part, but the numerical result looks similar to the case of overlapping above. That is to say, in such cases, the method provides that the perturbation overlaps a component of the periodic background, but hard to determine whether or not there is a part of the overlapping component missing.

<figure><img src="image_17.png" /><figcaption>Figure 9: Example 4: The case of “missing and overlapping”. (a) - the true background in 3D; (b) - the trace of the true background cut by a plane passing through the center of the ball and perpendicular with $Oy$. (c) - the reconstruction of the local perturbation using the indicator function $I_α(z)$ plotted in 3D; (d) - the trace of the reconstruction cut by a plane passing through the center of the ball and perpendicular with $Oy$.</figcaption></figure>

References

Bilo-thesis

[1] T. Arens. *Scattering by biperiodic layered media: The integral equation approach*. Habilitation theses, Karlsruhe Institute of Technology, January 2010.

Aud

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

BaoCow01

[3] G. Bao, L. Cowsar, and W. Masters. *Mathematical Modeling in Optical Science*. SIAM, Philadelphia, 2001.

BaoDob1995

[4] G. Bao, D. C. Dobson, and J. A. Cox. Mathematical studies in rigorous grating theory. *J. Opt. Soc. Amer. A*, 12:1029–1042, 1995.
BonSta1994

BFliss

CBMS

Thi-Phong1

Simon-peter

ChaZha1999

ChaEls2010

coat2011

a-Potthast

ElsSch1998

FliJol2016

comeon

Thi-Phong3

[5] A.-S. Bonnet-Bendhia and F. Starling. Guided waves by electromagnetic gratings and non-uniqueness examples for the diffraction problem. *Mathematical Methods in the Applied Sciences*, 17:305–338, 1994.

[6] 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.

[7] D. Cakoni, F. and Colton and H. Haddar. *Inverse Scattering Theory and Transmission Eigenvalues, second edition*, volume 98 of CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2022.

[8] F. Cakoni, H. Haddar, and T-P Nguyen. New interior transmission problem applied to a single Floquet-Bloch mode imaging of local perturbations in periodic media. *Inverse Problems*, 35(1):015009, 31, 2019.

[9] S. N. Chandler-Wilde and P. Monk. Existence, uniqueness, and variational methods for scattering by unbounded rough surfaces. *SIAM J. Math. Anal.*, 37(2):598–618, 2005.

[10] S. N. Chandler-Wilde and B. Zhang. Scattering of electromagnetic waves by rough interfaces and inhomogeneous layers. *SIAM J Math. Anal.*, 30:559–583, 1999.

[11] S.N. Chandler-Wilde and J. Elschner. Variational approach in weighted Sobolev spaces to scattering by unbounded rough surfaces. *SIAM J Math. Anal.*, 42:2554–2580, 2010.

[12] J. Coatléven. *Mathematical and numerical analysis of wave problems in periodic media*. Theses, Ecole Polytechnique X, November 2011.

[13] D. Colton, M. Piana, and R. Potthast. A simple method using morozov's discrepancy principle for solving inverse scattering problems. *Inverse Problems*, 13(6):1477–1493, dec 1997.

[14] J. Elschner and G. Schmidt. Diffraction in periodic structures and optimal design of binary gratings. part I: Direct problems and gradient formulas. *Math. Meth. Appl. Sci.*, 21:1297–1342, 1998.

[15] S. Fliss and P. Joly. Solutions of the time-harmonic wave equation in periodic waveguides: Asymptotic behaviour and radiation condition. *Arch. Rational Mech. Anal.*, 219:349–386, 2016.

[16] H. Haddar and A. Konschin. Factorization method for imaging a local perturbation in inhomogeneous periodic layers from far field measurements. *Inverse Probl. Imaging*, 14(1):133–152, 2020.

[17] H. Haddar and T-P. Nguyen. Sampling methods for reconstructing the geometry of a local perturbation in unknown periodic layers. *Comput. Math. Appl.*, 74(11):2831–2855, 2017.
[18] H. Haddar and T-P. Nguyen. A volume integral method for solving scattering problems from locally perturbed infinite periodic layers. *Appl. Anal.*, 96(1):130–158, 2017.

[19] A. Kirsch. Diffraction by periodic structures. In L. Pävarinta and E. Somersalo, editors, *Proc. Lapland Conf. on Inverse Problems*, pages 87–102. Springer, 1993.

[20] A. Kirsch and A. Lechleiter. The limiting absorption principle and a radiation condition for the scattering by a periodic layer. *SIAM J. Math. Anal.*, 50(3):2536–2565, 2018.

[21] A. Lechleiter and D-L. Nguyen. Volume integral equations for scattering from anisotropic diffraction gratings. *Math. Methods Appl. Sci*, 36(3):262–274, 2013.

[22] A. Lechleiter and S. Ritterbusch. A variational method for wave scattering from penetrable rough layers. *IMA J. Appl. Math.*, 75(3):366–391, 2010.

[23] A. Lechleiter and R. Zhang. Reconstruction of local perturbations in periodic surfaces. *Inverse Problems*, 34(3):035006, 2018.

[24] D-L Nguyen. *Spectral Methods for Direct and Inverse Scattering from Periodic Structures*. PhD thesis, Ecole Polytechnique X, 2012.

[25] T-P Nguyen. *Direct and inverse solvers for scattering problems from locally perturbed infinite periodic layers*. Theses, Ecole Polytechnique X, January 2017.

[26] T-P. Nguyen. Differential imaging of local perturbations in anisotropic periodic media. *Inverse Problems*, 36(3):034004, 30, 2020.

[27] R. Petit, editor. *Electromagnetic theory of gratings*. Springer, 1980.