The factorization method for the Drude-Born-Fedorov model for periodic chiral structures

We consider the electromagnetic inverse scattering problem for the 
Drude-Born-Fedorov model for periodic chiral structures known as chiral gratings 
both in $\mathbb{R}^2$ and $\mathbb{R}^3$. The Factorization method 
is studied as an analytical as well as a numerical tool for solving this inverse problem. 
The method constructs a simple criterion for characterizing shape of the periodic scatterer 
which leads to a fast imaging algorithm. This criterion is necessary and sufficient which 
 gives a uniqueness result in shape reconstruction of the scatterer. 
The required data consists of certain components of Rayleigh sequences of 
(measured) scattered fields caused by plane incident electromagnetic waves. 
We propose in this electromagnetic plane wave setting a rigorous analysis for the Factorization method. 
Numerical examples in two and three dimensions are also presented for showing the efficiency of the method.

1. Introduction. We consider inverse scattering of time-harmonic electromagnetic waves from penetrable periodic chiral structures known as chiral gratings both in R 2 and R 3 . For the sake of presentation, the grating in two-dimensional setting is supposed to be periodic in the x 1 -direction, invariant in the x 3 -direction, and bounded in the x 2 -direction, while in R 3 it is assumed to be periodic in the x 1 -and x 2 -direction, and bounded in the x 3 -direction. The inverse problem of our interest is the shape reconstruction of chiral gratings from measured data consisting, in principle, of scattered electromagnetic waves, when plane electromagnetic waves are used as incident waves. There has been recently a considerable interest on periodic chiral structures from the engineering and applied physics community, see for instance [11,15,17,18,20,24,25,27,38], thanks to their applications in photonics, such as optical filters, fiber optic sensors, optical switch and polarization mode converters.
In the applied math community, the forward problem for chiral media has excited interest during the past years, see for example [1,2,5,12,13,26,[34][35][36]. In particular, we refer to [37] for recent mathematical developments on chiral media. However, to the knowledge of the author, the number of works related to the inverse problem for chiral media is relatively small. Results on uniqueness can be found in [6,7,16,[31][32][33]41]. The paper [16] has further presented an iterative method of Newton-type for identifying shape of two-dimensional chiral scatterers. We refer to [9,10] for results on quantitative reconstruction of material parameters for the Helmholtz problem. Until recently, the shape identification problem for the case of electromagnetic waves has been studied in [19]. More precisely, the author has analyzed the Factorization method, introduced by A. Kirsch [21], as a tool for imaging the chiral scatterer. This Factorization method belongs to class of fast non-iterative methods, recently developed for solving inverse scattering problems. With its simple implementation and rigorous mathematical basis, the Factorization method has appeared as an ideal tool for imaging problems, see [23] for recent developments on this method.
We note that the result in [19] has been illustrated by numerical examples in 2D. Furthermore, the author in the latter work also stated in the conclusion that "A rigorous proof of the Factorization method for chiral periodic structures would be a good starting point to investigate scattering from chiral metamaterials. In contrast to the reconstruction of the scatterer, it would be interesting how the additional degree of freedom β influences the possibilities for cloaking". Motivated by this conclusion and applications of non-destructive evaluation, we aim in the present work to study the Factorization method as a tool for reconstructing shape of a chiral grating from measurements of scattered electromagnetic waves. More specifically, these measurements are the Rayleigh coefficients of evanescent and propagating modes of the scattered fields caused by incident electromagnetic plane waves. We prove that the Factorization method is able to solve this shape identification problem and that the method provides an efficient tool to image the structure via two-and threedimensional numerical experiments. Our analysis basically extends approaches for inverse scattering problems for bounded chiral media in [19] and periodic achiral media in [29] to the case of periodic chiral structures. It is worth to mention that we exploit the special plane incident fields introduced in [3,29] which allows us to suitably factorize the measurement operator. Further, a modified version of the method studied in [28] treats the case that the imaginary part of the middle operator in the factorization is just semidefinite. As mentioned above, the numerical examples in R 3 presented in this work are, to the best of our knowledge, the first numerical examples for this method in a chiral electromagnetic setting. Now let us describe some setting for Maxwell's equations in periodic chiral media. We consider scattering of time-harmonic electromagnetic waves at positive frequency ω by a periodic grating consisting of chiral dielectric materials. We assume that the electric permittivity ε, the magnetic permeability µ and the chiral admittance β are real-valued scalar functions in L ∞ (R 3 ). Further, there exists positive constants ε 0 and µ 0 such that ε = ε 0 , µ = µ 0 and β = 0 outside of the grating.
We denote by E the electric field, H the magnetic field while B and D are respectively the magnetic induction and the electric displacement. We assume transmission conditions on the boundary of the grating, see [19] for a detailed derivation. The fields are described by Maxwell's equations (with no free charge and current density) where constitutive relations are described by the Drude-Born-Fedorov model as 2. The Helmholtz case.
2.1. Problem setting. As being said in the introduction we consider in this case the periodic structure or the grating which is 2π-periodic in the x 1 -direction, invariant in the x 3 -direction and bounded in the x 2 -direction. In this case, we assume that the material parameters ε r , µ r and β are all 2π-periodic in x 1 , and that they do not depend on x 3 . Moreover, the incident field is also assumed to be invariant in x 3 and that its incident vector is orthogonal to the invariant axis of the grating. We then find that the system 3 is independent of x 3 . From a straightforward calculation (see [19]) we obtain, under the assumption 1 − k 2 ε r µ r β 2 = 0, the third components of the electric and magnetic fields satisfy Notice that the operation div in 4 is understood as The grating is illuminated by an electromagnetic plane wave with wave vector d = The polarizations p, s ∈ C 3 of the incident wave satisfy p · d = 0 and s = 1/(ωε 0 )(p × d). With these definitions, the incident plane waves E i and H i are given by Therefore if we define then u i : R 2 → C 2 is an α-quasi-periodic incident field in x 1 , i.e., for α = d 1 , n ∈ Z, and further Here we note that ∆ acts componentwise on vectorial functions. Therefore with the incident field u i we now reformulate 4 in terms of the scattered field u, defined by where P , Q are the contrasts defined by Here I n is the n × n identity matrix. Define for n ∈ Z, It is well-known for periodic scattering problems that (see, e.g., [4,8]) the scattered field u is required to be also α-quasi-periodic in x 1 and that it admits a Rayleigh expansion radiation condition of the form where h > sup{|x 2 | : (x 1 , x 2 ) ∈ supp(Q) ∪ supp(P )}, and the Rayleigh sequenceŝ u ± n ∈ C 2 are given bŷ From now, a function which satisfies 8 is called to be radiating. Note that we require that the series in 8 converges uniformly on compact subsets of {|x 2 | > h}. Further, note that only a finite number of terms in 8 are propagating plane waves which are called propagating modes, the rest are evanescent modes which correspond to exponentially decaying terms. For a technical reason we assume in the following that k is not a Wood's anomaly (or, equivalently, that the frequency ω is not a Rayleigh frequency), i.e., (9) β n = 0 for all n ∈ Z.
The technical reason behind this assumption is basically that the representation of the Green function that we introduce later on (see 27) is not well-defined at a Wood's anomaly. Now we define The problem 7-8 can be reduced to one period Ω due to its periodicity. We now consider a more general problem as follows: Given f ∈ L 2 (D, C 4 ), g ∈ L 2 (D, C 2 ). Find u : and additionally the radiation condition 8 where the matrices A and B are given in 4.
Well-posedness of the forward problem for periodic chiral structures has been done in [1] for the full Maxwell's equations. Variational solution theory for the scattering problem 8-10 can be done similarly as the achiral case, see for example [8].
Fredholm theory implies that existence of solution for problem 11 follows from uniqueness of solution. It is also well known that if −Im A, Im B are positive definite, and either A, B are globally smooth or Im B is strictly positive definite on D (absorbing medium), then one can obtain uniqueness of solution, see for instance [19]. We also refer to [8] for strategies to establish uniqueness of solution for all but possibly a discrete set of wave numbers k when the coefficients are not globally smooth and the grating is non-absorbing. In this work we assume that the wave number k is such that uniqueness of solution holds.

2.2.
Factorization of the near field operator. In this section we first set up the inverse problem relying on the near field operator. Secondly we factorize the near field operator which is an important step for the analysis of the factorization method.
To obtain data for the factorization method we use the incident fields as α-quasiperiodic plane waves where the polarizations p 1,j = 0. (13) Note that this special class of incident plane waves has been used and discussed in [3,29]. For the analysis of the factorization method we need the following assumption. Assumption 1. We assume that the periodic support D is a Lipschitz domain. We assume further that Ω \ D consists of at most two connected components and that each connected component of Ω \ D is unbounded. For all ξ ∈ C 4 , z ∈ C 2 and almost all x ∈ D, it is assumed that Im (Q(x)ξ · ξ) ≤ 0 and Im (P (x)z · z) ≥ 0, and that there exists positive constants C 1 and C 2 such that Re (Q(x)ξ · ξ) ≥ C 1 |ξ| 2 and Re (P (x)z · z) ≥ C 2 |z| 2 . Furthermore, the well-defined inverse (Re Q) −1 belong to L ∞ (D, C 4×4 ).
Thanks to well-posedness of the forward problem we can define the solution operator G : . Due to the linearity of problem 14, a linear combination of several incident fields will lead to a corresponding linear combination of resulting scattered fields. We obtain such linear combination using sequences and define the corresponding operator H : Here we divide by β j w ± j to simplify the computation. In our inverse problem the data measured are the Rayleigh sequences 8. It is known in previous works [3,29] that the factorization method for periodic structures needs near field measurements. We define the near field operator N : 2 (Z, C 4 ) → 2 (Z, C 4 ) to map a sequence (a j ) j∈Z to the Rayleigh sequences of the scattered field generated by the virtual incident field H(a j ) defined in 15, i.e. (16) [N (a j )] n := (û n ) n∈Z , is the radiating solution to 14 for the source (f, g) = H(a j ). Here we remark that H(a j ) is not the physical incident field used in our inverse problem. It is introduced just for analysis of the factorization method. If we compare to the original problem 7, then the solution to 14, where (f, g) = H(a j ), is exactly the scattered field generated by a linear combination of incident fields in 12. Now from the definition of the solution operator G we further have The inverse scattering problem is now to reconstruct the support D when the near field operator N is given. We solve this inverse problem using the Factorization method. For our analysis, we need the following auxiliary operators: We define the where j ) j∈Z are bounded sequences, it is easy to check that W is bounded and bijective. Indeed, the determinant of the matrix in 18 equals 4w * + j w * − j (p 1,j ) 2 , which is nonzero thanks to 13. Now we define the operator (19) E : . E is also bounded due to the well-posedness of 20. We study in the next lemma some properties of operator H and its adjoint H * . .
The problem 20 is uniquely solvable for all wave numbers k, see [8]. Now assume that u solves 20, we have from the radiation condition that u|Γ ±h = m∈Zû Thus we find that . By similar computations we also obtain that The latter indentities show that H * is given by 21. Note that H * is bounded since W and E are bounded operators. We next show the compactness of H * . Standard elliptic regularity results imply that u is which allows us to obtain the compactness of E. This implies that H * is a compact operator and so is H.
To prove the injectivity of H, we show that the range of H * is dense. Since W is bijective, it is sufficient so show that the sequences (δ nm , 0, 0, 0) m∈Z , (0, δ nm , 0, 0) m∈Z , (0, 0, δ nm , 0) m∈Z and (0, 0, 0, δ nm ) m∈Z belong to the range of E (by definition, the Kronecker symbol δ nm equals one for n = m and zero otherwise). Now we show that (δ nm , 0, 0, 0) m∈Z belong to the range of E (the other cases are similar). To this end, we first note that (1,0) Then the Rayleigh sequences of Φ n and ϕ n are equal. Due to Lax-Milgram theorem, there exists a unique solution v ∈ H 1 (D, C 2 ) to the equation Thus Green's first identity shows that . Adding the two last equations we obtain Now we show a factorization of the near field operator N in the following theorem.
Theorem 2.2. Assume that P and Q satisfy Assumption 1 and recall the operator W from 18. The operator T : is the solution to 14. Then near field operator N satisfies solves 14 which can also be rewritten as Comparing to the definition of the operator E in 19 we can deduce that G(f, g) = −ET (f, g). Now due to the fact that N = GH we have Additionally we know from 21 that H * = −W E which completes the proof.

2.3.
Characterization of the periodic support. It is well-known for the Factorization method that a characterization of the periodic support is from an application of the Theorem A.1. We thus need to study properties of the middle operator T in the factorization of Theorem 2.2.
Lemma 2.3. Suppose that P and Q satisfy Assumption 1. Let T be the operator defined as in Theorem 3.2. Let T 0 : Proof. (a) We show the injectivity of T by assuming that T (f, g) = 0 or where u is a radiating variational solution to the homogeneous problem ∆u+k 2 u = 0 in Ω h and ∂ x2 u = T ± k (u) on Γ ±h . The latter problem has only the trivial solution which implies that u = 0 in Ω h . Thus, f = g = 0 or T is injective.

Therefore it follows
Im (b) From the definitions of T and T 0 we note that where u,ũ ∈ H 1 α (Ω h , C 2 ) are respectively the solutions of 14 and 24. Consider now the sequence (f n , g n ) which converges weakly to zero in L 2 (D, C 4 ) × L 2 (D, C 2 ). It is sufficient to prove that (T − T 0 )(f n , g n ) L 2 → 0. Due to the boundedness of the solution operator from L 2 (D, we imply that u n andũ n converge weakly to zero in H 1 α (Ω h , C 2 ). That means that w n := u n −ũ n converge strongly to zero in L 2 (Ω h , C 2 ) because of the compact embedding . Now considering 14, 24 for u n ,ũ n , respectively, with (f n , g n ) in the right hand side, and making a subtraction we have Choosing ψ = w n and taking the real part of the equation we obtain We know that w n L 2 → 0 as n tends to ∞, and thatũ n is a bounded sequence. The latter facts allow us to conclude that ∇w n L 2 converges to zero. Thus we have (T − T 0 )(f n , g n ) L 2 → 0.
(c) We return to 25 forũ instead of u. Then taking the real part implies that since −Re Γ ±hũ · T ± k (ũ) ds ≥ 0. Now assume that there is no such a constant γ > 0 for the statement of (c), then we can find a sequence {(f n , g n )} such that (f n , g n ) L 2 = 1 and Re T 0 (f n , g n ), (f n , g n ) L 2 → 0. Due to 26, we imply that (Re Q) −1 f + ∇ũ n → 0 in L 2 (D, C 4 ) and g +ũ n → 0 in L 2 (D, C 2 ), whereũ n denotes the solution of 24 for f and g replaced by f n and g n , respectively. Also, we have from the variational formulation forũ n that which is a contradiction to (f n , g n ) L 2 = 1. Therefore Re T 0 is coercive.
Next, by using a special test sequence, we show a characterization of a point z belonging to the support D. We need the α-quasi-periodic Green's function. It is well known that the function G k (x, y) given by (27) G is the α-quasi-periodic Green's function of the Helmholtz equation in two dimensions. Now we denote by Rg(X) range of the operator X. The following lemma is important for the main theorem of the Factorization method.
Lemma 2.4. Recall the operators H * and W from 21 and 18, respectively. For z ∈ Ω we denote by (r ± z,n ) n∈Z the Rayleigh sequences of the Green's function G(·, z). Then z belongs to D if and only if W (r + z,n , 0, r − z,n , 0) ∈ Rg(H * ). Proof. From 27 we see that which implies that the Rayleigh sequences of G k (·, z) are given (28) r ± z,n = i 4πβ n e −i(αnz1±βn(z2∓h)) .
Recall the operator E defined in 19 and the fact that H * = −W E. It is sufficient to prove that (r + z,n , 0, r − z,n , 0) n∈Z ∈ Rg(E). Now we assume that z is not in D and (r + z,n , 0, r − z,n , 0) n∈Z ∈ Rg(E). Then there exists u ∈ H 1 α (Ω h , C 2 ) solving 20 with the source function (f, g) ∈ L 2 (D, C 4 ) × L 2 (D, C 2 ) in the right hand side. Further u n = (r + z,n , 0, r − z,n , 0) for n ∈ Z. Since the Rayleigh sequences of (1, 0) G k (·, z) and u are equal, both functions coincide in Ω \ (−h, h). Due to the analyticity of u and (1, 0) G k (·, z) in Ω \ D and Ω \ {z}, respectively, and the analytic continuation we conclude that u = (1, 0) G k (·, z) in Ω \ (D ∪ {z}). However from [4] we know that G k (·, z) has a logarithmic singularity at z. This is hence a contradiction since u ∈ H 1 (B, C 2 ) for some neighborhood B of z but (1, 0) G k (·, z) / ∈ H 1 (B, C 2 ) due to the singularity at z.
To show that for z ∈ D there exists (f, g) ∈ L 2 (D, C 4 ) × L 2 (D, C 2 ) such that E(f, g) = (r + z,n , 0, r − z,n , 0), we just apply the proof of injectivity of H in Lemma 3.1.
Theorem 2.5. Suppose that Q and P satisfy Assumption 1, and thatr z,j = (r + z,j , 0, r − z,j , 0) j∈Z is the test sequence from 28. For j ∈ Z, denote by (λ n , ψ j,n ) n∈N an orthonormal eigensystem of (W N ) = |Re (W N )| − Im (W N ). Then a point z ∈ Ω belongs to D if and only if Proof. As we assumed in the theorem, let (λ 1/2 n , ψ j,n ) n∈N be an orthonormal eigensystem of (W N ) where the indices j, n in each subblock belong both to Z M1,M2 . Here,û ± (1,2),n are the two components of the Rayleigh coefficients of the scattered field defined in 8. Moreover, for l = 1, 2, the notation ( · ) We compute another eigenvalue decomposition of (W N M1,M2 ) = ΦΛΦ −1 with a diagonal matrix Λ that contains the 4(M 1 +M 2 +1) eigenvalues λ n of (W N M1,M2 ) and an orthogonal matrix Φ = (φ j,n ) containing the eigenvectors. Then (31) ( From the discretization of the near-field operator, we denote by (r z,j ) the discretization of test sequences (r z,j ) j∈Z ∈ 2 (Z) 4 . Then the criterion 59 is numerically exploited for imaging by plotting the function where A n (z) =

4(M1+M2+1) j=1
W (r z,j )φ j,n . If the series in 32 approximates the true value of the exact Picard series in 29, then P M1,M2 should be very small outside of D and considerably larger inside D.
To show the performance of the method with noisy data, we perturb our synthetic data by artificial noise. More precisely, we add a complex-valued noise matrix X containing random numbers that are uniformly distributed on (−1, 1) to the data matrix N M1,M2 . Denoting by δ the noise level, the noisy data matrix (N M1,M2 ) δ is then given by where · 2 is the matrix 2-norm. Obviously, for such noisy discrete data, the eigenvalue decomposition in 31 has to be replaced by a singular value decomposition, that we will not detail here. Following the traditional way of implementing sampling methods, we apply Tikhonov regularization [14], that is, instead of implementing 32 we consider where λ δ n are the singular values of the perturbed data matrix (W N M1,M2 ) 1/2 ,δ and A δ n (z) is the corresponding perturbation of the function A n from 32 (that has now to be defined using the left singular vectors of the perturbed data matrix). The parameter γ is chosen according to Morozov's discrepancy principle by approximately solving the non-linear scalar equation at each sampling point z.
(iii) A periodic structure of rectangles with piecewise constant material parameters, For all experiments below we set that the wave number k equals 3.5, α equals cos(π/6) in the first example, cos(π/2) in the second and cos(π/3) in the last example. We plot 3 periods of the obtained images to better illustrate the periodicity although we merely reconstruct one period of the structure. Figure 1 contains the reconstructions of the piecewise linear continuous grating defined in 34. It shows that the evanescent modes are important for obtaining useful  reconstructions. We further see that the reconstructions still remain reasonable when the measurements are perturbed by 2 percent artificial noise, that means δ equals 0.02. It also can be typically seen that the top of the grating is easier to recover than the "convex" part. Figure 2 is devoted to the reconstructions of the periodic structure of ellipsoids defined in 36. Reasonable reconstructions can still be obtained for the case of the material parameters varying within D. Again evanescent modes turn out to be necessary for the reconstructions. The results are almost the same for the case of periodic structures of rectangles 38 in Figure 3 when the material parameters are piecewise constant.
3.1. Problem setting. For the case of Maxwell's equations, we assume that the material parameters ε r , µ r and β are all 2π-periodic in x 1 and x 2 . We now eliminate the electric field E from 3 to obtain that (40) Similar to the two-dimensional case, we define that a function u : R 3 → C 3 is called α-quasi-periodic for α := (α 1 , α 2 , 0) ∈ R 3 if u(x 1 + 2πn 1 , x 2 + 2πn 2 , x 3 ) = e 2πi α·n u(x 1 , x 2 , x 3 ) for all n = (n 1 , n 2 , 0) ∈ Z 3 and for all x = (x 1 , x 2 , x 3 ) ∈ R 3 . Assume that the biperiodic structure is illuminated by α-quasi-periodic incident electric and magnetic fields E i and H i , respectively, satisfying  Simple examples for such α-quasi-periodic fields are certain plane waves that we introduce below. We now reformulate 40 in terms of the scattered field u, defined by u := H − H i . Since by construction, curl curl H i − k 2 H i = 0, subtracting the latter equation from 40 we obtain that where the contrasts q and p are given by q := 1 − ε −1 r , p := µ r − 1. As in the Helmholtz case, the α-quasi-periodic scattered field u admits a Rayleigh expansion radiation condition of the form Here α n and β n in three-dimensional setting are given by α n := (α 1,n , α 2,n , 0) := (α 1 + n 1 , α 2 + n 2 , 0) ∈ R 3 , and β n := k 2 − |α n | 2 , k 2 ≥ |α n | 2 , i |α n | 2 − k 2 , k 2 < |α n | 2 .
Again we also assume that β n = 0 for all n ∈ Z 2 . We denote by Ω := (0, 2π) 2 × R the unit cell and by D the support of the contrasts q and p in Ω, that is, Considering a more general source term on the right hand side of 41, we have the following direct problem: Given f, g ∈ L 2 (D, C 3 ), find a radiating solution u : Ω → C 3 in a suitable function space to Obviously, if u is a solution of 41 then u solves 43 for the right-hand side f = curl H i , g = H i . For a variational formulation of 43, we define H α,loc (curl, Ω) = {u ∈ H loc (curl, Ω) : u = U | Ω for some and with boundaries Γ ±h = (0, 2π) 2 × {±h}. The variational formulation to the problem 43 is to find a radiating solution u ∈ H α,loc (curl, Ω) such that for all ψ ∈ H α (curl, Ω) with compact support. Existence and uniqueness of this problem can be obtained for all but possibly a discrete set of exceptional positive wave numbers, see [1].
In the sequel we assume that the wave number k > 0 is chosen such that 44 is uniquely solvable for all f, g ∈ L 2 (D, C 3 ). Then we can define a linear and bounded solution operator G mapping the source (f, g) to the Rayleigh sequences of the first two components of u ∈ H α,loc (curl, Ω), solution to 44, (45) G : L 2 (D, C 3 ) 2 → 2 (Z 2 , C 4 ), (f, g) → (û + 1,j ,û − 1,j ,û + 2,j ,û − 2,j ) j∈Z 2 ∈ 2 (Z 2 , C 4 ). For l = 1, 2, the Rayleigh sequences û + l,j j∈Z 2 of the first two components u l of u = (u 1 , u 2 , u 3 ) above the biperiodic electromagnetic structure are explicitly given by (46) where l = 1, 2. Analogously, the Rayleigh sequences û − l,j j∈Z 2 for l = 1, 2 can be computed by replacing h by −h in 46. Now we introduce the notationb = (b 1 , b 2 , −b 3 ) for any vector b = (b 1 , b 2 , b 3 ) . We consider the following α-quasiperiodic plane waves where the polarizations p Notice that these incident plane waves have been used and discussed in [29]. Now as in 15 using sequences (a j ) j∈Z 2 = (a where the role of the coefficients is to reduce technical difficulties in some of the later computations. Given a sequence (a j ) ∈ 2 (Z 2 , C 4 ) we can build incident electromagnetic waves of the form H(a j ) (see 48) and measure the tangential components of the Rayleigh sequences (see 46) of the scattered wave. As in the Helmholtz case the near field operator N is defined by where u ∈ H α,loc (curl, Ω) is the radiating solution to 44 for the source (f, g) = H(a j ). From the definition of the solution operator G it is now clear that (at least formally, since H has not yet been proven to be bounded) We aim to reconstruct the support D when the near-field operator N is given as data.
As in the Helmholtz case, we also need the sequences (w * ± j ) j∈Z 2 , defined by We now introduce the three-dimensional version of the operators W and E in 18-19. Let be given by 18, where p (l) j and (w * ± j ) j∈Z 2 are defined in 47 and 50, respectively, and −4π is replaced by 8π 2 . We set E : where u is the radiating variational solution to curl 2 u−k 2 u = curl f + g in Ω. Since the latter problem is well-posed for all frequencies (see [40]), E is well-defined and bounded. Next is the result corresponding to Lemma 2.1 of the Helmholtz case. 3,j ) , j ∈ Z 2 , l = 1, 2, defined in 48, the operator H : 2 (Z 2 , C 4 ) → L 2 (D, C 3 ) 2 is compact and injective, and its adjoint H * : L 2 (D, C 3 ) 2 → 2 (Z 2 , C 4 ) satisfies Proof. (a) We show the injectivity of T by assuming that T φ = 0 for some φ = (f, g) ∈ X. By definition of T this means that v ∈ H α,loc (curl, Ω) is a radiating variational solution to the homogeneous problem curl 2 v − k 2 v = 0 in Ω. The latter problem has only the trivial solution (see [39]) which implies that v = 0 in Ω. Thus, f = g = 0 and hence T is injective.
For φ = (f, g) ∈ X, we set w 2 = f + curl v, w 1 = g + v, and w = (w 1 , w 2 ) . Then T (f, g) = w, and we have We consider a smooth cut-off function χ ∈ C ∞ (R) such that χ(t) = 1 for |t| < h and χ(t) = 0 for |t| > 2h. Then x → χ(x 3 )v(x) belongs to H α,loc (curl, Ω) with compact support in Ω 3h . Exploiting that v ∈ H α,loc (curl, Ω) is the radiating solution to 57 and that v solves the Helmholtz equation in Ω \ Ω h , we have from Green's theorems that Recall that, for |x 3 | > h, v satisfies the radiating Rayleigh condition. Therefore, if we replace h in the last equation by r ≥ h, all the terms corresponding to evanescent modes tend to zero as r tends to infinity. A straightforward computation shows that −Im Since Im (M w · w) = 0 in D, this implies that As in the Helmholtz case, to explicitly characterize the periodic scatterer we use special test sequences. For fixed y ∈ Ω, the α-quasi-periodic Green's tensor G k (x, y) ∈ C 3×3 for the Maxwell problem is given by is the α-quasi-periodic Green's function for the 3D Helmholtz problem. Here, the div of a matrix and the ∇ are meant to be taken columnwise and componentwise, respectively. Note that G k satisfies the Rayleigh expansion condition. We refer to [29] for a proof of the next lemma.
Lemma 3.4. Recall the operators H † and W from 53 and 18, respectively. For any z ∈ Ω and fixed nonzero p = (p 1 , p 2 , p 3 ) ∈ C 3 , we denote by (Ψ z,j ) j∈Z 2 = (Ψ + z,j ,Ψ − z,j ) j∈Z 2 ∈ 2 (Z 2 , C 4 ) the upper and lower Rayleigh coefficients of the first two components of Ψ z (x) := k 2 G k (x, z)p where (Ĝ ± k,j (z)) j∈Z 2 are the upper and lower Rayleigh sequences of the α-quasiperiodic Green's function G k (·, z) given bŷ Then z belongs to D if and only if W (Ψ j,z ) ∈ Rg(H † ).
Our main theorem on the Factorization method is now a simple corollary of the previous chapters combined with the range identity from Theorem A.1.
Theorem 3.5. Suppose that Assumption 2 holds true and that the direct scattering problem 44 is uniquely solvable. Denote by λ n , (φ n,j ) j∈Z 2 n∈N the orthonormal eigensystem of (W N ) = |Re (W N )|+Im (W N ) and by (Ψ z,j ) j∈Z 2 the test sequence from Lemma 3.4. Then a point z ∈ Ω belongs to the domain D if and only if  [30]. However, for the case of noisy data, we simply add a realvalued noise matrix X containing random numbers that are uniformly distributed on (0, 1) to the data matrix (W N M1,M2 ) 1/2 instead of N M1,M2 . This simple model of artificial noise enables reasonable reconstructions with 5 percent noise perturbation. All the following experiments rely on three different biperiodic structures that are defined in one period Ω = (−π, π) 2 × R in terms of the periodic support D as follows: (i) A biperiodic binary grating with constant material parameters, q(x) = 0.5, p(x) = 0.5, β(x) = 0.001 in D.
(ii) A biperiodic structure of cubes with variable material parameters, q(x) = (x 3 + 1)(sin(x 1 ) 2 sin(x 2 ) 2 + 0.5)/3 in D, p(x) = (x 3 + 1)(cos(x 1 ) 2 cos(x 2 ) 2 + 0.5)/3 in D, (iii) A biperiodic structure of cylinders with piecewise constant material parameters, For all experiments below we set that the wave number k equals 2.5, α 1 = α 2 = k/2. Similar to Helmholtz case we plot 3 × 3 periods of the obtained images to better illustrate the periodicity. We further note that all reconstructions have been smoothed using the command smooth3 in Matlab. The isovalue for plotting the images is chosen by hand such that the size of the reconstruction roughly remains the same for all examples. At this point, a-priori knowledge about the size of the volume of the structure is, at least implicitly, used. Figure 4 presents the reconstructions of the biperiodic binary grating defined in 60. As in the Helmholtz case we see that the reconstructions crucially depend on a sufficient number of near-field measurements. If one only measures the propagating modes of the scattered field, then the resulting images of the biperiodically aligned cubes are not useful. The reconstructions still provide reasonable information on the biperiodic grating when the measurements are perturbed by 5 percent artificial noise. Figures 5 and 6 respectively show the same imaging experiment for biperiodically aligned cubes and cylinders defined in 61 and 62. Although the materials parameters vary smoothly or jump within D, the conclusion from the first experiment remains essentially the same. Further, the stability of the reconstructions for the case of artificial noise is partly due to the three dimensional setting, and due the relatively high dimension of the measurement operator.
Appendix A. The range identity theorem. This appendix presents an abstract result on range identities which is necessary to characterize the support D, see for instance [23] for a complete proof. To state the result, we introduce real and imaginary part of a bounded linear operator T : X * → X as follows Re (T ) := 1 2 (T + T * ), Im (T ) := 1 2i (T − T * ).
Theorem A.1. Let X ⊂ U ⊂ X * be a Gelfand triple with Hilbert space U and reflexive Banach space X. Furthermore, let V be a second Hilbert space and F : V → V , H : V → X and T : X → X * be linear and bounded operators with We make the following assumptions: a) H is compact and injective. b) There exists t ∈ [0, 2π] such that Re (e it T ) has the form Re (e it T ) = T 0 + T 1 with some positive definite selfadjoint operator T 0 and some compact operator T 1 : X → X * . c) Im T is non positive on X, i.e., Im T φ, φ ≤ 0 for all φ ∈ X. Moreover, we assume that one of the two following conditions is fulfilled: d) T is injective and t from b) does not equal π/2 or 3π/2.