Profile reconstruction of a local defect in a groove structure and the theoretical limit under the vector diffraction theory

Jun-ichiro Sugisaka; Takashi Yasui; Koichi Hirayama

doi:10.1364/OE.404067

1. Introduction

Advances in nanofabrication technology have led to the miniaturization of semiconductor devices. Dust adhesion and contamination during the fabrication process cause defects in the products. Even the slightest fabrication error can cause severe damage to device operation in extremely integrated circuits. Defect detection in products, semiconductor wafers in each fabrication step, and photomasks in the lithography system are necessary to enhance yields. An electron microscope and an atomic-force microscope provide sufficient resolution for defect detection. However, a substantial amount of time is necessary for scanning the entire device. In contrast, while an optical microscope can measure a wide range of wafer surfaces, the resolution of an image is restricted by the diffraction limit. Under the scalar diffraction theory, the amplitude and phase distribution of a light wave just over the sample are proportional to the reflectance and the surface profile of the sample, respectively. The limit of measurement is simply described as the highest spatial frequency in the sample that can be imaged. For example, Abbe’s criterion provides a limit resolution of $\lambda /2$ with the maximum numerical aperture (NA) of a dry objective lens.

Although a microscope cannot produce a sample image, the information on the sample structure can be captured as scattered or diffracted waves. The problem of reconstructing a sample from a scattered wave is called the inverse scattering problem. So far, various methods have been proposed to solve this problem. Optical diffraction tomography can provide a cross-sectional image of the sample with low-refractive-index contrast. In this method, the relationship between sample structure and scattered field is represented as an approximation, such as the Born and Rytov approximation [1–4]. When the sample contains high-refractive-index interfaces, an iterative method is employed. In this approach, we first define an objective function that represents the error between the measured scattered waves from the sample and the computed scattered waves for the estimated structure. The sample structure is then reconstructed to minimize the object function. Some studies have reconstructed a sample profile beyond the diffraction limit. In acoustic engineering, a particle shape has been reconstructed from the phaseless data of a scattered wave [5,6]. In this method, the integral equations provide a scattered field, and a minimum of the object function is found by the Newton method. This approach has been applied to the inverse electromagnetic-field-scattering problem [7,8]. The uniqueness of the results reconstructed from phaseless far-field data has also been discussed [9,10].

Reconstruction of the profiles of a structure on a planar substrate has also been studied. This technique is known as optical diffraction microscopy or optical profilometry. For example, Kress et al. reconstructed the profile of a local protrusion on a perfectly conductive plane with boundary-integral equations of the Helmholtz equation and Newton method [11]. Garnaes et al. determined the grating parameters from diffracted waves and measured them using an atomic-force microscope [12]. Hansen et al. applied a renormalized scalar diffraction theory to determine the grating parameter [13]. Karamehmedović et al. estimated the defect in a finite-sized grating. The scattered wave was computed using the auxiliary source method [14]. They also reconstructed the profile of a platinum submicron wire on a Si substrate using a scattered field expression with an approximated Green’s function [15]. Afifi et al. reconstructed the sample profiles with a rough surface whose height was much less than the illuminating wavelength [16]. Ishikawa et al. identified the nanostructure on a substrate by matching the observed data with the template data containing scattered fields of various structures [17]. Maire et al. reconstructed the shapes of submicron rods on a flat substrate [18]. The scattered fields were computed using a volume-integral equation. This method was also experimentally demonstrated [19]. Zang et al. expanded this approach to three-dimensional reconstruction by introducing the polarization characteristics of the scattered field [20]. Zhang et al. proposed an algorithm for imaging a locally rough surface with boundary-integral equations [21,22].

In this study, we reconstruct a more complicated structure, the profile of a local defect in an infinite number of periodically aligned groove structures. The defect is a deformed groove that can be seen in bio-fabricated devices. Because the groove period and depth are comparable to the illuminating wavelength, strong multiple scattering occurs over the entire device surface. Approximations such as scalar-wave approximation and the Born and Rytov approximations cannot be employed for computing the scattered and diffracted fields and for evaluating the objective function. We have previously proposed the difference-field boundary element method (DFBEM) [23]. The DFBEM computes the scattered wave from the complicated surface relief structure with low computational resources. We arrange the DFBEM for the inverse scattering problem and reconstruct the defect profile of the sample.

Unlike the resolution theory under scalar diffraction approximation, such as Abbe’s criterion and Rayleigh’s criterion, the resolution limit of the iterative method is less discussed. As the studies introduced above have demonstrated, an iterative method can reconstruct a sample profile beyond the resolution limit of an optical imaging system. We also discuss the profile shape that can be reconstructed under the given observation system and noise level. This discussion can be used as a generalized criterion for the reconstruction limit of the iterative method that can be used to estimate if the profile can be reconstructed for various sample shapes.

In the following section, we derive the integral equation of DFBEM for the inverse problem. We also present an optical system to acquire the scattered wave, and its modeling for numerical analysis. In Section 3, we numerically simulate the defect reconstruction. We show the reconstruction results as a relation between the reconstruction accuracy and the observation system conditions such as NA and noise level. Finally, we discuss the limit and accuracy estimation of profile reconstruction.

2. Principle of profile reconstruction

In the reconstruction process, we set an object function that evaluates the error between experimentally observed data of the sample and numerically computed data of the estimated structure. The estimated sample structure is iteratively modified to minimize the objective function. The observed data is acquired by the optical system depicted in Fig. 1. Electromagnetic fields and all structures, including the sample and lenses, are invariant along the $z$-axis. Therefore, we consider only two-dimensional field distribution in the $x$-$y$ plane. The sample is placed so that the surface is on $y = 0$. This optical system is equivalent to an optical microscope. The field at $y = {y_\textrm{o}}{\; }({ > 0} )$ is imaged on the observation line with the objective lens ${\textrm{L}_1}$ and imaging lens ${\textrm{L}_2}$, whose focal lengths are ${F_1}$ and ${F_2}$, respectively. Note that the defect in the sample cannot be imaged clearly because the sample surface is out of focus and the defect is less than the illuminating wavelength, $\lambda $. The defect groove is located on the optical axis, $x = 0$. This alignment is possible by referring to the defocused image of the sample.

Fig. 1. Optical system to acquire observed data ${\psi _e}({\boldsymbol r} )$ for defect reconstruction. The objective lens ${\textrm{L}_1}$ and the imaging lens ${\textrm{L}_2}$ image the field distribution over the sample grating on the observation line.

Download Full Size | PDF

The sample is a one-dimensional diffraction grating. All grooves have a rectangular profile with a period T. The groove shape and the period are known parameters. In contrast, the unknown defect profile $y(x )$ is expressed as

(1)$$y(x )= \mathop \sum \limits_{n = 1}^N {u_n}{U_n}(x ),$$

where ${\boldsymbol u} = {\left[ {\begin{array}{cccc} {{u_1}}&{{u_2}}& \cdots &{{u_N}} \end{array}} \right]}^{{\textrm{t}}}$ is a structural vector (unknown). The symbol $\textrm{t}$ denotes the vector transposition. ${U_n}(x )$ $({n = 1,2, \cdots N} )$ are the basis functions (known). The estimation of the defect profile is then reduced to a problem of determining the unknown vector ${\boldsymbol u}$.

We define the object function as

(2)$$h({\boldsymbol u} )\equiv \frac{1}{2}\mathop \sum \limits_i {|{{\psi_g}({{\boldsymbol u},{{\boldsymbol r}_i}} )- {\psi_e}({{{\boldsymbol r}_i}} )} |^2},$$

where ${\psi _e}({{{\boldsymbol r}_i}} )$ is the experimentally measured scattered field distribution on the observation line. The data ${\psi _e}({{{\boldsymbol r}_i}} )$ is a complex value, assuming that it is obtained by an interferometer or any phase retrieval method. The function ${\psi _g}({{\boldsymbol u},{{\boldsymbol r}_i}} )$ is the scattered field distribution computed with an estimated profile ${\boldsymbol u}$. These scattered fields are components radiated from the defect. It is not necessary to consider the field component diffracted from the periodic grooves because it does not depend on ${\boldsymbol u}$, and is always eliminated by subtraction ${\psi _g}({{\boldsymbol u},{{\boldsymbol r}_i}} )- {\psi _e}({{{\boldsymbol r}_i}} )$. The vector ${{\boldsymbol r}_i}$ is a sampled position on the observation line. We assume that ${\boldsymbol u}$ that minimizes $h({\boldsymbol u})$ is the exact sample profile. An optimization algorithm (the conjugate gradient method) is employed to find the solution.

In the conjugate gradient method, the evaluation of $h({\boldsymbol u} )$ and its gradient

(3)$${\nabla _{\boldsymbol u}}h({\boldsymbol u})\equiv \left[ \begin{array}{cccc}{\frac{{\partial h({\boldsymbol u} )}}{{\partial {u_1}}}}&{\frac{{\partial h({\boldsymbol u})}}{{\partial {u_2}}}}& \cdots &{\frac{{\partial h({\boldsymbol u} )}}{{\partial {u_N}}}}\end{array}\right]^\textrm{t}$$

is necessary. They are performed with the DFBEM using the following procedure. As illustrated in Fig. 2, we define the boundaries on the sample ${C_0}$, ${C_1}$, and ${C_2}$, which are the in tegral paths of the boundary-integral equations. The infinite boundary ${C_0}$ is the sample surface, except for the defect. The boundary ${C_1}$ is the original sample surface that disappears because of the defect. The boundary ${C_2}$ is the defect surface whose shape is expressed by Eq. (1). The region inside the sample that is surrounded by ${C_0}$ and ${C_1}$ is denoted as ${S_1}$. The region outside of the sample surfaces ${C_0}$ and ${C_2}$ is denoted as ${S_2}$. The defect region is surrounded by ${C_1}$ and ${C_2}$, and is denoted as ${S_3}$. The refractive indices inside ${S_p}$ are denoted by ${n_p}$ ($p = 1,\; 2,\; 3$). The electric field of the illumination is oriented in the $z$-direction ($s$-polarized). In the DFBEM integral equations, we treat only the $z$-component of the electric field, and it is denoted by a variable f. The incident wave is a plane wave that illuminates the sample from the ${S_2}$ side with an incident angle of $\theta $. The incident wave contains a single frequency, and the temporal factor $\textrm{exp} ({j\omega t} )$ that is common to the incident and scattered fields is not described in the following equations, where j, $\omega $, and t are an imaginary unit, an angular frequency, and a time variable, respectively. The fields in the regions ${S_1}$, ${S_2}$, and ${S_3}$ (including those on ${C_1}$ and ${C_2}$) are denoted by ${f_1}$, ${f_2}$, and ${f_3}$. ${f_1}$ and ${f_2}$ are decomposed into ${f_1} = {f_{01}} + \Delta {f_1}$ and ${f_2} = {f_{02}} + \Delta {f_2}$. The components ${f_{01}}$ and ${f_{02}}$ are fields of the defect-free (perfectly periodic) sample. The other components $\Delta {f_1}$ and $\Delta {f_2}$ are fields of difference between the perfectly periodic and defective structures, which correspond to the scattered field radiated from the defect. The components ${f_{0p}}$ ($p = 1,2$) satisfy the Bloch (Floquet) theorem:

(4)$${f_{0p}}({x + T,y} )= {f_{0p}}({x,y} )\exp [{ - j{n_2}k\sin (\theta )T} ],$$

where k is the wavenumber of the incident wave in vacuum. The computation of ${f_{0p}}$ can be performed by various methods such as the boundary element method, finite-difference time-domain method, finite element method, and rigorous coupled-wave analysis (Fourier modal method).

Fig. 2. Boundaries on the defective grating: ${C_0}$ is the dielectric interface of the correct groove shape, ${C_1}$ is the dielectric interface that disappears because of the defect, and ${C_2}$ is the dielectric interface on the defect.

Download Full Size | PDF

After obtaining ${f_{0p}}$, we can compute $\varDelta {f_1}$ and its normal derivatives $\partial \varDelta {f_1}/\partial n$ on ${C_0}$, ${f_3}$ and $\partial {f_3}/\partial n$ on ${C_1}$, and ${f_3}$ and $\partial {f_3}/\partial n$ on ${C_2}$ by solving the following integral equations.

(5)$$\frac{1}{2}\varDelta {f_1}({\boldsymbol r} )={-} A_{1,2}^{C_0}\varDelta {\phi _1} - A_{1,3}^{\prime C_1}{\phi _3} + A_{1,2}^{{C_1}}{\phi _1}{\; }({{\boldsymbol r} \in {C_0}} ),$$

(6)$$\frac{1}{2}{f_3} - \frac{1}{2}{f_1}({\boldsymbol r} )={-} A_{1,2}^{{C_0}}\varDelta {\phi _1} - A_{1,3}^{\prime C_1}{\phi _3} + A_{1,2}^{{C_1}}{\phi _1}{\; }({{\boldsymbol r} \in {C_1}} ),$$

(7)$$\frac{1}{2}\varDelta {f_1}({\boldsymbol r} )= A_{2,1}^{\prime C_0}\Delta{\phi _1} + A_{2,3}^{\prime C_2}{\phi_3} - A_{2,1}^{\prime C_1}{\phi _1}{\; }({{\boldsymbol r} \in {C_0}} ),$$

(8)$$\frac{1}{2}{f_3}(r )- {f_2}({\boldsymbol r} )= A_{2,1}^{\prime C_0}\Delta{\phi _1} + A_{2,3}^{\prime C_2}{\phi _3} - A_{2,1}^{\prime C_1}{\phi _1}{\; }({{\boldsymbol r} \in {C_2}} ),$$

(9)$$\frac{1}{2}{f_3}({\boldsymbol r} )= A_{3,1}^{{C_1}}{\phi _3} - A_{3,2}^{{C_2}}{\phi _3}{\; }({{\boldsymbol r} \in {C_1}} ),$$

(10)$$\frac{1}{2}{f_3}({\boldsymbol r} )= A_{3,1}^{{C_1}}{\phi _3} - A_{3,2}^{{C_2}}{\phi _3}{\; }({{\boldsymbol r} \in {C_2}} ),$$

where the integral terms $A_{p,q}^C{\phi _p}$ and $A_{p,q}^{\prime C}{\phi _q}$ are defined as

(11)$$\begin{array}{c} {A_{p,q}^C{\phi _p} \equiv \mathop \int \nolimits_C \left\{ {{G_p}({{\boldsymbol r};{\boldsymbol r^{\prime}}} )\frac{{\partial {f_p}({{\boldsymbol r^{\prime}}} )}}{{\partial n^{\prime}}} - {f_p}({{\boldsymbol r^{\prime}}} )\frac{{\partial {G_p}({{\boldsymbol r};{\boldsymbol r^{\prime}}} )}}{{\partial n^{\prime}}}} \right\}dl^{\prime}} \end{array},$$

(12)$$\begin{array}{c} {A_{p,q}^{\prime C}{\phi _q} \equiv \mathop \int \nolimits_C \left\{ {{G_p}({{\boldsymbol r};{\boldsymbol r^{\prime}}} )\frac{{n_p^2}}{{n_q^2}}\frac{{\partial {f_q}({{\boldsymbol r^{\prime}}} )}}{{\partial n^{\prime}}} - {f_q}({{\boldsymbol r^{\prime}}} )\frac{{\partial {G_p}({{\boldsymbol r};{\boldsymbol r^{\prime}}} )}}{{\partial n^{\prime}}}} \right\}dl^{\prime},} \end{array}$$

(13)$${\phi _p} \equiv \left[ {{f_p}({{\boldsymbol r^{\prime}}} ),\frac{{\partial {f_p}({{\boldsymbol r^{\prime}}} )}}{{\partial n^{\prime}}}} \right].$$

The variables $({n,l,z} )$ are components of the local coordinate system on each boundary, l is the tangential direction (indicated by arrows on the boundaries in Fig. 2), and n is in the normal direction. In Eqs. (11) and (12), ${G_p}$ is a Green function defined as

(14)$${G_p}({{\boldsymbol r^{\prime}},{\boldsymbol r}} )={-} \frac{j}{4}H_0^{(2 )}({{n_p}k|{{\boldsymbol r^{\prime}} - {\boldsymbol r}} |} ),$$

where $H_0^{(2 )}$ is the zero-order Hankel function of the second kind.

In the numerical computation, the integral path C and fields ${\phi _p}$ are discretized with boundary elements. For example, the integral path C is discretized with points ${{\boldsymbol r}_i}$ on $C$. The boundary elements are short straight lines between ${{\boldsymbol r}_i}$ and ${{\boldsymbol r}_{i + 1}}$. The position on the element ${\boldsymbol r^{\prime}}$ is represented by a parameter $s^{\prime}$ ($- 1 \le s^{\prime} \le 1$).

(15)$${\boldsymbol r}({s^{\prime}} )= \frac{{s^{\prime} + 1}}{2}{{\boldsymbol r}_{i + 1}} - \frac{{s^{\prime} - 1}}{2}{{\boldsymbol r}_i},$$

and the fields f and $\partial f/\partial n$ on the element are also expressed by

(16)$$f({{\boldsymbol r}({s^{\prime}} )} )= \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}{f_{i,1}} - \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}{f_{i,2}},$$

(17)$$\frac{{\partial f({{\boldsymbol r}({s^{\prime}} )} )}}{{\partial n}} = \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}{g_{i,1}} - \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}{g_{i,2}},$$

where ${s_1}$ and ${s_2}$ are ${\pm} 0.57735$, respectively [24]. Finally, $A_{p,q}^C{\phi _p}$ is a linear combination of ${f_{i,1}}$, ${f_{i,2}}$, ${g_{i,1}}$, and ${g_{i,2}}$ for all elements in C. Substituting ${\boldsymbol r}({{s_m}} )$ $({m = 1,2} )$ for each element ${\boldsymbol r}$ in Eqs. (5)–(10), the integral term is expressed as a product of the matrix $A_{p,q}^C$, which consists of the integral of the Green function, and the vector ${\phi _p}$, which consists of the field parameters ${f_{i,1}}$, ${f_{i,2}}$, ${g_{i,1}}$, and. ${g_{i,2}}$

(18)$$\begin{aligned} A_{p,q}^C{\phi _p}& \equiv \mathop \sum \limits_i \left[ {\frac{{{L_i}}}{2}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}{G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )ds^{\prime}} \right]{g_{i,1}}\\& - \mathop \sum \limits_i \left[ {\frac{{{L_i}}}{2}\mathop \int \limits_{ - 1}^1 \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}{G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )ds^{\prime}} \right]{g_{i,2}}\\& - \mathop \sum \limits_i \left[ {\frac{{{L_i}}}{2}\mathop \int \limits_{ - 1}^1 \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}\frac{{\partial {G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial n^{\prime}}}ds^{\prime}} \right]{f_{i,1}}\\& + \mathop \sum \limits_i \left[ {\frac{{{L_i}}}{2}\mathop \int \limits_{ - 1}^1 \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}\frac{{\partial {G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial n^{\prime}}}ds^{\prime}} \right]{f_{i,2}}, \end{aligned}$$

where ${L_i}$ is the length of the boundary elements ${L_i} = |{{{\boldsymbol r}_i} - {{\boldsymbol r}_{i + 1}}} |$. Equations (5)–(8) contain integral terms with an infinite integral path ${C_0}$. In this term, the integrand $\varDelta \phi $ rapidly decreases to zero as the distance from the defect increases. Hence, the integral path ${C_0}$ can be truncated near the defect. After solving Eqs. (5)–(10), we can compute the difference-field at the object plane ($y = {y_\textrm{o}}$) as

(19)$$\varDelta {f_2}({{\boldsymbol u},{\boldsymbol r}} )= A_{2,1}^{\prime C_0}\Delta{\phi _1} + A_{2,3}^{\prime C_2}{\phi _3} - A_{2,1}^{\prime C_1}{\phi _1}\; \; ({{\boldsymbol r} \in {\textrm{S}_2}} ).$$

Next, the field distribution at the back focal plane of the objective lens $g({{\boldsymbol u},{\xi_n}} )$ is computed.

(20)$$g({{\boldsymbol u},{{\boldsymbol \rho }_n}} )= \mathop \sum \limits_{i ={-} {N_\textrm{i}}}^{{N_\textrm{i}}} \varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_i}} )\exp \left( {jk\frac{{{\xi_n}}}{{\sqrt {F_1^2 + \xi_n^2} }}{x_i}} \right),$$

where ${{\boldsymbol r}_i}\; ({i ={-} {N_\textrm{i}}, - {N_\textrm{i}} + 1,\; \cdots ,{N_\textrm{i}}} )$ is the sampling position on $y = {y_\textrm{o}}$, and ${x_i}$ is the $x$-coordinate of ${{\boldsymbol r}_i}$. Equation (20) holds if we use an aberration-free lens. The vector$\; {{\boldsymbol \rho }_n}\; ({n ={-} {N_f}, - {N_f} + 1,\; \cdots ,{N_f}} )$ is a sampling position on the focal plane and ${\xi _n}$ is the $x$-coordinate of ${{\boldsymbol \rho }_n}$. The NA of the objective lens $\nu $ can be considered as a pupil function

(21)$$p({{\xi_n}} )= \left\{ {\begin{array}{c} {1\; \left( {\frac{{{\xi_n}}}{{\sqrt {F_1^2 + \xi_n^2} }} < \nu } \right)}\\ {0\; ({\textrm{otherwise}} )} \end{array}.} \right.$$

This function removes the high-spatial-frequency and evanescent field, which cannot be propagated through the observation system. The field distribution on the image plane ${\psi _g}({{\boldsymbol u},{\boldsymbol r}_m^{\boldsymbol{\prime}}} )$ $({m ={-} {N_\textrm{o}}, - {N_\textrm{o}} + 1,\; \cdots ,{N_\textrm{o}}} )$ is computed as

(22)$${\psi _g}({{\boldsymbol u},{\boldsymbol r}_m^{\boldsymbol{\prime}}} )= \mathop \sum \limits_{n ={-} {N_f}}^{{N_f}} g({{\boldsymbol u},{\xi_n}} )p({{\xi_n}} )\exp \left( {jk\frac{{x_m^{\prime}}}{{\sqrt {F_2^{2} + x^{\prime}{2}_m} }}{\xi_n}} \right).$$

The computation of ${\psi _g}({{\boldsymbol u},{\boldsymbol r}_m^{\prime}} )$ from $\varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_i}} )$ is a linear operation and can be expressed as

(23)$${{\boldsymbol \psi }_g}({\boldsymbol u} )= B\varDelta {{\boldsymbol f}_2}({\boldsymbol u} ),$$

where B is a matrix, and ${{\boldsymbol \psi }_g}$ and $\varDelta {{\boldsymbol f}_2}$ are vectors defined as

(24)$${{\boldsymbol \psi }_g}({\boldsymbol u} )= \left[ {\begin{array}{cccc} {{\psi_g}({{\boldsymbol u},{\boldsymbol r}_{ - {N_\textrm{o}}}^{\boldsymbol{\prime}}} )}&{{\psi_g}({{\boldsymbol u},{\boldsymbol r}_{ - {N_\textrm{o}} + 1}^{\boldsymbol{\prime}}} )}& \cdots &{{\psi_g}({{\boldsymbol u},{\boldsymbol r}_{{N_\textrm{o}}}^{\boldsymbol{\prime}}} )} \end{array}} \right],$$

(25)$$\varDelta {{\boldsymbol f}_2}({\boldsymbol u} )= \left[ {\begin{array}{cccc} {\varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_{ - {N_\textrm{i}}}}} )}&{\varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_{ - {N_\textrm{i}} + 1}}} )}& \cdots &{\varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_{{N_\textrm{i}}}}} )} \end{array}} \right].$$

The next step computes the gradient of the object function ${\nabla _{\boldsymbol u}}h({\boldsymbol u} )$. The adjoint-variable method [25,26] and a direct derivative can be used for this purpose. In this study, we apply the latter method for the analysis described in Section 4. From Eq. (2), the $l$-th component of ${\nabla _{\boldsymbol u}}h({\boldsymbol u} )$ is

(26)$$\frac{{\partial h({\boldsymbol u} )}}{{\partial {u_l}}} = \mathop \sum \limits_i \textrm{Re}\left\{ {[{\psi_g^\ast ({{\boldsymbol u},{{\boldsymbol r}_i}} )- \psi_e^\ast ({{{\boldsymbol r}_i}} )} ]\frac{{\partial {\psi_g}({{\boldsymbol u},{{\boldsymbol r}_i}} )}}{{\partial {u_l}}}} \right\}.$$

Because the matrix B does not depend on ${\boldsymbol u}$, the derivative of ${\psi _g}({{\boldsymbol u},{r_i}} )$ is given by the derivative of $\varDelta {{\boldsymbol f}_2}({\boldsymbol u} )$

(27)$$\frac{{\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )}}{{\partial {u_l}}} = B\frac{{\partial \varDelta {{\boldsymbol f}_2}({\boldsymbol u} )}}{{\partial {u_l}}},$$

where $\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )/\partial {u_l}$ and $\partial {{\boldsymbol f}_2}({\boldsymbol u} )/\partial {u_l}$ are vectors with $i$-th component $\partial {\psi _g}({{\boldsymbol u},{{\boldsymbol r}_i}} )/\partial {u_l}$ and $\partial {f_2}({{\boldsymbol u},{{\boldsymbol r}_i}} )/\partial {u_l}$, respectively. Based on Eq. (19), the derivative of $\varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_i}} )$ is

(28)$$\frac{{\partial \varDelta {f_2}({{\boldsymbol u},{{\boldsymbol r}_i}} )}}{{\partial {u_l}}} = A^{\prime C_0}_{2,1}\frac{{\partial \varDelta {\phi _1}}}{{\partial {u_l}}} + \frac{{\partial A^{\prime C_2}_{2,3}}}{{\partial {u_l}}}{\phi _3} + A^{\prime C_2}_{2,3}\frac{{\partial {\phi _3}}}{{\partial {u_l}}}.$$

According to Eq. (18), the derivative of the integral operator for ${C_2}$ is

(29)$$\begin{aligned} \frac{{\partial A_{p,q}^{{C_2}}}}{{\partial {u_l}}}{\phi _p} &\equiv \mathop \sum \limits_i \left[ {\frac{1}{2}\frac{{\partial {L_i}}}{{\partial {u_l}}}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}{G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )ds^{\prime}} \right.\\ &\left. { + \frac{{{L_i}}}{2}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}\frac{{\partial {G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial {u_l}}}ds^{\prime}} \right]{g_{i,1}}\\ &- \mathop \sum \limits_i \left[ {\frac{1}{2}\frac{{\partial {L_i}}}{{\partial {u_l}}}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}{G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )ds^{\prime}} \right.\\ &\left. { + \frac{{{L_i}}}{2}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}\frac{{\partial {G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial {u_l}}}ds^{\prime}} \right]{g_{i,2}}\\ &- \mathop \sum \limits_i \left[ {\frac{1}{2}\frac{{\partial {L_i}}}{{\partial {u_l}}}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}\frac{{\partial {G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial n^{\prime}}}ds^{\prime}} \right.\\ &\left. { + \frac{{{L_i}}}{2}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_2}}}{{{s_1} - {s_2}}}\frac{{{\partial^2}{G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial {u_l}\partial n^{\prime}}}ds^{\prime}} \right]{f_{i,1}}\\ &+ \mathop \sum \limits_i \left[ {\frac{1}{2}\frac{{\partial {L_i}}}{{\partial {u_l}}}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}\frac{{\partial {G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial n^{\prime}}}ds^{\prime}} \right.\\ &\left. { + \frac{{{L_i}}}{2}\mathop \int \nolimits_{ - 1}^1 \frac{{s^{\prime} - {s_1}}}{{{s_2} - {s_1}}}\frac{{{\partial^2}{G_p}({{\boldsymbol r}({s^{\prime}} );{\boldsymbol r}({{s_m}} )} )}}{{\partial {u_l}\partial n^{\prime}}}ds^{\prime}} \right]{f_{i,2}}. \end{aligned}$$

Variation in ${u_l}$ changes the defect shape $y(x )$, which leads to the displacement of the y-coordinate of the element nodes on ${C_2}$, as well as the position in the boundary element ${\boldsymbol r}({s^{\prime}} )$ and the element length ${L_i}$. The specific expressions for the derivatives of ${L_i}$, ${G_p}$, and $\partial {G_p}/\partial n^{\prime}$ are described in Ref. [8]. The derivatives of $\varDelta {\phi _1}$ and ${\phi _3}$ are also required to compute Eq. (28). These parameters are provided by solving the derivatives of Eqs. (5)–(10):

(30)$$\frac{1}{2}\frac{{\partial \varDelta {f_1}(r )}}{{\partial {u_l}}} + A_{1,2}^{{C_0}}\frac{{\partial \varDelta {\phi _1}}}{{\partial {u_l}}} + A^{\prime C_1}_{1,3}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} = 0\; ({r \in {C_0}} ),$$

(31)$$\frac{1}{2}\frac{{\partial {f_3}}}{{\partial {u_l}}} + A_{1,2}^{{C_0}}\frac{{\partial \varDelta {\phi _1}}}{{\partial {u_l}}} + A^{\prime C_1}_{1,3}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} = 0\; ({r \in {C_1}} ),$$

(32)$$\frac{1}{2}\frac{{\partial \varDelta {f_1}(r )}}{{\partial {u_l}}} - A^{\prime C_0}_{2,1}\frac{{\partial \varDelta {\phi _1}}}{{\partial {u_l}}} - A^{\prime C_2}_{2,3}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} = \frac{{\partial A^{\prime C_2}_{2,3}}}{{\partial {u_l}}}{\phi _3}\; \; ({r \in {C_0}} ),$$

(33)$$\frac{1}{2}\frac{{\partial {f_3}(r )}}{{\partial {u_l}}} - A^{\prime C_0}_{2,1}\frac{{\partial \varDelta {\phi _1}}}{{\partial {u_l}}} - A^{\prime C_2}_{2,3}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} = \frac{{\partial A^{\prime C_0}_{2,1}}}{{\partial {u_l}}}\varDelta {\phi _1} + \frac{{\partial A^{\prime C_2}_{2,3}}}{{\partial {u_l}}}{\phi _3} - \frac{{\partial A^{\prime C_1}_{2,1}}}{{\partial {u_l}}}{\phi _1}\; \; ({r \in {C_2}} ),$$

(34)$$\frac{1}{2}\frac{{\partial {f_3}(r )}}{{\partial {u_l}}} - A_{3,1}^{{C_1}}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} + A_{3,2}^{{C_2}}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} ={-} \frac{{\partial A_{3,2}^{{C_2}}}}{{\partial {u_l}}}{\phi _3}\; \; ({r \in {C_1}} ),$$

(35)$$\frac{1}{2}\frac{{\partial {f_3}(r )}}{{\partial {u_l}}} - A_{3,1}^{{C_1}}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} + A_{3,2}^{{C_2}}\frac{{\partial {\phi _3}}}{{\partial {u_l}}} = \frac{{\partial A_{3,1}^{{C_1}}}}{{\partial {u_l}}}{\phi _3} - \frac{{\partial A_{3,2}^{{C_2}}}}{{\partial {u_l}}}{\phi _3}\; \; ({r \in {C_2}} ).$$

The terms on the left-hand sides of Eqs. (30)–(35) are the product of the coefficients and unknown parameters $\partial \phi /\partial {u_l}$; the right-hand sides are constant terms. When solving as a set of simultaneous equations numerically, the coefficient matrix does not depend on l. Hence, $\partial \phi /\partial {u_l}$ for each l can be efficiently solved by the LU decomposition of the coefficient matrix.

3. Reconstruction simulation

The reconstruction process was validated through numerical simulation. We prepared the observation data ${\psi _e}({\boldsymbol r} )$ by numerical simulation, and the defect profile was reconstructed from ${\psi _e}({\boldsymbol r} )$ according to the process described in Section 2. The sample structure is Fig. 3(a) Profile of the sample grating. (b) Profile of the initial guess.illustrated in Fig. 3(a). Hereafter, all parameters of length, such as the dimensions of the sample shape and the sampling interval of the field, are normalized by the illuminating wavelength, $\lambda $. As long as the ratios of these parameters to $\lambda $ and the refractive indices are unchanged, the results also remain unchanged, regardless of the wavelength band. The period of the groove $T = 1.2\lambda $, and the width and depth of the groove are $0.5T$ and $0.3T$, respectively. The structural vector of the sample is represented by a four-dimensional vector:

(36)$${{\boldsymbol u}_0} \equiv {\left[ {\begin{array}{cccc} {{u_1}}&{{u_2}}&{{u_3}}&{{u_4}} \end{array}} \right]^\textrm{t}} = {\left[ {\begin{array}{ccccc} { - 0.23T}&0&{ - 0.06T}&0 \end{array}} \right]^\textrm{t}},$$

and the basis function

(37)$${U_n}(x )= \sin \left[ {\frac{{2n\pi }}{T}\left( {x + \frac{T}{4}} \right)} \right].$$

Fig. 3. (a) Profile of the sample grating. (b) Profile of the initial guess.

Download Full Size | PDF

Since we know that the defect position is on the optical axis ($x = 0$), the defect region is $- T/4 \le x \le T/4$.

In the computation of ${\psi _e}({\boldsymbol r} )$ and the reconstruction process, the lengths of the boundary elements on ${C_0}$ and ${C_1}$ are set to be less than $T/10$. On ${C_2}$, it is set to $T/40$ such that the fine shape can be accurately reflected in the computed result. The refractive indices ${n_1}$, ${n_2}$, and ${n_3}$ are set to 1.5, 1.0, and 1.5, respectively. The incident angle of the illumination is set to $\theta = 0^\circ $. The object plane ${y_\textrm{o}}$ of the optical imaging system is set to ${y_\textrm{o}} = T$. The field distribution on the object plane is discretized with an interval of $0.05T$. The focal lengths ${F_1}$ and ${F_2}$ are set to $9.0 \times {10^3}T$ and $3.6 \times {10^5}T$, respectively. The NA is varied from 0.05–0.95. The field distribution on the image plane is discretized with an interval of $10T.$ Assuming that the observed data ${\psi _e}({{{\boldsymbol r}_i}} )$ contains noise in the actual measurement, normal random numbers with the variance of ${\sigma ^2}$ are added to the real and imaginary parts of ${\psi _e}({{{\boldsymbol r}_i}} )$. For example, when we set a noise level of $x$%, $\sigma $ for the real part is $0.01x\textrm{Re}[{{\psi_e}({{{\boldsymbol r}_i}} )} ]$, and that for the imaginary part is $0.01x\textrm{Im}[{{\psi_e}({{{\boldsymbol r}_i}} )} ]$.

The reconstruction process was performed by the conjugate gradient method with an initial guess of ${\boldsymbol u} = \mathbf{0}$, as illustrated in Fig. 3(b). After computing the conjugate gradient from ${\nabla _{\boldsymbol u}}h({\boldsymbol u} )$, we found an ${\boldsymbol u}$ that provides the minimum $h({\boldsymbol u} )$ in the direction of the conjugate gradient with a golden-section search. This process was iteratively performed in 20 steps. The reconstruction process was performed ten times, while changing the seed of the random numbers. The results were evaluated by field error ${E_f}$ and structural error ${E_s}$ as follows.

(38)$${E_f} \equiv 2h({\boldsymbol u} )= \mathop \sum \limits_i {|{{\psi_g}({{\boldsymbol u},{{\boldsymbol r}_i}} )- {\psi_e}({{{\boldsymbol r}_i}} )} |^2},$$

(39)$${E_s} \equiv {|{{\boldsymbol u} - {{\boldsymbol u}_0}} |^2}.$$

The former is an error in the field on the image plane that can be evaluated in practice. The latter is the difference between the sample and the reconstructed defect profile.

The results for each NA and noise level are plotted in Fig. 4. The field error was normalized to that of the initial guess, $E_f^{(0 )}$. The reconstructed defect profiles for some parameters are illustrated in Fig. 5. The exact defect profile is reconstructed for a noise level of 0%, regardless of the NA. When the noise level is more than 0%, the reconstruction accuracy decreases with an decrease in the NA. For example, the noise level must be less than 0.5% and the NA must be larger than 0.1 to achieve ${E_s}$ less than 10%. ${E_s}$ less than 1% can be achieved when the noise level is not more than 0.1%.

Fig. 4. (a) Field error and (b) Structural error for each NA and noise level.

Download Full Size | PDF

Fig. 5. Sample structures (gray region) and reconstructed profiles (solid line) of the defective grating under various conditions with NA and noise level.

Download Full Size | PDF

4. Discussion on the limit of reconstruction

We discuss the limits of the defect shape that can be reconstructed by the procedure in this study, so as to estimate the potential of reconstruction for various types of defect shapes. An accurate reconstruction result can be achieved when the information on the defect shape is transmitted to the observation line through the optical system. For example, when the defect is deformed to a certain shape, but the field on the observation line varies less than the noise level, the deformation cannot be detected. To discuss quantitatively, we derive a relationship between the variation in defect shape and the variation in field on the observation line. We approximate the small variation in the field ${{\boldsymbol \psi }_g}({\boldsymbol u} )$ because of the variation in defect shape $\varDelta {\boldsymbol u}$ as

(40)$${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )\simeq {{\boldsymbol \psi }_g}({\boldsymbol u} )+ \frac{{\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )}}{{\partial {\boldsymbol u}}}\varDelta {\boldsymbol u}.$$

The coefficient of $\mathrm{\varDelta }{\boldsymbol u}\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )/\partial {\boldsymbol u}$ is a matrix consisting of the column vectors $\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )/\partial {u_i}$

(41)$$\frac{{\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )}}{{\partial {\boldsymbol u}}} = \left[ {\begin{array}{cccc} {\frac{{\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )}}{{\partial {u_1}}}}&{\frac{{\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )}}{{\partial {u_2}}}}& \cdots &{\frac{{\partial {{\boldsymbol \psi }_g}({\boldsymbol u} )}}{{\partial {u_N}}}} \end{array}} \right].$$

From Eq. (27), ${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )$ can be expanded as

(42)$${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )\simeq {{\boldsymbol \psi }_g}({\boldsymbol u} )+ B\frac{{\partial \varDelta {{\boldsymbol f}_2}({\boldsymbol u} )}}{{\partial {\boldsymbol u}}}\varDelta {\boldsymbol u}$$

The coefficient of $\varDelta {\boldsymbol u}$ is represented by a matrix $B^{\prime}$ as

(43)$${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )\simeq {{\boldsymbol \psi }_g}({\boldsymbol u} )+ B^{\prime}\varDelta {\boldsymbol u}.$$

The components of vector ${{\boldsymbol \psi }_g}$ and matrix $B^{\prime}$ in Eq. (43) are complex numbers, while those of $\varDelta {\boldsymbol u}$ are real numbers. In the following equations, the components of ${{\boldsymbol \psi }_g}$ and $B^{\prime}$ are rearranged so that they consist of real numbers, by separately aligning the real and imaginary parts of the components in the odd and even rows, respectively. The real matrix $B^{\prime}$ is factorized by singular value decomposition as

(44)$${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )\simeq {{\boldsymbol \psi }_g}({\boldsymbol u} )+ WS{V^\textrm{t}}\varDelta {\boldsymbol u},$$

(45)$$W = \left[ {\begin{array}{ccc} {{{\boldsymbol w}_1}}&{{{\boldsymbol w}_2}}& \cdots \end{array}} \right],$$

(46)$$V = \left[ {\begin{array}{ccc} {{{\boldsymbol v}_1}}&{{{\boldsymbol v}_2}}& \cdots \end{array}} \right].$$

The matrices W and V are real orthogonal matrices, that is, ${\boldsymbol w}_n^\textrm{t}{{\boldsymbol w}_m} = {\boldsymbol v}_n^\textrm{t}{{\boldsymbol v}_m} = {\delta _{nm}}$, where ${\delta _{nm}}$ is a Kronecker delta. S is a real diagonal matrix whose $({n,n} )$ components ${s_n}$ (singular value) are arranged in descending order $({{s_n} \ge {s_{n + 1}}} )$. We define another structural vector ${\boldsymbol u^{\prime}}$ and the base functions $U_n^{\prime}(x )$ as

(47)$${\boldsymbol u^{\prime}} \equiv {V^\textrm{t}}{\boldsymbol u},$$

(48)$$U_n^{\prime}(x )\equiv \mathop \sum \limits_m {V_{mn}}{U_m}(x ),$$

where ${V_{mn}}$ is the $({m,n} )$ component of V. A small variation in the defect shape $\varDelta y(x )$ is expressed using the components of $\varDelta {\boldsymbol u^{\prime}}({ \equiv {V^\textrm{t}}\varDelta {\boldsymbol u}} )$ as

(49)$$\varDelta y(x )= \mathop \sum \limits_n \varDelta u_n^{\prime}U_n^{\prime}(x ).$$

Finally, ${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )$ is represented as

(50)$${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )\simeq {{\boldsymbol \psi }_g}({\boldsymbol u} )+ \mathop \sum \limits_n {{\boldsymbol w}_n}{s_n}\varDelta u_n^{\prime}.$$

The contribution of the shape variation$\; \varDelta u_n^{\prime}$ to ${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )$ is determined by ${{\boldsymbol w}_n}{s_n}$. The vector ${{\boldsymbol w}_n}$ represents the variation in field distribution on the observation line because of the variation $\varDelta u_n^{\prime}$. Because the vectors ${{\boldsymbol w}_n}$ $({n = 1,2, \cdots } )$ are orthogonal to each other, the difference in field distribution ${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )- {{\boldsymbol \psi }_g}({\boldsymbol u} )$ is also orthogonal when different components of ${\boldsymbol u^{\prime}}$ are varied. Hence, by comparing ${{\boldsymbol \psi }_g}({\boldsymbol u} )$ and ${{\boldsymbol \psi }_g}({{\boldsymbol u} + \varDelta {\boldsymbol u}} )$, we can uniquely identify $\varDelta {\boldsymbol u}$. Because the norm on ${{\boldsymbol w}_n}$ is unit, the magnitude of the variation in ${{\boldsymbol \psi }_g}$ because of the $n$-th component of ${\boldsymbol u^{\prime}}$ is determined by ${s_n}$. If one ${s_n}$ is very small than the other ${s_n}$, the variation in ${{\boldsymbol \psi }_g}({\boldsymbol u} )$ is easily buried in the noise floor. Examples of $U_n^{\prime}(x )$ and ${s_n}$ at ${\boldsymbol u}_0^{\boldsymbol{\prime}} \equiv {V^t}{{\boldsymbol u}_0}$ (sample structure) and ${\boldsymbol u^{\prime}} = \mathbf{0}$ are depicted in Fig. 6. As the singular value decreases, the corresponding $U_n^{\prime}(x )$ tends to contain a higher spatial frequency; finally, ${s_4}$ is more than 30 times smaller than ${s_1}$ at ${\boldsymbol u}_0^{\boldsymbol{\prime}}$. To identify the cause of this relation, we computed the field difference around the defect before and after each component of ${\boldsymbol u^{\prime}}$ was increased by $0.05T$ from ${\boldsymbol u^{\prime}} = \mathbf{0}$. As plotted in Fig. 7, each difference component is distributed as a spherical wave. The negative wave source distributed in the region changed from ${\textrm{S}_3}$ to ${\textrm{S}_2}$ by the variation in ${\boldsymbol u^{\prime}}$, and the positive source distributed in the region changed from ${\textrm{S}_2}$ to ${\textrm{S}_3}$. When $U_n^{\prime}(x )$ contains only a lower spatial frequency, the wave source with the same sign is distributed, resulting in a stronger wave source. In contrast, when $U_n^{\prime}(x )$ contains a higher frequency, the wave sources with opposite signs are alternately distributed. In this case, the waves radiated with different signs destructively interfere at the observation line, resulting in a smaller contribution to ${{\boldsymbol \psi }_g}$. However, the contribution cannot be zero because there is a displacement between the positive and negative sources. In Fig. 7(d), the positive sources are distributed at $y < 0$, and the negative ones are distributed at $y > 0$. At a distance far from $y = 0$, these scattered waves cannot be perfectly anti-phase and cannot be canceled. In contrast, as the scalar diffraction theory derives, if all wave sources with the opposite sign are infinitely distributed with such a small period on $y = 0$, the radiated waves behave like an evanescent wave that cannot arrive at the observation line. Such a case does not exist as long as we consider the surface deformation of the defect profile, and the deformation is restricted in a finite region.

Fig. 6. (a) Original base functions of the defect shape ${U_n}(x )$ $({n = 1,2,3,4} )$, transformed base functions, $U_n^{\prime}(x )$, and the corresponding singular values, for ${s_n}$ $({n = 1,2,3,4} )$ at (b) the initial guess ${\boldsymbol u^{\prime}} = 0$ and (c) the sample shape ${\boldsymbol u}_0^{\prime}$.

Download Full Size | PDF

Fig. 7. Difference-field distribution before and after varying the components of (a) ${\boldsymbol u}_1^{\prime}$, (b) ${\boldsymbol u}_2^{\prime}$, (c) ${\boldsymbol u}_3^{\prime}$, and (d) ${\boldsymbol u}_4^{\prime}$ by $0.05T$ from ${\boldsymbol u^{\prime}} = \mathbf{0}$.

Download Full Size | PDF

Based on the above discussion, we estimate the reconstruction accuracy from the singular values and the noise level of the observed data, ${\psi _e}({{{\boldsymbol r}_i}} )$. In the following discussion, we do not consider the case of converging the local minimum in the searching algorithm on the conjugate gradient method. We can avoid such a case by restarting the reconstruction process by changing the initial guess of ${\boldsymbol u}$. Here, we assume that the object function is minimal at ${{\boldsymbol u}_0} + \varDelta {\boldsymbol u}$, which is slightly shifted by $\varDelta {\boldsymbol u}$ from the sample shape ${{\boldsymbol u}_0}$ owing to the noise in the observed data. From Eqs. (2) and (50), the object function $h({{{\boldsymbol u}_0} + \varDelta {\boldsymbol u}} )$ is expressed as

(51)$$h({{{\boldsymbol u}_0} + \varDelta {\boldsymbol u}} )= \frac{1}{2}\mathop \sum \limits_i {\left|{{\psi_g}({{{\boldsymbol u}_0},{{\boldsymbol r}_i}} )+ \mathop \sum \limits_n {w_{n,i}}{s_n}\varDelta u_n^{\prime} - {\psi_e}({{{\boldsymbol r}_i}} )} \right|^2},$$

where ${w_{n,i}}$ is the $i$-th component of ${{\boldsymbol w}_n}$. Because the derivative of the object function with $\varDelta u_m^{\prime}$ ($m = 1,2, \cdots ,N$) is zero at ${{\boldsymbol u}_0} + \varDelta {\boldsymbol u}$

\frac{{\partial h({{f_1},{f_2} \cdots } )}}{{\partial \varDelta u_m^{\prime}}} = \mathop \sum \limits_i \left[ {\mathop \sum \limits_n {s_n}\varDelta u_n^{\prime}{w_{n,i}} - ({{\psi_g}({{{\boldsymbol u}_0},{{\boldsymbol r}_i}} )- {\psi_e}({{{\boldsymbol r}_i}} )} )} \right]{s_m}{w_{m,i}}

(52)$$= {s_m}\mathop \sum \limits_n \left[ {\mathop \sum \limits_i {w_{m,i}}{w_{n,i}}} \right]{s_n}\varDelta u_n^{\prime} - {s_m}\mathop \sum \limits_i {w_{m,i}}({{\psi_g}({{{\boldsymbol u}_0},{r_i}} )- {\psi_e}({{r_i}} )} )= 0.$$

If we assume ${s_m} > 0\; $($m = 1,2, \cdots ,N$), we obtain a set of simultaneous equations

(53)$$\mathop \sum \limits_n \left[ {\mathop \sum \limits_i {w_{m,i}}{w_{n,i}}} \right]{s_n}\varDelta u_n^{\prime} = \mathop \sum \limits_i {w_{m,i}}({{\psi_g}({{{\boldsymbol u}_0},{r_i}} )- {\psi_e}({{r_i}} )} ).$$

Because ${w_{n,i}}$ is a component of the orthogonal matrix W,

\mathop \sum \limits_n {\delta _{mn}}{s_n}\varDelta u_n^{\prime} = \mathop \sum \limits_i {w_{m,i}}({{\psi_g}({{{\boldsymbol u}_0},{r_i}} )- {\psi_e}({{r_i}} )} )

(54)$$\varDelta u_m^{\prime} = \frac{1}{{{s_m}}}\mathop \sum \limits_i {w_{m,i}}({{\psi_g}({{{\boldsymbol u}_0},{r_i}} )- {\psi_e}({{r_i}} )} ).$$

In general, when arguments ${x_n}$ ($n = 1,2, \cdots $) of the function

(55)$$y = {F_n}({{x_1},{x_2}, \cdots } )$$

vary with the variation in $\sigma _{x,n}^2$, the variance of y, $\sigma _y^2$ is given by

(56)$$\sigma _y^2 = \mathop \sum \limits_n {\left[ {\frac{{\partial F({{x_1},{x_2}, \cdots } )}}{{\partial {x_n}}}} \right]^2}\sigma _{x,n}^2.$$

Applying Eq. (54) to Eq. (56), we obtain the variance of. $\varDelta u_m^{\prime}$ as

(57)$$\sigma _{u,m}^2 = \mathop \sum \limits_n {\left[ {\frac{{{w_{m,n}}}}{{{s_m}}}} \right]^2}\sigma _{\psi ,n}^2,$$

where $\sigma _{\psi ,n}^2$ is the variance of ${\psi _g}({{{\boldsymbol u}_0},{{\boldsymbol r}_n}} )- {\psi _e}({{{\boldsymbol r}_n}} )$. This variance ${\sigma _{u,m}}$ is equivalent to the variance of the error of the $m$-th element of ${\boldsymbol u}$, $|{u_m^{\prime} - u_{0m}^{\prime}} |$. Because $\sigma _{\psi ,n}^2$ is from the noise in the observed data, Eq. (57) represents the range of error in $\varDelta u_m^{\prime}$, which is caused by the statistical characteristics of the noise in the observed data.

We computed ${\sigma _{u,m}}$ for each component of $\varDelta {\boldsymbol u^{\prime}}$ at ${\boldsymbol u}_0^{\boldsymbol{\prime}}$. Because the noise in ${\psi _e}({{{\boldsymbol r}_n}} )$ (noise level of $x$%) is generated from a normal random number such that the variance of ${\psi _e}({{{\boldsymbol r}_n}} )$ is given by

(58)$$\sigma _{\psi ,n}^2 = {[{0.01x{\psi_g}({{{\boldsymbol u}_0},{{\boldsymbol r}_n}} )} ]^2}.$$

The variance ${\sigma _{u,m}}$ calculated from Eqs. (57) and (58) is plotted in Fig. 8. We also plot the actual error $|{u_m^{\prime} - u_{0m}^{\prime}} |$ obtained by the reconstruction process described in Section 3. The ${\sigma _{u,m}}$ is located near the center of the distribution of the reconstruction results, and is valid as an estimated accuracy of $|{u_m^{\prime} - u_{0m}^{\prime}} |$.

Fig. 8. Error $|{u_m^{\prime} - u_{0m}^{\prime}} |$ for each element m, NA, and noise level ((a) 0.1% and (b) 1.0%) estimated from ${\sigma _{u,m}}$ and evaluated from the reconstruction results.

Download Full Size | PDF

As an application, we evaluated ${\sigma _{u,m}}$ for other sample structures (different refractive-index contrast and different aspect ratio of the surface relief structure). First, the refractive index ${n_1}({ = {n_3}} )$ was varied from 1.5 to 2.5 or 3.5, while ${n_2}({ = 1.0} )$ remained unchanged. The NA and the noise level were set to 0.6 and 1.0%, respectively. As listed in Table 1, ${\sigma _{u,3}}$ and ${\sigma _{u,4}}$ decrease with increasing the contrast, ${n_1}/{n_2}$. As Eq. (57) indicates, the smaller ${\sigma _{u,m}}$ is due to the greater singular value, ${s_m}$. The singular value, ${s_m}$ becomes greater when the intensity of the field difference for structural variation of $u_m^{\prime}$ is larger. In general, the intensity of the field difference for $u_3^{\prime}$ and $u_4^{\prime}$ is small because of destructive interference between neighboring negative and positive sources as plotted in Figs. 7(c) and 7(d). However, in a higher refractive-index material, the wavelength becomes shorter, and the interval between the negative and positive source is relatively extended. As a result, the destructive interference diminishes, and the field difference becomes larger.

Table 1. Reconstruction accuracy (${{\boldsymbol \sigma }_{{\boldsymbol u},{\boldsymbol m}}}$) for different refractive index of the sample material.

View Table

Second, we evaluated ${\sigma _{u,m}}$ for different aspect ratio of the sample. The groove depth d was varied from $0.1T$ to $2.0T$, while the groove width ($0.5T$) remained unchanged (the aspect ratio ranging from 0.2 to 4.0). The defect depth was also varied in proportion to d;

(59)$${{\boldsymbol u}_0} = {\left[ {\begin{array}{cccc} { - 0.23T \times \frac{d}{{0.3T}}}&0&{ - 0.06T \times \frac{d}{{0.3T}}}&0 \end{array}} \right]^\textrm{t}}\; .$$

The NA and the noise level were set to 0.6 and 1.0%, respectively. As plotted in Fig. 9, ${\sigma _{u,m}}$ becomes the largest at $d = 1.0T$. This characteristic is because the field difference for a small variation of ${{\boldsymbol u}_0}$ becomes smaller at $d = 1.0T$. As depicted in Fig. 10, the field difference radiates from the bottom of the defect, and it propagates through two different paths ${L_1}$ and ${L_2}$. If we approximate the path is parallel to the $y$-axis, there exists a phase difference $2\pi ({1.5 - 1.0} )l/\lambda $ on the sample surface, where l is the defect depth, $l \simeq 0.23T \times d/({0.3T} )$. When $d = 1.09T$, the phase difference becomes $\pi $, and the scattered wave destructively interferes, and the field difference becomes minimum on the observation line, resulting in a small singular value. At that time, ${\sigma _{u,m}}$ becomes larger.

Fig. 9. Estimated reconstruction accuracy, ${\sigma _{u,m}}$ that corresponds to $|{u_m^{\prime} - u_{0m}^{\prime}} |$ for each groove depth, h. The defect depth is also varied in proportion to d. The groove width, the NA, and the noise level are constant, which are $0.5T$, $0.6$, and $1.0\%$, respectively.

Download Full Size | PDF

Fig. 10. Propagation paths, ${L_1}$ and ${L_2}$ of the field difference for small variation of the defect shape (gray area).

Download Full Size | PDF

Finally, we extended the dimension of the structural vector from four to twelve to investigate the limit of reconstruction. The basis function is defined by Eq. (37). Because the structural vector can express a finer (higher spatial frequency) profile, the defect surface (boundary ${C_2}$) was discretized by 40 elements. The base function $U_n^{\prime}(x )\; ({ - T/4 \le x \le T/4} )$ and the singular values are plotted in Fig. 11. The base functions $U_1^{\prime}(x )$ and $U_2^{\prime}(x )$, whose singular values are more than 75, have two peaks near the edges $({x ={\pm} T/4} )$. The shapes are mirror-symmetric for $U_1^{\prime}(x )$ and anti-symmetric for $U_2^{\prime}(x )$ with respect to a line $x = 0$. Because these functions represent a structural variation at two points apart from each other, ranging near the Abbe’s limit $T/2$, they can be detected as a strong signal as a variation of ${{\boldsymbol \psi }_e}$. $U_4^{\prime}(x )$ represents a variation whose depth in the defect is uniformly changed. As depicted in Fig. 7(a), this component can also be detected as a strong signal. $U_5^{\prime}(x )$ and $U_6^{\prime}(x )$ are similar to the triangular functions. Considering that the incident wave reflects on their linear slope, the reflected wave is detected as a strong signal on the observed line. The estimated ${\sigma _{u,m}}$ values are plotted in Fig. 12. The estimated error exponentially increases with m, and finally, ${\sigma _{u,12}}$ is ${10^4}$ times larger than the first component ${\sigma _{u,1}}$. As for $m > 10$ with a noise level of 0.1%, and for $m > 7$ with a noise level of 5%, ${\sigma _{u,m}}$ is larger than the groove depth, $0.3T$. Reconstruction of these components is virtually impossible under these conditions.

Fig. 11. Transformed base functions, $U_n^{\prime}(x )$ for 12-dimensional structural vector ${\boldsymbol u}$ and the corresponding singular values with ${s_n}$ at the initial guess ${\boldsymbol u} = 0$.

Download Full Size | PDF

Fig. 12. Standard variation ${\sigma _{um}}$ for each element of 12-dimensional structural vector ${\boldsymbol u}$ for NA of 0.5 and noise level of 0.1%, 0.5%, 1.0%, and 5.0%, estimated from the singular value.

Download Full Size | PDF

The reconstruction algorithm can be applied to other optical measurement systems if we can compute the system output, ${\psi _g}({\boldsymbol r} )$, from the scattered wave $\varDelta {f_2}({\boldsymbol r} )$. If this relationship is expressed by a matrix as Eq. (23), we can also estimate the reconstruction accuracy. For example, the optical system depicted in Fig. 1 is replaced with a scanning near-field optical microscope, which provides high-resolution images of samples by acquiring the evanescent field on the sample surface by a scanning probe. For simplicity, the microscope system is approximated by a spatial filter. The matrix B in Eq. (23) was configured by Fourier and inverse Fourier transform and a low-pass filter whose cut-off spatial frequency is $NA/\lambda $, where $NA$ corresponds to the numerical aperture of the objective lens. When $NA > 1$, ${\psi _g}({{\boldsymbol u},{\boldsymbol r}} )$ and ${\psi _e}({{{\boldsymbol r}_i}} )$ contain the evanescent field. The reconstruction accuracy, ${\sigma _{u,m}}$ at ${\boldsymbol u}_0^{\boldsymbol{\prime}}$ for $NA < 10$ is plotted in Fig. 13. The estimated ${\sigma _{u,m}}$ slowly decreases with increasing NA. Because this sample, which consists of a transparent dielectric material, produces little evanescent field (see the defect surface in Fig. 7(d)), the effect on improving the accuracy is small.

Fig. 13. Estimated accuracy of the defect shape, ${\sigma _{u,m}}$, reconstructed from the field distribution on the sample surface. The lateral axis is the cut-off spatial frequency truncated by the measurement system, which is represented by the corresponding NA of the objective lens. When $NA > 1$ (gray area), the measured field distribution contains the evanescent field.

Download Full Size | PDF

5. Conclusion

We developed a reconstruction process for local defect in a periodic surface relief structure. The reconstruction process applies an optimization algorithm using the conjugate gradient method. We arranged the DFBEM boundary-integral equations for efficient computation of the objective function and its gradient. Because discretization in the region using mesh and grid is not necessary, mesh regeneration and grid resampling are also not necessary when the estimated defect shape is updated in the reconstruction process. The integral equations also represent a direct relationship between the deformation of the defect profile and the observed field data. This relationship leads to the precise analysis of possible reconstruction of a fine defect profile.

Regarding the numerical simulation of defect reconstruction, we found that the reconstruction for a fine defect profile, whose spatial frequency is beyond the Abbe’s limit, can be achieved with the existence of noise in the observed data. For a general case of profile reconstruction, we classified a profile shape that can be easily or hardly reconstructed using singular value decomposition, and also quantified the degree of difficulty of reconstruction to estimate the reconstruction accuracy. The easily reconstructed profiles are those whose depth is uniformly varied, which are deformed at two positions apart from each other, and whose slope is changed. The hardly reconstructed profile fluctuates in a narrow region. The latter hardly contributes to the field distribution on the observation line, and the signal is easily buried in the noise floor. Unlike the scalar diffraction theory, it is not simply classified by the spatial frequency of the sample. Even when the spatial frequency of the defect profile is beyond the diffraction limit, the scattered wave can propagate and be detected through an optical imaging system.

The observed data that is referred to in the reconstruction contains phase information of the scattered field. For more practical cases, reconstruction from phaseless data that consists of the intensity or amplitude of a field can also be achieved by only arranging the objective function. We will analyze the limit of reconstruction under this condition of observation is our future work.

Funding

Japan Society for the Promotion of Science (19K12012).

Disclosures

The authors declare no conflicts of interest.

References

1. A. Devaney, “A filtered backpropagation algorithm for diffraction tomography,” Ultrason. Imaging 4(4), 336–350 (1982). [CrossRef]

2. Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. Dasari, and M. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express 17(1), 266–277 (2009). [CrossRef]

3. V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope,” J. Microsc. 205(2), 165–176 (2002). [CrossRef]

4. K. Kim, J. Yoon, and Y. Park, “Simultaneous 3D visualization and position tracking of optically trapped particles using optical diffraction tomography,” Optica 2(4), 343–346 (2015). [CrossRef]

5. O. Ivanyshyn and R. Kress, “Identification of sound-soft 3D obstacles from phaseless data,” Inverse Prob. Imaging 4(1), 131–149 (2010). [CrossRef]

6. T. Johansson and B. Sleeman, “Reconstruction of an acoustically sound-soft obstacle from one incident field and the far-field pattern,” IMA J. Appl. Math. 72(1), 96–112 (2007). [CrossRef]

7. O. I. Yaman and F. L. Louër, “Material derivatives of boundary integral operators in electromagnetism and application to inverse scattering problems,” Inverse Prob. 32(9), 095003 (2016). [CrossRef]

8. J. Sugisaka, T. Yasui, and K. Hirayama, “Reconstruction of scatterer shape from relative intensity of scatteted field by using linearized boundary element method,” IEICE Trans. Electron. E103.C(2), 30–38 (2020). [CrossRef]

9. D. Zhang and Y. Guo, “Uniqueness results on phaseless inverse acoustic scattering with a reference ball,” Inverse Problems 34(8), 085002 (2018). [CrossRef]

10. X. Xu, B. Zhang, and H. Zhang, “Uniqueness in inverse scattering problems with phaseless far-field data at a fixed frequency,” SIAM J. Appl. Math. 78(3), 1737–1753 (2018). [CrossRef]

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

12. J. Garnaes, P.-E. Hansen, N. Agersnap, J. Holm, F. Borsetto, and A. Kühle, “Profiles of a high-aspect-ratio grating determined by spectroscopic scatterometry and atomic-force microscopy,” Appl. Opt. 45(14), 3201–3212 (2006). [CrossRef]

13. P.-E. Hansen and L. Nielsen, “Combined optimization and hybrid scalar–vector diffraction method for grating topography parameters determination,” Mater. Sci. Eng. B 165(3), 165–168 (2009). [CrossRef]

14. M. Karamehmedović, M.-P. Sørensen, P.-E. Hansen, and A. Lavrinenko, “Application of the method of auxiliary sources to a defect-detection inverse problem of optical diffraction microscopy,” J. European Opt. Soc.-Rapid publications 5, 10021 (2010). [CrossRef]

15. M. Karamehmedović, P.-E. Hansen, K. Dirscherl, E. Karamehmedović, and T. Wriedt, “Profile estimation for Pt submicron wire on rough Si substrate from experimental data,” Opt. Express 20(19), 21678–21686 (2012). [CrossRef]

16. S. Afifi and M. Diaf, “Scattering by random rough surfaces: Study of direct and inverse problem,” Opt. Commun. 265(1), 11–17 (2006). [CrossRef]

17. S. Ishikawa and Y. Hayasaki, “Super-resolution complex amplitude reconstruction of nanostructured binary data using an interference microscope with pattern matching,” Opt. Express 21(15), 18424–18433 (2013). [CrossRef]

18. K. Belkebir and A. Sentenac, “High-resolution optical diffraction microscopy,” J. Opt. Soc. Am. A 20(7), 1223–1229 (2003). [CrossRef]

19. G. Maire, F. Drsek, J. Girard, H. Giovannini, A. Talneau, D. Konan, K. Belkebir, P. Chaumet, and A. Sentenac, “Experimental demonstration of quantitative imaging beyond Abbe’s limit with optical diffraction tomography,” Phys. Rev. Lett. 102(21), 213905 (2009). [CrossRef]

20. T. Zhang, Y. Ruan, G. Maire, D. Sentenac, A. Talneau, K. Belkebir, P. Chaumet, and A. Sentenac, “Full-polarized tomographic diffraction microscopy achieves a resolution about one-fourth of the wavelength,” Phys. Rev. Lett. 111(24), 243904 (2013). [CrossRef]

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

22. B. Zhang and H. Zhang, “Imaging of locally rough surfaces from intensity-only far-field or near-field data,” Inverse Prob. 33(5), 055001 (2017). [CrossRef]

23. J. Sugisaka, T. Yasui, and K. Hirayama, “Expansion of the difference-field boundary element method for numerical analyses of various local defects in periodic surface-relief structures,” J. Opt. Soc. Am. A 32(5), 751–763 (2015). [CrossRef]

24. A. Tadeu and J. António, “Use of constant, linear and quadratic boundary elements in 3D wave diffraction analysis,” Eng. Anal. Bound. Elem. 24(2), 131–144 (2000). [CrossRef]

25. N. Georgieva, S. Glavic, M. Bakr, and J. Bandler, “Feasible adjoint sensitivity technique for EM design optimization,” IEEE Trans. Microwave Theory Tech. 50(12), 2751–2758 (2002). [CrossRef]

26. S. Ali, N. Nikolova, and M. Bakr, “Recent advances in sensitivity analysis with frequency-domain full-wave EM solvers,” Appl. Comput. Electromag. Soc. J. 19(3), 147–154 (2004).

$m$	$σ_{u, m} (n_{1} = 1.5)$	$σ_{u, m} (n_{1} = 2.5)$	$σ_{u, m} (n_{1} = 3.5)$
1	$0.412 \times 10^{- 2}$	$0.368 \times 10^{- 2}$	$0.441 \times 10^{- 2}$
2	$0.551 \times 10^{- 2}$	$0.893 \times 10^{- 2}$	$0.576 \times 10^{- 2}$
3	$2.12 \times 10^{- 2}$	$1.51 \times 10^{- 2}$	$0.647 \times 10^{- 2}$
4	$2.31 \times 10^{- 2}$	$1.78 \times 10^{- 2}$	$1.23 \times 10^{- 2}$

Profile reconstruction of a local defect in a groove structure and the theoretical limit under the vector diffraction theory

Abstract

1. Introduction

2. Principle of profile reconstruction

3. Reconstruction simulation

4. Discussion on the limit of reconstruction

5. Conclusion

Funding

Disclosures

References

Cited By

Figures (13)

Tables (1)

Equations (61)

Optics Express