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

We discuss a novel approach for imaging local faults inside an infinite bi-periodic layered medium in 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.


Formulation of the Problem
We consider nondestructive testing of an infinite bi-periodic penetrable layer in 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 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 its indicator function computable from scattering data, leading to computationally efficient non-iterative imaging method.The connection between scattering data and the support of local perturbations is made through a rigorous mathematical analysis of the scattering problem by the bi-periodic layer with and without local perturbations.To be more specific and set the notations, let x := (x 1 , x 2 , x 3 ) ∈ R 3 , and assume that the scattering media is a penetrable infinite layer periodic in the x 1 and x 2 variables with period L 1 and L 2 , respectively.Given a 2-dimensional vector L = (L 1 , L 2 ), we call a function w defined in R 3 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 ∈ L ∞ (R 3 ) is L-periodic such that Re (n p ) > 0 and Im (n p ) ≥ 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) is governed by ∆u + k 2 n p u = 0 in R 3 , u is L-periodic (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 structures such as gratings, inhomogeneous layers or waveguides, is a research topic that has received enormous attention due to contemporary applications in material science.Among vast literature in the topic, we refer the reader to [1,3,4,5,12,14,15,19,20,21,24,25,27] (the list is not exclusive by any means).For 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, where an application of well-known techniques from analysis such as variational methods and analytic Fredholm theory [14], provides existence and uniqueness of the 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.Our main interest is in the case when local perturbations are present inside the bi-periodic layer.We denote by ω the support of perturbations, such that ω is a compact set with connected complement R 3 \ ω.We may assume without loss of generality (up to a possible rearrangement of the cell), that ω 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 ∈ L ∞ (R 3 ), and it satisfies Re (n) ≥ n 0 > 0, Im (n) ≥ 0. Note that ω 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 ω by using the measured scattered fields outside the layer due to appropriate incident fields.The challenging task however is to resolve ω without an explicit knowledge of n p (which in practice can have a complicated form) nor reconstructing it, but just using the fact that n p is L-periodic.This problem was first considered in [25], where the GLSM with differential measurements introduced in [2] was modified to the current periodic configuration.In this context, as oppose to [2], 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.The analysis in [25] was further developed in [8] and [26], where numerical examples were also presented only for the 2-dimensional case.However, the imaging method was justified under some restrictive assumptions on the location of local perturbations.In the current work, we remove these restrictions and complete the justification of the imaging method for a general setting of local perturbations.In particular, our analysis includes for the first time the case where components in some periodic cells are missing, or where the perturbation is entirely inside a component.In addition, this paper presents numerical examples for local perturbations of a bi-periodic layer in R 3 , which are much more challenging and closer to real applications.
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 ) ∈ N 2 sufficiently large and then extend it as M L-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 M L-periodic extension is needed only for the analysis of derivation of the indicator function of the set ω, and 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 assumption in the analysis, but all of them rely an explicit knowledge of the Green's function of the periodic layered background, which is not the case for our method.Although it is desirable to remove this technical step, one could see it as a 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 : where u i is the probing incident field and u s is the scattered field.Here M = (M 1 , M 2 ) ∈ 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 M L-periodicity, solving (2) is equivalent to solving it in the period , ℓ = 1, 2}, with ⌊•⌋ denoting the floor function and Z the set of integers.Note that in (2) we still call n the M L-periodic extension of n| Θ and without loss of generality assume that the defective period is 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 where for each mode j = (j We remark that in terms of (17), 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: where the Rayleigh coefficients u s (ℓ) are given by Recall that the M − L , M + L is a rectangle on x 1 x 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 We shall use the notations together with the Rayleigh radiation condition (5).We remark that the solution w ∈ H 1 # (Θ h ) of ( 7) can be extended to a function in Θ satisfying ∆w + k 2 nw = k 2 (1 − n)f , using the Rayleigh expansion (5), and hence by M L-periodicity to a solution in the entire R 3 .Note that the scattering problem for L-periodic layer (1) is equivalent to (7) where w := u s , f := u i | h Θ and n := n p , whereas the scattering problem for M L-periodic layer (2) is equivalent to (7) where w := u s and f := u i | h Θ .Throughout the paper we make the following assumption: 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 ∈ L 2 (Θ h ).
For sufficient conditions that guarantee Assumption 1 we refer the reader to [26], [20], [25] and the references therein.If Φ(n p ; •) is the fundamental solution to and the Rayleigh radiation condition (5).(8) then the solution w of (7) has the representation as Finally, as it becomes clear latter in the paper, our inversion algorithm is well-suited to the case when a period of healthy L-periodic layer consists of several compactly supported inhomogeneities sitting in homogeneous structure (see Figure 1), which is the case in many applications.For sake of presentation we assume that the homogenous base structure of the L-period media has refractive index one.With this in mind, we introduced the following notations which will be used throughout the paper.To this end, let us first denote by ab the element wise multiplication of two generic vectors a = (a 1 , a 2 ) and b Recalling that Ω 0 is the period containing ω, we denote by and assume that the exterior of each of D p , D and ω is connected.Next we denote by Let ν m := (mL, 0) ∈ R 3 be the translation vector Ω 0 → Ω m (m-th period) for m ∈ Z 2 and denote by which are L-periodic copies of the respective aforementioned regions.Finally we denote by ω mis := ω \ D (possibly part of ) missing components of D p and ω mis p its L-periodic copies.Note that D p is periodic and We refer the reader to Figures 2-3 for an illustration of type of defects ω we consider in this paper and for an illustration of the notation introduced above.
Remark 1.In the special case when the defect ω consists only of missing components of D p in Ω 0 , that in when n = 1, we have D p = D p and Λ = O = ω mis ∩ Ω 0 .
Remark 2. All the results presented here can be readily 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 | ≤ h.

The domain
The domain D = Supp(n − 1) The domain ω = Supp(n p − n) The imaging method discussed in this paper falls into the recently developed qualitative approach to inverse scattering or otherwise referred to as non-iterative methods [7].The underlying 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 its standard presentation this approach requires an expression for the Green's function of unperturbed background.However, in many applications an accurate modeling of the background is difficult to compute, hence 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, were mathematically analyzed to arrive at an indicator function for the support of perturbation to the healthy based structure.This idea was latter adapted to periodic layers in [25] (see also [17]).The advantage of applying differential measurements to the inverse problem for periodic layer like the one considered here, is that instead of comparing the measured data against measurements corresponding to the healthy structure, one makes use its periodicity.More precisely, the measurement operator restricted to single Floquet-Bloch modes plays the role of data operator corresponding to the healthy background.As a result, the support ω of perturbations is obtained without knowing or reconstructing the constitutive material properties of the periodic background.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 interior transmission problems related to the problems considered here, are studied in details in [8,26].In this paper we provide a comprehensive presentation of this novel imaging method for locally perturbed bi-periodic layer in R 3 in acoustic scattering.We include in our discussion new configurations of defective areas, including missing components or part of components of the healthy cell.In addition, we present the first 3D numerical implementation of the imaging algorithm.The paper is organized as follows.In the next section, we develop the analytical tools of our inversion method, in particular the properties of the data (near field) operator corresponding the faulty periodic layer, and its restriction to single Floquet-Bloch modes.Based on this analysis, in Section 3 we design an imaging algorithm for the support ω of the defective region.In addition, we implement numerically this algorithm, and present a large collection of 3D numerical examples for a various type of defects.

The Inverse Problem
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 (that is 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 (that is the transmitters are above the layer whereas the receivers are under the layer).The inverse problem reads: from a knowledge of Rayleigh sequences u s + (ℓ; j) due to all incident waves u i,+ (x; j) for j ∈ Z 2 (or Rayleigh sequences u s − (ℓ; j) due to all incident waves u i,− (x; j) for j ∈ Z 2 ) determine ω = Supp(n − n p ) without knowing n and 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, from now on we consider the Rayleigh sequences u s + (ℓ; j) ℓ,j∈Z 2 ×Z 2 due to incident waves u i,+ (x; j), and to simplify the notation we let u s (ℓ; j) := u s + (ℓ; j) and u i (x; j) := u i,+ (x; j) for all ℓ ∈ Z 2 and j ∈ Z 2 .This scattering data defines the so-called near field (or data) operator N : (recall that ℓ 2 (Z 2 ) is the Hilbert space of square summable sequences in Z 2 ).This operator is one of the main objects of our imaging method.We show that its properties are connected to the reconstruction D (i.e. the support on n − 1).To this end, we define the Herglotz operator H : ℓ 2 (Z 2 ) → L 2 (D) by Its adjoint ) is given by Let us denote by H inc (D) the closure of the range of H in L 2 (D).We then consider the compact operator G : where { w(ℓ)} ℓ∈Z 2 is the Rayleigh sequence of the solution w ∈ H 1 # (Θ h ) of (7).Then by linearity N : The properties of G and H are crucial to our inversion method.To state them, we must recall the standard interior transmission problem: for given (φ, ψ) ∈ H 3/2 (∂D) × H 1/2 (∂D) where ν denotes the outward normal on ∂D.
The wave number k is called a (standard) transmission eigenvalue if the corresponding homogeneous problem, i.e. (15), with φ = 0 and ψ = 0, has non-trivial solutions.Up-todate 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 ∂D∩∂Ω 0 = ∅ where Ω 0 is given by (3).If the boundary of D intersects the vertical sides of the boundary Ω 0 , then the previous interior transmission problem should be augmented with periodicity conditions on ∂D ∩ ∂Θ (intersection of D with horizontal sides of the boundary Ω 0 causes no problems).For sake of presentation, since this case does not affect the assumptions on the solvability of the interior transmission problem (in H 2 (D) with periodic conditions on ∂D ∩ ∂Θ), we omitting it in our discussion.In the sequel we make the following assumption.
Assumption 2. ∂D ∩ ∂Ω 0 = ∅, and the refractive index n and the wave number k > 0 are such that (15) has a unique solution.
In particular, if Re (n − 1) > 0 or −1 < Re (n − 1) < 0 uniformly in a neighborhood of ∂D inside D the interior transmission problem (15) satisfies the Fredholm alternative, and the set of real standard transmission eigenvalues is discrete (possibly empty) without interior accumulation points.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 H and G in relation to the operator N which is available from the measurements.
Finally we introduce the solution operator T : with w being the solution of (7).By construction we have the following relation Gf = H * Tf which leads to the following symmetric factorization of N 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 where Re (N) := 1 2 (N + N * ), Im (N) := 1 2i (N − N * ), and N * :  [17]).Under Assumptions 1 and 2 we have that where 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. In particular that Picard series for the range of N ♯ converges if and only if the sampling point z is in D. However, for complicated periodic structures and relatively small defects or defects with peculiar location such as ω ⊂ D p , reconstructing D is unsatisfactory.Our aim is to design an indicator function of the set ω without knowing or recovering D p in the same spirit as the criterion in Proposition 3. To achieve this, the idea of using the data (near field) operator corresponding to one Floquet-Bloch mode and perform the same type of analysis to it as for N was first introduced in [25].This operator plays the role of the data operator corresponding to unperturbed 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 ∈ Z 2 M we consider only q + ℓM Rayleigh coefficients of scattered waves u s (q+ℓM ; q+jM ) generated by incident waves u i (•; q+jM ).Thus, single Floquet-Bloch mode data operator: The operators N q and N are related through the projection operator given by with its adjoint I * q : ℓ 2 (Z 2 ) → ℓ 2 (Z 2 ) given by (I * q b)(j) = b(q + jM ).
Hence we have (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 α q −quasi-periodic fields with period L, where α q = 2π M L q, since the sequence N q a corresponds to the Fourier coefficients of α q −quasi-periodic component of the scattered field.Recall that a function u is called α q −quasi-periodic fields with period L = (L 1 , L 2 ) if u(x + jL, x 3 ) = e iαq•(jL) u(x, x 3 ), for all j ∈ Z 2 .
To understand the operator N q we need the α q −quasi-periodic fundamental solution Φ q (•) that satisfies ∆Φ q (• − z) for z ∈ R d .Its Rayleigh coefficients Φ q (•; z) of Φ q (• − z) are given by Φ q (j; z) Similarly to the Herglotz operator H, the single Floquet-Bloch mode Herglotz operator H q : ℓ 2 (Z 2 ) → L 2 (D) is defined by where we note that H q a restricted to D p is also α q −quasi-periodic function with period L (see Remark 3 if missing (parts of) components are present).Note that R(H q ) is characterized in Proposition 4. Then the factorization (19) along with (24) immediately implies where T is defined by (18).The role of G given by ( 13) with respect to N q is now played by G q : R(H q ) → ℓ 2 (Z 2 ) by G q = H * q T| R(Hq) .
Observing that φ(j; is a Fourier basic of M L periodic function in L 2 (Ω M ), we have that any w ∈ L 2 (Ω M ) which is M L periodic, has the expansion w(x) = j∈Z w(j, x 3 )φ(j; x), where w(j, x 3 ) : (here the line over φ denotes the conjugation whereas x = (x 1 , x 2 )).Splitting j by module M component by component (i.e.splitting j ℓ by module M ℓ , ℓ = 1, 2) we can arrange the expansion of w as where φ(q + M ℓ; x) is α q −quasi-periodic with period L. We also have that where w q is α q −quasi-periodic with period L.Moreover, by the orthogonality of the Fourier basis {φ(j; •)} j∈Z , we have that w q (j) = 0 if j ̸ = q + M ℓ, ℓ ∈ Z 2 and w(q + M ℓ) = w q (q + M ℓ) where w q (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 w(j) corresponding to the indices j = q + M ℓ, ℓ ∈ Z 2 , where w is solution of (7).
The above discussion is helpful to proving the following properties for 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. The operator H q : ℓ 2 (Z 2 ) → L 2 (D) is injective.Furthermore the closure of its range is given by R(H q ) = H q inc (D) where We remind the reader that our discussion in the following heavily relies on Notation 1.Before we prove the preposition, we make some clarification in the remark below.
Remark 3. If ω contains missing (part of) components of D p in (34) we mean that v has an extension as α q −quasi-period in D p .More specifically this extension is defined as (35) Note also that, in the case D p ⊆ D (this is the case studied in [8]), i.e. there are no missing components, the definition in (34) becomes Proof of Proposition 4. Injectivity of H q follows from injectivity of the operators H and I q .Hence it suffices to show that H * q is injective in H q inc (D).The case when ω ⊂ D corresponds to Lemma 3.1 in [8].We sketch here this proof to confirm that it works also in the case when ω ̸ ⊂ D. In particular, ω is considered to be the entire component O or a part of the component O.To this end, let φ ∈ H q inc (D) and assume H * q (φ) = 0. Then we define where Φ q has the Rayleigh coefficients given by (26).By the definition of u and using (26) we see that the Rayleigh coefficients of u are given by u(j) = 0 for all j ̸ = q + M ℓ and u(q + M ℓ) = (H * (φ)) (q + M ℓ) = (H * q (φ))(ℓ) = 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 y ∈ D, ∆Φ q (•; y) + k 2 Φ q (•; y) = 0 in the complement of D p .This implies that Using a unique continuation argument we infer that u = 0 in Ω M \ D p .Therefore, u ∈ H 2 0 ( D p ) by the regularity of volume potentials.We remark that ω may include components that are part of 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 fo Lemma 3.1].For this reason, here we focus on the case where the defect occupies a region in O that is missing (or partially missing) in D. Specifically in this case, ω = D p \ D, and n = 1 in ω.Therefore D p = D p and D ∩ ω = ∅.We then obtain u ∈ H 2 0 (D p ), and it rewrite it as where φ is α q −quasi-periodic in D p .From the definition of u(x), we first observe that u(x) is α q −quasi-periodic.Indeed for all m ∈ Z 2 M .We now consider x ∈ D p ∩ Ω 0 .By the α q −quasi-periodicity of Φ q and φ with period L, we have and hence Next, we define φ(x) for all x ∈ D p ∩ Ω 0 by Evidently, extending φ as α q −quasi-periodic to D p we have Furthermore, linearity of the Laplace operator and the fact that φ(x) satisfies equation ∆φ + Since u ∈ H 2 0 (D p ) we then have , then the definition of φ implies φ = 0 in D p ∩ Ω = D, which ends the proof.Otherwise, when M is such that M 1 M 2 ≥ 2 we have φ = 0 in D p ∩ Ω 0 , and therefore by the quasi-periodicity φ = 0 in D p ⊃ D. This proves the injectivity of (H ± ) * on H q inc (D) and hence proves the Lemma.
In the case of missing component, we assume that f ∈ H q inc (D) then f (defined by the same way as ṽ in (35)) is α q −quasi-periodic in D p .Furthermore, since in this case, n = 1 in ω, solution w to equation (7) Using decomposition (32) for the solution w of ( 7) along with the facts that n p is periodic, f is α q −quasi-periodic and that n − n p is compactly supported in one period ω ⊂ Ω 0 , then (7) in terms of the coefficients w q in (32) for this w takes the form However, since f = f in the support of n − 1, (39) becomes Therefore, operator G q : R(H q ) = H q inc (D) → ℓ 2 (Z 2 ) can be equivalently defined by G q (f ) : where w q solves (40), w is solution to (7), and { w q (ℓ)} ℓ∈Z 2 is the Rayleigh sequence of w q .
We need to introduce one more interior transmission problem that is related to the characterization of the range of G q . where Φ(n p ; •) is the M L-periodic fundamental solution given by ( 8) and ν denotes the unit normal on ∂Λ outward to Λ.
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.
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 in addition when (n, D) is replaced by (n p , D p ). Then the operator G q : H q inc (D) → ℓ 2 (Z 2 ) is injective with dense range.
Proof.As in the previous proof, we only consider the case where the defect is constituted by components in D p that are missing in D, i.e. ω ⊆ O and ω ∩ D = ∅.The other cases can be treated as in [8,Theorem 3.2].Assume that f ∈ H q inc (D) such that G q (f ) = 0 and f ∈ L 2 (D p ) is an extension of f given in (35).Let w be solution of ( 7) with data f .We have that G q (f ) := I * q { w q (ℓ)} ℓ∈Z 2 , where and w q is solution to Note that ∆w q + k 2 w q = 0 in Ω M \ D p .Using a similar unique continuation argument as in the beginning of proof of Proposition 4, we deduce that In other words, w q = 0 outside D p .From now to the end of the proof, we consider only the period Ω 0 .Since Supp(n − n p ) ∩ Supp(1 − n) = ω ∩ D = ∅, and n = n p in D, (44) can be split into two equations where w ∈ H 1 loc (Θ) is solution to (7) and Observe that w satisfies ∆w + k 2 w = 0 in ω.Therefore, from (45 Assumption 2 with (n, D) replaced by (n p , D p ) implies that (47) has only the trivial solution and therefore w q = −w = 0 in ω.
We again see that, in the domain Assumption 2 with (n, D) replaced by (n p , D p ) again implies that (48) has only the trivial solution, and therefore This proves the lemma.
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 plays an important role in the design of our imaging method.
Theorem 2. Suppose that Assumptions 1, 2 and 3 hold.Then, I * q Φ q (•; z) ∈ R(G q ) if and only if z ∈ D p .
We now have all the ingredients to arrive at our imaging method in the next section.

A Differential Imaging Algorithm
We now apply theoretical results of Section 2 to design an algorithm that provides us with an indicator function of the support of the defect ω 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 in Section 2. Roughly speaking we can construct three appropriate sequences a α,z , a α,z q and ãα,z q (to become precise latter), as nearby solutions to as α → 0, where Φ(•; z) are the Rayleigh coefficients of Φ(1, z) (i.e.Φ(n p ; z) defined by (8) with with n p = 1) given by (17), and Φ q (•; z) are the Rayleigh coefficients of Φ q (• − z) given by (26).Then we show that these nearby solutions satisfy: The domain D = Supp(n − 1) where (u z , v z ) is the solution of problem (15) with φ = Φ q (• − z) and ψ = ∂Φ q (• − z)/∂ν on ∂D and ( u z , v z ) is α q −quasi-periodic extension of the solution (u, f ) of the interior transmission problem in Definition 1 with φ = Φ q (• − z) and ψ = ∂Φ q (• − z)/∂ν on ∂Λ.
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, where for a and b in ℓ Thus we reconstruct the support of defects that do not intersect healthy components together with periodic copies of the components that intersect the defective region as illustrated below.
The domain D \ O c p Proof.This theorem complements Theorem 5.2 in [8] where the configuration with D \ O c p = Λ p was investigated only.Thus the same arguments as in the proof of Theorem 5.2 in [8] show that lim α→0 I α (z) = 0 for z / ∈ D and z ∈ O c p .Therefore we only need to show that Note that D \ O c p ⊂ Λ p (see the figure above and Figure 3 for an illustration).By Lemma 1, factorization of N q,♯ and identity N q,♯ = I * q N ♯ I q we obtain (N ♯ I q ãα,z q , I q ãα,z p ) \ D p (parts of the defect that do not intersect healthy components of the periodic layer).By Lemma 1(ii) (N ♯ a α,z q , a α,z q ) → +∞ as α → 0. This implies, D(a α,z q , ãα,z q ) ≥ (N ♯ a α,z q , a α,z q ) − (N ♯ I q ãα,z q , I q ãα,z q ) → +∞.
We then conclude that lim α→0 I α (z) > 0. Case 2: We assume that z ∈ corresponds to the case where ω mis ̸ = ∅.In this domain n = 1 and n p ̸ = 1.If z is not in ω mis p (i.e. one of the periodic copies of ω mis ) then the same arguments as in [8] apply and show that D(a α,z q , ãα,z q ) remains bounded away from 0 as α → 0. This is due to the fact that h z ̸ = v z where v z and h z are defined in Lemma 1 (ii) and (iii), respectively.This implies lim α→0 I α (z) > 0. The case where z ∈ ω mis p \ ω mis can treated similarly to Case 1.By Lemma 1 (ii) and (iii), N ♯ a α,z q , a α,z q → ∞ while N ♯ a α,z q , a α,z q remains bounded as α → 0. Consequently D(a α,z q , ãα,z q ) → ∞ as α → 0. Lemma 1 (i) indicates that ⟨N ♯ a α,z , a α,z ⟩ remains bounded as α → 0. The last two items show that lim α→0 I α (z) > 0, which ends the proof of the theorem.

Numerical Studies
In this section we present some numerical examples using simulated data in R 3 showing viability of our imaging method.In our examples, the probing region is Ω M containing M 1 × M 2 periods including the defective cell.Data is generated by a family of incident plane waves as described in (4) with the indices j = (j 1 , j 2 ) in the set Here, q ∈ Z 2 M is fixed and N min , N max ∈ Z 2 + are given.The scattered wave associated with each individual incident wave is computed numerically by implementing the Floquet-Bloch transform and volume integral method described in [18].The collection of Rayleigh coefficients corresponding to ℓ ∈ Z 2 inc 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 u s (ℓ; j) ℓ,j∈Z Since in practice the measured data is always corrupted by noisy, we simulate the noise in our computed data.In particular, if δ > 0 is the noise level, the noisy near-field operator N δ is computed by adding random noise to N, that is where A = [A(j, ℓ)] N inc ×N inc is a matrix of uniform complex random variables with real and imaginary parts in [−1, 1] 2 .We define the functionals J δ α and J δ α,q associated with the noisy data, by We then consider a α,z δ , a α,z q,δ and ãα,z q,δ in ℓ 2 (Z 2 ) as minimizing sequences of J δ α ( Φ(•; z), a), J δ α ( Φ q (•; z), a), and J δ α,q ( Φ q (•; z), a), z) are positive for z inside the perturbation and in the periodic copies of the background that intersects the perturbation.However, we do not have quantitative information about the values of I α (z) at each point in the probing region.A careful observation of iso-surfaces in the reconstructions presented in Figure 5 suggests that the value of I(z) at the points z in the periodic copies of the part of the cube intersecting the defect is much smaller than the values at the points in the perturbed region.One can argue that this difference clearly distinguishes the true defect from the periodic artifacts.The presence of periodic copies of intersection subsets also indicates whether the defective region overlaps with components of the healthy periodic layer or not (compare to Example 1).Example 4..We end our numerical investigations by showing three more numerical examples illustrating a case where a whole component of the periodic domain D p is missing (see Figure 7), a case where only a part of a component of D p is missing (see Figure 8), and a case where both scenarios are present (see Figure 9).The 3D reconstruction obtained in Figure 7 confirms the result of our theory as we obtain periodic copies of the missing component except in Ω 0 .
In the case of a partially missing component presented in Figure 8, the 3D reconstruction resembles the case of an entire component missing: only the missing part is repeated periodically except in Ω 0 .In the last case presented in Figure 9

Figure 1 :
Figure 1: Sketch of the geometry for the L−periodic problem.The healthy bi-periodic layer consists of a homogenous layer occupying −h ≤ x 3 ≤ h with periodically distributed inhomogeneities indicated by red balls.Blue balls indicate the compactly supported perturbations ω located in one period denoted by Ω 0 .
) is defined by (52).Let us split domain D \ O c p into (D \ O c p ) \ D p and (D \ O c p ) ∩ D p , and treat each domain separately.Case 1: We assume that z ∈ (D \ O c

Figure 4 :
Figure 4: (a) exact geometry of the domain D with a defect being a ball (see details in the description of Example 1), (b) a slice of D in the y direction passing through the center of the ball, (c) a 3D view of the reconstruction given by (60), (d) a slice view of the indicator function I δ α that corresponds to (b).

Example 2 :
. The periodic background in this example is the same in Example 1.The perturbation ω is a ball that intersects the smaller cube (see Figures5(a) and 5(b)).In this example n = 4 in the part of the ball inside the cube and n = 3 in the rest of the ball (see Figure 5(b)).The 3D reconstruction using the indicator function (60) is displayed in Figure 5(c) and Figure 5(d) displays the indicator function at the slice plane.Theorem 3 indicates that the values of I δ α

Figure 5 :
Figure 5: (a) exact geometry of the domain D with a defect being a ball that intersects the periodic domain D p (see details in the description of Example 2), (b) a slice of D in the y direction passing through the center of the ball, (c) a 3D view of the reconstruction given by (60), (d) a slice view of the indicator function I δ α

Example 3 :
. In this example, we change the configuration of the periodic background which now includes two cubes and one ball in each period (see Figure6(a)) with refractive index n b = 2.We consider a local perturbation of the refractive index that changes the value of n b inside the ball in Ω 0 (see Figure 6(b)).We set n = 4 in this defective component.According to the theoretical result in Theorem 3 the indicator function I α (z) should visualize periodic copies of this ball.This is what is observed in Figures 6(c) and 6(d)).We also observe that the values of the indicator function are larger in the ball inside Ω 0 which may also be used as an indicator for the location of the defect.

Figure 6 :
Figure 6: (a) exact geometry of the domain D with a defect being a ball that coincides with one component of the periodic domain D p (see details in the description of Example 3), (b) a slice of D in the y direction passing through the center of the ball, (c) a 3D view of the reconstruction given by (60), (d) a slice view of the indicator function I δ α that corresponds to (b).
, we have a partially missing component and an additional defect in the form of ball intersecting the component that has a missing part.This is indeed a complication over the previous case and the obtained reconstructions are compatible with what one would expect from the theory, i.e. the indicator function shows the defective component repeated periodically.

Figure 7 :Figure 8 :
Figure 7: This example illustrates the case of a defect consisting of a missing component of D p as shown in (a).(b) displays a 3D view of the reconstruction given by (60).

Figure 9 :
Figure 9: (a) exact geometry of the domain D with a defect being partially missing component of the periodic domain D p and a ball intersecting this defective component, (b) a slice of D in the y direction passing through the middle of the defective component, (c) a 3D view of the reconstruction given by (60), (d) a slice view of the indicator function I δ α that corresponds to (b) h}.We denote by H 1 # (Θ h ) the restrictions to Θ h of functions that are in the Sobolev space H 1 loc (|x 3 | ≤ h) and are M L-periodic.The space H