Inverse Diffraction Grating of Maxwell's Equations in Biperiodic Structures References and Links

Consider a time-harmonic electromagnetic plane wave incident on a perfectly conducting biperiodic surface (crossed grating). The diffrac-tion is modeled as a boundary value problem for the three-dimensional Maxwell equation. The surface is assumed to be a small and smooth deformation of a planar surface. In this paper, a novel approach is developed to solve the inverse diffraction grating problem in the near-field regime, which is to reconstruct the surface with resolution beyond Rayleigh's criterion. The method requires only a single incident field with one polarization, one frequency , and one incident direction, and is realized by using the fast Fourier transform. Numerical results show that the method is simple, efficient, and stable to reconstruct biperiodic surfaces with subwavelength resolution. An adaptive finite element method with perfectly matched absorbing layers for the wave scattering by periodic structures, " SIAM J. Integral equation methods in a quasi-periodic diffraction problem for the time-harmonic Maxwell's equations, " SIAM J. Analyzing diffraction gratings by a boundary integral equation Neumann-to-Dirichlet map method, " J. An inverse problem in diffractive optics: conditional stability, " In-A computational inverse diffraction grating problem, " J. Numerical solution of an inverse diffraction grating problem from phasless data, " J. A two-step algorithm for the reconstruction of perfectly reflecting periodic profiles , Grating profile reconstruction based on finite elements and optimization techniques, " SIAM J. Operator expansions and constrained quadratic optimization for interface reconstruction: Impenetrable periodic acoustic media, " Wave Motion, to appear. 20. D. Dobson, " Optimal design of periodic antireflective structures for the Helmholtz equation, " Eur. Numerical solution of optimal design problems for binary gratings, " J. Comput. An improved operator expansion algorithm for direct and inverse scattering computations, " Waves Random Media 9, 441–457 (1999). 26. J. A. DeSanto and R. J. Wombell, " The reconstruction of shallow rough-surface profiles from scattered field data, " Inverse Probl. An adaptive edge element method with perfectly matched absorbing layers for wave scattering by biperiodic structures, " Math. Unique determination of periodic polyhedral structures by scattered electromagnetic fields, " Trans. An inverse problem for scattering by a doubly periodic structure, " Trans. An inverse electromagnetic scattering problem for a bi-periodic inhomogeneous layer on a perfectly conducting plate, " Appl. The linear sampling method for inverse electromagnetic scattering by a partially coated bi-periodic structures, " Math. Inverse electromagnetic scattering problems by a doubly periodic structure, " Math. Appl. The factorization method for …


Introduction
Consider the scattering of a time-harmonic electromagnetic plane wave by a biperiodic structure, known as crossed grating or two-dimensional grating.Scattering theory in periodic structures has many applications in micro-optics including the design and fabrication of optical elements such as corrective lenses, anti-reflective interfaces, beam splitters, and sensors.Depending on the direction and polarization of the incident plane wave, the governing mathematical model can be simplified from the full three-dimensional Maxwell equations to two fundamental polarizations: the transverse electric polarization and the transverse magnetic polarization, known as linear grating or one-dimensional grating.In both polarizations, the scalar components of electromagnetic waves satisfy the two-dimensional Helmholtz equation.We refer to the monograph [1] for a good introduction to the problems of electromagnetic diffraction.
Recently, the scattering problems in periodic structures have been studied extensively on both mathematical and numerical aspects.We refer to [2] and references therein for the mathematical studies of existence and uniqueness of the diffraction grating problems.Numerical methods can be found in [3][4][5] for either an integral equation approach or a variational approach.A comprehensive review can be found in [6] on diffractive optics technology and its mathematical modeling as well as computational methods.
We consider an inverse problem, which is to reconstruct the grating surface from a measured data field at a constant height above the surface.The mathematical questions on uniqueness and stability of the inverse problem have been studied by many researchers for the one-dimensional grating [7][8][9][10][11].Computationally, a number of methods have been developed for the reconstruction of perfectly conducting grating surfaces in the transverse electric polarization [12][13][14][15][16][17][18][19].We refer to [20][21][22][23] for related optimal design problems in diffractive optics, which are to design grating structures to obtain some specified diffraction patterns, and [24][25][26][27] for general inverse surface scattering problems.These work addressed conventional far-field imaging, where the scattering data is taken at distances which are greater than the wavelength of the incident field.The role of evanescent wave components were ignored and the resolution of reconstructions was limited by Rayleigh's criterion, approximately half of the incident wavelength, also known as the diffraction limit.We refer to [28][29][30][31] for the existence, uniqueness, and numerical approximations of solutions for the direct two-dimensional grating problems.Mathematical studies can be found in [32][33][34][35][36][37] on the uniqueness results for detecting biperiodic grating surfaces.In contrast, numerical results are very rare due to nonlinearity and ill-posedness of the inverse problem plus complexity of the three-dimensional Maxwell equations.A qualitative imaging method can be found in [38] for solving an inverse scattering problem from penetrable biperiodic structures.More recent reviews can be found in [39,40] on the direct and inverse scattering problems in periodic media.
In this work, we develop an efficient and stable computational method to solve the inverse problem.The grating surface is assumed to be a small and smooth deformation of a planar surface.Based on the transformed field expansion, the method reduces the boundary value problem into a successive sequence of two-point boundary value problems.For transformed field expansion and related boundary perturbation method, we refer to [41][42][43][44][45][46][47] for solving the direct and inverse diffraction grating problems.An explicit reconstruction formula is derived for the linearized inverse problem by dropping higher order terms in the expansions.A spectral cutoff regularization is adopted to suppress the exponential growth of the noise in the evanescent wave components, which carry high spatial frequency of the scattering surface and contribute to the super resolution in the near-field regime.The method requires only a single incident field with one polarization, one frequency, and one incident direction, and is realized by using the fast Fourier transform.The numerical results are computed by using synthetic scattering data provided by an adaptive edge element method with a perfectly matched absorbing layer [29].
Two numerical examples, one smooth surface and one non-smooth surface, are presented to demonstrate the effectivenss of the proposed method.The influence is carefully investigated on the reconstructions for such parameters as surface deformation, measurement distance, and noise level.The numerical results show that the method is simple, efficient, stable to reconstruct biperiodic grating surfaces with subwavelength resolution.
This paper significantly extends our previous work on near-field imaging of one-dimensional surfaces [48][49][50], where the two-dimensional scalar Helmholtz equation was considered, to two-dimensional grating surfaces.Apparently, the techniques differ greatly from existing work because of the complicated model problem of three-dimensional Maxwell's equations.To the best of our knowledge, we develop the first quantitative method for solving the inverse diffraction grating problem of Maxwell's equations in biperiodic structures and provide numerical examples of reconstruction with super resolved resolution.We point out a closely related work on the inverse surface scattering in near-field imaging [51], where the scattering surface is assumed to be a small and local perturbation of a planar surface.Other related work may be found in [52,53] for solving an inverse medium scattering problem in near-field optical imaging and in [54,55] for resolution and stability analysis of conductivity imaging.

Model problem
In this section, we define some notations and introduce a boundary value problem for the diffraction by a biperiodic structure.

Maxwell's equations
Let us first specify the diffraction grating problem geometry.Denote (ρ, z) ∈ R 3 , where ρ = (x, y) ∈ R 2 .As seen in Fig. 1, the problem may be restricted to a single period of Λ = (Λ 1 , Λ 2 ) in ρ due to the periodicity of the structure.Let the surface in one period be described by The grating surface function φ is assumed to be in the form φ (ρ) = εψ(ρ), where ψ ∈ C 2 (R 2 ) is a biperiodic function with period Λ and is called the grating profile, and ε is sufficiently small and is called the surface deformation parameter.

Denote by
the space above S, which is filled with some homogeneous medium characterized by a positive constant wavenumber κ.
be the domain bounded below by S and bounded above by the plane surface Let (E inc , H inc ) be the incoming plane waves that are incident upon the grating surface from above, where ( Here α = (α 1 , α 2 ), α 1 = sin θ 1 cos θ 2 , α 2 = sin θ 1 sin θ 2 , and β = cos θ 1 , where θ 1 and θ 2 are the latitudinal and longitudinal incident angles, which satisfy 0 ≤ θ 1 < π/2, 0 ≤ θ 2 < 2π.Denote by d = (α 1 , α 2 , −β ) the unit propagation direction vector.The unit polarization vectors t = (t 1 ,t 2 ,t 3 ), s = (s 1 , s 2 , s 3 ) satisfy t • d = 0 and s = d × t, which gives explicitly For normal incident, i.e., θ 1 = 0, we have α 1 = 0, α 2 = 0, β = 1, and Hence we get from |t| = |s| = 1 that t 2 1 + t 2 2 = 1,t 3 = 0.For the sake of simplicity, we focus on the case of normal incidence from now on since our method requires only a single incident wave for solving the inverse problem.We mention that the method works for general non-normal incidence with obvious modifications. Denote . Under the normal incidence, the incoming plane waves (2) reduce to It can be verified that the incident electromagnetic waves satisfy the three-dimensional timeharmonic Maxwell equation: The scattering of time-harmonic electromagnetic waves follows Maxwell's equations in the space above the grating surface: where E is the total electric field and H is the total magnetic field.Due to the homogeneous medium, the electromagnetic fields satisfy the divergence free condition: We consider the perfect electric conductor condition: where ν S = (ν 1 , ν 2 , ν 3 ) ∈ R 3 is the unit normal vector on S, given explicitly as Here φ x = ∂ x φ (x, y) and φ y = ∂ y φ (x, y) are the partial derivatives.

Transparent boundary condition
To reduced the diffraction grating problem from an unbounded domain Ω S into a bounded domain Ω, a transparent boundary condition needs to be imposed on Γ.
For any tangential vector u(ρ, h) = (u 1 (ρ, h), u 2 (ρ, h), 0) on Γ, where u j is a biperiodic function in ρ with period Λ, we define a boundary operator T : where v j is also a biperiodic function in ρ with the same period Λ.Here u j and v j have the following Fourier expansions and the Fourier coefficients u jn and v jn satisfy Using the boundary operator ( 9), we may derive the transparent boundary condition [29]: where Recalling the incident fields (3) and using the boundary operator (9), we have explicitly that

Reduced model problem
Taking curl on both sides of (5), we may eliminate the magnetic field and deduce a boundary value problem for the electric field: It is convenient introducing an equivalent scalar form of the problem (11) in order to apply the transformed field expansion.Denote E = (E 1 , E 2 , E 3 ).Noting the divergence condition (6), we may reduce Maxwell's equations to the Helmholtz equation for E j : The divergence free condition ( 6) can be explicitly written as Substituting ( 8) into ( 7) yields The transparent boundary condition (10) becomes where the Fourier coefficients of H 1 and H 2 are given by Here E 1n (h) and E 2n (h) are the Fourier coefficients of E 1 (ρ, h) and E 2 (ρ, h), respectively.Given the incident field, the direct problem is to solve the boundary value problem ( 12)-( 15) for the known surface function φ (ρ).The inverse problem is to reconstruct the function φ (ρ) from the tangential trace of the total field measured at , where δ is the noise level.In particular, we are interested in the inverse problem in the near-field regime where the measurement distance h is much smaller than the wavelength λ = 2π/κ.

Transformed field expansion
In this section, we introduce a transformed field expansion to find a power series solution for the direct problem ( 12)-(15).

Change of variables
The transformed field expansion method begins with the change of variables: which maps the domain Ω into a rectangular slab Introduce a new function Ẽ = ( Ẽ1 , Ẽ2 , Ẽ3 ) and let Ẽ j ( x, ỹ, z) = E j (x, y, z) under the transformation.After tedious but straightforward calculations, it can be verified from ( 12) that the total electric field, upon dropping the tilde, satisfies the equation where x + φ 2 y ) .The divergence free condition (13) becomes The perfect electric conductor condition ( 14) is The transparent boundary condition (15) on z = h reduces to

Power series solution
Consider a formal expansion of E j in a power series of ε: Substituting φ = εψ into c j and inserting (20) into ( 16), we may derive where Here ψ x = ∂ x ψ(x, y) and ψ y = ∂ y ψ(x, y) are the partial derivatives.Substituting (20) into the divergence free condition (17) yields where .
The perfect electric conductor boundary condition ( 18) can be written as where .
Inserting (20) into the transparent boundary condition (19), we get where , and the Fourier coefficients of Here 1 (ρ, h), respectively.The divergence free condition (22) on z = h reduces to where v .
Clearly, the problem ( 21)-( 26) for j , which depend only on previous two terms of . Thus, the transformed problem ( 21)- (26) indeed can be solved efficiently in a recursive manner starting from k = 0.

Zeroth order term
Recalling the recurrence relation (21) and letting k = 0, we have The divergence free condition (22) reduces to The perfect electric conductor boundary condition ( 23) is Using ( 28) and ( 29), we have The transparent boundary condition (25) becomes In addition, the divergence free condition (28) gives one more boundary condition: Since E (0) j (ρ, z) and f j are periodic functions of ρ, they have the Fourier expansions where f j0 = −2iκt j exp(−iκh) and f jn = 0 for n = 0. Substituting (33) into (27), we derive a second order ordinary differential equation for the Fourier coefficient E (0) jn (z): Similarly, we have the boundary conditions at z = 0: and the boundary conditions at z = h: Simple calculations yield that the solution of the two-point bounary value problem (34)-( 36) is

Reconstruction formula
Assume that the noisy data takes the form where E j (ρ, h), j = 1, 2 is the exact data and δ is the noise level.
Evaluating the power series (20) at z = h and replacing E j (ρ, h) with E δ j (ρ, h), we have Rearranging (45), and dropping O(ε 2 ) and O(δ ) yield εE (1) which is the linearization of the nonlinear inverse problem and enables us to find an explicit reconstruction formula for the linearized inverse problem.Noting φ = εψ and thus φ n = εψ n , where φ n is the Fourier coefficient of φ .Substituting ( 44) into (46), we deduce that where E δ jn (h) is the Fourier coefficient of the noisy data E δ j (x, h) and Here δ 0n the Kronecker's delta function.
It follows from the definition of β n and ( 47) that it is well-posed to reconstruct those Fourier coefficients φ n with |α n | < κ, since the small variations of the measured data will not be amplified and lead to large errors in the reconstruction, but the resolution of the reconstructed function φ is restricted by the given wavenumber κ.In contrast, it is severely ill-posed to reconstruct those Fourier coefficients φ n with |α n | > κ, since the small variations in the data will be exponentially enlarged and lead to huge errors in the reconstruction, but they contribute to the super resolution of the reconstructed function φ .To obtain a stable and super-resolved reconstruction, we may adopt a regularization to suppress the exponential growth of the reconstruction errors.
For fixed h, the cut-off wavenumber κ c is chosen in such a way that exp (κ 2 c − κ 2 ) 1/2 h = SNR, which implies that the spatial frequency will be cut-off for those below the noise level.More explicitly, we have which indicates κ c > κ as long as SNR > 0 and super resolution may be achieved.
Taking into account the frequency cut-off, we have a regularized reconstruction formulation where the characteristic function Once φ n are computed, the grating surface function can be approximated by Hence, the method requires only two fast Fourier transforms: one is done for the data to obtain E δ jn (h) and another is done to obtain the approximated function φ .

Numerical experiment
In this section, we discuss the algorithmic implementation for the direct and inverse problems, present two numerical examples to illustrate the effectiveness of the proposed method, and examine influence of the parameters ε, h, and δ on the reconstruction results.As seen in Fig. 2, we consider two types of grating profiles: one is a smooth function with finitely many Fourier modes and another is a non-smooth function with infinitely many Fourier modes.Although the method requires that the grating profile function ψ(ρ) is C 2 (R 2 ), it is still applicable to non-smooth functions numerically.The first-order Nédélec edge element is used for solving the direct problem and obtaining the synthetic scattering data.Uniaxial perfect matched layer (PML) boundary condition is imposed on z direction so that no artificial wave reflection occurs to ruin the wave field inside the domain.Adaptive refinement technique [29] is used to achieve the solution having a specified accuracy in an optimal fashion.Our implementation is based on parallel hierarchical grid (PHG) [56], which is a toolbox for developing parallel adaptive finite element programs on unstructured tetrahedral meshes and it is under active development at the State Key Laboratory of Scientific and Engineering Computing.The finite elements implemented in PHG are the Largrange element, hierarchical H 1 and H(curl) element.In order to generate the tetrahedral mesh with a biperiodic structure, we generate a non-uniform hexahedral mesh firstly and divide each hexahedron into six tetrahedrons.The linear system resulted from finite element discretization is solve by the multifrontal massively parallel sparse direct solver [57,58].
In the following two examples, the incident wave is taken as E inc = (1, 0, 0) exp(−iκz), i.e., t 1 = 1 and t 2 = t 3 = 0, and only the first component of the electric field, E 1 (ρ, h), needs to be measured.The wavenumber is κ = 2π, which corresponds to the wavenlength λ = 1.Define by R the unit rectangular domain, i.e., R = [0, 1.0λ ] × [0, 1.0λ ].The computational domain is R × [φ , 1.0λ ] with the PML region R × [0.5λ , 1.0λ ].The scattering data E 1 (ρ, h) is obtained by interpolation into the uniform 256 × 256 grid points on the measurment plane z = h.In all the figures, the plots are rescaled with respect to the wavelength λ to clearly show the relative size, and the meshes are done in 32 × 32 instead of 256 × 256 grid points in order to reduce the display sizes.To test the stability of the method, some relative random noise is added to the scattering data, i.e., the scattering data takes the form where rand stands for uniformly distributed random numbers in [−1, 1].The relative L 2 (R) error is defined by where φ is the exact surface function and φ δ ,ε is the reconstructed surface function.Example 1.This example illustrates the reconstruction results of a smooth grating profile with finitely many Fourier modes.The exact grating surface function is given by φ (ρ) = εψ(ρ), where the grating profile function ψ(x, y) = 0.6 sin(2πx) sin(2πy) + sin(4πx) sin(4πy). ( First, consider the surface deviation parameter ε.The measurement is taken at h = 0.4λ and no additional random noise is added to the scattering data, i.e., δ = 0.This test is to investigate the influence of surface deviation parameter on the reconstructions.In (46), higher order terms of ε are dropped in the power series to linearize the inverse problem and to obtain the explicit reconstruction formula.As expected, the smaller the surface deviation ε is, the more accurate is the approximation of the linearized model to the original nonlinear model problem.Table 1 shows the relative L 2 (R) error of the reconstructions with four different surface deformation parameter ε = 0.2λ , 0.1λ , 0.05λ , 0.025λ for fixed measurement distance h = 0.4λ .It is clear to note that the error decreases from 85.0% to 9.86% as ε decreases from 0.2λ to 0.025λ .Next, consider the noise level δ and the measurement distance h.In practice, the scattering data always contains certain level of noise.To test the stability and super resolving capability of the method, we add an amount of 5% random noise to the scattering data.Table 2 reports the relative L 2 (R) error of the reconstructions with four different measurement distance h = 0.4λ , 0.3λ , 0.2λ , 0.1λ for fixed ε = 0.025λ .Comparing the results for the same ε = 0.025λ and h = 0.4λ in Tables 1 and 2, we can see that the relative error increases dramatically from 9.86% by using noise free data to 86.3% by using 5% noise data.The reason is that a smaller cut-off should be chosen to suppress the expotentially increasing noise in the data and thus the Fourier modes of the exact grating surface function can not be recovered for those higher than the cutoff frequency, which leads to a large error and poor resolution in the reconstruction.A smaller measurement distance is desirable in order to have a large cut-off frequency, which enhances the resolution and reduces the error.As can be seen in Table 2, the reconstruction error decreases from 86.3% by using h = 0.4λ to as low as 12.0% by using h = 0.1λ even for 5% noise data.Figure 3 plots the reconstructed surfaces by using h = 0.4λ , 0.3λ , 0.2λ , 0.1λ .Comparing the exact surface profile in Fig. 2(a) and the reconstructed surface in Fig. 3(d), we can see that the reconstruction is nearly perfect and the difference is really minor by carefully checking the contour plots.
Consider the influence of ε by using noise-free data.The measurement is taken at h = 0.2λ .Table 3 presents the relative L 2 (R) error of the reconstructions with four different surface deformation parameter ε = 0.1λ , 0.05λ , 0.025λ , 0.0125λ .The error decreases from 72.9% to 15.2% as ε decreases from 0.1λ to 0.0125λ .Based on these results, the following observation can be made: a smaller deformation parameter ε yields a better reconstruction; smaller ε and h are required in order to obtain comparable error with that in Table 2 for Example 1 due to the non-smooth nature of the grating surface function of Example 2.  Next, consider the influence of the noise level δ and the measurement distance h.An amount of 5% random noise is added to the scattering data.Table 4 reports the relative L 2 (R) error of the reconstructions with four different measurement distance h = 0.2λ , 0.1λ , 0.05λ , 0.025λ for fixed ε = 0.0125λ .Comparing the results for the same ε = 0.0125λ and h = 0.2λ in Tables 3  and 4, we can see that the relative error is more than doubled from 15.2% by using noise-free data to 39.9% by using 5% noise data.Again, the reason is that a smaller cut-off is chosen to suppress the expotentially increasing noise in the data and thus higher Fourier modes of the exact grating surface function can not be recovered.A smaller measurement distance helps to enhance the resolution and reduce the error.In Table 4, the reconstruction error decreases from 39.9% by using h = 0.2λ to as low as 11.9% by using h = 0.025λ .Figure 4 shows the reconstructed surfaces by using h = 0.2λ , 0.1λ , 0.05λ , 0.025λ .Comparing the exact surface profile in Fig. 2(b) and the reconstructed surface in Fig. 4(d), we can see that a good reconstruction can still be possible when using a small measurement distance.

Conclusion
We have presented a simple, stable, and effective computational method for reconstructing biperiodic grating surfaces and achieved subwavelength resolution.Using the transformed field and Fourier expansions, we deduced a power series solution for the direct problem.By dropping higher order terms in power series, we linearized the nonlinear inverse problem and obtained an explicit reconstruction formula, which was implemented by using the fast Fourier transform.We considered two examples, one of which has finite Fourier modes and another one has infinite Fourier modes, and investigated how the parameters influence the reconstructions.The results show that super resolution may be achieved by using small measurement distance.As for future work, we will study a more complicated transmission problem where the surface is penetrable, and investigate the mathematical issues such as uniqueness, stability, and resolution.

Table 1 .
Example 1: Relative error of the reconstructions by using different ε with h = 0.4λ and δ = 0.0.

Table 2 .
Example 1: Relative error of the reconstructions by using different h with ε = 0.025λ and δ = 5%.

Table 4 .
Example 2:Relative error of the reconstructions by using different h with ε = 0.0125λ and δ = 5%.