Inverse scattering of periodic surfaces with a superlens

We propose a scheme for imaging periodic surfaces using a superlens. By employing an inverse scattering model and the transformed field expansion method, we derive an approximate reconstruction formula for the surface profile, assuming small amplitude. This formula suggests that unlimited resolution can be achieved for the linearized inverse problem with perfectly matched parameters. Our method requires only a single incident wave at a fixed frequency and can be efficiently implemented using fast Fourier transform. Through numerical experiments, we demonstrate that our method achieves resolution significantly surpassing the resolution limit for both smooth and non-smooth surface profiles with either perfect or marginally imperfect parameters.


Introduction
Resolution is a crucial factor in any wave imaging system.In traditional imaging systems, such as optical microscopes, the resolution is limited by approximately half of the wavelength, a rule known as the Rayleigh criterion or the Abbe diffraction limit.Decreasing the wavelength may not always be feasible, especially in source imaging problems, and high-frequency incident waves may cause side effects such as sample damage.Therefore, it is desirable to surpass the diffraction limit, enabling high resolution to be achieved with low-frequency waves.
Typically, the diffraction limit arises from the inability to capture evanescent waves, which convey fine details of the imaged object.However, these waves are confined near the object's surface, and their amplitude diminishes rapidly as the distance from the object increases.Consequently, a viable strategy for overcoming the diffraction limit involves recovering evanescent waves through measurements taken in close proximity to the object.This approach has been employed in near-field optical microscopy techniques, such as scanning near-field optical microscopy [1,2,3,4] and photon scanning tunneling microscopy [5].
Optical microscopy techniques generate images of samples directly on the measuring device, often resulting in the loss of depth information or its recovery at a lower resolution.To reconstruct the full profile of a sample, it is desirable to solve inverse scattering problems based on the underlying wave equations.For infinite surfaces, the inverse scattering problem involves reconstructing the surface profile using scattered waves measured on a surface.One prevalent approach is to construct an objective functional derived from the discrepancy between measured and computed data for a given surface profile.The surface profile is then identified as a minimizer of the objective functional through optimization algorithms [6,7,8,9].While these methods offer versatility in terms of objective functionals and optimization algorithms, they are computationally expensive and prone to issues related to local minima.
An alternative approach, known as sampling methods, involves designing an indicator function with high contrast values across the surface profile.The surface can then be visually identified by plotting the indicator function within a sampling region [10,11,12,13].While sampling methods offer the advantages of speed and minimal reliance on prior knowledge, they typically necessitate a substantial volume of measurement data, which may not be readily available in practical scenarios.Furthermore, these methods are intrinsically qualitative, leading to potential limitations in the quantitative precision of the resulting images.
In specific practical situations, surfaces may display a low profile with respect to period and wavelength.In such instances, techniques based on the transformed field expansion (TFE) can be utilized.The TFE was devised in [14,15] to address direct rough surface scattering problems.It was initially employed in [16] to solve inverse scattering problems on infinite rough surfaces and has since been expanded to various contexts, encompassing periodic surfaces [17], obstacles [18,19], interior cavities [20], penetrable surfaces [21], elastic waves [22,23], and three-dimensional problems [24,25,26].These methods leverage the TFE and linearization to derive explicit and accessible reconstruction formulas.They are not only straightforward but also effective and efficient, necessitating merely a single incident field.Another advantage of these techniques is their capacity to attain convergence results and error estimates [27,28], attributable to the explicit and accessible reconstruction formulas.
In this paper, we re-examine the near-field inverse scattering problem for periodic surfaces, aiming to further improve imaging resolution by integrating a superlens into the model.An ideal superlens is a slab of material characterized by a relative permittivity ε = −1 and relative permeability µ = −1.The notion of materials with negative permittivity and permeability was initially proposed by Veselago in [29], who contended that such materials could exhibit atypical optical properties, including a negative refractive index.Subsequently, Pendry introduced the concept of employing a slab of this material as an imaging lens and illustrated the potential for eliminating the resolution limit, coining the term "superlens" in [30].
Following this, a multitude of numerical and experimental demonstrations of superlenses for electromagnetic waves have been showcased in the literature [31,32,33].In more recent developments, the concept has been expanded to include acoustic waves and elastic waves.Nevertheless, the general imaging approach in these studies remains akin to that of a traditional microscope, where light from a source passes through the lens, ultimately forming a two-dimensional projection of the imaged profile on the image plane.In contrast, our imaging scheme relies on inverse scattering, providing a comprehensive reconstruction of the entire profile.
In this study, we propose an imaging model that comprises an impenetrable periodic surface, a slab of negative index material situated above the surface, a measurement surface atop the slab, and a single incident field.By employing the low-profile assumption and the TFE method, we derive an explicit reconstruction formula.We illustrate that, in the ideal case, unlimited resolution can be achieved, while high resolution can be obtained with slightly non-ideal parameters.To the best of our knowledge, this research represents the first instance of integrating negative index material into inverse scattering problems, thereby unlocking the full-profile imaging resolution enhancement potential offered by the superlens.
This work can be viewed as an extension of our previous study [34], wherein we attained enhanced resolution utilizing a slab with a high refractive index.The enhancement factor in [34] is approximately proportional to the slab's refractive index; however, it tends toward infinity in this paper as the parameters approach ideal values.Furthermore, we apply the TFE differently in this work by treating the entire domain as a single boundary value problem.As a result, the reconstruction formula becomes substantially simpler, facilitating preliminary resolution analysis.Moreover, we commence with the most general settings for the slab's material parameters to derive a comprehensive formula, which encompasses the no-slab, high-index slab, and negative-index slab as special cases.
The remainder of the paper is organized as follows.In Section 2, we establish the physical model and derive the boundary value problem for the entire domain using interface conditions and a transparent boundary condition.Section 3 presents the TFE, which reduces the problem into a recursive system of one-dimensional boundary value problems.In Section 4, we obtain the analytical solutions for the leading and linear terms.By solving the linearized inverse problem, we derive the reconstruction formula and conduct elementary resolution analysis in Section 5. Numerical experiments are presented in Section 6 to demonstrate the effectiveness of the proposed imaging scheme.Finally, Section 7 offers a conclusion and discusses future research directions.

The model problem
Consider the electromagnetic scattering by a periodic surface defined as where f is a periodic function with period Λ.We assume that f can be expressed as where 0 < δ 1 is a small constant, referred to as the surface deformation parameter.The surface Γ f is assumed to be perfectly electrically conducting (PEC), and the half space above it is assumed to be a vacuum.
We introduce a slab of double negative index material above the scattering surface, with its two surfaces defined by Γ a = {(x, a) : x ∈ R} and Γ b = {(x, b) : x ∈ R}.The domain bounded between Γ f and Γ a is denoted by Ω = {(x, y) ∈ R 2 : f (x) < y < a}, and the domain between Γ a and Γ b is represented by D = {(x, y) ∈ R 2 : a < y < b}.The problem geometry is illustrated in Fig. 1.The propagation of electromagnetic waves is described by the time-harmonic Maxwell's equations: where E and H are the electric and magnetic fields, respectively, ω represents the angular frequency, ε and µ are the permittivity and permeability, respectively.In this work, we consider the transverse electric (TE) mode, assuming all variables are independent of the z coordinate, and . Then, Eq. ( 2) reduces to The relative permittivity and permeability are denoted by ε r = ε/ε 0 and µ r = µ/µ 0 , respectively, where ε 0 and µ 0 are the permittivity and permeability of vacuum, respectively.Define the free-space wavenumber κ = ω √ ε 0 µ 0 , then Eq. ( 3) can be rewritten as where η = κ √ ε r µ r .The positive x-axis is taken as the branch cut for any square root in this paper.The term √ ε r µ r is commonly referred to as the refractive index.
A plane wave with an incident angle θ ∈ (−π/2, π/2) is described by u in (x, y) = e i(αx−βy) , where α = κ sin θ and β = κ cos θ.However, we will focus on normal incidence, i.e., θ = 0.In this case, the incident field simplifies to where both u in and u are periodic in x with a period of Λ.
Given a scattering surface and incident field, the direct scattering problem aims to determine the total field u(x, y).The near-field imaging problem we consider in this paper is an inverse scattering problem, seeking to determine the scattering surface function f (x) from the measured data of the total field at Γ b .
The incident field given by Eq. ( 4) clearly satisfies the Helmholtz equation: Let u sc = u − u in denote the scattered field, which satisfies the same Helmholtz equation: Assuming that u sc consists only of upward propagating and evanescent waves, we can obtain the Rayleigh expansion from Eq. ( 5) as follows: sc (b) denotes the Fourier coefficient of u sc (x, b), and Taking ∂ y in Eq. ( 6) and evaluating at Γ b , we get where the boundary operator T is defined by and the "+" sign indicates taking partial derivative or limit from above.It is straightforward to verify where Adding Eq. ( 9) to Eq. ( 8) yields the transparent boundary condition: From the interface conditions ν where the "−" sign indicates taking partial derivative or limit from below.Combining Eqs.(11) and ( 12) leads to Similarly, we have For simplicity of notations, we will denote the relative permittivity and permeability by and µ, respectively, in the following.To summarize, we obtain the following boundary-interface value problem for the total field u in D ∪ Ω: 3

Transformed field expansion
We apply the TFE method, along with power series and Fourier series expansions, to reduce the governing equations ( 13) to a recursive system of ordinary differential equations.This approach has been previously employed in [16] for the purpose of near-field imaging of infinite rough surfaces.
Consider the change of variables which transforms the domain Ω to the rectangle Ω 0 = (0, Λ) × (0, a), the surface Γ a to itself, the surface Γ f to the plane Γ 0 = {(x, 0) : x ∈ R}, and the domain D to itself.It is important to note that the change of variables is defined in the whole computational domain, as opposed to focusing solely on the domain Ω as in [34].This approach allows us to derive significantly simpler solution forms and reconstruction formulas.
Applying the chain rule leads to the following differentiation rules in Ω: Substituting the differentiation rules to the boundary-interface value problem Eqs. ( 13) and dropping the tilde over the variables for simplicity of notations, we obtain the transformed problem: where Based on the smallness assumption (1), we consider the power series expansion Substituting Eqs. ( 1) and ( 15) into Eqs.( 14) yields the recursive system of equations where where δ i,j is the Kronecker delta.The differential operators D 1 , D 2 , T are given by and ρ is given by Eq. (10).
As the variables u m , v m , τ m , and ρ m are periodic in the x direction with a period of Λ, the Fourier series expansion can be applied to the governing equations ( 16) to obtain a recursive system of boundary-interface value problems: where β n are given in Eq. ( 7) and 4 Analytical solutions Solving Eqs. ( 18) by the variation of parameters leads to the general solution where a m,n , b m,n , c m,n , d m,n are coefficients to be determined, and the integral kernel Applying the boundary and interface conditions (18a), (18c), and (18e), we obtain the linear system where A direct calculation shows that the determinant of the matrix in Eq. ( 21) is given by
Solving Eq. ( 21) by Cramer's rule yields the coefficients Substituting Eqs.(24) into Eq.( 20), we obtain the leading term of the solution It should be noted that u 0 is independent of x and corresponds to the total field when the scattering surface is flat, i.e., when the surface height deviation δ is zero.

First order term
For m = 1 we have ρ 1 = 0 and hence r 1,n = 0. Solving Eq. ( 21) by Cramer's rule for a 1,n , b 1,n yields Substituting Eq. ( 22) with m = 1 into the above equations and using the difference formula for the sine function, we obtain where By substituting Eq. ( 25) into Eq.( 17) and taking the Fourier transform, we obtain the following expressions for v By substituting Eqs.(28) into Eq.( 27) and using the identity κ 2 = α 2 n + β 2 n , we derive the expression for ψ n through a lengthy but straightforward calculation: Substituting Eq. ( 29) into Eq.( 26) yields the following expressions for a 1,n and b 1,n : Finally, substituting Eqs.(30) into Eq.( 20) provides the Fourier coefficients of the linear term: Evaluating Eq. ( 31) at y = b results in the fundamental identity of this paper: where the scaling factor is given by: Remark 1 For the sake of conciseness, we omit specific cases where either β n = 0 or γ n = 0 from this paper.However, it is important to note that these cases can be addressed separately if necessary.

Reconstruction Formula
By retaining only the zeroth and first order terms in the power series expansion given by Eq. ( 15), we arrive at the approximation: Evaluating Eq. ( 34) at y = b and applying the Fourier transform with respect to x gives Substituting Eqs. ( 32) and (1) into Eq.( 35), we deduce This implies that the Fourier coefficients of the scattering surface function f (x) can be computed directly from the Fourier coefficients of the total field u(x, b).The scattering surface function is then reconstructed as follows: where N ∈ Z represents the cut-off frequency.By examining Eq. ( 36), we observe that the coefficients Υ n play a significant role in characterizing the stability and resolution of the inverse scattering problem.These coefficients typically have complex forms.In the following section, we explore some special cases.Referring to the Abbe diffraction limit, we define the minimum feature size d n for the n-th Fourier mode of f as half of the period of e iαnx , which is expressed as d n = Λ/(2|n|).Additionally, we denote the wavelength of the incident field in free space as λ = 2π/κ.

No-slab Case
First, we consider the case in which the slab is absent, corresponding to the parameter values = 1 and µ = 1.It follows that η = κ and γ n = β n .Substituting these parameters into Eq.( 23), we obtain Plugging Eq. (37) into Eq.( 33) yields which is consistent with the result presented in [17].We distinguish two cases: 1.If |α n | < κ, i.e., |n| < Λ/λ, then β n is a real number, and Therefore, we can achieve stable reconstruction of the corresponding Fourier modes f (n) .Moreover, since the minimum feature size d n > λ/2, the resolution of these modes falls within the limits set by the Abbe diffraction limit.

If |α
and is approximately an exponential function of |n|.Thus, the reconstruction of the corresponding Fourier modes f (n) becomes increasingly unstable as |n| increases.
A comprehensive and precise analysis of errors in this scenario can be found in [27].

Ordinary Dielectric Material
Let us assume that the slab is made of an ordinary dielectric material with parameter values > 1 and µ = 1.This case was previously analyzed by [34].The novel reconstruction formula presented in this paper, Eq. ( 36), establishes a one-to-one correspondence between the Fourier coefficients.This is in contrast to the convolutional type found in [34].As a result, the formula is significantly simpler and more convenient to analysis.

Double Negative Metamaterial
In this section, we focus on the main topic of our paper, which is the examination of slabs that have negative permittivity ε (or negative real part) and negative permeability µ.

Perfectly Matched Parameters
We start by assuming = µ = −1, referred to as perfectly matched parameters, in line with the original concept of the superlens proposed by [30].As a result, we have η = κ and γ n = β n .By substituting these parameters into Eq.( 23), we obtain Substituting Eq. (38) into Eq.( 33) gives ) for all n ∈ Z.Therefore, we can reconstruct all frequency modes with stability, leading to unlimited resolution!(ii) If b < 2a, then Consequently, this scenario is equivalent to the case without a slab but with a measurement distance of 2a−b.Furthermore, since 2a−b < b, we can expect an improvement in resolution in this situation.

Imperfect Parameters
Ideal parameters are often unattainable in practical applications due to physical and engineering limitations.In particular, the lens constituent material is typically lossy, which leads to an effective permittivity that is a complex-valued quantity.To address this issue, we introduce the parameters where ε and µ are small-amplitude, real or complex numbers.Using Equation ( 19), we obtain Let n be a fixed integer from Z.As ι → 0, we can easily confirm that |e iβn(b−2a) |, which is the scaling factor for perfectly matched parameters.This implies that the resolution limit approaches that of the perfectly matched parameters as ε → 0 and µ → 0.

Numerical Experiments
In the following numerical experiments, we keep the period fixed at Λ = 1 and the wavelength at λ = 1.1.To obtain the total field u for a given scattering surface f , we solve the corresponding direct scattering problem using the finite element method along with the perfectly matched layer Figure 2: Sketch of the computational domain for the finite element method.A PML is placed above Γ b to absorb the scattered field, allowing u = 0 to be used as the boundary condition on the outer boundary of the PML.Periodic boundary conditions are enforced on the left and right boundaries, while the usual PEC boundary condition u = 0 is imposed on Γ f .(PML) technique to truncate the infinite computational domain to a finite one.Fig. 2 shows a sketch of the computational domain.
The field data is then sampled at points uniformly distributed on Γ b , and random values are added to the samples to simulate measurement noise.Specifically, the measurement data is given by where M = 100, and r m are independent random numbers drawn from the uniform distribution in [−0.05, 0.05], i.e., the relative noise level is 5%.Throughout the experiment, the lower boundary of the slab is fixed at a = 0.1 while the upper boundary is fixed at b = 0.2.

Smooth profile function
The scattering surface function is defined as f (x) = δg(x), where δ = 0.01, and the function represents frequency modes with |n| = 1, 3, 10.The clear cut-off of the Fourier modes is used to test the proposed method's capabilities.
To recover the first, first two, or all three Fourier modes of the profile, the cut-off frequency N must be set to at least 1, 3, or 10, respectively.The reconstructions for the case without the slab are shown in Row 1 of Figure 3.When N = 1, only the first Fourier mode, 0.4 cos(2πt), is recovered.To recover the surface profile up to the 3rd Fourier mode, N should be set to at least 3.However, using N = 3 leads to a solution with a significant error.Attempting to recover the 10th mode with N = 10 proves unsuccessful.Better resolution can be achieved by decreasing the measurement distance b, but there is a lower limit for b, at least on top of the surface, and the error is not a monotonic function of b.A detailed error estimate in this case can be found in [27].
We can increase the resolution by introducing a dielectric slab with a high refractive index.Let us set ε = 16 and µ = 1, which correspond to a refractive index of 4. According to the analysis presented in [34], the resolution can be increased by approximately a factor of 4. Therefore, the profile can be reconstructed stably up to the 3rd Fourier mode, but not the 10th mode, as confirmed by the numerical results shown in Fig. 3 (Row 2).
We now replace the dielectric slab with a double negative metamaterial having perfectly matched parameters, i.e., ε = µ = −1.According to the analysis presented in Section 5.3, we should be able to recover all the frequency modes of the profile.This claim is verified by the numerical results shown in Figure 3 (Row 3), where the reconstructed profile remains unpolluted even when N = 10.
Lastly, we consider the case of imperfect parameters, taking ε = −1 + 0.05i and µ = −0.97 as an example.The corresponding reconstructions are shown in Fig. 3 (Row 4).As expected, the results deteriorate with increasing deviation from the ideal case, particularly for N = 10, but they are still relatively close to the ideal case.The reconstruction accuracy will naturally    33) is plotted on a logarithmic scale for various values of ε and µ in the numerical experiments presented in Fig. 3.

Nonsmooth Profiles
Although the derivation of the reconstruction formula requires the existence of g and g , our method can still be applied to nonsmooth profiles due to the approximation properties of Fourier series.We first consider a nonsmooth profile defined by f (x) = δg(x), where δ = 0.01 and Note that this profile contains infinitely many nonzero Fourier modes.The results of the experiment are displayed in Fig. 5.For the case without a slab (ε = 1.0, µ = 1.0), we begin with the cut-off frequency N = 1.The reconstructed profile is smooth, but the two humps are indistinguishable.With N = 2, the reconstruction starts capturing the humps, but the accuracy remains low.When we attempt to increase the accuracy by raising N to 4, the reconstruction becomes significantly contaminated by noise.
Next, we insert a slab with a high refractive index (ε = 16.0,µ = 1.0).The reconstruction now becomes stable and more accurate with N = 4.However, increasing N to 5 leads to unstable results again.
Finally, we replace the slab with a superlens (ε = −1.0,µ = −1.0).The reconstruction is now stable and accurate up to at least N = 8.
We conducted another experiment using a setup similar to the previous one and present the results in Fig. 6.In the case without a slab, with ε = 1.0 and µ = 1.0, we did not obtain satisfactory results for any cut-off frequency, N .For the slab with a high refractive index characterized by ε = 16.0 and µ = 1.0, the optimal reconstruction is achieved with N = 2, while the quality deteriorates for N 3.In the superlens case with ε = −1.0 and µ = −1.0, the reconstruction remains stable up to at least N = 20.However, the Gibbs phenomenon becomes more pronounced for larger values of N .

Conclusion
In this paper, we investigated the superresolution effect of double negative index materials using an inverse scattering problem approach for diffraction gratings with low amplitude profiles.By applying the transformed field expansion method, we derived a simple reconstruction formula for the linear approximation of the scattering problem.Our findings suggest that it is feasible to achieve unlimited resolution with perfectly matched values of permittivity and permeability.To demonstrate the effectiveness of our reconstruction method, we conducted numerical experiments for both smooth and nonsmooth profiles with perfect or imperfect parameters.
There are several directions for future research on this topic.Firstly, the direct scattering problem with sign-changing coefficients for periodic structures requires further exploration to establish its well-posedness.The uniqueness of the inverse problem can also be investigated, and the convergence and error estimates of our proposed method can be analyzed.Moreover, it is possible to extend our method to other forms of waves and to the three-dimensional scenarios.For more general profile or inverse medium scattering problems, optimization or sampling methods may prove effective.Finally, to apply our method in real-world scenarios, it is crucial to incorporate the metamaterial's constitution and measuring device into the model.

Figure 1 :
Figure 1: Geometry of the model problem.The slab is positioned with its lower boundary Γ h and upper boundary Γ b above the scattering surface Γ f , and the measurements are taken on Γ b .

Figure 4 :
Figure 4: The scaling factor |Υ n | in Eq. (33) is plotted on a logarithmic scale for various values of ε and µ in the numerical experiments presented in Fig.3.