Modelling of coherence scanning interferometry for complex surfaces based on a boundary element method

Coherence scanning interferometry (CSI) is a well-established technique for measuring surface topography based on the coherence envelope and phase of interference fringes. The most commonly used surface reconstruction methods, i.e. frequency domain analysis, the envelope detection method, and the correlogram correlation method, obtain the phase of the measured field for each pixel and, from this obtain the surface height, by assuming the two are directly proportional. For surfaces with minor deviations from a plane, it is straightforward to show that the scattered field’s phase is a linear function of surface height. An alternative approach known as the “foil model” gives more generally the scattered field as the result of a linear filtering process operating on a “foil” representation of the surface. This model assumes that the surface slowly varies on the optical scale and that there is no multiple scattering. However, for surfaces that are rough at the optical scale or have coherent features (e.g. vee-grooves), the effect of multiple scattering cannot be neglected and remains a problem for reconstruction methods. Linear reconstruction methods cannot provide accurate surface topographies for complex surfaces, since for such surfaces, the measurement process of CSI is fundamentally non-linear. To develop an advanced reconstruction method for CSI, an accurate model of the imaging process is required. In this paper, 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 image formation as back-propagation of the scattered field, combined with the reflected reference field. Through this approach, the optical response of a CSI system can be predicted rigorously for almost any arbitrary surface geometry. Future work will include a comprehensive experimental verification of this model, and development of the non-linear surface reconstruction algorithm.


INTRODUCTION
Coherence scanning interferometry (CSI) offers high precision measurement of surface topography, capable of achieving a sub-nanometre noise level for planar and relatively smooth surfaces 1 . By using a broadband illumination source and a reference mirror, fringes that are localised to the surface topography are formed as the instrument scans, and the coherence envelope and phase of these fringes can be obtained and used to estimate the surface topography 2 . These reconstruction methods include among others the envelope detection method 3,4 , which obtains the centroid of the coherence envelope to estimate surface height; the frequency domain analysis method 5,6 , which uses both envelope and phase information for a more refined surface estimation; and the correlogram correlation method 7 , which identifies the coherence peak through correlation to a reference signal.
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 (EFO), 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 8 . Nonetheless, CSI models based on this simple assumption can still predict the main features of an interference signal 9,10 , and reconstruction methods that assume this are effective 6 . 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 two-dimensional (2D) point spread function of the imaging system 11,12 . A model that is based on the Debye approximation 13,14 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 15 ).
Some of the advanced approximate surface scattering models are based on the Kirchhoff-tangent (or physical optics) approximation [16][17][18] (KA); an approximation which requires that surfaces vary slowly on the optical scale. CSI models reliant on the KA include one used to model rectangular grating structures with periods much larger than the grating heights 14 ; and one based on linear systems theory called the foil model 19,20 . The foil model can provide the CSI signal from a surface so long as the KA is satisfied, and the model is currently being used to predict the effects of reference mirror defocus 21 , provide new methods of instrument calibration 22 , and characterise lateral resolution from measurements of micro-spheres 23 .
The foil model also assumes that multiple scattering is negligible 20 . 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 24,25 . For such complex surfaces, the CSI measurement process is fundamentally non-linear, 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 time-domain, finite element methods (FEMs), boundary element methods (BEMs) (also known as Method of Moments), and rigorous coupled-wave analysis (RCWA) 26,27 . A RCWA CSI model for high numerical aperture (NA) objectives has been developed for parameter determination of unresolvable etched grating structures 28 , 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 discretisation, 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 29 , the theory being developed earlier by Maradudin et al 30 . This method finds along a surface the field and its surface normal derivative by taking advantage of the Ewald-Oseen extinction theorem 31,32 , and solves the subsequent set of inhomogeneous integral equations through conversion to matrix equations by appropriate spatial discretisation of the integrals. Such an approach is formally exact, and accounts for both surface plasmons, polarisation effects, and structures which contain overhangs and other complex re-entrant features. This BEM model has been used to calculate the far field scattering of a sinusoidal surface and the result agrees with that of laser scatterometry measurement 33 . The current algorithm restricts the BEM model to surfaces that only scatter within the plane of incidence, 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. However, a full three-dimensional (3D) algorithm that can handle any arbitrary complex surfaces is under development.
The CSI model presented in this paper generates fringes based on the interference between the reference field and the scattered field, which is calculated through a BEM monochromatic scattering model for different illumination conditions that together comprise a broadband incoherent source. The image formation can be considered as a filtering and demodulation operation applied to the scattered field 20 ; more details are given in the following section.

MODELLING 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 illustrated in Figure 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 polarisation of illumination is selected, between either the transverse electric (TE) or transverse magnetic (TM) polarisation (i.e. s-or p-polarisations). 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 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 29,31,34 . 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.

Generating the measured field for broadband illumination
In order to reconstruct the scattered field along the surface, the measured scattered field in the far field ( ), for position vector , is demodulated by multiplying by the reference field that is reflected from the reference mirror in a real system, given by the conjugate of the illuminating field ( ) * . The reconstructed field is given by 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 reference field to calculate the values of ̃( ) = FT{ ( )}. The spatial frequency components of the scattered field ̃( ) can only be found on a spherical shell in k-space with a radius of 1/ 35 . The reference field is a spherical shell in k-space with the same radius. Both spherical shells are truncated due to the finite NA 19 . As shown in Figure 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 measured field ̃( ) for each illumination wavevector inc must be found over a range of observation vectors obs , and the set of measured fields for each inc summed. For each possible pairing of inc and obs (where | inc | = | obs | = 0 = 1/ ), the complex scattered field value is iteratively added at the position obs − 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 ( 0 ). The fringe signal ̃( ) 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 ( ) = FT −1 {̃( )}.

METHOD
In order to verify the BEM CSI model, qualitative comparisons were made between results from the model and those from experimental measurements. A range of prismatic surfaces were measured using a Zygo Nexview™ NX2 CSI instrument, see Table 1. In each case, a 10 µm scan along the optical axis was performed using a 50× objective lens, which has a NA of 0.55 (acceptance angle of ~33°), a Sparrow criteria optical resolution of 0.52 µm, a field of view (FOV) of (0.17 × 0.17) mm when using the 1.0× zoom lens, and from the 1000 × 1000 pixel FOV, a spatial sampling of 0.174 µm per pixel. The signal data measured and recorded by the instrument, i.e. the intensity measured at each pixel for each scan position, is exported as a 3D array of integers. The k-space fringe data is then obtained through use of the 3D fast Fourier transform (3D FFT), and a band-pass filter (BPF) applied to isolate the high spatial frequency fringe components. A subsequent inverse 3D FFT of the filtered signal provides the real-space experimental fringe data without DC 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. It was assumed that the light is incident upon a perfect conductor and arbitrarily use TE polarisation illumination. The spectral density as a function of wavenumber is modelled as a Gaussian with a mean of 1.72 cycles/µm, and FWHM of 0.24 cycles/µm (corresponding to a mean of 0.58 µm and a FWHM of 0.08 µm) and approximately matching the corresponding parameters of the instrument. The real-space fringe signal ( ) 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 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 cycles/µm and 0.0673 cycles/µm for the lateral and axial directions respectively. The spectrum is sampled fifteen times over three standard deviations of the total spectrum (i.e. from 1.42 cycles/µm to 2.03 cycles/µm), and the incident angles sampled eighteen 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. Each surface in Table 1, with the exception of the vee-groove, has a length of 170 µm along the lateral direction, i.e. xdirection, matching the FOV of the experiment, and the surface geometry was sampled with a spacing of 0.099 µm. Under these conditions, the CSI signal simulation took around forty-five minutes 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. Additionally, it is expected that more computationally efficient approaches, e.g. parallelisation, would reduce this time considerably. Step height with sharp edge NPL ACG-2.1 XP01 Step height: 2.1 µm

RESULTS AND DISCUSSIONS
First the experimental results when measuring an optical flat are compared with results from the model, as seen in Figure  3 and Figure 4. Good 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 given that the instrument's source spectrum is not exactly Gaussian, with similar issues seen elsewhere 22 .
In Figure 5 and Figure 6, a high spatial frequency, low amplitude sinusoidal grating is used, as specified in Table 1. As expected from EFO, the diffraction orders produced are spaced relatively far apart due to the surface wavelength of 2.54 µm. The agreement between experiment and simulation is again good, with the signals in both the real space and kspace matching. 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 apodisation due to the reference mirror in its Mirau objective 23 .
In Figure 7 and Figure 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 pattern in k-space to be closer together, and the increased amplitude gives the more complex pattern seen here. Good 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 Figure 9 and Figure 10, and as specified in Table 1. Agreement between experimental results and the modelled results in Figure 9 is good in the flat regions, but differs significantly 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 3D measurement; this can in part be seen through examination of the fringes after filtering of the out-of-plane k-space signal (not shown here). Polarisation effects introduced by the instrument's optical elements, and in particular the Mirau interferometric objective 36 , are not considered in our model, which could contribute to this difference. In addition to these effects, the tilted fringes near the corner and the vertical wall of the modelled 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.
In addition to comparisons to experimental measurements, in Figure 11 the model's results of a 10 µm deep vee-groove with 70° 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 37,38 ; and the relationship, described elsewhere 37 , 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 25 . This result, therefore, presents evidence that the model can account for multiple scatter.  . The k-space signal magnitude has been normalised to +1 in each case. Figure 5. Cross-sectional CSI signal of a sinusoidal grating after BPF filtering a), and corresponding simulation b). Note that the fringes were measured a) and generated b) over the same FOV as in Figure 3, but the display window has been shrunk for better visual comparison.  Proc. of SPIE Vol. 11057 1105713-8 Figure 9. Cross-sectional CSI signal of a step height obtained from cross grating sample after BPF filtering a), and corresponding simulation for a step height, assuming a step inclination of 90° b). Figure 10. Cross section of the magnitude of the k-space CSI signal from a step height found on a cross grating sample a), and corresponding simulation for a step height b). Figure 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 modelled. b) The magnitude of the k-space CSI signal.

CONCLUSIONS
In this paper, a rigorous model of CSI based on BEM is presented as an approach for generating fringes for arbitrarily complex surfaces. Current CSI models cannot accurately predict interferometric fringes for surfaces with complex geometries that cause multiple scattering. Here, a rigorous CSI model that accounts for multiple scattering is presented. Evidence of the model's validity is provided by comparison to experimental measurements from a commercial CSI instrument for a number of surfaces, giving good qualitative agreement.
In future work, it is planned to demonstrate clearly that this rigorous CSI model can provide accurate results for structures where multiple scatter visibly dominates, and results depart significantly from that predicted by traditional linear models.
In particular, it would be valuable to show that fringes can be accurately predicted for surfaces that are difficult to reconstruct using current surface reconstruction methods. Development of a three-dimensional version of the BEM model is already in progress, and once completed, a full 3D CSI fringe model can be constructed and results for various surfaces compared to the current model.