Inverse diffraction grating of Maxwell’s equations in biperiodic structures

Gang Bao; Tao Cui; Peijun Li

doi:10.1364/OE.22.004799

1. 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–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–11]. Computationally, a number of methods have been developed for the reconstruction of perfectly conducting grating surfaces in the transverse electric polarization [12–19]. We refer to [20–23] for related optimal design problems in diffractive optics, which are to design grating structures to obtain some specified diffraction patterns, and [24–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–31] for the existence, uniqueness, and numerical approximations of solutions for the direct two-dimensional grating problems. Mathematical studies can be found in [32–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–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–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.

2. Model problem

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

2.1. Maxwell’s equations

Let us first specify the diffraction grating problem geometry. Denote (ρ, z) ∈ ℝ³, where ρ = (x, y) ∈ ℝ². As seen in Fig. 1, the problem may be restricted to a single period of Λ = (Λ₁, Λ₂) in ρ due to the periodicity of the structure. Let the surface in one period be described by S = {(ρ, z) ∈ ℝ³ : z = ϕ(ρ), 0 < x < Λ₁, 0 < y < Λ₂}, where ϕ ∈ C²(ℝ²) is a biperiodic function satisfying ϕ(ρ + Λ, z) = ϕ(ρ, z). The grating surface function ϕ is assumed to be in the form

ϕ (ρ) = ε ψ (ρ),

where ψ ∈ C²(ℝ²) is a biperiodic function with period Λ and is called the grating profile, and ε is sufficiently small and is called the surface deformation parameter.

Fig. 1 Geometry of the diffraction grating problem.

Download Full Size | PDF

Denote by Ω_S = {(ρ, z) ∈ ℝ³ : z > ϕ(ρ), 0 < x < Λ₁, 0 < y < Λ₂} the space above S, which is filled with some homogeneous medium characterized by a positive constant wavenumber κ. Denote Ω = {(ρ, z) ∈ ℝ³ : ϕ(ρ) < z < h, 0 < x < Λ₁, 0 < y < Λ₂} be the domain bounded below by S and bounded above by the plane surface Γ = {(ρ, z) ∈ ℝ³ : z = h, 0 < x < Λ₁, 0 < y < Λ₂}, where h > max_{0 < x < Λ₁,0 < y < Λ₂}ϕ(ρ).

Let (E^inc, H^inc) be the incoming plane waves that are incident upon the grating surface from above, where

E^{inc} = t exp (i κ (α \cdot ρ - β z)), H^{inc} = s exp (i κ (α \cdot ρ - β z)) .

Here α = (α₁, α₂), α₁ = sinθ₁ cosθ₂, α₂ = sinθ₁ sinθ₂, and β = cosθ₁, where θ₁ and θ₂ are the latitudinal and longitudinal incident angles, which satisfy 0 ≤ θ₁ < π/2, 0 ≤ θ₂ < 2π. Denote by d = (α₁, α₂, −β) the unit propagation direction vector. The unit polarization vectors t = (t₁, t₂, t₃), s = (s₁, s₂, s₃) satisfy t · d = 0 and s = d × t, which gives explicitly

s_{1} = α_{2} t_{3} + β t_{2}, s_{2} = - (α_{1} t_{3} + β t_{1}), s_{3} = α_{1} t_{2} - α_{2} t_{1} .

For normal incident, i.e., θ₁ = 0, we have α₁ = 0, α₂ = 0, β = 1, and s₁ = t₂, s₂ = −t₁, s₃ = 0. Hence we get from |t| = |s| = 1 that

t_{1}^{2} + t_{2}^{2} = 1

, t₃ = 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 $E^{inc} = (E_{1}^{inc}, E_{2}^{inc}, E_{3}^{inc})$ and $H^{inc} = (H_{1}^{inc}, H_{2}^{inc}, H_{3}^{inc})$ . Under the normal incidence, the incoming plane waves (2) reduce to

E_{j}^{inc} = t_{j} exp (- i κ z), H_{j}^{inc} = s_{j} exp (- i κ z) .

It can be verified that the incident electromagnetic waves satisfy the three-dimensional time-harmonic Maxwell equation:

\nabla \times E^{inc} - i κ H^{inc} = 0, \nabla \times H^{inc} + i κ E^{inc} = 0, in ℝ^{3} .

The scattering of time-harmonic electromagnetic waves follows Maxwell’s equations in the space above the grating surface:

\nabla \times E - i κ H = 0, \nabla \times H + i κ E = 0, in Ω_{S},

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:

\nabla \cdot E = 0 and \nabla \cdot H = 0 in Ω_{S} .

We consider the perfect electric conductor condition:

E \times ν_{S} = 0 on S,

where ν_S = (ν₁, ν₂, ν₃) ∈ ℝ³ is the unit normal vector on S, given explicitly as

ν_{1} = \frac{ϕ_{x}}{{(1 + ϕ_{x}^{2} + ϕ_{y}^{2})}^{1 / 2}}, ν_{2} = \frac{ϕ_{y}}{{(1 + ϕ_{x}^{2} + ϕ_{y}^{2})}^{1 / 2}}, ν_{3} = \frac{- 1}{{(1 + ϕ_{x}^{2} + ϕ_{y}^{2})}^{1 / 2}} .

Here ϕ_x = ∂_xϕ(x, y) and ϕ_y = ∂_yϕ(x, y) are the partial derivatives.

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

Let n = (n₁, n₂) ∈ 𝕑² and denote α_n = (α₁_n, α₂_n), where α₁_n = 2πn₁/Λ₁ and α₂_n = 2πn₂/Λ₂. For a biperiodic function u(ρ) with period Λ in ρ, it has the Fourier series expansion

u (ρ) = \sum_{n \in ℤ^{2}} u_{n} exp (i α_{n} \cdot ρ), u_{n} = Λ_{1}^{- 1} Λ_{2}^{- 1} \int_{0}^{Λ_{1}} \int_{0}^{Λ_{2}} u (ρ) exp (- i α_{n} \cdot ρ) d ρ .

For any vector field u = (u₁, u₂, u₃), denote its tangential component on Γ by

u_{Γ} = ν_{Γ} \times (u \times ν_{Γ}) = (u_{1} (ρ, h), u_{2} (ρ, h), 0),

where ν_Γ = (0, 0, 1) is the unit normal vector on Γ.

For any tangential vector u(ρ, h) = (u₁(ρ, h), u₂(ρ, h), 0) on Γ, where u_j is a biperiodic function in ρ with period Λ, we define a boundary operator T :

T u = (v_{1} (ρ, h), v_{2} (ρ, h), 0),

where v_j is also a biperiodic function in ρ with the same period Λ. Here u_j and v_j have the following Fourier expansions

u_{j} (ρ, h) = \sum_{n \in ℤ^{2}} u_{j n} (h) exp (i α_{n} \cdot ρ), v_{j} (ρ, h) = \sum_{n \in ℤ^{2}} v_{j n} (h) exp (i α_{n} \cdot ρ),

and the Fourier coefficients u_jn and v_jn satisfy

{\begin{matrix} v_{1 n} (h) = \frac{1}{κ β_{n}} [(κ^{2} - α_{2 n}^{2}) u_{1 n} (h) + α_{1 n} α_{2 n} u_{2 n} (h)], \\ v_{2 n} (h) = \frac{1}{κ β_{n}} [(κ^{2} - α_{1 n}^{2}) u_{2 n} (h) + α_{1 n} α_{2 n} u_{1 n} (h))] . \end{matrix}

Using the boundary operator (9), we may derive the transparent boundary condition [29]:

(\nabla \times E) \times ν_{Γ} = i κ T E_{Γ} + f, on Γ,

where

f = i κ (H^{inc} \times ν_{Γ} - T E_{Γ}^{inc}) = (f_{1}, f_{2}, 0) .

Recalling the incident fields (3) and using the boundary operator (9), we have explicitly that

f_{1} = - 2 i κ t_{1} exp (- i κ h) and f_{2} = - 2 i κ t_{2} exp (- i κ h) .

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

{\begin{matrix} \nabla \times (\nabla \times E) - κ^{2} E = 0 in Ω, \\ E \times ν_{S} = 0 on S, \\ (\nabla \times E) \times ν_{Γ} - i κ T E_{Γ} = f on Γ . \end{matrix}

It is convenient introducing an equivalent scalar form of the problem (11) in order to apply the transformed field expansion. Denote E = (E₁, E₂, E₃). Noting the divergence condition (6), we may reduce Maxwell’s equations to the Helmholtz equation for E_j:

Δ E_{j} + κ^{2} E_{j} = 0 in Ω .

The divergence free condition (6) can be explicitly written as

\partial_{x} E_{1} + \partial_{y} E_{2} + \partial_{z} E_{3} = 0 in Ω .

Substituting (8) into (7) yields

E_{2} + ϕ_{y} E_{3} = 0, E_{1} + ϕ_{x} E_{3} = 0, ϕ_{y} E_{1} - ϕ_{x} E_{2} = 0, on z = ϕ (ρ) .

The transparent boundary condition (10) becomes

\partial_{z} E_{1} - \partial_{x} E_{3} = i κ H_{1} + f_{1}, \partial_{z} E_{2} - \partial_{y} E_{3} = i κ H_{2} + f_{2}, on z = h,

where the Fourier coefficients of H₁ and H₂ are given by

{\begin{matrix} H_{1 n} (h) = \frac{1}{κ β_{n}} [(κ^{2} - α_{2 n}^{2}) E_{1 n} (h) + α_{1 n} α_{2 n} E_{2 n} (h)], \\ H_{2 n} (h) = \frac{1}{κ β_{n}} [(κ^{2} - α_{1 n}^{2}) E_{2 n} (h) + α_{1 n} α_{2 n} E_{1 n} (h)] . \end{matrix}

Here E_1n(h) and E₂_n(h) are the Fourier coefficients of E₁(ρ, h) and E₂(ρ, 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 Γ, i.e., $E^{δ} (ρ, h) \times ν_{Γ} = (E_{2}^{δ} (ρ, h), - E_{1}^{δ} (ρ, h), 0)$ , 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π/κ.

3. Transformed field expansion

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

3.1. Change of variables

The transformed field expansion method begins with the change of variables:

\tilde{x} = x, \tilde{y} = y, \tilde{z} = h (\frac{z - ϕ}{h - ϕ}),

which maps the domain Ω into a rectangular slab

D = {(\tilde{x}, \tilde{y}, \tilde{z}) \in ℝ^{3} : 0 < \tilde{z} < h} = ℝ^{2} \times (0, h) .

Introduce a new function Ẽ = (Ẽ₁, Ẽ₂, Ẽ₃) 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

c_{1} \frac{\partial^{2} E_{j}}{\partial x^{2}} + c_{1} \frac{\partial^{2} E_{j}}{\partial y^{2}} + c_{2} \frac{\partial^{2} E_{j}}{\partial z^{2}} - c_{3} \frac{\partial^{2} E_{j}}{\partial x \partial z} - c_{4} \frac{\partial^{2} E_{j}}{\partial y \partial z} - c_{5} \frac{\partial E_{j}}{\partial z} + κ^{2} c_{1} E_{j} = 0 in D,

where

\begin{array}{l} c_{1} = {(h - ϕ)}^{2}, \\ c_{2} = (ϕ_{x}^{2} + ϕ_{y}^{2}) {(h - z)}^{2} + h^{2}, \\ c_{3} = 2 ϕ_{x} (h - z) (h - ϕ), \\ c_{4} = 2 ϕ_{y} (h - z) (h - ϕ), \\ c_{5} = (h - z) [(ϕ_{x x} + ϕ_{y y}) (h - ϕ) + 2 (ϕ_{x}^{2} + ϕ_{y}^{2})] . \end{array}

The divergence free condition (13) becomes

\partial_{x} E_{1} + \partial_{y} E_{2} - (\frac{h - z}{h - ϕ}) (ϕ_{x} \partial_{z} E_{1} + ϕ_{y} \partial_{z} E_{2}) + (\frac{h}{h - ϕ}) \partial_{z} E_{3} = 0 in D .

The perfect electric conductor condition (14) is

E_{2} + ϕ_{y} E_{3} = 0, E_{1} + ϕ_{x} E_{3} = 0, ϕ_{y} E_{1} - ϕ_{x} E_{2} = 0, z = 0 .

The transparent boundary condition (15) on z = h reduces to

(\frac{h}{h - ϕ}) \partial_{z} E_{1} - \partial_{x} E_{3} = i κ H_{1} + f_{1}, (\frac{h}{h - ϕ}) \partial_{z} E_{2} - \partial_{y} E_{3} = i κ H_{2} + f_{2} .

3.2. Power series solution

Consider a formal expansion of E_j in a power series of ε:

E_{j} (ρ, z; ε) = \sum_{k = 0}^{\infty} E_{j}^{(k)} (ρ, z) ε^{k} .

Substituting ϕ = εψ into c_j and inserting (20) into (16), we may derive

Δ E_{j}^{(k)} + κ^{2} E_{j}^{(k)} = F_{j}^{(k)} in D,

where

\begin{array}{l} F_{j}^{(k)} & = \frac{2 ψ}{h} \frac{\partial^{2} E_{j}^{(k - 1)}}{\partial x^{2}} + \frac{2 ψ}{h} \frac{\partial^{2} E_{j}^{(k - 1)}}{\partial y^{2}} + \frac{2 (h - z) ψ_{x}}{h} \frac{\partial^{2} E_{j}^{(k - 1)}}{\partial x \partial x} + \frac{2 (h - z) ψ_{y}}{h} \frac{\partial^{2} E_{j}^{(k - 1)}}{\partial y \partial z} \\ + \frac{(h - z) (ψ_{x x} + ψ_{y y})}{h} \frac{\partial E_{j}^{(k - 1)}}{\partial z} + \frac{2 κ^{2} ψ}{h} E_{j}^{(k - 1)} - \frac{ψ^{2}}{h^{2}} \frac{\partial^{2} E_{j}^{(k - 2)}}{\partial x^{2}} - \frac{ψ^{2}}{h^{2}} \frac{\partial^{2} E_{j}^{(k - 2)}}{\partial y^{2}} \\ - \frac{{(h - z)}^{2} (ψ_{x}^{2} + ψ_{y}^{2})}{h^{2}} \frac{\partial^{2} E_{j}^{(k - 2)}}{\partial z^{2}} - \frac{2 ψ ψ_{x} (h - z)}{h^{2}} \frac{\partial^{2} E_{j}^{(k - 2)}}{\partial x \partial z} - \frac{2 ψ ψ_{y} (h - z)}{h^{2}} \frac{\partial^{2} E_{j}^{(k - 2)}}{\partial y \partial z} \\ + \frac{(h - z) [2 (ψ_{x}^{2} + ψ_{y}^{2}) - ψ (ψ_{x x} + ψ_{y y})]}{h^{2}} \frac{\partial E_{j}^{(k - 2)}}{\partial z} - \frac{κ^{2} ψ^{2}}{h^{2}} E_{j}^{(k - 2)} . \end{array}

Here ψ_x = ∂_xψ(x, y) and ψ_y = ∂_yψ(x, y) are the partial derivatives.

Substituting (20) into the divergence free condition (17) yields

\partial_{x} E_{1}^{(k)} + \partial_{y} E_{2}^{(k)} + \partial_{z} E_{3}^{(k)} = w^{(k)} in D,

where

w^{(k)} = \frac{ψ}{h} (\partial_{x} E_{1}^{(k - 1)} + \partial_{y} E_{2}^{(k - 1)}) + (\frac{h - z}{h}) (ψ_{x} \partial_{z} E_{1}^{(k - 1)} + ψ_{y} \partial_{z} E_{2}^{(k - 1)}) .

The perfect electric conductor boundary condition (18) can be written as

E_{1}^{(k)} = u_{1}^{(k)}, E_{2}^{(k)} = u_{2}^{(k)},

where

u_{1}^{(k)} (ρ) = - ψ_{x} E_{3}^{(k - 1)}, u_{2}^{(k)} = - ψ_{y} E_{3}^{(k - 1)} .

Evaluating (22) at z = 0, we have

\partial_{x} E_{1}^{(k)} + \partial_{y} E_{2}^{(k)} + \partial_{z} E_{3}^{(k)} = u_{3}^{(k)}, z = 0,

where

u_{3}^{(k)} (ρ) = \frac{ψ}{h} (\partial_{x} E_{1}^{(k - 1)} + \partial_{y} E_{2}^{(k - 1)}) + (ψ_{x} \partial_{z} E_{1}^{(k - 1)} + ψ_{y} \partial_{z} E_{2}^{(k - 1)}) .

Inserting (20) into the transparent boundary condition (19), we get

{\begin{matrix} \partial_{z} E_{1}^{(k)} - \partial_{x} E_{3}^{(k)} = i κ H_{1}^{(k)} + v_{1}^{(k)}, \\ \partial_{z} E_{2}^{(k)} - \partial_{y} E_{3}^{(k)} = i κ H_{2}^{(k)} + v_{2}^{(k)}, \end{matrix}

where

\begin{array}{l} v_{1}^{(0)} = f_{1}, v_{1}^{(1)} = - \frac{ψ}{h} \partial_{z} E_{1}^{(0)}, v_{1}^{(k)} = - \frac{ψ}{h} (\partial_{x} E_{3}^{(k - 1)} + i κ H_{1}^{(k - 1)}), \\ v_{2}^{(0)} = f_{2}, v_{2}^{(1)} = - \frac{ψ}{h} \partial_{z} E_{2}^{(0)}, v_{2}^{(k)} = - \frac{ψ}{h} (\partial_{y} E_{3}^{(k - 1)} + i κ H_{2}^{(k - 1)}), \end{array}

and the Fourier coefficients of

H_{1}^{(k)} (ρ, h)

and

H_{2}^{(k)} (ρ, h)

are

{\begin{matrix} H_{1 n}^{(k)} (h) = \frac{1}{κ β_{n}} [(κ^{2} - α_{2 n}^{2}) E_{1 n}^{(k)} (h) + α_{1 n} α_{2 n} E_{2 n}^{(k)} (h)], \\ H_{2 n}^{(k)} (h) = \frac{1}{κ β_{n}} [(κ^{2} - α_{1 n}^{2}) E_{2 n}^{(k)} (h) + α_{1 n} α_{2 n} E_{1 n}^{(k)} (h)] . \end{matrix}

Here

E_{1 n}^{(k)} (h)

and

E_{2 n}^{(k)} (h)

are the Fourier coefficients of

E_{1}^{(k)} (ρ, h)

and

E_{1}^{(k)} (ρ, h)

, respectively. The divergence free condition (22) on z = h reduces to

\partial_{x} E_{1}^{(k)} + \partial_{y} E_{2}^{(k)} + \partial_{z} E_{3}^{(k)} = v_{3}^{(k)},

where

v_{3}^{(k)} (ρ) = \frac{ψ}{h} (\partial_{x} E_{1}^{(k - 1)} + \partial_{y} E_{2}^{(k - 1)}) .

Clearly, the problem (21)–(26) for $E_{j}^{(k)}$ involves $F_{j}^{(k)}$ , $u_{j}^{(k)}$ , $v_{j}^{(k)}$ , which depend only on previous two terms of $E_{j}^{(k - 1)}$ and $E_{j}^{(k - 2)}$ . Thus, the transformed problem (21)–(26) indeed can be solved efficiently in a recursive manner starting from k = 0.

3.3. Zeroth order term

Recalling the recurrence relation (21) and letting k = 0, we have

Δ E_{j}^{(0)} + κ^{2} E_{j}^{(0)} = 0 in D .

The divergence free condition (22) reduces to

\partial_{x} E_{1}^{(0)} + \partial_{y} E_{2}^{(0)} + \partial_{z} E_{3}^{(0)} = 0 in D .

The perfect electric conductor boundary condition (23) is

E_{1}^{(0)} (ρ, 0) = 0, E_{2}^{(0)} (ρ, 0) = 0 .

Using (28) and (29), we have

\partial_{z} E_{3}^{(0)} (ρ, 0) = - \partial_{x} E_{1}^{(0)} (ρ, 0) - \partial_{y} E_{2}^{(0)} (ρ, 0) = 0 .

The transparent boundary condition (25) becomes

{\begin{matrix} \partial_{z} E_{1}^{(0)} (ρ, h) - \partial_{x} E_{3}^{(0)} (ρ, h) = i κ H_{1}^{(0)} (ρ, h) + f_{1} (ρ), \\ \partial_{z} E_{2}^{(0)} (ρ, h) - \partial_{y} E_{3}^{(0)} (ρ, h) = i κ H_{2}^{(0)} (ρ, h) + f_{2} (ρ) . \end{matrix}

In addition, the divergence free condition (28) gives one more boundary condition:

\partial_{x} E_{1}^{(0)} (ρ, h) + \partial_{y} E_{2}^{(0)} (ρ, h) + \partial_{z} E_{3}^{(0)} (ρ, h) = 0 .

Since $E_{j}^{(0)} (ρ, z)$ and f_j are periodic functions of ρ, they have the Fourier expansions

E_{j}^{(0)} (ρ, z) = \sum_{n \in ℤ^{2}} E_{j n}^{(0)} (z) exp (i α_{n} \cdot ρ), f_{j} = \sum_{n \in ℤ^{2}} f_{j n} exp (i α_{n} \cdot ρ),

where f_j₀ = −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_{j n}^{(0)} (z)$ :

\frac{d^{2} E_{j n}^{(0)} (z)}{d z^{2}} + (κ^{2} - {| α_{n} |}^{2}) E_{j n}^{(0)} (z) = 0, 0 < z < h .

Similarly, we have the boundary conditions at z = 0:

E_{1 n}^{(0)} = 0, E_{2 n}^{(0)} = 0, E_{3 n}^{{(0)}^{'}} = 0,

and the boundary conditions at z = h:

{\begin{array}{l} E_{1 n}^{{(0)}^{'}} - i α_{1 n} E_{3 n}^{(0)} = \frac{i}{β_{n}} [(κ^{2} - α_{2 n}^{2}) E_{1 n}^{(0)} + α_{1 n} α_{2 n} E_{2 n}^{(0)}] + f_{1 n}, \\ E_{2 n}^{{(0)}^{'}} - i α_{2 n} E_{3 n}^{(0)} = \frac{i}{β_{n}} [(κ^{2} - α_{1 n}^{2}) E_{2 n}^{(0)} + α_{1 n} α_{2 n} E_{1 n}^{(0)}] + f_{2 n}, \\ E_{3 n}^{{(0)}^{'}} + i α_{1 n} E_{1 n}^{(0)} + i α_{2 n} E_{2 n}^{(0)} = 0 . \end{array}

Simple calculations yield that the solution of the two-point bounary value problem (34)–(36) is

E_{j}^{0} (ρ, z) = t_{j} (exp (- i κ z) - exp (i κ z)) .

3.4. First order term

Recalling (21) and letting k = 1, we have

Δ E_{j}^{(1)} + κ^{2} E_{j}^{(1)} = F_{j}^{(1)} in D,

where

\begin{array}{l} F_{j}^{(1)} = & \frac{2 ψ}{h} \frac{\partial^{2} E_{j}^{(0)}}{\partial x^{2}} + \frac{2 ψ}{h} \frac{\partial^{2} E_{j}^{(0)}}{\partial y^{2}} + \frac{2 (h - z) ψ_{x}}{h} \frac{\partial^{2} E_{j}^{(0)}}{\partial x \partial z} + \frac{2 (h - z) ψ_{y}}{h} \frac{\partial^{2} E_{j}^{(0)}}{\partial y \partial z} \\ + \frac{(h - z) (ψ_{x x} + ψ_{y y})}{h} \frac{\partial E_{j}^{(0)}}{\partial z} + \frac{2 κ^{2} ψ}{h} E_{j}^{(0)} . \end{array}

It follows from (37) that we have explicitly

\begin{array}{l} F_{j}^{(1)} (ρ, z) = & \frac{2 κ^{2} t_{j}}{h} ψ (exp (- i κ z) - exp (i κ z)) \\ - \frac{i κ t_{j} (h - z)}{h} (ψ_{x x} + ψ_{y y}) (exp (- i κ z) + exp (i κ z)) . \end{array}

The divergence free condition (22) reduces to

\partial_{x} E_{1}^{(1)} + \partial_{y} E_{2}^{(1)} + \partial_{z} E_{3}^{(1)} = w^{(1)} in D,

where

\begin{array}{l} w^{(1)} (ρ, z) & = \frac{ψ}{h} (\partial_{x} E_{1}^{(0)} + \partial_{y} E_{2}^{(0)}) + (\frac{h - z}{h}) (ψ_{x} \partial_{z} E_{1}^{(0)} + ψ_{y} \partial_{z} E_{2}^{(0)}) \\ = - \frac{i κ (h - z)}{h} (t_{1} ψ_{x} + t_{2} ψ_{y}) (exp (- i κ z) + exp (i κ z)) . \end{array}

The perfectly conducting boundary condition (23) on z = 0 is

\begin{array}{l} E_{1}^{(1)} (ρ, 0) = u_{1}^{(1)} (ρ) = - ψ_{x} (ρ) E_{3}^{(0)} (ρ, 0) = 0, \\ E_{2}^{(1)} (ρ, 0) = u_{2}^{(1)} (ρ) = - ψ_{y} (ρ) E_{3}^{(0)} (ρ, 0) = 0 . \end{array}

Evaluating (39) at z = 0 gives

\partial_{z} E_{3}^{(1)} (ρ, 0) = w^{(1)} (ρ, 0) - \partial_{x} E_{1}^{(1)} (ρ, 0) - \partial_{y} E_{2}^{(1)} (ρ, 0) = - 2 i κ (t_{1} ψ_{x} + t_{2} ψ_{y}) .

The transparent boundary condition (25) becomes

\partial_{z} E_{1}^{(1)} - \partial_{x} E_{3}^{(1)} = i κ H_{1}^{(1)} + v_{1}^{(1)}, \partial_{z} E_{2}^{(1)} - \partial_{y} E_{3}^{(1)} = i κ H_{2}^{(1)} + v_{2}^{(1)},

where

v_{j}^{(1)} (ρ) = - \frac{ψ}{h} \partial_{z} E_{j}^{(0)} (ρ, h) = \frac{i κ t_{j}}{h} (exp (- i κ h) + exp (i κ h)) ψ .

The divergence free condition (26) gives one more boundary condition on z = h:

\partial_{x} E_{1}^{(1)} (ρ, h) + \partial_{y} E_{2}^{(1)} (ρ, z) + \partial_{z} E_{3}^{(1)} (ρ, z) = 0 .

Consider the Fourier expansions for periodic functions of

E_{j}^{(1)} (ρ, z)

,

F_{j}^{(1)} (ρ, z)

, and ψ(ρ):

\begin{matrix} ψ (ρ) = \sum_{n \in ℤ^{2}} ψ_{n} exp (i α_{n} \cdot ρ), \\ E_{j}^{(1)} (ρ, z) = \sum_{n \in ℤ^{2}} E_{j n}^{(1)} (z) exp (i α_{n} \cdot ρ), \\ F_{j}^{(1)} (ρ, z) = \sum_{n \in ℤ^{2}} F_{j n}^{(1)} (z) exp (i α_{n} \cdot ρ), \end{matrix}

where

\begin{array}{l} F_{j n}^{(1)} (z) = & [\frac{2 κ^{2} t_{j}}{h} (exp (- i κ z) - exp (i κ z)) \\ + \frac{i κ t_{j} (h - z)}{h} (α_{1 n}^{2} + α_{2 n}^{2}) (exp (- i κ z) + exp (i κ z))] ψ_{n} . \end{array}

Substituting the above Fourier expansions into (38), we derive an equation for $E_{j n}^{(1)} (z)$ :

\frac{d^{2} E_{j n}^{(1)} (z)}{d z^{2}} + (κ^{2} - {| α_{n} |}^{2}) E_{j n}^{(1)} (z) = F_{j n}^{(1)} (z), 0 < z < h,

together with boundary conditions at z = 0:

E_{1 n}^{(1)} = 0, E_{2 n}^{(1)} = 0, E_{3 n}^{{(1)}^{'}} (0) = 2 κ (α_{1 n} + α_{2 n}) ψ_{n},

and boundary conditions at z = h:

{\begin{array}{l} E_{1 n}^{{(1)}^{'}} - i α_{1 n} E_{3 n}^{(1)} = \frac{i}{β_{n}} [(κ^{2} - α_{2 n}^{2}) E_{1 n}^{(1)} + α_{1 n} α_{2 n} E_{2 n}^{(1)}] + v_{1 n}^{(1)}, \\ E_{2 n}^{{(1)}^{'}} - i α_{2 n} E_{3 n}^{(1)} = \frac{i}{β_{n}} [(κ^{2} - α_{1 n}^{2}) E_{2 n}^{(1)} + α_{1 n} α_{2 n} E_{1 n}^{(1)}] + v_{2 n}^{(1)}, \\ E_{3 n}^{{(1)}^{'}} + i α_{1 n} E_{1 n}^{(1)} + i α_{2 n} E_{2 n}^{(1)} = 0, \end{array}

where

v_{1 n}^{(1)}

and

v_{1 n}^{(2)}

are the Fourier coefficients of

v_{1}^{(1)} (ρ)

and

v_{2}^{(1)} (ρ)

. Explicitly, we have

v_{j n}^{(1)} = \frac{i κ t_{j}}{h} (exp (- i κ h) + exp (i κ h)) ψ_{n} .

Following straightforward but tedious calculations, we may solve the two-point boundary value problem (41)–(43) and obtain elegant equations

E_{1 n}^{(1)} (h) = 2 i κ t_{1} exp (i β_{n} h) ψ_{n} and E_{2 n}^{(1)} (h) = 2 i κ t_{2} exp (i β_{n} h) ψ_{n},

which relate the Fourier coefficients of order one terms

E_{j}^{(1)} (ρ, h)

with the Fourier coefficient of the grating profile function ψ(ρ).

4. Reconstruction formula

Assume that the noisy data takes the form

E_{j}^{δ} (ρ, h) = E_{j} (ρ, h) + 𝒪 (δ),

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

E_{j} (ρ, h) = E_{j}^{(0)} (ρ, h) + ε E_{j}^{(1)} (ρ, h) + 𝒪 (ε^{2}) + 𝒪 (δ) .

Rearranging (45), and dropping 𝒪(ε²) and 𝒪(δ) yield

ε E_{j}^{(1)} (ρ, h) = E_{j}^{δ} (ρ, h) - E_{j}^{(0)} (ρ, h)

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

ϕ_{n} = {(2 i κ t_{j})}^{- 1} [E_{j n}^{δ} (h) - E_{j n}^{(0)} (h)] exp (- i β_{n} h),

where

E_{j n}^{δ} (h)

is the Fourier coefficient of the noisy data

E_{j}^{δ} (x, h)

and

E_{j n}^{(0)} (h)

is the Fourier coefficient of

E_{j}^{(0)} (x, h)

given as

E_{j n}^{(0)} (h) = t_{j} (exp (- i κ h) - exp (i κ h)) δ_{0 n} .

Here δ₀_n 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.

Following [51], we consider the spectral cut-off regularization. Define the signal-to-noise ratio (SNR) by

SNR = min {ε^{- 2}, δ^{- 1}} .

For fixed h, the cut-off wavenumber κ_c is chosen in such a way that

exp ({(κ_{c}^{2} - κ^{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

\frac{κ_{c}}{κ} = {[1 + {(\frac{log SNR}{κ h})}^{2}]}^{1 / 2},

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

ϕ_{n} = {(2 i κ t_{j})}^{- 1} [E_{j n}^{δ} (h) - E_{j n}^{(0)} (h)] exp (- i β_{n} h) χ_{n},

where the characteristic function

χ_{n} = {\begin{array}{l} 1 & for & | α_{n} | \leq κ_{c}, \\ 0 & for & | α_{n} | > κ_{c} . \end{array}

Once ϕ_n are computed, the grating surface function can be approximated by

\begin{array}{l} ϕ (ρ) & \approx \sum_{n \in ℤ} ϕ_{n} exp (i α_{n} \cdot ρ) = \sum_{| α_{n} | \leq κ_{c}} {(2 i κ t_{j})}^{- 1} [E_{j n}^{δ} (h) - E_{j n}^{(0)} (h)] exp (i (α_{n} \cdot ρ - β_{n} h)) \\ = \sum_{| α_{n} | \leq κ_{c}} {(2 i κ t_{j})}^{- 1} E_{j n}^{δ} (h) exp (i (α_{n} \cdot ρ - β_{n} h)) + {(2 i κ)}^{- 1} (1 - exp (- 2 i κ h)) . \end{array}

Hence, the method requires only two fast Fourier transforms: one is done for the data to obtain

E_{j n}^{δ} (h)

and another is done to obtain the approximated function ϕ.

5. 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²(ℝ²), it is still applicable to non-smooth functions numerically.

Fig. 2 The exact grating profile ψ. (a) Example 1: smooth grating profile with finite Fourier modes; (b) Example 2: non-smooth grating profile with infinite Fourier modes.

Download Full Size | PDF

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¹ 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 and t₂ = t₃ = 0, and only the first component of the electric field, E₁(ρ, 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₁(ρ, 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

E_{1}^{δ} (ρ, h) = E_{1} (ρ, h) (1 + δ rand),

where rand stands for uniformly distributed random numbers in [−1, 1]. The relative L²(R) error is defined by

e = \frac{{‖ ϕ - ϕ_{δ, ε} ‖}_{0, R}}{{‖ ϕ ‖}_{0, R}},

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²(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λ.

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

View Table | View all tables in this article

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

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

View Table | View all tables in this article

Fig. 3 Example 1: Reconstructed grating surfaces by using different h with ε = 0.025λ and δ = 5%. (a) h = 0.4λ; (b) h = 0.3λ; (c) h = 0.2λ; (d) h = 0.1λ.

Download Full Size | PDF

Example 2. This example illustrates the reconstruction results of a non-smooth grating profile with infinitely many Fourier modes, as seen in Fig. 2(b). The exact grating surface function is given by ϕ(ρ) = εψ(ρ), where the grating profile function

ψ (x, y) = | sin (2 π x) sin (2 π y) | - | cos (2 π x) cos (2 π y) | .

Clearly, the profile function (51) is nondifferentiable and its Fourier coefficients decay slowly. Comparing with the grating profile (50), it is more challenging to obtain as good reconstructions as those in Example 1 since a much higher cutoff frequency is desirable to recover as many Fourier modes as possible for (51).

Consider the influence of ε by using noise-free data. The measurement is taken at h = 0.2λ. Table 3 presents the relative L²(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.

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

View Table | View all tables in this article

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

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

View Table | View all tables in this article

Fig. 4 Example 2: Reconstructed grating surfaces by using different h with ε = 0.0125λ and δ = 5%. (a) h = 0.2λ; (b) h = 0.1λ; (c) h = 0.05λ; (d) h = 0.025λ.

Download Full Size | PDF

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

Acknowledgments

The research of GB was supported in part by the NSF grants DMS-0968360 and DMS-1211292, the ONR grant N00014-12-1-0319, a Key Project of the Major Research Plan of NSFC (No. 91130004), and a special research grant from Zhejiang University. The research of TC was partially supported by the National Basic Research Project Grant 2011CB309700, the National High Technology Research and Development Program of China Grant 2012AA01A309, and NSFC Grants 11101417 and 11171334. The research of PL was supported in part by the NSF grant DMS-1151308.

References and links

1. R. Petit, ed., Electromagnetic Theory of Gratings (Springer, 1980). [CrossRef]

2. G. Bao, D. Dobson, and J. A. Cox, “Mathematical studies in rigorous grating theory,” J. Opt. Soc. Am. A 12, 1029–1042 (1995). [CrossRef]

3. Z. Chen and H. Wu, “An adaptive finite element method with perfectly matched absorbing layers for the wave scattering by periodic structures,” SIAM J. Numer. Anal. 41, 799–826 (2003). [CrossRef]

4. J. C. Nédélec and F. Starling, “Integral equation methods in a quasi-periodic diffraction problem for the time-harmonic Maxwell’s equations,” SIAM J. Math. Anal. 22, 1679–1701 (1991). [CrossRef]

5. Y. Wu and Y. Y. Lu, “Analyzing diffraction gratings by a boundary integral equation Neumann-to-Dirichlet map method,” J. Opt. Soc. Am. A 26, 2444–2451 (2009). [CrossRef]

6. G. Bao, L. Cowsar, and W. Masters, Mathematical Modeling in Optical Science, Vol. 22 of Frontiers in Applied Mathematics (SIAM, 2001). [CrossRef]

7. G. Bao, “A unique theorem for an inverse problem in periodic diffractive optics,” Inverse Probl. 10, 335–340 (1994). [CrossRef]

8. G. Bao and A. Friedman, “Inverse problems for scattering by periodic structure,” Arch. Rational Mech. Anal. 132, 49–72 (1995). [CrossRef]

9. G. Bruckner, J. Cheng, and M. Yamamoto, “An inverse problem in diffractive optics: conditional stability,” Inverse Probl. 18, 415–433 (2002). [CrossRef]

10. F. Hettlich and A. Kirsch, “Schiffer’s theorem in inverse scattering theory for periodic structures,” Inverse Probl. 13, 351–361 (1997). [CrossRef]

11. A. Kirsch, “Uniqueness theorems in inverse scattering theory for periodic structures,” Inverse Probl. 10, 145–152 (1994). [CrossRef]

12. T. Arens and A. Kirsch, “The factorization method in inverse scattering from periodic structures,” Inverse Probl. 19, 1195–1211 (2003). [CrossRef]

13. G. Bao, P. Li, and H. Wu, “A computational inverse diffraction grating problem,” J. Opt. Soc. Am. A 29, 394–399 (2012). [CrossRef]

14. G. Bao, P. Li, and J. Lv, “Numerical solution of an inverse diffraction grating problem from phasless data,” J. Opt. Soc. Am. A 30, 293–299 (2013). [CrossRef]

15. G. Bruckner and J. Elschner, “A two-step algorithm for the reconstruction of perfectly reflecting periodic profiles,” Inverse Probl. 19, 315–329 (2003). [CrossRef]

16. J. Elschner, G. Hsiao, and A. Rathsfeld, “Grating profile reconstruction based on finite elements and optimization techniques,” SIAM J. Appl. Math. 64, 525–545 (2003).

17. F. Hettlich, “Iterative regularization schemes in inverse scattering by periodic structures,” Inverse Probl. 18, 701–714 (2002). [CrossRef]

18. K. Ito and F. Reitich, “A high-order perturbation approach to profile reconstruction: I. Perfectly conducting gratings,” Inverse Probl. 15, 1067–1085 (1999). [CrossRef]

19. A. Malcolm and D. P. Nicholls, “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. J. Appl. Math. 4, 321–340 (1993). [CrossRef]

21. D. Dobson, “Optimal shape design of blazed diffraction grating,” J. Appl. Math. Optim. 40, 61–78 (1999). [CrossRef]

22. J. Elschner and G. Schmidt, “Diffraction in periodic structures and optimal design of binary gratings: I. Direct problems and gradient formulas,” Math. Methods Appl. Sci. 21, 1297–1342 (1998). [CrossRef]

23. J. Elschner and G. Schmidt, “Numerical solution of optimal design problems for binary gratings,” J. Comput. Phys. 146, 603–626 (1998). [CrossRef]

24. I. Akduman, R. Kress, and A. Yapar, “Iterative reconstruction of dielectric rough surface profiles at fixed frequency,” Inverse Probl. 22, 939–954 (2006). [CrossRef]

25. R. Coifman, M. Goldberg, T. Hrycak, M. Israeli, and V. Rokhlin, “An improved operator expansion algorithm for direct and inverse scattering computations,” Waves Random Media 9, 441–457 (1999). [CrossRef]

26. J. A. DeSanto and R. J. Wombell, “The reconstruction of shallow rough-surface profiles from scattered field data,” Inverse Probl. 7, L7–L12 (1991). [CrossRef]

27. R. Kress and T. Tran, “Inverse scattering for a locally perturbed half-plane,” Inverse Probl. 16, 1541–1559 (2000). [CrossRef]

28. G. Bao, “Variational approximation of Maxwell’s equations in biperiodic structures,” SIAM J. Appl. Math. 57, 364–381 (1997). [CrossRef]

29. G. Bao, P. Li, and H. Wu, “An adaptive edge element method with perfectly matched absorbing layers for wave scattering by biperiodic structures,” Math. Comput. 79, 1–34 (2009). [CrossRef]

30. D. Dobson, “A variational method for electromagnetic diffraction in biperiodic structures,” Math. Model. Numer. Anal. 28, 419–439 (1994).

31. A. Lechleiter and D. L. Nguyen, “On uniqueness in electromagnetic scattering from biperiodic structures,” ESAIM Math. Modell. Numer. Anal. 47, 1167–1184 (2013). [CrossRef]

32. H. Ammari, “Uniqueness theorems for an inverse problem in a doubly periodic structure,” Inverse Probl. 11, 823–833 (1995). [CrossRef]

33. G. Bao, H. Zhang, and J. Zou, “Unique determination of periodic polyhedral structures by scattered electromagnetic fields,” Trans. Am. Math. Soc. 363, 4527–4551 (2011). [CrossRef]

34. G. Bao and Z. Zhou, “An inverse problem for scattering by a doubly periodic structure,” Trans. Am. Math. Soc. 350, 4089–4103 (1998). [CrossRef]

35. G. Hu, J. Yang, and B. Zhang, “An inverse electromagnetic scattering problem for a bi-periodic inhomogeneous layer on a perfectly conducting plate,” Appl. Anal. 90, 317–333 (2011). [CrossRef]

36. G. Hu and B. Zhang, “The linear sampling method for inverse electromagnetic scattering by a partially coated bi-periodic structures,” Math. Methods Appl. Sci. 34, 509–519 (2011). [CrossRef]

37. J. Yang and B. Zhang, “Inverse electromagnetic scattering problems by a doubly periodic structure,” Math. Appl. Anal. 18, 111–126 (2011).

38. A. Lechleiter and D. L. Nguyen, “Factorization method for electromagnetic inverse scattering from biperiodic structures,” SIAM J. Imaging Sci. 6, 1111–1139 (2013). [CrossRef]

39. D. L. Nguyen, Spectral Methods for Direct and Inverse Scattering from Periodic Structures (PhD thesis, Ecole Polytechnique, Palaiseau, France, 2012).

40. K. Sandfort, The factorization method for inverse scattering from periodic inhomogeneous media (PhD thesis, Karlsruher Institut für Technologie, 2010).

41. O. Bruno and F. Reitich, “Numerical solution of diffraction problems: a method of variation of boundaries,” J. Opt. Soc. Am. A 10, 1168–1175 (1993). [CrossRef]

42. O. Bruno and F. Reitich, “Numerical solution of diffraction problems: a method of variation of boundaries. III. Doubly periodic gratings,” J. Opt. Soc. Am. A 10, 2551–2562 (1993). [CrossRef]

43. Y. He, D. P. Nicholls, and J. Shen, “An efficient and stable spectral method for electromagnetic scattering from a layered periodic struture,” J. Comput. Phys. 231, 3007–3022 (2012). [CrossRef]

44. A. Malcolm and D. P. Nicholls, “A field expansions method for scattering by periodic multilayered media,” J. Acoust. Soc. Am. 129, 1783–1793 (2011). [CrossRef] [PubMed]

45. A. Malcolm and D. P. Nicholls, “A boundary perturbation method for recovering interface shapes in layered media,” Inverse Probl. 27, 095009 (2011). [CrossRef]

46. D. P. Nicholls and F. Reitich, “Shape deformations in rough surface scattering: cancellations, conditioning, and convergence,” J. Opt. Soc. Am. A 21, 590–605 (2004). [CrossRef]

47. D. P. Nicholls and F. Reitich, “Shape deformations in rough surface scattering: improved algorithms,” J. Opt. Soc. Am. A 21, 606–621 (2004). [CrossRef]

48. G. Bao and P. Li, “Near-field imaging of infinite rough surfaces,” SIAM J. Appl. Math. 73, 2162–2187 (2013). [CrossRef]

49. G. Bao and P. Li, “Near-field imaging of infinite rough surfaces in dielectric media,” SIAJ J. Imaging Sci., to appear.

50. T. Cheng, P. Li, and Y. Wang, “Near-field imaging of perfectly conducting grating surfaces,” J. Opt. Soc. Am. A 30, 2473–2481 (2013). [CrossRef]

51. G. Bao and J. Lin, “Near-field imaging of the surface displacement on an infinite ground plane,” Inverse Probl. Imag. 7, 377–396 (2013). [CrossRef]

52. S. Carney and J. Schotland, “Inverse scattering for near-field microscopy,” App. Phys. Lett. 77, 2798–2800 (2000). [CrossRef]

53. S. Carney and J. Schotland, “Near-field tomography,” MSRI Ser. Math. Appl. 47, 133–168 (2003).

54. H. Ammari, J. Garnier, and K. Solna, “Resolution and stability analysis in full-aperture, linearized conductivity and wave imaging,” Proc. Am. Math. Soc. 141, 3431–3446 (2013). [CrossRef]

55. H. Ammari, J. Garnier, and K. Solna, “Partial data resolving power of conductivity imaging from boundary measurements,” SIAM J. Math. Anal. 45, 1704–1722 (2013). [CrossRef]

56. PHG (Parallel Hierarchical Grid), http://lsec.cc.ac.cn/phg/.

57. P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L’Excellent, “A fully asynchronous multifrontal solver using distributed dynamic scheduling,” SIAM J. Matrix Anal. Appl. 23, 15–41 (2001). [CrossRef]

58. P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, and S. Pralet, “Hybrid scheduling for the parallel solution of linear systems,” Parallel Comput. 32, 136–156 (2006). [CrossRef]

Inverse diffraction grating of Maxwell’s equations in biperiodic structures

Abstract

1. Introduction

2. Model problem

2.1. Maxwell’s equations

2.2. Transparent boundary condition

2.3. Reduced model problem

3. Transformed field expansion

3.1. Change of variables

3.2. Power series solution

3.3. Zeroth order term

3.4. First order term

4. Reconstruction formula

5. Numerical experiment

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Tables (4)

Equations (87)

Optics Express

ε	0.2λ	0.1λ	0.05λ	0.025λ
e	8.50 × 10⁻¹	5.27 × 10⁻¹	1.72 × 10⁻¹	9.86 × 10⁻²

h	0.4λ	0.3λ	0.2λ	0.1λ
e	8.63 × 10⁻¹	8.59 × 10⁻¹	1.88 × 10⁻¹	1.20 × 10⁻¹