1. 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–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–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–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–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 $\varepsilon = - 1$ and relative permeability $\mu = - 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–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. When it comes to quantifying imaging performance, a succession of studies from the literature has delved into the geometrical parameters of the background medium and the configuration of the measurement. The works in [35,36] offer a quantitative analysis of how inhomogeneous scenarios influence performance. Meanwhile, [37,38] distinctly establish the impact of the evanescent contribution and the confines of the measurement domain, respectively.

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.

2. Model problem

Consider the electromagnetic scattering by a periodic surface defined as

We introduce a slab of double negative index material above the scattering surface, with its two surfaces defined by $\Gamma _a = \{ (x, a): x \in \mathbb {R} \}$ and $\Gamma _b = \{ (x, b): x \in \mathbb {R} \}$. The domain bounded between $\Gamma _f$ and $\Gamma _a$ is denoted by $\Omega = \{ (x, y) \in \mathbb {R}^2: f (x) < y < a \}$, and the domain between $\Gamma _a$ and $\Gamma _b$ is represented by $D = \{ (x, y) \in \mathbb {R}^2: a < y < b \}$. The problem geometry is illustrated in Fig. 1.

 

Fig. 1. Geometry of the model problem. The slab is positioned with its lower boundary $\Gamma _a$ and upper boundary $\Gamma _b$ above the scattering surface $\Gamma _f$, and the measurements are taken on $\Gamma _b$.

Download Full Size | PDF

The propagation of electromagnetic waves is described by the time-harmonic Maxwell’s equations:

The relative permittivity and permeability are denoted by $\varepsilon _r = \varepsilon / \varepsilon _0$ and $\mu _r = \mu / \mu _0$, respectively, where $\varepsilon _0$ and $\mu _0$ are the permittivity and permeability of vacuum, respectively. Define the free-space wavenumber $\kappa = \omega \sqrt {\varepsilon _0 \mu _0}$, then Eq. (3) can be rewritten as

A plane wave with an incident angle $\theta \in (- \pi / 2, \pi / 2)$ is described by $u_{\mathrm {in}} (x, y) = e^{i (\alpha x - \beta y)}$, where $\alpha = \kappa \sin \theta$ and $\beta = \kappa \cos \theta$. However, we will focus on normal incidence, i.e., $\theta = 0$. In this case, the incident field simplifies to

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 $\Gamma _b$.

The incident field given by Eq. (4) clearly satisfies the Helmholtz equation:

Let $u_{\mathrm {sc}} = u - u_{\mathrm {in}}$ denote the scattered field, which satisfies the same Helmholtz equation:

Assuming that $u_{\mathrm {sc}}$ consists only of upward propagating and evanescent waves, we can obtain the Rayleigh expansion from Eq. (5) as follows:

Taking $\partial _y$ in Eq. (6) and evaluating at $\Gamma _b$, we get

Adding Eq. (9) to Eq. (8) yields the transparent boundary condition:

From the interface conditions $\nu \times \mathbf {E}^+=\nu \times \mathbf {E}^-$ and $\nu \times \mathbf {H}^+=\nu \times \mathbf {H}^-$ on $\Gamma _b$ with $\nu =[0, 0, 1]^\top$, we obtain

Similarly, we have

For simplicity of notations, we will denote the relative permittivity and permeability by $\epsilon$ and $\mu$, respectively, in the following. To summarize, we obtain the following boundary-interface value problem for the total field $u$ in $\overline {D \cup \Omega }$:

Remark 1 It is worth noting that our model presumes an idealized infinite transverse dimension for the surface. When dealing with periodic surfaces with bounded transverse dimensions, it might be beneficial to incorporate approximation methodologies, as studied in [39]. However, such intricacies fall outside the scope of this particular study.

3. Transformed field expansion

We apply the TFE method, along with power series and Fourier series expansions, to reduce the governing Eqs. (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

Applying the chain rule leads to the following differentiation rules in $\Omega$:

Let

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:

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

As the variables $u_m, v_m, \tau _m$, and $\rho _m$ are periodic in the $x$ direction with a period of $\Lambda$, the Fourier series expansion can be applied to the governing Eqs. (16) to obtain a recursive system of boundary-interface value problems:

4. Analytical solutions

Solving Eqs. (18) by the variation of parameters leads to the general solution

Applying the boundary and interface conditions (18a), (18c), and (18e), we obtain the linear system

A direct calculation shows that the determinant of the matrix in Eq. (21) is given by

4.1 Zeroth order term

For $m = 0$ it follows from Eq. (17) that

The Fourier coefficients are 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 $\delta$ is zero.

4.2 First order term

For $m = 1$ we have $\rho _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

By substituting Eq. (25) into Eq. (17) and taking the Fourier transform, we obtain the following expressions for $v_1^{(n)}(y)$ and $\tau _1^{(n)}$:

By substituting Eqs. (28) into Eq. (27) and using the identity $\kappa ^2 = \alpha _n^2 + \beta _n^2$, we derive the expression for $\psi _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:

Remark 2 For the sake of conciseness, we omit specific cases where either $\beta _n = 0$ or $\gamma _n = 0$ from this paper. However, it is important to note that these cases can be addressed separately if necessary.

5. 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:

By examining Eq. (36), we observe that the coefficients $\Upsilon _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 \alpha _n x}$, which is expressed as $d_n = \Lambda / (2|n|)$. Additionally, we denote the wavelength of the incident field in free space as $\lambda = 2 \pi / \kappa$.

5.1 No-slab case

First, we consider the case in which the slab is absent, corresponding to the parameter values $\epsilon = 1$ and $\mu = 1$. It follows that $\eta = \kappa$ and $\gamma _n = \beta _n$. Substituting these parameters into Eq. (23), we obtain

Plugging Eq. (37) into Eq. (33) yields

5.2 Ordinary dielectric material

Let us assume that the slab is made of an ordinary dielectric material with parameter values $\epsilon > 1$ and $\mu = 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.

5.3 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 $\varepsilon$ (or negative real part) and negative permeability $\mu$.

5.3.1 Perfectly matched parameters

We start by assuming $\epsilon = \mu = - 1$, referred to as perfectly matched parameters, in line with the original concept of the superlens proposed by [30]. As a result, we have $\eta = \kappa$ and $\gamma _n = \beta _n$. By substituting these parameters into Eq. (23), we obtain

Substituting Eq. (38) into Eq. (33) gives

5.3.2 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

Remark 3 The proposed method presumes prior knowledge of the unknown surface’s period $\Lambda$ , which may not be known in real-world scenarios. One potential approach to address this problem is establishing an inverse problem wherein the period acts as an unknown parameter, set to be reconstructed separately or together with the profile function.

Remark 4 In this study, we focused solely on periodic surfaces. The derivation of our method carries through to non-periodic surfaces by replacing Fourier series with Fourier transform. Since the scattered field is not periodic anymore, the inverse problem then requires the measurement data on an infinite plane above the surface, which is infeasible in practice. In fact, a more suitable incident wave for such surfaces would be a tapered wave instead of a plane wave. We refer to [7] for an iterative method developed to reconstruct multiscale rough surfaces (non-periodic) with tapered waves of multiple frequencies.

6. Numerical experiments

In the following numerical experiments, we keep the period fixed at $\Lambda = 1$ and the wavelength at $\lambda = 1.1$. The lower boundary of the slab is fixed at $a = 0.1$ while the upper boundary is fixed at $b = 0.2$. 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 (PML) technique to truncate the infinite computational domain to a finite one. Figure 2 shows a sketch of the computational domain.

 

Fig. 2. Sketch of the computational domain for the finite element method. A PML is placed above $\Gamma _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 $\Gamma _f$.

Download Full Size | PDF

The field data undergoes sampling at uniformly distributed points on $\Gamma _b$. To simulate measurement noise, random values are incorporated into these samples. The measurement data is represented as:

6.1 Smooth profile function

The scattering surface function is defined as $f (x) = \delta g (x)$, where $\delta = 0.01$, and the function

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 Fig. 3. When $N=1$, only the first Fourier mode, $0.4 \cos (2 \pi x)$, 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].

 

Fig. 3. Numerical experiments for the smooth profile function depicted in Eq. (39). The red solid line represents the true profile, while the blue dashed line illustrates the reconstructed profile. The relative permittivity and permeability of the slab are specified as follows: Row 1: $\varepsilon = 1.0$ and $\mu = 1.0$; Row 2: $\varepsilon = 16.0$ and $\mu = 1.0$; Row 3: $\varepsilon = -1.0$ and $\mu = -1.0$; Row 4: $\varepsilon = -1 + 0.05i$ and $\mu = -0.97$; Row 5: $\varepsilon = -1 + 0.1i$ and $\mu = -1.06$. The cut-off frequency is specified as follows: Column 1: $N = 1$; Column 2: $N = 3$; Column 3: $N = 10$.

Download Full Size | PDF

We can increase the resolution by introducing a dielectric slab with a high refractive index. Let us set $\varepsilon = 16$ and $\mu = 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., $\varepsilon = \mu = -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 Fig. 3 (Row 3), where the reconstructed profile remains unpolluted even when $N = 10$.

Next, we consider the case of imperfect parameters, taking $\varepsilon = - 1 + 0.05 i$ and $\mu = - 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 decrease as $\varepsilon$ and $\mu$ move further away from $-1$. Figure 3 (Row 5) displays the reconstruction with $\varepsilon = - 1 + 0.1 i$ and $\mu = - 1.06$. Note that the first and third modes are still well recovered, while the 10th mode is not.

Lastly, we fix $\varepsilon = -1.0$, $\mu = -1$ and $N=10$ to investigate the effect of the sampling rate $M$ and the noise level $\sigma$. With $\sigma$ fixed at $0.05$, Fig. 4 (Row 1) shows the reconstructions with $M=25, 50, 100$ respectively. Note that the sampling rate $M=25$ is barely above the threshold specified by the Nyquist-Shannon sampling theorem. Consequently, the outcome is noticeably less optimal compared to scenarios with $M=50$ or $M=100$. On the other hand, the results from $M=50$ and $M=100$ are virtually indistinguishable, suggesting that $M=50$ suffices for this instance.

 

Fig. 4. Numerical experiments for the smooth profile function depicted in Eq. (39). The red solid line represents the true profile, while the blue dashed line illustrates the reconstructed profile. In all results, $\varepsilon = -1.0$, $\mu = -1$ and $N=10$. Row 1: $\sigma =0.05$; Column 1–3: $M=25,50,100$. Row 2: $\sigma =1.0$; Column 1–3: $M=50,100,200$.

Download Full Size | PDF

To observe a noticeable impact from the noise, the noise level must be raised significantly beyond $0.05$. When setting $\sigma$ to $1.0$, Fig. 4 (Row 2) displays the reconstructions for $M=50, 100, 200$, respectively. Observe that the outcomes at $M=100$ or $M=200$ are only marginally less optimal than those achieved at $\sigma =0.05$, indicating a high tolerance to noise.

In Fig. 5, we present the absolute value of the scaling factor $|\Upsilon _n|$ defined in Eq. (33) plotted against the mode number $n$ for the various values of $\varepsilon$ and $\mu$ used in the previous numerical experiments. We make the following observations:

 

Fig. 5. The scaling factor $|\Upsilon _n|$ in Eq. (33) is plotted on a logarithmic scale for various values of $\varepsilon$ and $\mu$ in the numerical experiments presented in Fig. 3.

Download Full Size | PDF

6.2 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) = \delta g(x)$, where $\delta = 0.01$ and

Note that this profile contains infinitely many nonzero Fourier modes. The results of the experiment are displayed in Fig. 6. For the case without a slab ($\varepsilon = 1.0$, $\mu = 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.

 

Fig. 6. Numerical experiments for the nonsmooth profile function (40). The red solid line represents the true profile, while the blue dashed line corresponds to the reconstructed profile. Row 1: $\varepsilon = 1.0$, $\mu = 1.0$, $N = 1, 2, 4$; Row 2: $\varepsilon = 16.0$, $\mu = 1.0$, $N = 1, 4, 5$; Row 3: $\varepsilon = -1.0$, $\mu = -1.0$, $N = 1, 4, 8$.

Download Full Size | PDF

Next, we insert a slab with a high refractive index ($\varepsilon = 16.0$, $\mu = 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 ($\varepsilon = - 1.0$, $\mu = - 1.0$). The reconstruction is now stable and accurate up to at least $N = 8$.

Next, we consider a discontinuous profile given by $f (x) = \delta g (x)$, where $\delta = 0.001$ and

 

Fig. 7. Numerical experiments for the nonsmooth profile function (41). The red solid line represents the true profile, while the blue dashed line corresponds to the reconstructed profile. Row 1: $\varepsilon = 1.0, \mu = 1.0, N = 1, 2, 3$; Row 2: $\varepsilon = 16.0, \mu = 1.0, N = 1, 2, 3$; Row 3: $\varepsilon = - 1.0, \mu = - 1.0, N = 1, 2, 20$.

Download Full Size | PDF

7. 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.