Open Access
Add to BibSonomy
Abstract
A reconstruction algorithm for the defect profile beneath a dielectric rough surface using scattered far fields is proposed. This is a fundamental technique for optical measurements, such as diffraction tomography, for defect inspection of microfabricated devices. In general, the profile reconstruction process is considerably time consuming. We propose a domain-boundary integral hybrid method to reduce the number of operations while maintaining the rigor of scattering; the polarization properties, scattering, and multiscattering in the sample are considered. We present reconstruction simulations to validate the proposed algorithm. The reconstruction limit as well as its dependency on sample illumination and rough surface shape are also discussed.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. INTRODUCTION
The determination of the dielectric distribution of an object using a scattered electromagnetic field is called the inverse scattering problem. This is a fundamental problem in various non-destructive measurement techniques using electromagnetic waves, such as radio waves, visible light, and x rays.
Diffraction tomography is a measurement technique based on the inverse scattering problem. This technique reconstructs a cross-sectional image of the target sample from the scattered fields. In general, the relationship between the sample profile and the scattered field is considerably complicated, and many algorithms for the efficient and accurate reconstruction of target samples have been proposed.
When the refractive-index contrast in the sample is small, the first Born approximation can be applied [1–3]. In such a low-contrast sample, multiple scattering in the sample can be ignored, and the relationship between the scattered field and the refractive index distribution is expressed by a simple integral equation. It consists of the domain integral of a free-space Green’s function, contrast function (sample profile), and incident field. It is also a linear equation with respect to the contrast function. The Born iterative method (BIM) is an improved method for achieving a high noise tolerance [4].
Alternative iteration algorithms have also been proposed for high-contrast samples. In the distorted BIM (DBIM) [5], a Green’s function containing multiscattering in the sample is employed. In the iteration process, both the Green’s function and contrast function are simultaneously updated. The DBIM has been applied to finite-sized scatterers [5–7]. Ultrasonic waves have also been used to reconstruct three-dimensional scatterers [8]. A reconstruction algorithm using a scattered field of multilevel frequencies has been proposed [9,10]. The multilevel fast multipole algorithm, which is known as the rapid solution algorithm for scattering simulations, was applied to the DBIM [11]. The BIM-interpolated DBIM is a method that takes advantage of high noise tolerance (BIM) and fast convergence (DBIM) [12]. Remis and Van den Berg introduced two types of regularization processes to the DBIM to obtain a high edge-preserving effect [13], which was later modified to eliminate the redundancy of the two regularization processes and the necessity of artificial determination of the regularization parameters [14].
Nonlinear optimization algorithms have also been applied to solve the inverse scattering problem [15–18]. In these methods, a reconstruction process is performed to minimize the cost function. The Newton–Kantorovitch method, which is a nonlinear optimization algorithm, is known to be equivalent to DBIM [13]. In the contrast-source method [19–21], the product of two unknown parameters, the dielectric contrast and field in the sample, is considered as an unknown parameter.
The target problem in our study was to reconstruct small defects buried beneath the surfaces of microfabricated devices, for example, a particle buried under a layer such as metal, oxide, or other epitaxial layers on an electric device, and metal or dielectric coating on a diffractive optical device. An irregular small area whose physical or chemical properties are changed by processes such as ion implantation and oxidation in a device, such as solid-state memory, processors, lasers, and diffused channel waveguides, is also a reconstruction target. As in general optical measurements, the incident wave illuminates the sample in a wide region, and scattered fields are observed at far distances through any optical system. We assume that the dielectric constants inside and outside the sample, including the rough surface shape, are known, and that only the inhomogeneous dielectric constant distribution in the defect is unknown.
Similar problems have been solved for some applications, such as the exploration of objects buried in the Earth. In some studies, the DBIM has been applied to reconstruct a buried object under a flat surface [22–24] and under a rough surface [25–27]. Reconstruction using the contrast-source method has also been employed [28,29]. However, a drawback of these methods is that a large amount of computational resources is required; in the integral equation for obtaining the sample profile, a Green’s function that contains the scattering effect on the rough surface is required. The Green’s function is given by the domain integral of the field near the rough surface, and this field is obtained by solving a large number of simultaneous equations. Moreover, these simultaneous equations consist of another Green’s function, which includes reflection and refraction on a flat dielectric interface. An infinite integral in the wavenumber domain is necessary to compute the Green’s function [30].
In this study, we propose a domain-boundary integral hybrid method for reconstructing a target defect. In this method, clear dielectric interfaces bounded by homogeneous dielectric regions on a rough surface are discretized with boundary (one-dimensional) elements, and inhomogeneous regions are discretized with two-dimensional elements. The region where two-dimensional discretization is necessary is only a small region near the defect. Thus, computation resources are reduced compared with the numerical methods, such as the finite-difference time-domain and finite element methods. The hybrid method also overcomes the drawback of the standard boundary-integral method {boundary element method (BEM) [31] or method of moment [32,33]}, which cannot treat an inhomogeneous scatterer. Green’s functions that include scattering off the rough surface and reflection/refraction on a flat surface are not necessary. Because the discretization of the field near the rough surface is reduced from two-dimensional to one-dimensional, memory resources in the computation process can be saved. It is easy to analyze the limit of the defect profile reconstruction because the equations of the proposed method simply represent the relationship between the dielectric distribution and the scattered field.
The remainder of this paper is organized as follows. Section 2 describes the domain-boundary integral hybrid method and its integral equations for the forward scattering problem. In Section 3, the domain-boundary integral hybrid method is applied to the inverse scattering problem. Section 4 presents the numerical results of the reconstruction simulations. In Section 5, we discuss the reconstruction limit determined by factors such as illumination and rough surface shape of the sample. Finally, we present conclusions in Section 6.
2. DOMAIN-BOUNDARY INTEGRAL HYBRID METHOD
In this section, the numerical process for the forward scattering computation using the domain-boundary integral hybrid method is described. First, we consider a finite-sized scatterer that contains a small inhomogeneous region, and then we describe a semi-infinite sample with a buried defect. Throughout this paper, we assume that the sample structure and incident wave are constant along the $z$ axis. Thus, we consider only two-dimensional field and dielectric constant distributions on the $x-y$ plane. The position on this plane is denoted by the vector ${\boldsymbol r}$. The field distribution can be decomposed into $s$ polarization [transverse magnetic (TM) mode] and $p$ polarization [transverse electric (TE) mode]. We separate these components so that we can describe the field distribution with only the $z$ component of the electric field, ${E_z}$, for $s$ polarization, or the magnetic field, ${H_z}$, for $p$ polarization. The incident wave has a single angular frequency $\omega$, and the time dependency of the field is expressed by $\exp (j\omega t)$, where $j$ is the imaginary unit, $\sqrt {- 1}$.
A. Finite-Sized Scatterer
Fields ${E_z}$ and ${H_z}$ obey the following wave equations:
Integral expressions and equations for $p$ polarization are derived using a similar process. Equation (2) is first arranged as
When numerically computing the integral expressions and integral equations, boundary ${C}$ and region ${{S}_3}$ are discretized with boundary elements and triangle elements, respectively. For each boundary element, two sampling points ${A}$ and ${B}$ are placed, as shown in Fig. 2(a). When an $s$ coordinate is placed such that $s = \pm 1$ at both element terminals, $s = - 0.57735$ at ${A}$ and $s = 0.57735$ at ${B}$ [34]. The field distribution on the element is linearly interpolated and is represented by a linear combination of ${E_z}$ at ${A}$ and ${B}$ and two basis functions, $u({\boldsymbol r})$. Consequently, the field on ${C}$ can be written as
Fig. 2. (a) Boundary and (b) triangular elements and basis functions for discretizing the field and scattering-potential distributions.
Region ${{S}_3}$ is also discretized using triangular elements, as in the finite element method. The field distribution in ${{S}_3}$ is approximated by a linear combination of the basis functions. The basis function ${v_i}({\boldsymbol r})$ is a linear function that is unity at the $i$th node and zero at the other nodes, as depicted in Fig. 2(b). The field distribution in ${{S}_3}$ is represented as
where ${N_S}$ is the number of nodes. The discretized operators are expressed as a product of vectors $\bar {\cal J}_{n,S}^{({S})}({\boldsymbol r})$ and ${{\boldsymbol E}_z}$:Fig. 3. (a) Profile of a sample with a rough surface and an inhomogeneous defect beneath the rough surface. (b) Profile of the sample without the rough surface and defect.
B. Semi-infinite Sample
The cross section of the target sample is illustrated in Fig. 3(a). Regions ${{S}_1}$, ${{S}_2}$, and ${{S}_3}$ are denoted according to the scatterer illustrated in Fig. 1. Region ${{S}_4}$ is a projection on the surface. We also define boundaries ${{C}_0}$ (between ${{S}_1}$ and ${{S}_2}$), ${{C}_1}$ (between ${{S}_1}$ and ${{S}_4}$), and ${{C}_2}$ (between ${{S}_2}$ and ${{S}_4}$). The integral expressions for ${\boldsymbol r} \in {S}$ are as follows:
3. INVERSE SCATTERING ALGORITHM
It is assumed that the scattered far fields of the sample, ${f^{{(s)}}}({\boldsymbol \rho})$ $(= {f_0}({\boldsymbol \rho}) + f_{S}^{{(s)}}({\boldsymbol \rho}))$, at ${{\boldsymbol r}_n}$ $({{\boldsymbol r}_n} \in {{S}_2},n = 1,2, \ldots ,{N_{o}})$ are experimentally measured. The unit vector ${\boldsymbol \rho}$ represents the propagation direction of the far field. If we have a scattering potential distribution ${F^{(1)}}$ close to that of the sample ${F^{{(s)}}}$, the following linear equation holds for all $n$:
The coefficient Eqs. (44)–(50) are a set of simultaneous equations. For different $m$, the coefficients of the unknowns $\partial {f^{(s)}}/\partial\! {F_m}$ and ${\partial ^2}{f^{(s)}}/\partial\! {F_m}\partial n$ are invariant. Thus, we can obtain the solution for all $m$ efficiently using $LU$ decomposition of the coefficient matrix.
After obtaining $\partial {f^{(1)}}/\partial\! {F_m}$, Eq. (42) for all incident waves is solved simultaneously with respect to $\Delta F$. The defect profile ${F^{(s)}}$ is given by ${F^{(1)}} + \Delta F$. When the initial guess ${F^{(1)}}$ is not close to ${F^{(s)}}$, the result ${F^{(1)}} + \Delta F$ is often inaccurate. An accurate solution can be obtained by letting ${F^{(1)}} + \Delta F$ equal ${F^{(1)}}$ and iterating the above process until ${F^{(1)}}$ converges.
When the numbers of boundary elements and triangle elements are $N$ and $M$, respectively, the necessary number of operations for the reconstruction process is $O((N + M{)^3})$. It is determined by solving Eqs. (44)–(50), requiring the maximum number of operations in the reconstruction process. If we applied an efficient algorithm such as the fast multipole method [33] and hierarchical matrices [39], the number of operations is reduced to $O((N + M)\log (N + M))$.
4. RECONSTRUCTION SIMULATION
A. Simulation Parameters
The geometry of the target sample is shown in Fig. 4. For simulations, the coordinate system was chosen such that the surface of the sample was at $y = 0$. The center, height, and width of the surface projection were $x = 0$, $0.1\lambda$, and $1.0\lambda$, respectively, where $\lambda$ is the illumination wavelength in a vacuum. The defect region ${{S}_3}$ was a square of size $0.6\lambda$. The center position $({x_c},{y_c}) = (0, - 0.65\lambda)$. The incident wave was a Gaussian beam with a beam waist of $40\lambda$, and the center was located at the origin. The propagating directions of the incident waves in ${{S}_2}$ and ${\theta _{ i}}$ were ${-}{15^ \circ}$, ${-}{90^ \circ}$, and ${-}{165^ \circ}$. For each incident wave, the scattered wave was observed at a far distance in the direction of a unit vector ${\boldsymbol \rho}$, which was oriented from 30° to 150° sampled at equal angular intervals. The number of samples, ${N_{o}}$, was 144. The relative permittivities in ${{S}_1}$, ${{S}_2}$, ${{S}_4}$ were 1.5, 1.0, and 1.5, respectively.
The boundaries ${{C}_0}$, ${{C}_1}$, and ${{C}_3}$ were discretized with boundary elements whose lengths were less than $\lambda /12$. Region ${{S}_3}$ was discretized with 144 isosceles right triangles as illustrated in Fig. 4 (The number of nodes was 169).
Two samples with defects were prepared. Sample A was a circular defect whose center position $(x,y) = ({x_c},{y_c})$, and the relative permittivity in ${{S}_3}$ was
B. Validation of Forward Scattering Simulation
Before reconstruction simulation, we validate the domain-boundary integral hybrid method for forward scattering simulation. When $F({\boldsymbol r})$ is set to be constant (dielectric permittivity is homogeneous in ${{S}_3}$), the scattered wave can be computed by the standard BEM [32]. If the sample is illuminated only around ${{S}_4}$, the infinite integral path on the sample surface can be truncated near ${{S}_4}$.
We prepared the sample data as illustrated in Fig. 4. The scattering potential was set to a constant ${(1.7^2} - {1.5^2})$. The incident beam is a Gaussian beam, whose width is $4\lambda$ or $10\lambda$, and it perpendicularly illuminates the sample surface.
The scattered far fields in ${{S}_2}$ are plotted in Fig. 5. The incident beam width is $4\lambda$, and the length of the truncated integral paths (${C_0}$) for the left and right sides is $10\lambda$. The field distribution agreed well, and the relative error, defined by
Next, we varied the length of the truncated path ${C_0}$, and examined the convergence of the far field to evaluate the necessary computation resources. We assumed that the far field computed with the ${C_0}$ length of $60\lambda$ is a true value and evaluated the relative error. The results for the incident beam width of $4\lambda$ and $10\lambda$ are plotted in Figs. 6(a) and 6(b), respectively.
Fig. 6. Convergence property with respect to the ${C_0}$ (sample surface) length for the incident beam width of (a) $4\lambda$ and (b) $10\lambda$.
The convergence speed for the standard BEM is decreased with expanding beam width. As the field on ${C_0}$ is widely distributed as the incident beam width is expanded, sufficiently longer ${C_0}$ than the incident beam width is necessary. By contrast, the hybrid method computes ${f_0}$ and ${f_s}$ separately. The ${f_s}$ component localizes near the ${{S}_4}$, independent of the incident beam width. The ${f_0}$ component is a simple reflected beam on the flat surface, and an integral equation is not necessary. Consequently, it rapidly converges with increasing ${C_0}$ length. This convergence property holds for a complicated random rough surface if we use the algorithm presented in Ref. [38].
C. Inverse Scattering Simulation
During the reconstruction process, we defined two types of errors to evaluate the accuracy and convergence of the reconstruction:
Fig. 8. Reconstructed scattering-potential distribution (${F^{({S})}}$) after one, five, and 29 iteration steps: (a) $s$ polarization, sample A; (b) $s$ polarization, sample B; (c) $p$ polarization, sample A; (d) $p$ polarization, sample B.
Next, the scattering potential was reconstructed from noisy data. The noisy data were prepared by adding $w({\boldsymbol r}){e^{j\theta ({\boldsymbol r})}}$ to ${f^{(s)}}({\boldsymbol r})$, where $w({\boldsymbol r})$ is a normal random number, and $\theta ({\boldsymbol r})$ is a uniform random number, ${-}\pi \le \theta ({\boldsymbol r}) \lt \pi$. The amount of noise is determined by $\sigma ({\boldsymbol r}) \equiv {n_f}|{f^{(s)}}({\boldsymbol r})|$, which is the standard deviation of $w({\boldsymbol r})$. The regularization parameter $\chi$ was determined using Morozov’s discrepancy principle; $\chi$ for each ${n_f}$ is listed in Table 1. The structural errors after 29 iterations are plotted as a function of ${n_f}$ in Fig. 9. The scattering potential distributions (${F^{(1)}}$) for ${n_f}= 10^{- 4}$, ${10^{- 3}}$, and ${10^{- 2}}$ are shown in Fig. 10.
Table 1. Regularization Parameters for Reconstruction from Noisy Data
Fig. 9. Relation between the amount of noise (${n_f}$) and the structural error, ${e_s}$, after 29 iteration steps.
Fig. 10. Scattering-potential distribution reconstructed from noisy reference data obtained after 29 iteration steps: (a) $s$ polarization, sample A; (b) $s$ polarization, sample B; (c) $p$ polarization, sample A; (d) $p$ polarization, sample B.
We also reconstructed the defect profile for different heights of the surface projections, $0.1\lambda$, $0.2\lambda$, $0.4\lambda$, and $0.6\lambda$. As shown in Fig. 11, the reconstructed defect profiles after 29 iteration steps are comparable; there is little dependence of the reconstruction result on the rough-surface height. However, it is necessary to set an accurate height when reconstructing the defect. For example, when the projection height value in the reconstruction process is varied by only $0.01\lambda$ from the sample height, the structural error is increased from 0.04608 to 33.55. This is because the refractive-index contrast on the sample surface is much larger than that of the buried defect, and a slight variation of the surface height significantly changes the scattered field compared with that from the buried defect. The relation between the surface shape and the reconstruction property is discussed in Section 5.
Fig. 11. Reconstructed profiles of buried defects for different heights of surface projections, $0.1\lambda$, $0.2\lambda$, $0.4\lambda$, and $0.6\lambda$. The incident light is $s$ polarized, and the noise level ${n_f}$ is ${10^{- 3}}$.
Sample B was reconstructed from different initial guesses. The incident beam was $s$ polarized, and the reference data were a noisy with ${n_f}= 10^{- 3}$. As shown in Fig. 12, each initial guess provided comparable reconstruction results.
Fig. 12. Scattering-potential distributions reconstructed from different initial guesses. The defect profile ${F^{(s)}}$ is sample B, and the incident light is $s$ polarized.
5. LIMIT OF RECONSTRUCTION
Equation (42) can be rewritten in matrix form as follows:
The singular values ${\sigma _i}$ and the corresponding ${{\boldsymbol v}_i}$ for $F({\boldsymbol r}) = 0$ are shown in Figs. 13 and 14, respectively, in decreasing order of the singular values ${\sigma _i}$. For both $s$ and $p$ polarizations, the contrast of the singular value ${\sigma _j}/{\sigma _i}$ and the distribution ${{\boldsymbol v}_i}$ are comparable. This means that the reconstruction accuracies for $s$ and $p$ polarizations are similar.
Fig. 13. Singular value of the matrix $A$ in Eq. (60) for $s$ and $p$ polarizations. The singular value ${\sigma _i}$ is sorted in decreasing order, and is zero for $i \geqslant 85$.
Fig. 14. Real part of ${{\boldsymbol v}_i}$ $(1 \le i \le 5)$ and the corresponding singular values for (a) $s$ polarization and (b) $p$ polarization.
The singular value ${\sigma _i}$ $(i \gt 40)$ is less than ${10^{- 16}}$ times as small as ${\sigma _1}$. Thus, the corresponding ${{\boldsymbol v}_i}$ $(i \gt 40)$, which is a granular distribution, is unstably amplified and causes a large error and multiplexity of the result. Tikhonov’s regularization eliminates such instabilities and multiplexity and provides a solution such that the result contains little components ${{\boldsymbol v}_i}$ for the small singular values. This is the reason that we obtained comparable results for different initial guesses (see Fig. 12).
The scattering potential for the largest singular value, ${{\boldsymbol v}_i}$, is a layered distribution, as shown in Fig. 14. The layer thickness was approximately ${\lambda _1}/4$, which corresponds to the dimensions of a dielectric multilayer mirror. This structure reflects the incident wave (${\theta _{i}} = - {90^ \circ}$) toward 90° such that $|{\boldsymbol f}|$ $({=} {\sigma _1}|{{\boldsymbol u}_1}|)$ becomes large.
To discuss the factor that increases ${\sigma _j}/{\sigma _i}$ in detail, we analyze ${\boldsymbol v}$ in the wavenumber space. We consider vector ${{\boldsymbol v}_1}$ to be the product of
and a square window function $S({\boldsymbol r})$, which is unity in ${\boldsymbol r} \in {{S}_3}$ and zero otherwise. In wavenumber space, ${\cal F}{{\boldsymbol v}_1}$ is expressed as the convolution of ${\cal F}\exp (- j{{\boldsymbol K}_1}{\boldsymbol r}) = \delta ({\boldsymbol k} - {{\boldsymbol K}_1})$ and $[{\cal F}S({\boldsymbol r})]({\boldsymbol k})$, namely, $[{\cal F}S({\boldsymbol r})]({\boldsymbol k} - {{\boldsymbol K}_1})$, where ${\cal F}$ represents the two-dimensional Fourier transform.As is often discussed in Bragg’s diffraction and x-ray diffraction in crystals, the wavenumber vector of the diffracted wave ${{\boldsymbol k}_s}$ is given by ${{\boldsymbol k}_s} = {{\boldsymbol k}_i} + {\boldsymbol K}$ as depicted in Fig. 15(a). The vector ${{\boldsymbol k}_i}$ is the wavenumber of the incident wave in region ${{S}_3}$, and the circle with the dashed line is Ewald’s circle. The amplitude of the scattered wave is determined by that of the incident wave and ${C_{{{\boldsymbol K}_1}}}$ in Eq. (68). Considering refraction on the sample surface, the wavenumber vector ${{\boldsymbol k}_s}$ must be oriented between 54.74° and 125.3° in ${{S}_3}$ such that the diffracted wave arrives in the observation range of 30°–150° in ${{S}_2}$. Arc ${{D}_1}$ in Fig. 15(b) is a set of ${\boldsymbol K}$ such that ${{\boldsymbol k}_s}$ satisfies this condition for ${\theta _{ i}} = - {90^ \circ}$. For ${\theta _{ i}} = - {15^ \circ}$ (${-}{49.91^ \circ}$ in ${{S}_3}$), as shown in Fig. 15(c), the set of ${\boldsymbol K}$ is shifted, and is shown as ${{ D}_2}$ in Fig. 15(d).
Fig. 15. Relation between the wavenumber vector of the incident wave (${k_{ i}}$) and the scattered wave (${k_{S}}$) for various directions of ${{\boldsymbol k}_i}$ in ${{S}_1}$ (a) ${-}{90^ \circ}$ and (b) ${-}{49.91^ \circ}$ (${-}{15^ \circ}$ in ${{S}_2}$). ${{D}_1}$ in (c) and ${{D}_2}$ in (d) are sets of ${\boldsymbol K}$ such that the scattered wave arrives within the observation range.
Figure 16(a) shows the distribution of ${\cal F}{{\boldsymbol v}_1}$ for $s$ polarization. Because the maximum value of ${\cal F}{{\boldsymbol v}_1}$ is located on ${{D}_1}$, ${C_{{{\boldsymbol K}_1}}}$ in Eq. (68) becomes large, and ${{\boldsymbol v}_1}$ produces a strong scattered wave. This results in a large ${\sigma _1}$. As plotted in Fig. 16(b), the maximum of ${\cal F}{{\boldsymbol v}_2}$ is located on ${{\rm D}_2}$ and ${{D}_3}$, where ${{D}_3}$ is the set of ${\boldsymbol K}$ such that the scattered wave arrives within the observation range for the incident direction of ${-}{165^ \circ}$. However, the peak in ${\cal F}{{\boldsymbol v}_2}$ is split into two peaks: ${{\boldsymbol K}_{11}}$ and ${{\boldsymbol K}_{12}}$. Under the constraint $|{\cal F}{{\boldsymbol v}_1}| = |{\cal F}{{\boldsymbol v}_2}| = 1$, both coefficients ${C_{{{\boldsymbol K}_{11}}}}$ and ${C_{{{\boldsymbol K}_{12}}}}$ for the two peaks are smaller than ${C_{{{\boldsymbol K}_1}}}$. Thus, ${\sigma _2}$ is smaller than ${\sigma _1}$. For ${{\boldsymbol v}_3}$, the maximum of ${\cal F}{{\boldsymbol v}_3}$ deviates from ${{D}_1}$, ${{D}_2}$, and ${{D}_3}$ such that ${\cal F}{{\boldsymbol v}_3}$ is determined to be orthogonal to both ${\cal F}{{\boldsymbol v}_1}$ and ${\cal F}{{\boldsymbol v}_2}$. As a result, ${\sigma _3}$ decreases further than ${\sigma _2}$.
Fig. 16. Distributions of (a) ${\cal F}{{\boldsymbol v}_1}$, (b) ${\cal F}{{\boldsymbol v}_2}$, and (c) ${\cal F}{{\boldsymbol v}_3}$. The dashed circle is the Ewald’s circle, and ${{D}_1}$, ${{D}_2}$, and ${{D}_3}$ are the sets of ${\boldsymbol K}$ that produce scattered waves that propagate toward the observation range for incident directions of ${-}{90^ \circ}$, ${-}{15^ \circ}$, and ${-}{165^ \circ}$, respectively.
According to the above discussion, the secondary and subsequent singular values $({\sigma _2},{\sigma _3}, \ldots)$ increase upon performing multiple measurements using incident waves corresponding to various directions. For example, if the reference data are acquired using only ${\theta _{i}} = - {90^ \circ}$, ${{ D}_2}$ and ${{ D}_3}$ shown in Fig. 16 disappear, and scattered waves in the observation range for the scattering potentials ${{\boldsymbol v}_2}$ and ${{\boldsymbol v}_3}$ are considerably decreased. In fact, when the reference data contained scattered waves with only ${\theta _{i}} = - {90^ \circ}$, the contrast of the singular values ${\sigma _1}/{\sigma _2}$ increased from 3.363 to 7.766.
To decrease ${\sigma _1}/{\sigma _i}$ $(i \gt 1)$, an effective method is to increase the difference in ${\theta _{i}}$. Then, the distance between ${{D}_1}$ and ${{D}_2}$ or ${{D}_1}$ and ${{D}_3}$ can be increased. If the difference in incident angles becomes small, ${{D}_2}$ and ${{D}_3}$ approach ${{D}_1}$, and we cannot find ${{\boldsymbol v}_2}$ whose maximum is on ${{D}_1}$, ${{D}_2}$, or ${{D}_3}$. Thus, ${\boldsymbol v}_1^{H}{{\boldsymbol v}_2} = 0$. This feature corresponds to resolution enhancement by microscopy with a high numerical aperture. However, in the target sample, the incident wave was refracted at the surface of the sample, and the difference in the direction of ${{\boldsymbol k}_i}$ in ${{S}_3}$ was reduced. If region ${S_2}$ is filled by a medium with a higher refractive index material, the angular difference in ${{\boldsymbol k}_i}$ is increased. For example, if ${n_2}$ is changed from 1.0 to 1.3, ${\sigma _1}/{\sigma _2}$ decreases from 3.363 to 1.990. This feature corresponds to resolution enhancement by immersion microscopy.
Using incident light with a higher frequency is also effective to decrease ${\sigma _1}/{\sigma _i}$ $(i \gt 1)$. The singular values when the frequency of incident light is doubled are plotted in Fig. 17. Comparing ${{\boldsymbol v}_i}$ whose corresponding $|{\sigma _i}|$ is more than 0.1, ${{\boldsymbol v}_i}$ for the doubled frequency contains about twice as high spatial frequency as that for the original frequency. In the wavenumber domain, the radius of the Ewald’s circle and $|{\boldsymbol K}|$ are doubled, and we can obtain a defect profile with a higher spatial frequency. This analysis indicates that the resolving power is proportionally increased with increasing the frequency of incident light.
Fig. 17. Singular values when the frequency of incident light is doubled. The distributions, ${{\boldsymbol v}_7}$ (doubled frequency) and ${{\boldsymbol v}_3}$ (the original frequency), are scattering potentials for $|{\sigma _i}| \gt 0.1$.
Finally, the roughness of the sample was considered. The incident wave in ${{S}_3}$ contained a scattered wave (${f_{S}}$) from the rough surface. As shown in Fig. 18(a), the scattered wave propagates from the projection ${{S}_4}$ as a cylindrical wave. The propagating direction in ${{S}_3}$ was close to that of the direct incident wave after refraction at the sample surface. The rough structure of the target sample had little effect on the singular value.
Fig. 18. Phase distributions of the scattered field component (${f_{S}}$) for (a) single-bump surface structure (Fig. 4) and (b) two-bump structure. The square with a dashed line indicates region ${{S}_3}$.
If the rough surface is changed to a two-bump structure, as illustrated in Fig. 18(b), a large variation is observed in the propagating direction of the scattered wave in ${{S}_3}$. The scattered wave from each bump was superimposed in ${{S}_3}$ with the same phase when the incident wave was ${-}{90^ \circ}$. In contrast, when the incident wave was ${-}{15^ \circ}$, the phase difference from each bump was close to $\pi /2$, and a significantly different field was produced in ${{S}_3}$. In this case, the contrast ${\sigma _1}/{\sigma _2}$ decreases from 3.363 to 2.505, which is 25.5% smaller than that for the single-bump structure. Analogous to coded-aperture imaging [40], if we employ any optimization process considering the surface structure, we can find incident wave distributions that maximize the difference in the propagating direction in ${{S}_3}$.
In the imaging optical system, the image resolution is proportional.
6. CONCLUSION
We presented an efficient algorithm for reconstructing the profile of a defect buried beneath a rough dielectric surface. The domain-boundary integral hybrid method does not require the time-consuming computation of Green’s function, which includes the scattering effect on the rough surface. Consequently, this method rapidly provides the scattered field during iterative reconstruction. The hybrid method has been validated by comparing with the standard BEM. We have also verified rapid convergence of the solution independent of the incident beamwidth. Our algorithm also provides a simple expression that relates the target defect profile to the scattered field. It enables us to easily analyze the limit of reconstruction using this expression and the singular-value decomposition. These analyses yielded the following conclusions: first, there was little polarization dependency in the reconstruction limit; second, multiple measurements obtained by varying the illumination direction are effective in enhancing the reconstruction limit; finally, we suggest that the incident wave distribution should be optimized to enhance the reconstruction limit.
For practical applications, another measurement of the sample surface with a scanning electron microscope or an atomic force microscope may be necessary, because the proposed reconstruction process must acquire the accurate surface shape on the sample. In a previous work, we proposed another inverse scattering algorithm to reconstruct the surface shape on a surface-relief structure [41]. By combining this algorithm and the proposed algorithm in this paper, we can simultaneously reconstruct the shapes of the rough surface, layered structure inside the sample, and the buried defect from the scattered field data.
Funding
Japan Society for the Promotion of Science (19K05642).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
REFERENCES
1. K. Iwata and R. Nagata, “Calculation of refractive index distribution from interferograms using the Born and Rytov’s approximation,” Jpn. J. Appl. Phys. 14, 379 (1975). [CrossRef]
2. A. Devaney, “Inversion formula for inverse scattering within the Born approximation,” Opt. Lett. 7, 111–112 (1982). [CrossRef]
3. B. Chen and J. J. Stamnes, “Validity of diffraction tomography based on the first Born and the first Rytov approximations,” Appl. Opt. 37, 2996–3006 (1998). [CrossRef]
4. M. Moghaddam and W. C. Chew, “Study of some practical issues in inversion with the Born iterative method using time-domain data,” IEEE Trans. Antennas Propag. 41, 177–184 (1993). [CrossRef]
5. W. C. Chew and Y.-M. Wang, “Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method,” IEEE Trans. Med. Imaging 9, 218–225 (1990). [CrossRef]
6. G. Leone, R. Persico, and R. Pierri, “Inverse scattering under the distorted Born approximation for cylindrical geometries,” J. Opt. Soc. Am. A 16, 1779–1787 (1999). [CrossRef]
7. R. Lavarello and M. Oelze, “A study on the reconstruction of moderate contrast targets using the distorted Born iterative method,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 55, 112–124 (2008). [CrossRef]
8. R. J. Lavarello and M. L. Oelze, “Tomographic reconstruction of three-dimensional volumes using the distorted Born iterative method,” IEEE Trans. Med. Imaging 28, 1643–1653 (2009). [CrossRef]
9. O. S. Haddadin and E. S. Ebbini, “Imaging strongly scattering media using a multiple frequency distorted Born iterative method,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 45, 1485–1496 (1998). [CrossRef]
10. A. G. Tijhuis, K. Belkebir, and A. Litman, “Multiple-frequency distorted-wave Born approach to 2D inverse profiling,” Inverse Probl. 17, 1635 (2001). [CrossRef]
11. A. J. Hesford and W. C. Chew, “Fast inverse scattering solutions using the distorted Born iterative method and the multilevel fast multipole algorithm,” J. Acoust. Soc. Am. 128, 679–690 (2010). [CrossRef]
12. T. Q. Huy, T. D. Tan, and N. Linh-Trung, “An improved distorted Born iterative method for reduced computational complexity and enhanced image reconstruction in ultrasound tomography,” in International Conference on Advanced Technologies for Communications (ATC) (IEEE, 2014), pp. 703–707.
13. R. F. Remis and P. Van den Berg, “On the equivalence of the Newton-Kantorovich and distorted Born methods,” Inverse Probl. 16, L1 (2000). [CrossRef]
14. H. Zheng, C. Wang, and E. Li, “Modification of enhanced distorted Born iterative method for the 2D inverse problem,” IET Microw. Antennas Propag. 10, 1036–1042 (2016). [CrossRef]
15. H. Harada, D. J. Wall, T. Takenaka, and M. Tanaka, “Conjugate gradient method applied to inverse scattering problem,” IEEE Trans. Antennas Propag. 43, 784–792 (1995). [CrossRef]
16. S. Barkeshli and R. Lautzenheiser, “An iterative method for inverse scattering problems based on an exact gradient search,” Radio Sci. 29, 1119–1130 (1994). [CrossRef]
17. P. Lobel, R. Kleinman, C. Pichot, L. Blanc-Feraud, and M. Barlaud, “Conjugate-gradient method for solving inverse scattering with experimental data,” IEEE Antennas Propag. Mag. 38(3), 48 (1996). [CrossRef]
18. A. Roger, “Newton-Kantorovitch algorithm applied to an electromagnetic inverse problem,” IEEE Trans. Antennas Propag. 29, 232–238 (1981). [CrossRef]
19. P. M. Van Den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Probl. 13, 1607 (1997). [CrossRef]
20. P. M. van den Berg, A. Van Broekhoven, and A. Abubakar, “Extended contrast source inversion,” Inverse Probl. 15, 1325 (1999). [CrossRef]
21. A. Abubakar, P. M. Van den Berg, and J. J. Mallorqui, “Imaging of biomedical data using a multiplicative regularized contrast source inversion method,” IEEE Trans. Microw. Theory Tech. 50, 1761–1771 (2002). [CrossRef]
22. R. Pierri, G. Leone, F. Soldovieri, and R. Persico, “Electromagnetic inversion for subsurface applications under the distorted Born approximation,” Nuovo Cimento C 24, 245–262 (2001).
23. T. J. Cui, W. C. Chew, A. A. Aydiner, and S. Chen, “Inverse scattering of two-dimensional dielectric objects buried in a lossy earth using the distorted Born iterative method,” IEEE Trans. Geosci. Remote Sens. 39, 339–346 (2001). [CrossRef]
24. G. Leone and F. Soldovieri, “Analysis of the distorted Born approximation for subsurface reconstruction: Truncation and uncertainties effects,” IEEE Trans. Geosci. Remote Sens. 41, 66–74 (2003). [CrossRef]
25. Y. Altuncu, F. Akleman, O. Semerci, and C. Ozlem, “Imaging of dielectric objects buried under a rough surface via distorted Born iterative method,” J. Phys. Conf. Ser. 135, 12006 (2008). [CrossRef]
26. F. Li, Q. H. Liu, and L.-P. Song, “Three-dimensional reconstruction of objects buried in layered media using Born and distorted Born iterative methods,” IEEE Geosci. Remote Sens. Lett. 1, 107–111 (2004). [CrossRef]
27. Y. Altuncu, T. Durukan, and R. E. Akdogan, “Reconstruction of two-dimensional objects buried into three-part space with locally rough interfaces via distorted Born iterative method,” Prog. Electromagnet. Res. 166, 23–41 (2019). [CrossRef]
28. E. Tetik and I. Akduman, “3D imaging of dielectric objects buried under a rough surface by using csi,” Int. J. Antennas Propag. 2015, 179304 (2015). [CrossRef]
29. T. U. Gürbüz and B. Aslanyürek, “A two-stage procedure for microwave imaging of a buried dielectric along with the randomly rough surface above it,” Int. J. Antennas Propag. 2015, 928450 (2015). [CrossRef]
30. Y. Altuncu, A. Yapar, and I. Akduman, “Numerical computation of the Green’s function of a layered media with rough interfaces,” Microw. Opt. Technol. Lett. 49, 1204–1209 (2007). [CrossRef]
31. N. Kumagai, N. Morita, and J. R. Mautz, Integral Equation Methods for Electromagnetics (Artech House Antenna Library, 1990).
32. C. Bourlier, N. Pinel, and G. Kubické, Method of Moments for 2D Scattering Problems, Basic Concepts and Applications (Wiley, 2013).
33. W. C. Gibson, The Method of Moments in Electromagnetics (Chapman & Hall/CRC Press, 2014).
34. A. Tadeu and J. Antonio, “Use of constant, linear and quadratic boundary elements in 3D wave diffraction analysis,” Eng. Anal. Bound. Elem. 24, 131–144 (2000). [CrossRef]
35. J.-I. 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, 751–763 (2015). [CrossRef]
36. J.-I. Sugisaka, T. Yasui, and K. Hirayama, “Efficient scattering analysis of arbitrarily shaped local defect in diffraction grating,” IEICE Trans. Electron. E99.C, 76–80 (2016). [CrossRef]
37. J.-I. Sugisaka, T. Yasui, and K. Hirayama, “Efficient analysis of diffraction grating with 10000 random grooves by difference-field boundary element method,” IEICE Trans. Electron. E100.C, 27–36 (2017). [CrossRef]
38. J.-I. Sugisaka, T. Yasui, and K. Hirayama, “Fast actual-size vectorial simulation of concave diffraction gratings with structural randomness,” J. Opt. Soc. Am. A 34, 2157–2164 (2017). [CrossRef]
39. M. Bebendorf, Hierarchical Matrices (Springer, 2008).
40. E. E. Fenimore and T. M. Cannon, “Coded aperture imaging with uniformly redundant arrays,” Appl. Opt. 17, 337–347 (1978). [CrossRef]
41. J. Sugisaka, T. Yasui, and K. Hirayama, “Profile reconstruction of a local defect in a groove structure and the theoretical limit under the vector diffraction theory,” Opt. Express 28, 30908–30927 (2020). [CrossRef]
