Modeling of interference microscopy beyond the linear regime

Abstract. Coherence scanning interferometry (CSI), a type of interference microscopy, has found broad applications in the advanced manufacturing industry, providing high-accuracy surface topography measurement. Enhancement of the metrological capability of CSI for complex surfaces, such as those featuring high slopes and spatial frequencies and high aspect-ratio structures, requires advances in modeling of CSI. However, current linear CSI models relying on approximate surface scattering models cannot accurately predict the instrument response for surfaces with complex geometries that cause multiple scattering. A boundary elements method is used as a rigorous scattering model to calculate the scattered field at a distant boundary. Then, the CSI signal is calculated by considering the holographic recording and reconstruction of the scattered field. Through this approach, the optical response of a CSI system can be predicted for almost any arbitrary surface geometry.


Introduction
Detailed information about a part's surface topography is valuable, with a surface's roughness, texture, and form determining functional macroscopic properties such as friction, adhesion, lubrication, and wear. [1][2][3][4] Coherence scanning interferometry (CSI) is a high precision surface topography measurement method, achieving surface measurement noise levels in the subnanometer range for planar and relatively smooth surfaces 5 while also accurately handling much rougher surfaces. 6,7 As a CSI instrument scans, the interference between the light from the surface and reference mirror causes fringes to form, and these fringes are localized to the surface topography due to the low-coherence illumination used; by obtaining the fringes' coherence envelope and phase, the surface topography can then be estimated. 8 Common reconstruction methods include the envelope detection method, 9,10 which estimates the surface height from the coherence envelope's center or peak; methods that rely on frequency-domain analysis, 11,12 which obtain a more refined surface estimation by acquiring and combining both envelope and phase information through Fourier transform of the fringes; and the correlogram correlation method, 13 which through correlation to a reference signal can identify the location of coherence profile.
Reconstruction methods must rely on an assumed relationship between the measured field and the true surface topography, and the aforementioned common reconstruction methods all assume that the phase of the measured field at a point is proportional to the surface height at that point. Under elementary Fourier optics, it is assumed that the measured field's phase is a linear function of surface height for surfaces close to a plane, beyond which the relationship no longer exactly holds. 14 Nonetheless, CSI models based on this simple assumption can still *Address all correspondence to Rong Su, E-mail: rong.su@nottingham.ac.uk predict the main features of an interference signal, 15,16 and reconstruction methods that assume this are effective. 12 In particular, a CSI model based on this assumption has been used to simulate batwing effects, accounting for diffraction effects by convolution of the theoretical twodimensional (2-D) point spread function of the imaging system. 17,18 A model that is based on the Debye approximation 19,20 can calculate the electromagnetic field near to the focus of an aplanatic optical system for polarized light, but it does not account for multiple scattering (as noted elsewhere 21 ).
Some approximate surface scattering models are based on the Kirchhoff (or physical optics) approximation, [22][23][24] an approximation that requires that surfaces vary slowly on the optical scale. CSI models reliant on the Kirchhoff approximation include one used to model rectangular grating structures with periods much larger than the grating heights 20 and one based on threedimensional (3-D) imaging theory called the foil model. 25,26 The foil model can provide the CSI signal from a surface so long as the Kirchhoff approximation is satisfied, and the model is currently being used to predict the effects of reference mirror defocus, 27 provide new methods of instrument calibration, 28 characterize resolution from measurements of microspheres, 29 correct measured transfer characteristics, 30 and allow for lens aberration compensation in CSI. 31 The foil model also assumes that multiple scattering is negligible. 26 Despite the mitigating effect of a finite spatial frequency bandwidth, the effect of multiple scattering and loss of diffraction orders cannot be neglected and remains a problem for surfaces that are rough at the optical scale, or when coherent features such as vee-grooves or sharp edges are present. 32,33 For such complex surfaces, the CSI measurement process is fundamentally nonlinear, and consequently the linear reconstruction methods cannot reconstruct accurate surface topographies. Only an advanced reconstruction method that accounts for these effects could provide an accurate surface topography estimate for these surfaces, and such a method must be based on a rigorous scattering model.
Rigorous models for optical scattering typically use numerical techniques to solve Maxwell's equations exactly. A review on the field of computational electromagnetics (CEM) is beyond the scope of this paper, but some of the rigorous CEM methods include finite-difference timedomain, finite element methods (FEMs), boundary element methods (BEMs) (also known as method of moments), and rigorous coupled-wave analysis (RCWA). 34,35 An RCWA CSI model for high numerical aperture (NA) objectives has been developed for parameter determination of unresolvable etched grating structures. 36 However, RCWA approaches are typically only appropriate for periodic structures.
In order to solve the scattering problem for arbitrarily complex surfaces efficiently, a rigorous BEM-based optical scattering model has been chosen. In contrast to FEM's volume discretization, BEM solves linear partial differential equations along only the boundaries and, therefore, is in principle faster for surface scattering. The method used in this work is based on that of Simonsen, 37 the theory being developed earlier by Maradudin et al. 38 This method finds the field and its surface normal derivative along a surface by taking advantage of the Ewald-Oseen extinction theorem 39,40 and solves the subsequent set of inhomogeneous integral equations through conversion to matrix equations by appropriate spatial discretization of the integrals. Such an approach is formally exact, and accounts for surface plasmons, polarization effects, and structures, which contain overhangs and other complex re-entrant features.
The current BEM algorithm is restricted to surfaces that only scatter within the plane of incidence (i.e., 2-D), i.e., surfaces fully described by lines on the plane of incidence (x-z plane) and infinitely extended along the third dimension (y direction) perpendicular to the plane of incidence. The algorithm has been used to calculate the far-field scattering of a sinusoidal surface, and the result agrees with that of a calibrated scattering instrument. 41 The objective of this paper is to demonstrate a new way of modeling CSI images using a BEM-based rigorous surface scattering model and synthesizing images in k-space. The model considers the effects of multiple scattering and works for surfaces with arbitrary geometries. Since the current BEM model is only 2-D but the experimental instrument generates 3-D images, comparison between simulation and experiment is intended only to be qualitative, not quantitative. A novel fully 3-D BEM-like model has recently been developed to overcome these limitations 42 and will be used in the future to extend the proposed CSI model to 3-D.
One potential application of the proposed CSI model is to measure complex surfaces that are beyond the linear regime. One possible approach is outlined in Ref. 33: with some a priori information about a surface, iterative improvement of the surface topography through minimization of the differences between measured and modeled CSI images can provide accurate surface topography, even in regions typically not measured. Another application of this model is to build a virtual CSI instrument, allowing for the prediction of measurement uncertainty. 43,44 2 Modeling CSI The procedure that allows a BEM-based CSI model to generate fringe signal data for surfaces can be broken into a number of steps, which are shown in Fig. 1.

Choice of Inputs
A surface's coordinates along the lateral axis x and the optical axis z can be generated by an analytical function or numerically specified. Naturally, the surface described by these coordinates defines the boundary between two homogeneous mediums of different refractive indices, and for each medium, the complex refractive index must be specified. Next, the optical parameters, such as the NA of the lens, are chosen, providing the range of angles that the incident illumination can take and the acceptance angle for filtering of the scattered field. The polarization of illumination is selected, between either the transverse electric (TE) or transverse magnetic polarization (i.e., s-or p-polarizations). The illumination's broadband spectrum can be defined by a Gaussian distribution with a given mean wavelength and full-width at half-maximum (FWHM).

BEM for Surface Scattering
Once the surface, optics, and illumination have been defined, the broadband spectrum and the angles of illumination are sampled, and for each possible pairing of wavelength and incident angle, the surface field values and far-field scatter are found using the BEM method, which provides the total field and its surface normal derivative along the surface. The BEM model takes advantage of the extinction theorem to form a pair of inhomogeneous surface integral equations for the two media, which are then coupled together by the boundary conditions that the field and its normal derivative must satisfy along the media's interface. The extinction theorem used here can be considered equivalent to Kirchhoff's integral equation, both consequences of applying Green's second integral identity to the Helmholtz equation. 37,39,45 These equations can be solved computationally to find the surface "source" fields, from which the scattered field at any point can be found. For these, far-field scatter calculations to be accurate, the surface must be resampled equidistantly before the surface field values are found, with the resampling distance typically set to λ∕5 or smaller for illumination wavelength λ. To ensure that the same surface coordinates are used for each wavelength of light sampled from the spectrum, the smallest wavelength sampled is chosen to determine the resampling distance.

Synthesis CSI Signal
The CSI fringe is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 3 1 OðrÞ ¼ E s ðrÞE r ðrÞ Ã ; (1) where the scattered field in the far field E s ðrÞ, for position vector r, is demodulated by multiplying by the conjugate of reference field E r ðrÞ that is reflected from the reference mirror in a real system. The demodulation can be carried out in the spatial frequency domain, i.e., k-space, through a convolution of the far-field scatter and the conjugate of the reference field to calculate the values ofÕðkÞ ¼ FTfOðrÞg. The spatial frequency components of the scattered fieldẼ s ðkÞ can only be found on a spherical shell in k-space with a radius of 1∕λ. 46 The reference field is a spherical shell in k-space with the same radius. Both spherical shells are truncated due to the finite NA. 25 As shown in Fig. 2, the two truncated spherical shells are convolved as a consequence of the demodulation process, i.e., the reflected reference field shifts the scattered field values to higher spatial frequencies.
For broadband illumination, the scattered fieldẼ s ðkÞ for each illumination wavevector k inc must be found over a range of observation vectors k obs , and the set of scattered fields for each k inc summed. For each possible pairing of k inc and k obs (where jk inc j ¼ jk obs j ¼ k 0 ¼ 1∕λÞ, the complex scattered field value is iteratively added at the position k obs − k inc to any existing value at that position. The calculation of the far-field scatter and the demodulation process are repeated for each spectral component of the light source, weighted by the spectral density, which is Fig. 2 (a) Construction of k -space fringe dataÕðkÞ is achieved by adding the computed scattered field values for each k inc and k obs at the k obs − k inc position, after suitable weighting. In (b), for a specific k 0 , the nonzero field values ofÕðkÞ for a limited set of k inc are shown in black, illustrating how scattered field values (red) are shifted by the range of k inc for a limited acceptance angle. Note that here each black arc originates from a different set of scattered field values due to a dependence on k inc . usually assumed to be Gaussian-distributed. The fringe signalÕðkÞ is obtained through the superposition of the signal for each wavelength and angle of incidence. The CSI fringe image in real space is then given by the real part of OðrÞ ¼ FT −1 fÕðkÞg.

Methods and Materials
In order to verify the BEM CSI model, we compare the results from the model with those from experimental measurements. A range of prismatic surfaces were measured using a Zygo Nexview™ NX2 CSI instrument, see Table 1. Results for similar surfaces from a linear model of CSI can be seen in Ref. 28. In each case, a 10-μm scan along the optical axis (z axis) was performed using a 50× objective lens, which has an NA of 0.55 (acceptance angle of ∼33 deg), a Sparrow criteria optical resolution of 0.52 μm, a field of view (FOV) of (0.174 × 0.174) mm when using the 1.0× zoom lens, and from the 1000 × 1000 pixel FOV, a spatial sampling of 0.174 μm per pixel. Here, the circular illumination aperture fills the whole NA of the lens; i.e., the system's illumination NA is equal to the NA of the objective lens. This is confirmed by the experimentally measured 3-D surface transfer function of the same instrument configuration seen in Ref. 31. Due to the design of the objective lens, the polarization of illumination incident on the surface is nominally circular, with a small radial component resulting from the high NA of the objective lens. The signal data measured and recorded by the instrument, i.e., the intensity measured at each pixel for each scan position, are exported as a 3-D array of integers. The k-space fringe data are then obtained through the use of the 3-D Fourier transform, and a band-pass filter (BPF) applied to isolate the high spatial frequency fringe components. A subsequent inverse Fourier transform of the filtered signal provides the real-space experimental fringe data without low spatial frequency components. The BEM CSI model is provided with the corresponding curves that specify the real surfaces, e.g., a sinusoid for a sinusoidal grating and a horizontal line for an optical flat. Each surface in Table 1, with the exception of the vee-groove, has a length of 170 μm along the lateral direction, i.e., x direction, matching the FOV of the experiment, and the surface geometry was sampled with a spacing of 0.099 μm.
It was assumed that the light is incident upon a perfect conductor and can be treated as linear TE polarization illumination. The spectral density as a function of wavenumber is modeled as a Gaussian with a mean of 1.72 μm −1 , and FWHM of 0.24 μm −1 (corresponding to a mean of 0.58 μm and an FWHM of 0.08 μm) and approximately matching the corresponding parameters of the instrument. The real-space fringe signal OðrÞ is determined at coordinates with lateral and axial spacing that match that of the real instrument's signal data, i.e., using a lateral spacing of Step height with sharp edge NPL ACG-2.1 XP01 Step height: 2.1 μm Vee-groove N/A (only modeled) Vee-groove dihedral angle: 70 deg Depth: 10 μm 0.174 μm and an axial spacing of 0.071 μm, with 1000 lateral points and 205 axial points. This corresponds to a k-space grid spacing of 0.0058 and 0.0673 μm −1 for the lateral and axial directions, respectively. The spectrum is sampled 15 times over three standard deviations of the total spectrum (i.e., from 1.42 to 2.03 μm −1 ), and the incident angles sampled 18 times over the angles within the acceptance angle for the NA. Over the same range of angles, 1113 observation angles are chosen, at which the far-field scatter is calculated. Under these conditions, the CSI signal simulation took around 45 min for each surface on a PC with Intel ® Xeon ® E5-1620 v4 @ 3.50 GHz CPU and 64 GB RAM. However, for these surfaces, a reduction to nine incident angles halves this time with very little effect on the generated fringes. In addition, it is expected that more computationally efficient approaches, e.g., parallelization, would reduce this time considerably.
These comparisons are intended to be qualitative, not quantitative, due to the limitations of the BEM-based CSI model in matching the 3-D nature of the experimental measurement, e.g., the circular aperture of the instrument. Due to the commercial nature of the system, we were unable to adjust this aperture to a slit. Therefore, to allow for qualitative comparison, this difference is mitigated through the measurement of grating or prismatic surfaces, due to their strong scattering characteristics into the plane in which BEM is limited to.
Note that the definition used for spatial frequency k x (along the x axis), as seen in Fig. 4 and all subsequent k-space figures, is not the angular frequency definition, but the linear frequency definition (an unscaled reciprocal of space). This is done as the grating pitch can more easily be connected to the k-space diffraction patterns produced.

Results and Discussions
First, the experimental results when measuring an optical flat are compared with results from the model, as seen in Figs. 3 and 4. Qualitative agreement is achieved, showing that the model generates fringes that match the experimental results. The coherence envelope of the fringes slightly differs between the experimental results and those from the CSI model, which is expected due to the 2-D limitation of the BEM modeling and because the instrument's source spectrum is not exactly Gaussian. Differences in spectral distribution causing the coherence envelope to differ are shown in Fig. 4 of Ref. 28, where a linear 3-D CSI model based on the Kirchhoff approximation and CSI experimental measurement data are compared. It is well known that a narrower spectrum will lead to a broader coherence envelope. In Figs. 5 and 6, a high spatial frequency, low amplitude sinusoidal grating is used, as specified in Table 1. As expected from elementary Fourier optics, the diffraction orders produced are spaced relatively far apart due to the surface wavelength of 2.54 μm. Qualitative agreement between experiment and simulation is again seen, with the signals in both the real space and k-space appearing to match. However, the amplitudes of the higher diffraction orders, relative to the zeroth-order, are weaker in the experimental measurement. This is partly because the current BEM algorithm only considers in-plane illumination, but in the experiment, the zeroth-order will have a number of contributions from off-axis illumination that increase its magnitude, i.e., this effect can be considered as the difference between measuring a grating using a spherical lens and a cylindrical lens. In addition, this effect is partially caused by the imperfect transfer function of the instrument due to optical aberration and apodization due to the reference mirror in its Mirau objective. 29,31 In Figs. 7 and 8, comparison is made using a sinusoidal grating with higher amplitude and longer wavelength, as specified in Table 1. The surface wavelength of 50 μm causes the resulting Fig. 5 Cross-sectional CSI signal of a sinusoidal grating Rubert 543E, (a) experimental result and (b) corresponding simulation. Note that the fringes were (a) measured and (b) generated over the same FOV as in Fig. 3, but the display window has been shrunk for better visual comparison. pattern in k-space to be closer together, and the increased amplitude gives the more complex pattern seen here. Qualitative agreement is again seen for both the real and k-space fringe data.
The final comparison between the model and experiment uses a step height, as shown in Figs. 9 and 10, and as specified in Table 1. The experimental results and the modeled results in Figs. 9(a) and 9(b), respectively, differ around the step itself. This discrepancy likely occurs due to the inherent difference between the CSI model, which is restricted to surfaces that only scatter within the plane of incidence, and a real 3-D measurement. This can in part be seen through examination of the fringes after filtering of the out-of-plane k-space signal in Fig. 9(c).
This filter was a simple binary mask where only a central slice of the 3-D Fourier transform of the signal data was taken, and the rest discarded. The data were then inverse Fourier transformed, the real part taken, and the central slice displayed as the fringe data. This filter does not work exactly as a physical slit aperture (located at the back focal plane of the objective lens) but can nonetheless remove the out-of-plane scatter of the in-plane illumination, as well as part of the scatter of the out-of-plane illumination. Therefore, the filter should moderately improve the comparison.  Polarization effects introduced by the instrument's optical elements, and in particular the Mirau interferometric objective, 47 are not considered in our model, which could also contribute to this difference. In addition to these effects, the tilted fringes near the corner and the vertical wall of the modeled step height is probably caused by the double reflection between the two orthogonal surfaces, which are likely less pronounced in the experiment because the texture of the vertical wall is higher compared to the smooth surface assumed in the simulation. This discrepancy will be investigated in future work with a full 3-D BEM-based CSI model.
In addition to comparisons to the experimental measurements, in Fig. 11, the model's results of a 10-μm deep vee-groove with 70-deg dihedral angle are presented, where the sampling in  wavenumber and incident angle has been increased. The inverted "v" fringe pattern seen at the pit of the vee-groove is understood to be virtual images of the two vee-groove walls, generated by multiple reflection; 48,49 and the relationship (described in Ref. 48) that relates the dihedral angle of the multiple reflection image to the vee-groove dihedral angle appears to be satisfied here. The result also visually matches that found elsewhere. 33 This result, therefore, presents evidence that the model can account for multiple scatter.

Summary
In this paper, a rigorous model of CSI based on BEM is presented. The model accounts for multiple scattering effects and is a promising approach for generating fringes for arbitrarily complex surfaces. Such a model will eventually provide opportunities to accurately measure complex surfaces, through, e.g., iterative improvement of surface topography through the minimization of the differences between measured and modeled CSI images. Evidence of the model's validity is provided by comparison to experimental measurements from a commercial CSI instrument for several surfaces, giving qualitative agreement; modeling of a vee-groove provides evidence toward the model's capacity to account for multiple scattering. Fig. 11 (a) Simulated CSI signal in real space for a vee-groove as described in Table 1. Note that the blue dashed line denotes the geometry of the vee-groove modeled. (b) The magnitude of the kspace CSI signal. Due to the recent completion of the 3-D BEM-like surface scattering model, 42 the current 2-D CSI model can be directly extended to a full 3-D CSI model.