Multi-wavelength multi-angle reflection tomography

TING ZHANG,1,2 KÉVIN UNGER,2 GUILLAUME MAIRE,2 PATRICK C. CHAUMET,2,* ANNE TALNEAU,3 CHARAN GODHAVARTI,2 HUGUES GIOVANNINI,2 KAMAL BELKEBIR,2 AND ANNE SENTENAC2 1College of Information Science and Electronic Engineering (ISEE), Zhejiang University, Zhejiang Provincial Key Laboratory of Information Processing, Communication and Networking (IPCAN), Hangzhou 310027, China 2Aix Marseille Université, CNRS, Centrale Marseille, Institut Fresnel, France 3CNRS, Lab. Photon & Nanostruct, 91460 Marcoussis, France *patrick.chaumet@fresnel.fr


Introduction
Computational tomographic diffraction microscopy consists in illuminating the sample under varying angles of incidence and processing the resulting images to reconstruct a super-resolved map of the object. The reconstruction schemes are usually based on a simple model of the wave-target interaction which links the image recorded for a given incident angle to a specific domain of the Fourier spectrum of the sample. This approach, either using field data (via an interferometric mounting) [1] or intensity data [2], is increasingly popular because it provides 3D quantitative images of both absorbing and phase objects [2] with a better resolution than analogical microscopes. Most computational tomographic microscopes work in transmission. In this configuration, the sample spectrum obtained at a given wavelength for a wide set of incident and observation angles encompasses the spectra accessible at larger wavelengths [3][4][5][6]. Thus, the angular diversity (of both the illumination and observation) using monochromatic illumination with the smallest wavelength is sufficient for retrieving all the possible information on the sample. In transmission, multi-wavelength illuminations are essentially useful for improving the signal to noise ratio or compensating a rough discretization of the incident angles and have seldom been implemented [5,7].
In the reflection configuration, on the contrary, angular diversity at the smallest wavelength possible does not permit the retrieval of all the sample information. Long wavelengths allows the scanning of low spatial frequency regions that are inaccessible at small wavelengths. These additional information permit to diminish the axial side lobes of the imager point-spreadfunction [8,9]. Hence, most analogical reflection microscopes, such as Optical Coherence Tomography or Optical Profilometry, take advantage of both angular and wavelength diversity. On the other hand, in computational microscopy, the few reflection configurations have either been implemented with angular diversity without wavelength diversity [10,11] or with wavelength diversity but without angular diversity [12,13].
In this work, we present the first, to our knowledge, reflection tomographic microscope taking advantage of both angular and wavelength diversity. We study a complex configuration, often encountered in the semi-conductor industry, in which a contrasted sample (resin in air), deposited on a silicon substrate, cannot be reconstructed using direct reconstruction schemes (backpropagation, Born or Rytov approximation) because of the substrate reflection and the multiple scattering within the sample. We consider an iterative reconstruction procedure [14,15] based on a rigorous model of the light-sample interaction which proved effective with monochromatic illumination [11,15] and investigate different implementations to account for the multi-wavelength data.

Experimental device
The multi-wavelength Tomographic Diffraction Microscope (TDM) is based on a synthetic aperture interferometric reflection microscope [16], as depicted in Fig. 1. The light source is a super continuum collimated laser (NKT Photonics SuperK Extreme EXW-12), and we consider p = 1, ..., P = 5 wavelengths, λ = 475 nm, 525 nm, 575 nm, 625 nm, and 675 nm, covering thus the visible spectrum with 50 nm step in our experiment. For each wavelength, the output light is linearly polarized and divided into an illumination and a reference path by the beam splitter (BS1).
The light shining the sample is a collimated beam directed towards a rotating mirror M (Newport FSM-300), which is optically conjugated with the object plane in front of the objective OL (Zeiss Epiplan-Apochromat 50× with numerical aperture NA= 0.95 in air). The sample is successively illuminated with l = 1, ..., L illumination angles by changing the orientation of the mirror. The reflected and diffracted field by the sample, is collected by the same objective and detected on a sCMOS camera (Andor Zyla), optically conjugated with the object plane, with a global magnification of about 300.
In the reference path, the beam is spatially filtered with a pinhole P and collimated to obtain an aberration free reference wave. This beam interferes with the sample reflected field thanks to the beam splitter BS 3 . Off-axis digital holography technique is used to retrieve the amplitude and phase of the diffracted field.
Two half-wave plates are placed in the illumination and the detection paths, (HW 1 and HW 2 ) to record the vectorial reflected field for both p and s incident polarizations [17]. The global phase and amplitude of the field recorded at each incidence and wavelength are adjusted so that its specularly reflected part agrees with the theoretical Fresnel reflection coefficient (one assumes that the sample does not distort the specular reflection).

Formulation of the forward scattering problem
The link between the scattered field and the sample is modeled using an electromagnetic forward scattering solver based on the Discrete Dipole Approximation or equivalently the finite volume method [18]). We consider a reference medium consisting in the microscope elements and the substrate holding the target and we introduce the reference monochromatic Green tensor for the p th wavelength, G p such that G p (r, r )s is the field radiated at r in the microscope environment by a source s placed at r [4,19]. The target is viewed as a perturbation of the reference medium. It is defined by a permittivity contrast χ(r) = ε(r) − 1 which is non zero only in a bounded domain Ω. In absence of the target, the field in the microscope for the l−th incident angle and the p−th wavelength, is noted E inc l, p . In Ω it is the sum of the incident collimated beam and its specular reflection by the planar substrate. In presence of the target, the field in the microscope for the l−th incident angle and the p−th wavelength satisfies the volume integral equation [4], From this self-consistent equation, one can obtain the field E l, p in Ω and the scattered field . In symbolic notations, these fields satisfy the equations, where A Ω p and B Γ p , are the Green operators acting from L 2 (Ω) into L 2 (Ω) and from L 2 (Ω) to L 2 (Γ), respectively. Notice that A Ω p and B Γ p do not depend on the angle of incidence of the impinging beam but depend on the wavelength.

Formulation of the inverse scattering algorithm
The reconstruction of the sample from the scattered field recorded by the TDM is performed using a non-linear iterative inversion algorithm which is able to account for multiple scattering. In this technique, known as the Hybrid Method (HM), both χ and E l, p inside Ω are sought simultaneously by minimizing a cost functional involving the distance between the measurements, f mes l, p and the simulated field f diff l, p . The HM has been detailed extensively in [14,15,20] for the monochromatic configuration.
We assume herein that the contrast χ is frequency independent. Two strategies can then be used to solve the inverse scattering problem in presence of wavelength diversity. The frequency-hopping approach approach (FH) [21,22] consists in processing monochromatic data, from the low to the high frequency, where the initial estimate of the contrast distribution at higher frequencies is given by the final result obtained at the previous sequence. Thus, P iterative monochromatic inversions are achieved separately. The multiple-frequencies inversion method (MF), on the contrary, minimizes a cost functional involving the entire data set [23,24]. Both approaches have the same unknown-over-data ratio, L × P + 1 unknowns in Ω for L × P data in Γ.
For self-consistency, we sketch briefly the MF inversion scheme, the frequency hopping approach being easily deduced by taking only one wavelength in the data set. The cost functional of the MF reads, where normalizing coefficients are defined as The subscripts Ω and Γ are included in the norm . and later in the inner product ., . to indicate the domain of integration. The residual errors h (1) l, p,n and h (2) l, p,n are defined as In the HM, two sequences related to the contrast and to the fields inside the investigating domain, χ n and E l, p,n , respectively, are built up according to the following recursive relations E l, p,n = E l, p,n−1 + α l, p,n v l, p,n + β l, p,n w l, p,n , where v l, p,n , w l, p,n and d n are updating directions with respect to the total field E l, p,n and to the contrast χ n , respectively, and α l, p,n , β l, p,n , and κ n are scalar coefficients. The updating directions v l, p,n and d n are chosen as the standard Polak-Ribière conjugate-gradient directions [25], while w l, p,n is given by w l, p,n =Ẽ l, p,n−1 − E l, p,n−1 , whereẼ l, p,n−1 represents the total field inside the investigating domain Ω, calculated from Eq. (1) with contrast χ n−1 . The scalar weight α l, p,n , β l, p,n and κ n are chosen at each iteration step n so as they minimize the normalized cost functional mentioned in Eq. (4).
The sample under test is assumed to be dielectric and non-magnetic. Subsequently, a priori information stating that the contrast permittivity distribution is real and positive. This is introduced by retrieving a real auxiliary functions ξ n in stead of χ n , such that χ n = ξ 2 n . The recursive relation with respect to the contrast χ n , Eq. (9), is changed to ξ n = ξ n−1 + β n;ξ d n;ξ .
As updating direction d n;ξ , we take d n;ξ = g n;ξ + γ n;ξ d n−1;ξ with γ n;ξ = g n;ξ , g n;ξ − g n−1;ξ Ω g n−1;ξ where g ξ is the gradient of the cost functional F(ξ, E ·, ·,n ) with respect to ξ, evaluated at the (n − 1)-th step, assuming that the total fields inside the investigating domain do not change. The gradient is given by where the over bar denotes the complex conjugate, and (A Ω p ) † and (B Γ p ) † are the adjoint operators of A Ω p and B Γ p , respectively. The updating directions v l, p,n for the total field inside the test domain is similar to that chosen for the contrast function ξ: v l, p,n = g l, p,n;E l, p + γ l, p,n;E l, p v l, p,n−1 , γ l, p,n;E l, p = g l, p,n;E l, p , g l, p,n;E l, p − g l, p,n−1;E l, p Ω g l, p,n−1;E l, p 2 Ω , where g l, p,n;E l, p is the gradient of the cost functional F(ξ, E l, p,n ) with respect to the field E l, p , evaluated at the (n − 1)-th step, assuming that ξ does not change:

Synthetic and experimental results
The object under test is a resin star with six branches of length 936 nm and height 130 nm radiating from a disk of diameter 208 nm and deposited on a silicon substrate. The resin, Hydrogen silsesquioxane, is not dispersive in the visible range [26] so that the contrast χ can be assumed to be constant whatever the wavelength. On the other hand, the silicon substrate is dispersive and the variation of its permittivity has been accounted for in the inversion scheme. The Scanning Electronic Microscope (SEM) image is given in Fig. 2(a). The image obtained through a bright field microscope is presented in Fig. 2(b). In the following, the investigation domain Ω will be restricted to a box of dimensions 6 µm× 6 µm and 600 nm (mesh size 50 nm). This a priori information is important as it limits the number of unknown and permits a reduction of the number of illumination and observation angles.
The object is successively illuminated by 40 plane waves with 20 illumination angles, polar angles ranging from θ inc = 16 • up to 59 • and azimuthal angles φ inc distributed within [0, 360 • ] and two different polarizations. The projections of the incident wavevectors on the transverse plane are plotted in Fig. 3. It is seen that they roughly fill the numerical aperture of the objective. Combined with the 31329 pixels of observation, they were found sufficient for imaging samples contained in the box Ω. Note that for more generality, the illumination transverse wavevectors do not follow the symmetry of the sample under study.

Synthetic results
To investigate the interest of multi-wavelength data in a reflection imaging configuration, without being disturbed by experimental noise issues, we first consider synthetic data. We simulated the field diffracted by the sample previously described in the same conditions as the experimental configuration (the mesh size for the direct and inverse problems are different for avoiding the inverse crime). Figure 4 displays transverse and axial cuts of the reconstructions obtained with HM from monochromatic data at two different wavelengths. At λ = 475 nm, the resin branches are well retrieved in the (x, y) plane, but strong ghost images in the (x, z) plane deteriorate the quality of the reconstruction. At λ = 675 nm the branches are poorly retrieved (their dimensions are significantly underestimated), but the ghost images are weaker and their distance from the actual target is increased. Both features can be explained with an analysis (under single scattering approximation) of the sample Fourier domain (support of the Optical Transfer Function) that is scanned with a monochromatic imaging system in reflection when the NA of the illumination coincides with that of the observation. The latter is the top-slice of the Ewald sphere of radius 2k 0 , where k 0 = 2π/λ, that does not include any axial frequency below 2k 0 1 − NA 2 [4], Fig. 5(a). Note that this domain is obtained for a reflection system without a highly reflecting substrate. In presence of a mirror, more Fourier components of the sample are theoretically accessible [27] but they do not change the conclusions.
First, it is seen that the high spatial transverse Fourier components of the sample is better recovered at small wavelength than at long wavelength. Using a conventional Fourier inverse algorithm, the branches reconstructed at λ = 675 nm appear fuzzier than that obtained at λ = 475 nm. our technique which includes a positivity prior, the lack of information at λ = 675 nm on the branch edges and tip (which generate the high spatial components), translates into an underestimation of the structure length and width.
The second important characteristic of the OTF support is the missing horizontal band. The latter has a major impact on the reconstruction in the axial direction as it prevents the distinction between the spectrum of a horizontal plane, which is constant along k z , and that of a z-oscillating structure, which is constant on the available segment of the OTF, k z ∈ [2k 0 1 − NA 2 , 2k 0 ] and null outside [23]. This indetermination explains the presence of the ghost images observed in Fig. 4  Combining the different monochromatic data is thus expected to improve significantly the reconstruction, especially in the (x, z) plane as it reduces the central missing band, see Fig. 5(b). In practice, with multi-wavelength data, the image of the actual target should be enhanced, as its position does not depend on the wavelength, while the ghost images should be dimmed, as they appear at different distances for the different wavelengths. Now the issue is to determine the best way to process the multi-wavelength data. In Fig. 6, we plot the reconstruction obtained with the FH and MF algorithm. Clearly, the reconstruction is better resolved and the ghost images better dimmed with MF than with FH.

Experimental results
We now turn to the experimental data. To give an idea of the noise, we display in Fig. 7 a comparison between the synthetic and experimental diffracted field for a given angle of illumination and wavelength. Actually, the noise level is strongly dependent on the illumination wavelength. The error between the experimental and synthetic diffracted field at the p−th wavelength defined as is given in table 1. Following the study conducted on synthetic data, we first provide in Fig. 8 the reconstruction obtained with the smallest available wavelength λ = 475 nm. As expected, we obtain an unsatisfactory reconstruction in the x, z plane with an important ghost image. We then display the reconstruction obtained with the frequency hopping [21,22] technique (FH) in Fig. 9(a). We observe that the reconstruction is very noisy and heterogeneous. This disappointing result can be understood by noting that the data at the largest wavelength λ = 675 nm are very noisy compared to the others. These data being used as the starting point in the hopping procedure, the whole chain of reconstruction is spoiled and trapped in false solution.
To verify this interpretation, we have performed the FH approach by accounting only for the 4 wavelengths presenting the same level of noise. The reconstruction shown in Fig. 9(a) is indeed much better than that of Fig. 9(b). Yet, the branches of the resin star are still strongly inhomogeneous and noisy-like. In Fig. 10, we present the reconstructions obtained with the Multiple Frequencies (MF) Inversion scheme without and with the largest wavelength λ = 675 nm. We observe that both reconstructions are satisfactory, all the branches being homogeneous and the reconstructed relative permittivity being close to the actual value of 2. Contrary to the FH approach, the noisy measurement at λ = 675 nm had no deleterious impact on the reconstruction profile. The reconstruction errors obtained using the MF or FH methods with four or five wavelengths defined as  They confirm that processing the multi-wavelength data set as a whole should always be preferred to successive monochromatic inversions, both for noise robustness and overall performance.

Conclusion
In conclusion, we have shown that using several wavelengths in a reflection tomographic microscope is an asset for 3D imaging as it permits a substantial decrease of the ghost objects in the reconstruction. Yet, the accuracy of the reconstruction depends strongly on the inversion scheme. Multiple Frequency Inversion algorithms which account for the whole multi-wavelength data set in the cost functional should be preferred to frequency hopping techniques that perform successive monochromatic inversions. In this work, we did not account for sample dispersion in the reconstruction scheme. An important perspective of this work would be to introduce a dispersive model on the permittivity contrast, which would add only a few scalar unknowns to the inverse problem. This extension would be particularly useful for studying absorptive samples, which are often more wavelength dependent than phase objects, and help characterizing natural or artificial samples [28][29][30].

Disclosures
The authors declare that there are no conflicts of interest related to this article.