Factorization Method for Imaging a Local Perturbation in Inhomogeneous Periodic Layers from Far Field Measurements

We analyze the Factorization method to reconstruct the geometry of a local defect in a periodic absorbing layer using almost only incident plane waves at a ﬁxed frequency. A crucial part of our analysis relies on the consideration of the range of a carefully designed far ﬁeld operator, which characterizes the geometry of the defect. We further provide some validating numerical results in a two dimensional setting.


Introduction
We consider the time harmonic inverse scattering problem from inhomogeneous infinite periodic layers that contain a local perturbation by developing the Factorization method [15] for reconstructing the support of the perturbation from measurements of scattered waves at a fixed frequency. Moreover, we indicate how other sampling methods can be applied (we refer to the monographs [11,8,9] for applications of these methods to other types of scattering problems).
The use of sampling methods to identify periodic media has been treated by several work in the literature and without being exhaustive, we refer to [19,20,2,1,21,12,24,14]. For the inverse problem in a locally perturbed periodic waveguide we refer to [23,7] and references therein. The specificity of our work comes from formulating the forward and inverse problems for incident plane waves and define far field data using measurements of the scattered field at finite distance surface. This setting allows us to provide a formulation of the sampling methods, and more specifically the Factorization method, that are similar to the one for bounded inhomogeneous inclusions in free space. Up to our knowledge, the sampling methods for locally perturbed infinite periodic layers and far field settings have not been treated in the literature. We refer to [13,10] where differential sampling methods have been designed to reconstruct defects in periodic backgrounds without knowledge of 1 the background using propagative and evanescent modes. However, the method applied in [13,10] has been only justified for the case of periodic defects with periodicity length that is the multiple of the background periodicity.
In this work the general problem of infinite periodic layer with a local perturbation is considered and we focus on the scalar problem. We make the assumption that the index of refraction of the background has at least a positive imaginary part on a non-empty open set. This prevents, for instance, guided modes and ensures well-posedness of the scattering problem from local defects and sources with bounded support in classical Sobolev spaces (see [16]). We explain how the results in [16] can be exploited to define scattered waves from incident plane waves, which furthermore allow us to define associated far fields in a similar manner as for the scattering in free space. We show that this definition coincides with scaled Fourier coefficients of the scattered field at a finite distance from the periodic layer.
The inverse problem is addressed by considering sets of incident plane waves that can be arbitrarily close to the set of propagative plane waves, where the associated far fields are defined as measurements. We provide in this setting the theoretical foundations of the Factorization method to reconstruct the support of the local perturbation from a well designed far field operator. However, our arguments does not cover the limiting case where only measurements of Fourier coefficients associated with propagating modes are used. The main difficulty comes from obtaining injectivity of the far field operator in this case. We refer to [5] where a closely related unique continuation principle is discussed for homogeneous layers (without absorption).
The absorption assumption for the background refractive index introduces difficulty in the design of the Factorization method since non-physical incident waves need to be introduced in order to obtain the desired factorization of the far field. Nevertheless, we explain in the numerical part, how one can build an approximation of this non-physical operator from the physical ones if the inclusion does not intersect the absorbing part of the background. The argument relies on analyzing the closure of the range of associated Herglotz operators. We make the observation that the reconstructions obtained by the approximate operator or the physical one are qualitatively the same. After recalling the Linear Sampling Method, the performance of the sampling algorithms is tested for two dimensional configurations and some toy problems including some numerical examples that mimic non-destructive testing of nano-grass structures. Although not covered by the theory, we use measurements that only corresponds to propagative plane waves for the numerical experiments.
The outline of the paper is as follows. In Section 2 the scattering problem for plane waves and associated far fields is introduced. This necessitates the definition of quasiperiodic scattering problems for the unperturbed periodic background. We then relate the far fields for the perturbed problem to the Fourier coefficients of the scattered field. Section 3 is dedicated to the analysis of the Factorization method. The theory is provided for the "non-physical" far field operator associated with incident fields that satisfies the Helmholtz equation with conjugated refractive index. In the last two section we show some numerical examples in the two dimensional setting, where in Section 4 the"non-physical" setting and in Section 5 the "physical" far field operator is considered.

Setting of the scattering problem
In this section, we introduce the scattering problem for a periodic refractive index with a local perturbation and define the far field of the scattered field. We first introduce the (unperturbed) quasi-periodic scattering problem for plane waves as incident fields, since we will consider the quasi-periodic total field as the incident field in the perturbed problem.
Let k > 0 be the wave number and fix some small η > 0 for the rest of the paper, if not said otherwise. For the numerical part, we will drop this requirement and use η = 0. We call a function u : R 2 → R 2 α-quasi-periodic for some α ∈ R, if the identity The square root is hereby extended by a branch cut at the negative imaginary axis. We can rewrite the plane waves as u inc (x, d) = e −ikd·x with some d ∈ R × C, such that d · d = 1, and |Im d 2 | ≤ η. We define the set of all directions as and consider α as the function depending on d by defining α : S → − k 2 + k 2 η 2 , k 2 + k 2 η 2 , α(d) := kd 1 . (1)

The quasi-periodic scattering problem
Let R ≥ R 0 > 0 and n p ∈ L ∞ (R 2 ) be a 2π-periodic refractive index w.r.t. first component x 1 , for which it holds n p (x) = 1 for |x 2 | ≥ R 0 . Since the parameter is periodic and the incident fields are α-quasi-periodic, we can reduce the problem to one periodic cell Ω R 0 := (−π, π) × (−R, R) with the upper boundary Γ R 0 := (−π, π) × {R} and the lower boundary Γ −R 0 := (−π, π) × {−R}. For some d ∈ S and the incident field u inc (·, d) the scattering problem is to find the α-quasi-periodic scattered field u s qp , which lies in the space H 1 α (Ω R 0 ) defined as the closure of regular α-quasi-periodic functions with respect to the H 1 (Ω R 0 )-norm. We analogously define the spaces H s α (Γ ±R 0 ) for s ∈ R. For a regular α-quasi-periodic function φ defined on Γ R 0 (and analogously for functions defined on Γ −R 0 ), we define the j-th Fourier coefficient Since the refractive index equals to one for x 2 ∈ (−R 0 , R 0 ), we can explicitly write down the solution of the form is the j-th Fourier coefficient of the trace of u s qp on Γ R 0 , as well as which is called Rayleigh radiation condition. The radiation condition allows us to define the α-quasi-periodic Dirichlet-to-Neumann operator T ±R α : H , which is given by . Hence, we seek the scattered field u s qp ∈ H 1 α (Ω R 0 ), which solves the variational problem related to Then, referring for instance to [17] one has the following theorem.
Theorem 1. Assume that the refractive index satisfies Im n p ≥ 0. Moreover, let k be outside the countable set of eigenvalues, or, let the set {Im n p > 0} contain an open ball. Then there exists a unique solution u s qp ∈ H 1 α (Ω R 0 ) to the quasi-periodic scattering problem (2).

The locally perturbed scattering problem
Having the total field of the periodic problem, we consider the scattering problem in the case of the locally perturbed refractive index n := n p + q ∈ L ∞ (R 2 ), where the perturbation is described by some q ∈ L ∞ (R 2 ) with the compact support D := supp(q) ⊆ Ω R 0 0 . In Figure 1 4 one can see an example of such a perturbed refractive index. The perturbation prevents the reduction of the problem to one periodic cell and we have to treat the problem as an unbounded domain problem. Define for R ≥ 0 the sets The perturbed scattering problem is described by the Helmholtz equation, where the total wave field u t : R 2 → C satisfies Instead of solving the problem for the total field, we consider for some d ∈ S the total field u t as a superposition of a scattered field u s and the total field of the quasi-periodic scattering problem u qp (·, d) as the incident field, and seek for the scattered field u s ∈ H 1 (Ω R ) for every R > R 0 . We can reformulate the scattered problem for the scattered field u s = u s (·, d) that satisfies the equation Analogously to the quasi-periodic scattering problem, we can prescribe the solution in R 2 \ Ω R with the so-called angular spectrum representation of the form where u s is the Fourier transform of the trace of u s on Γ ±R . Moreover, we obtain the exterior Dirichlet-to-Neumann map T ±R by which is a bounded linear operator from H 1 /2 (Γ ±R ) to H −1 /2 (Γ ±R ).
Assuming there exists a unique solution for the scattered field, the total field u t = u s + u qp (·, d) would solve the scattering problem (3). For the Factorization method, however, we shall modify the incident field by using the complex complemented function u qp instead of u qp on the right hand side of (4). In consequence, we consider the following scattering problem that we formulate for a general source term f ∈ L 2 (Ω R 0 ) that we extends by 0 outside Ω R 0 .
Problem 1. Find a function u s ∈ H 1 (Ω R ) for some R > R 0 as a solution to the variational problem related to In order to ensure that this problem is well-posed one needs additional assumptions on the index of refraction n p . Following [16, Theorem 1], we shall make the following assumption.
Assumption 1. We assume that {Im n p > 0} contains an open ball and the refractive index satisfies Im n p ≥ 0 and Re n p ≥ 0.
We then obtain the following result that can be proved similarly to [16, Theorem 1].
Theorem 2. If Assumption 1 holds and if Im q ≥ 0, then there exists a unique solution to Problem 1.

The far field of a solution
We proceed with the definition of the far field for a solution to Problem 1. We shall utilize the horizontal Bloch-Floquet transform J R , which is defined by for smooth functions with compact support φ ∈ C ∞ 0 (Ω R ). The Bloch-Floquet transform extends for s ∈ R to an isomorphism between H s (Ω R ) and L 2 (I; H s α (Ω R 0 )), as well as between H s (Γ ±R ) and L 2 (I; H s α (Γ ±R 0 )), where the index α indicates that the space depends on α ∈ I (see [18]). For some φ ∈ C ∞ 0 (Ω R ) and w := J R φ the inverse of the transform is given by Definition 3. Let u s ∈ H 1 (Ω R ) be the solution to Problem 1 for some right hand side f ∈ L 2 (Ω R 0 ). We consider β(d ) := kd 1 as a function of d ∈ S, like α in (1), and set the notation We define the far field u ∞ ∈ L 2 (S) of u s by The following lemma gives a first characterization of the far field.
Proof. Fix d ∈ S and set β = β(d ), then we obtain by the definition of the inverse of the Bloch-Floquet transform the identity The smooth and periodic function δ −β,M (α) := M m=−M e −i2πm(β+α) is the truncated Fourier series expansion of the 1-periodic delta distribution, since the j-th Fourier coefficient is given by δ −β,j = e −i2πjβ , the sum converges to δ −β with respect to H −s per (I) for some s > 1 /2 (see [22,Chapter 5]). Since the Bloch transformed right hand side is analytical in β, we know by [16, Theorem 8] that the scattered field is continuous in β ∈ I and the evaluation at point α is allowed, such the the evaluation of the distribution is well-defined. Hence, we conclude that for the zeroth coefficient. Therefore, we can reduce the representation to This implicates, in particular, that a vanishing far field is equivalent to the condition that the functions Proof. If the far field is zero everywhere, then the Fourier transform of the trace F u s Since the right hand side f is compactly supported in Ω R 0 , the Bloch-Floquet transformed function is given by f = J R f , which is independent of the artificial component β. In particular, the right hand side is analytical w.r.t. β and it follows by [16, Theorem 8], that the solution J R u s is analytical in β ∈ R \ {−k, k}. Therefore, the Fourier transform of the traces F u s Γ ±R are analytical in ξ ∈ R \ {−k, k} and vanishes everywhere in R, and in particular, the traces u s Γ ±R equal to zero on Γ ±R . Since the scattered field can be extended by the radiation condition, which is uniquely determined by the trace in Ω R \ Ω R , the scattered field is zero for some R > R, and there exists an open set, where the scattered field is vanishing. The right hand side f is only supported in D and hence, the theorem of unique continuation implies that the scattered field u s is vanishing everywhere outside of D.

The Factorization method
The aim of this section is to introduce the Factorization method to localize the support of the perturbation. The Factorization method gives a theoretically justified numerical scheme for the inverse problem of reconstructing the shape of the perturbation by knowing the far fields of the scattered waves induced by plane waves as the incident fields for some d ∈ S.
The key idea is to define a far field operator F and decompose it into three operators F = H * qp T H qp , where T is a linear and invertible operator. Having the decomposition, we will derive a binary criterion, whether some point z ∈ R 2 lies inside of D, or, outside of D.
Note that, if u inc (·, d) is the plane wave depending on d ∈ S, we write u s (·, d) and u qp (·, d) to show the dependency on d, u s (·, d, d ) for the Bloch-Floquet transformed scattered field, and u ∞ (d, d ) for the far field in d ∈ S. Before we define the Herglotz operator H qp , we need to define integral over S for functions defined on S. The set S can be decomposed into the two sets Therefore, we define the two functions φ ± : − 1 + η 2 , 1 + η 2 → C depending on d 1 for a function φ : S → C as the restrictions of φ to S ± . We then define the integral of We define the Herglotz operator H qp : L 2 (S) → L 2 (D) and the far field operator F : L 2 (S) → L 2 (S) by where u ∞ (d, ·) is the far field of the scattered field u s (·, d) to Problem 1 with the right hand side f = k 2 q u qp (·, d). We denote the adjoint operator of H qp w.r.t. the L 2 (D) scalar product by H * qp : L 2 (D) → L 2 (S). A simple calculation shows that the adjoint operator is given by We observe that by changing the definition of the incident field, we have changed the nature of the considered scattering problem. This change allows us to obtain the desired symmetric factorization of the far field as stated in the lemma below. However, this change does not correspond to what physically one can observe if the index of refraction n p has non-vanishing imaginary part. This is why we shall also discuss in the numerical section when and how one can approximate this non-physical operator from the one corresponding to f = k 2 q u qp (·, d).
In addition to the assumptions made to obtain a well-posed forward problem, we need the following assumption to show the decomposition.
Assumption 2. We assume that Im q ≥ 0 and Re q ≥ c > 0 on D and D c := R 2 \ D is connected.
Lemma 7. The operator F can be decomposed as F = H * qp T H qp with the invertible bounded operator T : where v s is the scattered field to the scattering problem 1 with the right hand side f = k 2 qv.
Proof. (a) It holds for β = β(d ) and φ ∈ H wherefrom we obtain In consequence, we can add this term to the far field to derive the identity Since the integrand is periodic w.r.t. x 1 , we can apply the Green formula to obtain If we multiply the function g ∈ L 2 (S) to the far field and integrate over S, we get the claimed decomposition by where U s (x) is the scattered field S u s (x, d)g(d) dS(d) and U inc = H qp (g).
(b) It is assumed in Assumption 2 that Re q ≥ c > 0 holds on D. If we call L : L 2 (D) → L 2 (D) the solution operator, which maps some right hand side f ∈ L 2 (D) to the solution of Problem 1 restricted to D, and divide T v by k 2 q for a v ∈ L 2 (D), we obtain which is a Fredholm operator of order zero. Thus, it remains to check the injectivity.
Let v ∈ L 2 (D) satisfy T v = 0 in D. Therefore, it holds k 2 qv s = −k 2 qv in D, and v s satisfies ∆v s + k 2 (n p + q)v s = k 2 qv s in Ω R .
Since this scattering problem for q = 0 is uniquely solvable by Theorem 2, v s has to be the zero function in Ω R and we conclude v = v s = 0 in D.
Fix z ∈ (−π, π)×R, and let Φ z,per be the periodic fundamental solution to the Helmholtz equation Let χ be some cut-off function for some δ > 0, where χ(x) = 0 for |x − z| > δ and χ(x) = 1 for |x − z| < δ /2 and Φ z the solution to Problem 1 with q = 0 and then the fundamental solution to Problem 1 with q = 0 can be defined as Φ z := χΦ z,per + Φ z .
Theorem 8. Assume that the domain D has a connected complement. Let φ ∞ z be the far field of the (fundamental) solution Φ z of Problem 1 with the right hand side f = δ z for some z ∈ Ω R 0 . Then it holds z ∈ D, if and only if φ ∞ z is in the range R(H * qp ).
Proof. Fix z ∈ D. Then there exists an open ball B ε (z) with radius ε > 0, such that the closure B ε (z) lies in D. We choose some cut-off function χ, which fulfills χ(x) = 1 if |x − z| ≥ ε and χ(x) = 0 if |x − z| ≤ ε /2 and consider the function w := χΦ z ∈ H 1 (Ω R ). Define moreover the function ψ := ∆w + k 2 n p w ∈ L 2 (Ω R ), which is supported in D. By Lemma 4 and the Green formula, it follows that and thus the Green formula and Lemma 4 imply the identity Consequently, the far field φ ∞ z of Φ z is in the range of H * qp . Assume otherwise, that φ ∞ z ∈ R(H * qp ), but z ∈ D. Since the far field is in the range, there exist a solution w s ∈ H 1 (Ω R ) to the scattering problem, for which the far field is given by φ ∞ z . By Lemma 6, we conclude that w s equals the fundamental solution Φ z on R 2 \ (D ∪ {z}). But, the analyticity of the solution contradicts the case, that z ∈ R 2 \ D. If z ∈ D \ D, then there is also a contradiction, since Φ z ∈ H 1 /2 (∂D) for z ∈ ∂D.
We showed that a point z ∈ Ω R 0 lies inside D if and only if the far field φ ∞ z is an element in the range of the operator H * qp . The next step is to derive a characterization of the range of H * qp in terms of F . For that, we will apply [15,Theorem 2.15] in the version of [6, Theorem 3.2], that we give in the following proposition. Proof. The total field of the quasi-periodic problem is continuous in d 1 and analytical outside of d 1 = ±1 (see [16, Theorem 7]). Since the total field is the kernel of the integral operator H qp the operator H qp is compact.
For the injectivity, set I := [−1, 1], Hg := (H qp g) and note that Hg ∈ L 2 (D) is the restriction of the function which solves the homogeneous Helmholtz equation ∆v g + k 2 n 2 p v g = 0 in R 2 . At first, we compute the Bloch-Floquet transform of v g . For that, we choose a smooth function g ± ∈ C ∞ 0 ( I), which we can extend by zero to R. Since all of the derivatives of g ± vanish in d 1 = ±1 and the quasi-periodic total field is analytical outside of d 1 = ±1, we conclude for Further, u qp (·, d) ∈ H 2 loc (R 2 ) for all d ∈ S, wherefrom we conclude by the Sobolev inequalities that u qp is continuous in x. Therefore, we obtain pointwise The smooth periodic function d 1 → M j=−M e 2πi(α 0 − d 1 )j converges in the distributional sense (i.e. in D (R)) to the periodic delta distribution m∈Z δ α 0 +m uniformly w.r.t. α 0 ∈ I (see [22,Chapter 5]). Hence, for every α 0 ∈ I and x ∈ Ω R 0 the limit is given by the finite sum The limit still holds w.r.t. the L 2 (I × Ω R 0 )-norm, as one concludes by the Fourier series representation of the limit in [22,Chapter 5].
To extend the formula to general g ∈ L 2 ( I), set u ± per (x, kd 1 ) := j∈Z u ± g (x, d 1 + j) with u ± g defined in (6). For the L 2 (Ω R )-norm, we see analogously to (7) that and therefore, the function v g lies in L 2 (Ω R ) and it holds Therefore, the equation (8) holds for every g ∈ L 2 (S) with supp g ± ⊆ I. Let g ∈ L 2 (S) be some function of the kernel satisfying Hg = 0 ∈ L 2 (D). The function v g solves the homogeneous Helmholtz equation and vanishes on D and the theorem of unique continuation implies v g = 0 in R 2 . The periodic scattered field is bounded by || u s qp (·, d)|| H 1 (Ω R 0 ) ≤ CR || u s qp (·, d)|| H 1 (Ω R 0 0 ) for every R > R 0 and the incident field u inc is exponentially growing in ∓x 2 direction for |d 1 | > 1 on S ± . From this it follows that supp g ± ⊆ I must hold. Hence, we can the identity (8) holds for g and we can apply the boundary operator (T ±R α 0 − ∂ /∂x 2 ) for some α 0 ∈ I to the trace of the total field. For Since the functions {e −i(α 0 +m)x 1 } m∈Z are linear independent in L 2 ([−π, π]), we finally have g ± = 0 in I, and in consequence, g = 0 in S.
For the characterization of the range of H qp , note that clearly R(H qp ) ⊆ H inc (D).
To prove the inverse inclusion, it is sufficient to prove injectivity of the operator H * qp on H inc (D). Let φ ∈ H inc (D) be some function with H * qp φ = 0 and consider the scattered wave w ∈ H 1 (Ω R ) as the solution to Problem 1 satisfying ∆w +k 2 n 2 p w = φ. Thus, the application of the Green formula gives for every d ∈ S Now, Lemma 6 implicates that w = 0 on D c , and in particular, w ∈ H 1 0 (D). Therefore, since φ ∈ H inc (D), we get and consequently, H qp is one-to-one on H inc (D).
The denseness of the range of H * qp is a direct implication of the injectivity of H qp .
Corollary 11. The operator H * qp is compact with dense range.
For given functions n, n p and q we define the homogeneous interior transmission problem as the problem for (w s , v) ∈ L 2 (D) 2 satisfying w s ∈ H 2 (D) and For the Factorization method we will assume that there is no non-trivial solution to this problem. We now prove the needed properties for the middle operator T .
Lemma 12. Assume that the assumptions 1 and 2 hold and assume further, that there is no non-trivial solution to (9). Then the operator (Re T ) can be decomposed into (Re T ) = C + K, where C is a self-adjoint and coercive operator and K is compact, and the operator (Im T ) is strictly positive on R(H qp ).
Proof. For every v ∈ L 2 (D), we call v s the solution to Problem 1 with the right hand side k 2 qv. Define C : L 2 (D) → L 2 (D) by Cv := −k 2 (Re q)v, which is a coercive and self-adjoint operator. For the remaining part K := (Re T ) − C : L 2 (D) → L 2 (D) and for w, v ∈ L 2 (D) we have wherefrom we conclude that K is a compact operator, since w s and v s ∈ H 1 (D). It remains to show the positiveness of (Im T ). At first, we see, that Hence, we will consider −Im (T v, v) L 2 (D) instead of (v, (Im T )v) L 2 (D) . Further, we can write the scalar product (T v, v) L 2 (D) as we conclude and therefore, (Im T ) is a semi-positive operator. Assume now, that for some v ∈ R(H qp ) it holds Im (T v, v) L 2 (D) = 0. The calculation (10) implicates that v s vanishes on {Im n p > 0} and the theorem of unique continuation gives v s = 0 on the connected set D c including {Im n p > 0} \ D. Since we assumed to have no transmission eigenfunctions, the scattered field v s has to vanish everywhere and in consequence v = 0 holds as claimed.
The last two results and Theorem 8 show that we fulfill the requirements to apply Proposition 9 and obtain the following characterization result as the theoretical justification for the numerical approach we will discuss in the next section.
Theorem 13. Under the same hypothesis as in Lemma 12 we have that the operator F # := |(Re F )| + (Im F ) : L 2 (S) → L 2 (S) is strictly positive. If {λ j , ψ j } ∞ j=1 is an eigensystem of F # , then we can characterize the support D by 4 Numerical examples related to the theoretical part In this section we present some numerical examples that validate the theoretical results of the previous section. For the numerics, we consider η = 0 instead of an arbitrary small η > 0 and discretize the interval We set d j,± 2 := ± 1 − |d j 1 | 2 and define the directions d j,± := (d j 1 , d j,± 2 ). To generate the data we take the plane waves u inc (·, d j,± ) as the right hand side and compute the quasiperiodic scattered fields by the Finite-Elements method. Having the approximated total field of the quasi-periodic scattering problem, we approximate the scattered field to the perturbed problem with the help of the Bloch-Floquet transform based algorithm described in [16]. Moreover, we approximate the fundamental solution for the unperturbed problem by solving the variational problem corresponding to Problem 1 with q = 0 and f = δ z by applying the Bloch-Floquet transform based algorithm. Afterwards, we can compute the far fields with the formula given in Remark 5 for the same points {d j 1 } N d j=1 . We define the discrete far field operator by where F ± j ∈ C 2N d , j = 1, . . . , N d , is the vector of values of the far field for the incident field u inc (·, d j,± ) evaluated at the points d j,+ in the first part, and at d j,− in the second part. If we call {( λ j , ψ j )} 2N d j=1 the discrete eigensystem of F # := |(Re F )| + (Im F ) and φ ∞ z the discrete far field of the fundamental solution, we compute for some z ∈ Ω R 0 the sum which is a discrete version of (11), and visualize the values 1 /Mz for some mesh points z ∈ [−π, π] × [0, 5].
For the examples we give the results with the far field operator without noise and with relative additive noise of 1% with unified distribution around zero, i.e., we use the date given by F ε = F + c with c ∈ C 2 N d , ||c|| 2 = 0.01 || F || 2 , where || · || 2 is the biggest singular value.
To generate the data, we set the cut-off of the Dirichlet-to-Neumann operator to |j| ≤ 300, discretize the domain Ω R 0 in 2 2·6 rectangular cells, and the interval I in 2 7 subintervals. The relative tolerance for GMRES is chosen as 10 −10 , and we use a variable transform for the discrete inverse Bloch-Floquet transform for better approximation (see [16] for the description of the parameter). We set N d = 64, which means, that we take 128 direction for the incident fields and the far field. For the right hand sides we choose 25 2 = 625 points {z j } 625 j=1 inside Ω R 0 , which are equally distributed in every direction.

Example 1
For the first example we look at a simple model of nano-grass and simulate a production failure where some rods are missing. In Figure 2  (e) Real part of second perturbed refractive index between −3π and 3π. Imaginary part stays the same.
(f) Result for second refractive index with 1% noise for k 2 = 3. In the next examples, we will analyze the influence of the reconstruction quality w.r.t. the wave number and the absorption area. For that, we look at the unperturbed refractive index where ν > 0 is some number we will choose later, and the perturbation Figure 3 for the parameter).
(a) Real part of refractive index Re np + q1 between −3π and 3π.
(b) Imaginary part of refractive index Im np between −3π and 3π.

Example 2
For the second example we look at the influence of the value of the absorption area. We choose a relatively small ν = 10 −3 and an even smaller ν = 10 −6 with k 2 = 0.4 and q = q 1 .
The results with and without noise are visualized in Figure 4 for (a) ν = 10 −3 without extra noise.

Example 4
This example shows the reconstruction results for the small wave number k 2 = 0.09 with ν = 10 −3 , as well as for the larger wave number of k 2 = 3 with ν = 1. The results with and without noise for q = q 1 are visualized in Figure 6 for Ω R 0 = [−π, π] × [0, 5]. As a conclusion, the binary criterion for the shape reconstruction of the perturbation, which is theoretically justified in Theorem 13, has been confirmed by the numerical examples although the quality of the reconstructions is less sharp than those observed in the free space.   The reconstruction quality of the Factorization method seems to be less dependent on the imaginary part of the refractive index ν (compare Figure 4). This differs from the Newtonlike methods, where a small ν decreases the reconstruction quality and increases the number of steps (compare remark in [16]). Nevertheless, the reconstruction quality depends on the adequate choice of the wave number for the size of the perturbation. If the wave number is too large, the quality of the reconstruction of the Factorization method decreases and a smaller ν has a negative impact on it, as we see for k 2 = 1.8 in Figure 5. For smaller objects a larger wave number gives better results (compare Figure 2).

Numerical examples for physical data
Up to now, we were using the non-physical discrete far field operator F for the numerical examples, for which the scattered waves satisfy Problem 1 with the right hand sides f = k 2 q u qp (·, d). In this section, we will approximate the non-physical discrete far field F ε approx operator from the data given by the physical discrete far field operator F ε p , which is derived from scattered waves satisfying Problem 1 with the right hand side f = k 2 q u qp (·, d). The physical data is numerically generated using the same finite element method as the one described in the previous section.
For this, note that in the case D∩{Im n p > 0} = ∅ we can prove analogously to Lemma 10 that the closure of the range of H qp : L 2 (S) → L 2 ( Ω), where Ω := Ω R 0 \ {Im n p > 0}, is the space H inc ( Ω). In the same way we can show that H inc ( Ω) is the closure of the range of H T qp : L 2 (S) → L 2 ( Ω), H T qp g := S u qp (·, d)g(d) dS(d) by observing that H T qp g = H qp g. Since the physical far field operator F p has the decomposition F p = H * qp T H T qp , and the theoretical can be decomposed in F = H * qp T H qp , we approximate the discrete theoretical far field operator F ε approx row by row from the data given by F ε p by approximating a solution g δ i to the discrete equation for some parameter δ > 0 and the canonical vectors e i ∈ R 2 N d , i = 1, . . . , 2N d , with Tikhonov as the regularization method, such that F e i ≈ F ε p g δ i =: F ε approx e i .
For the discretization of the Herglotz operator, we use a grid of 100 · 95 uniformly distributed points in each direction for [−π, π] × [ 1 /2, 5]. We apply the L-curve criteria plotting ||( to derive an estimation of 3% for the noise level and regularize with Tikhonov for a decreasing sequence of δ → 0 until holds for the first time. For comparison, we will furthermore include the results of the Linear Sampling Method (LSM) (see for instance [9, Chapter 2]), for which we approximate a solution to the equation F ε p g z ≈ φ ∞ z by minimizing the Tikhonov functional together with the Morozov principle, i.e., we compute a δ z > 0, such that || F ε p g δz z − φ ∞ z || 2 = 0.01 || F ε p || 2 ||g δz z || 2 . The image is obtained by plotting z → 1/||g δz z ||. An alternative method would be to use the generalized version of the Linear Sampling Method as proposed in [4,3]. We did not test the performances of that method and postpone it to a future work.
In Figure 7 we show the results for the approach of the approximated theoretical far field operator compared to the results for the theoretical far field operator and for the Linear Sampling Method. For every example we use k 2 = 3 and add unified distributed noise of 1% to F and to F p in the way we described in the last section. The relative residuum between the theoretical discrete far field operator F and the approximated theoretical discrete far field operator F ε approx is around 33%, nevertheless, the localization of the perturbation still works. Moreover, we added the result of applying the Factorization method directly to the physical far field operator for comparison. Although we do not have a justification for the Factorization method for the noisy physical data, the results in (a) and (c) of Figure 7 are very comparable and somewhat better than (b) and (d). For the other nano-grass simulation example we used in Figure 2 we observed the same quality of numerical results.