Surface Depth ‐ Mapping of Material via the Transport ‐ of ‐ Intensity Equation

: We present a new approach for a surface characterization based on the TIE method combined with the SEM. Experimental verification is carried out on the example of characterization of a crater on the surface of monocrystalline silicon (111). The approach is universal and can be used for any opaque object. It improves the robustness and stability of the quantitative phase re ‐ trieval process and has two important features. Firstly, it allows one to quantitatively retrieve the phase in a region of arbitrarily chosen dimensions. Secondly, phase retrieval process does not re ‐ quire the choice of boundary conditions.


Introduction
One of the leading areas in modern laser physics research is the study of the interaction of ultrashort laser pulses with the surface of various materials, since these processes are fundamental in modern medicine [1], femtochemistry [2], creation of new nano-and microstructures for electronics [3,4], photonics [5], etc. The obvious advances in such studies are linked to the development of solid-state physics and theoretical methods for its modeling [6]. Combined, these methods can be used to determine surface structure [7,8] and chemical composition [9], study the electronic structure of solids, acquire images with atomic resolution and study processes occurring on the surface. It is generally accepted that there is no universal experimental method, and a set of mutually complementary methods should be used to solve these problems. For example, through combining methods of scanning probe microscopy, scanning electronic microscopy (SEM), energy dispersive X-ray spectroscopy and Raman spectroscopy, one can obtain complete information about the surface structure and its characteristics [10][11][12][13].
Among the abovementioned methods a special place is occupied by SEM, a powerful method of physical surface analysis allowing the establishment of composition, structure and some other properties of near-surface layers of the sample in the form of changes of intensity distribution in the observation plane. In most cases the interpretation of such images can be explained by phase contrast effects [14]. In addition to the resolution SEM, the phase contrast method was historically established as a qualitative tool to improve the image detail of the object under study [15]. At the same time, there are methods relating to objects' phase retrieval which are based on the registration of several intensity distributions, acquired at displacement relative to the focus plane. They belong to the category of methods called inline holography [16][17][18], and one of them-based on the transport-of-intensity equation-will be discussed in this article.
The transport-of-intensity equation (TIE) is a well-established deterministic phase retrieval approach which can be used for quantitative measurements of bulk and surface characteristics of transparent media by analyzing phase distributions [19]. For instance, quantitative phase imaging of objects under study can be used to measure their refractive index [20,21], thickness, mass, birefringence and magnetic properties [22,23]. This paper proposes the use of TIE combined with SEM as an alternative to atomic force microscopy (AFM) [24,25] to address problem of surface characterization of solid opaque materials (on the example of monocrystalline silicon). The main result of this work is to provide a new approach to improve the robustness and stability of the phase retrieval process which can be used in case of uneven background and object illumination. It is shown that this approach solves the main difficulty when working with opaque objects, namely the need for a clear formulation of the boundary conditions.

Materials and Methods
The propagation of an electromagnetic wave with a wavelength , a complex field , , along the z axis in an isotropic medium is described by a paraxial differential equation , , -intensity of the wave, and , , -phase of the wave; ∆-two-dimensional Laplacian, k-wavenumber, / -longitudinal derivative.
The resulting equation is nothing but the real part of the Helmholtz equation and is called the transport-of-intensity equation. It determines the one-to-one relationship between intensity , its longitudinal derivative / and the phase of the electromagnetic wave.
There are many methods and approaches used to solve Equation (2). For instance, in [27] an approach based on the Helmholtz decomposition theorem, which allows the representation of the vector field as a scalar and vector potential, was proposed [29]. Since quantitative phase imaging via TIE belongs to the deterministic phase retrieval approach [30], it is necessary to restrict it only to cases where the field is free from phase singularities [31]. Under such constraints, the vector potential can be neglected and the TIE solution can be obtained using only the phase scalar potential without any generality restriction. Then the original Equation (2) is simplified to the standard Poisson equation written with respect to the phase.
One of the most popular methods for solving Equation (3) is a method based on expressing differential operators through the fast Fourier transform (FFT) [32]. In [33] it is shown that this approach bypasses the serious mathematical limitations in solving TIE which arise when using finite elements methods. Using the properties of the Fourier transform, the differential operators in Equation (3) can be written as [34]: where F и F −1 -operators of direct and inverse Fourier transforms, respectively; kx,y = 2πνx,y-frequency coefficients, νx,y-frequency grids.
In [35] it was proved that this transition preserves the correctness and uniqueness of the solution to an additive constant. It should be noted that this approach requires strict adherence to boundary conditions. Violation of boundary conditions leads to distortion of boundaries and overlapping of spectra in the whole selected area. Due to the periodic nature of the Fourier transform, periodic boundary conditions (PBC) must be satisfied, where the function value at the boundaries must cyclically repeat [36], before applying relations (4) and (5).
In general, the choice of the boundary conditions when solving TIE is a very difficult task, because the effect of the boundary conditions is independent of the chosen solution method. Due to the periodic nature of the FFT, the boundary conditions are implicit. However, this implication can be eliminated with an extended FFT method that uses a symmetrization rule [37]. This allows the use of special cases of Neumann boundary conditions (NBC, Equation (6)), Dirichlet boundary conditions (DBC, Equation (7)), PBC or any other mixed boundary conditions [38]. The correctness of using NBC, DBC and PBC when solving TIE using the extended FFT method has been proven by numerical analysis, which investigates the effect of boundary conditions on the accuracy of phase retrieval [36]. In this case, it is assumed that the intensity distribution turns to zero outside the selected area [35]. (6) where g is a smooth function on the boundary Ω; is the outward normal derivative.
where g is a smooth function on the boundary Ω. In all cases, the derivation of the boundary conditions follows from the law of energy conservation, which in a more general case is replaced by the law of intensity conservation which is assumed to be true for TIE [26,35]. According to [35,39], the law of intensity conservation will be satisfied only if the dissipation of energy flux through the selected circuit is negligible. In practice, however, this condition is met by increasing the ratio between the size of the selected area and the size of the object itself. This minimizes the change in overall illumination in the image plane in the selected area.

Experimental Demonstration
An experimental demonstration of the method was carried out by the characterization of an ablation crater obtained on the surface of single-crystal silicon with a crystallographic orientation (111). This crater was obtained by tight focusing of laser radiation with wavelength λ = 515 nm and pulse duration τ = 300 fs through a micro-objective with numerical aperture NA = 0.65 into a spot with radius R1/e ≈ 0.6 μm. Ablation was performed in a single-pulse mode, with the energy per pulse being E = 450 nJ.
Pre-characterization of the ablation crater was performed using Certus Standard V AFM («Nano Scan Technology» Ltd., Dolgoprudnyy, Russia), which has lateral spatial resolution <0.1 nm and vertical resolution <0.01 nm. AFM scanning was performed in the semi-contact mode using a NSG10 cantilever. Figure 1 shows a normalized 2D image of the crater under study with a depth profile. Image size is 6.05 × 6.05 μm and the spatial resolution is 108 nm. To demonstrate the applicability of TIE for three-dimensional measurement of the surface depth profile with atomic resolution, SEM images (Tescan Vega, Brno, Czech Republic) were used. A series of images was made with 3 keV accelerating voltage and 30 pA beam current. According to [40], this mode provides high spatial resolution and minimum depth of focus, and is optimal for this particular sample.
The series of images was formed by tuning the configuration of the electron column, namely by changing the working distance (WD) of the system (the distance between the lower pole tip of the lens and the focusing plane of the electron beam). The nominal WD value was ~11.9 mm. The focusing plane was displaced relative to the initial WD position by 100 μm with 5 μm steps, with a displacement accuracy of ~1 nm. The choice of the displacement range and steps was based on two limiting factors: to prevent high spatial frequency attenuation due to spatial coherence; and to ensure that the change in the intensity distribution was sufficient to calculate the longitudinal intensity derivative. Some

Results and Discussion
The AFM data has a non-normalized dynamic range and so before analysis it was normalized to a minimum (that is, the point with the greatest depth) in order to simplify the analysis. Figure 3 shows cross-sections of the depth profile along the 0X and 0Y axes through the zero-depth point. The horizontal dashed line indicates the silicon surface level, while vertical lines show the abstract intersection of the crater's inner surface with silicon surface level. The resulting crater depth is 432 ± 1 nm. Quantitative phase imaging requires a priori information about the wavelength of the radiation. For SEM the theoretical value of the wavelength will be equal to the de Broglie wavelength of the electron passing the accelerating potential where h-Planks constant, -mass of the electron, e-elementary charge, -acceleration voltage.
The boundary conditions and the size of the area under study are the parameters when solving a TIE. Since the boundary conditions of the problem to be solved are indefinite, we consider them in the framework of the Dirichlet-Neumann boundary problem. In turn, the choice of the size of the area under study is derived from the law of energy conservation, which requires a negligibly small value of energy flux dissipation through the selected contour. Since such a calculation is not possible in this case, we will empirically determine the optimal ratio between the area under study and the object size.  In particular, in Figures 4a-c the object is positioned so that it is almost in contact with the image boundaries, resulting in serious boundary artifacts that affect the accuracy of phase retrieval and the quantification of crater depth. Figures 4d-f show that the influence of the boundary conditions on the measuring of the crater depth and the contribution of heterogeneity (non-uniformity of the background value and object illumination) decreases as the area under study increases, at least visually. However, it will be numerically shown below that an increase in the area under study increases the error for some cases of boundary conditions while reconstructing the depth maps.
The accuracy of phase retrieval using TIE depends not only on the appropriately chosen boundary conditions, as demonstrated earlier, but also on the quasi-uniformity of the intensity distribution. It is worth noting that the background intensity values on the SEM-images have a normal Gaussian intensity distribution with a mean value of 45 and dispersion of 22.5. In practice, this condition is unfeasible, so it was ensured by methods of digital image processing. The use of such methods does not contradict the conditions of intensity singularity (presence of small and/or zero intensity values). Figure 5 shows 2D crater images obtained after SEM-images were processed digitally. The boundary conditions and TIE calculation are similar to those in Figure 4. It can be seen from Figure 5 that images with a uniform background resulted in a more stable TIE solution. As a result, the error of the depth profile reconstruction decreased. This finding is confirmed by the calculation of the root-mean-square error (RMSE) between depth maps obtained by AFM ℎ , and SEM ℎ , with TIE, according to The RMSE was calculated for all cases of boundary conditions and several values of the size of the area under study, which ranged from 4.875 to 15.6 μm. Its calculation was performed only within the area where the ablation crater was located, regardless of the size of the area under study used for the TIE. The resulting dependences of the RMSE on the size of the studied region used for TIE for different boundary conditions are shown in The results of the calculations show that before the digital processing of SEM images, the best result was achieved using Dirichlet boundary conditions. Figure 6c shows that the error value decreases as the size of the study area increases, while for no boundary (Figure 6a) and Neumann boundary (Figure 6b) conditions the error monotonically increases. If we speculate on the area size which could net the minimum RSME an extrapolation of the data from Figure 6c is needed. Assuming the dependence is linear the optimal size of the region under study for the Dirichlet boundary conditions would be ~50 μm. If we are considering the results obtained after digital processing, again Di-richlet boundary conditions give the best results. In this case, the error of the depth map measurement remains practically unchanged in the selected range of the areas under study (except for some points).
As mentioned earlier, AFM has lateral resolution <0.1 nm (X and Y directions), and longitudinal resolution <0.01 nm (Z direction). By comparison, the SEM has high lateral resolution (~2 nm [40]) and low longitudinal resolution [41]. Therefore, combining SEM-imaging with the TIE phase retrieval method proposed in this paper, we can obtain the depth map of objects with much higher longitudinal resolution. As seen in Figures  6a-c, the cumulative RSME of depth maps reconstruction was no more than 6% for all cases considered, and shows that the proposed method allows the carrying out of surface characterization compatible with that of AFM.

Conclusions
We presented a new approach for quantitative surface characterization based on applying TIE to SEM-images, and demonstrated its feasibility using the example of a crater on the surface of monocrystalline silicon. The proposed approach is universal and can be used for any opaque object, can increase the robustness and stability of the quantitative phase retrieval process and has two important features. Firstly, it allows one to retrieve the quantitative phase in an area with arbitrarily chosen dimensions. Secondly, no choice of boundary conditions is required for phase retrieval. We proved that if the intensity is positive in a given region of the object and converts to a constant value there is no violation of the intensity singularity condition. The solution of the TIE in such a case is singular to the additive constant in the absence of any boundary conditions. We have demonstrated that by applying this TIE solution to SEM-images it is possible to produce high-resolution surface characterization compatible with that of AFM. The cumulative error in depth map reconstruction was no more than 6% for all cases considered, but the best results were obtained by solving this problem using Dirichlet boundary conditions. Data Availability Statement: Data underlying the results presented in this paper are available from the corresponding author upon reasonable request.