Considerations on the proposed linear theory of surface measurement for coherence scanning interferometers

It was suggested in [Appl. Opt. 52, 3662 (2013)] that the result of a measurement via coherence scanning inter-ferometry could be viewed as the convolution of a point spread function of the instrument and an open surface in 3D space that lies at the air/material interface over a portion of the object ’ s surface. Further, it was suggested that by measuring certain objects, such as ones that are very close to spherical, and whose surface is known to a sufficient level of accuracy, that a point spread function for the instrument could be determined from the measurement result. We conclude that the approximations used in this calculation do not give sufficient accuracy to allow this to be achieved, and that the truncation of the surface function from the closed surface surrounding the object is not defined sufficiently well in order to give a unique solution to the problem. The physical justification for the truncation of the surface in this manner is also questioned.


INTRODUCTION
When an object is sufficiently weakly scattering so that the firstorder Born approximation and Fourier optics are applicable [1], then the imaging process of 3D optical microscopes can be described using a transfer function (TF) [1,2]. The TF is a complex valued function that describes how the imaging process modifies the relative amplitude and phase of each of the spatial frequency terms when mapping the spatial frequencies in the scattering potential to those in the image created, and will be equal to zero for spatial frequencies that are not passed by the instrument. Thus, if the scattering potential is known, it can be Fourier transformed (FT) to obtain the set of spatial frequencies it contains, and multiplied by the TF in order to get the result that would be obtained in spatial frequency space when imaging that object. An inverse FT of this result provides the measurement result in real space. Alternatively, this process can be viewed as a convolution of the point spread function (PSF), which is the inverse Fourier transform (IFT) of the TF, with the scattering potential to yield the measurement result in real space. This is a very quick and simple explanation of the process, and in reality such things as the coherence of the light used needs to be considered [1]; however, the essential points presented here are that there is a useful link between a property of the object being imaged, i.e., its scattering potential, and a function relating to the instrument alone, the PSF or TF, that allows the imaging process to be described. This is, however, only true in the weak scattering limit. Away from this limit, the relationship between the scattering potential and the image becomes far more complex.
In [3], the starting point of an incident plane wave being scattered from a strongly scattering, slowly varying, surface was considered. The scattered field on a surface far away from the object was calculated, reduced to those limited sections where a CSI would record the field, and then back-propagated in order to see what information is carried in this recorded field. The statement was made that this was then integrated over all the illumination angles, which appeared to yield the result that the measurement appeared to be equivalent to the convolution of the PSF of the instrument with a function that only takes a value at the air/material interface. This is equivalent to saying that the result of a measurement when viewed in spatial frequency space is given by the product of the TF of the instrument and the FT of the function defining the air/material interface. This is similar to the ideas for the weak scattering case described above, only with the surface function taking the role of the scattering potential. The surface function was, however, defined very loosely in [3], with a statement that a window function would be applied in physical space in order to limit it to a section of interest on the top surface. In [5], the surface function was broadened into a Gaussian form along the optical axis of the system in order to allow it to be represented by a discrete set of points in space to ease numerical calculations; this was then cropped to a region within a given angle from the optical axis, and then was multiplied by a Gaussian function perpendicular to the optical axis in order to remove the hard cutoff [22]. This loose definition leads to problems in separating the measurement result into the instrument-dependent and object-dependent parts.
In the following, we re-examine the claim that the result of a measurement of a strongly scattering object can be viewed as the convolution of the PSF and a function that corresponds to a region of the upper surface of the scattering object. In Section 2, we examine the work in [5], showing that the loose definition of the surface function as an arbitrarily windowed section of the full air/material interface cannot lead to an accurate answer, as slight variations of the surface retained significantly affect the determined phase of the PSF. In Section 3, we examine in more detail the mathematics in the paper by Coupland et al. [3] and draw different conclusions from those of the authors. It will be shown that if the surface function were to be filtered, that it should be done by angle with respect to the angle of incidence of the illuminating light, though this is still only a first approximation of the correct result. The section of the surface that information is recorded about therefore varies as the angle of the incident light, a result that corresponds to that which would be expected from geometrical optics. This angular dependence leads to the conclusion that the result is not the convolution of a PSF with part of the surface function, and in any case, all the results are deemed likely to be too approximate in order to allow the calibration of an instrument.

IMAGING LARGE STRONGLY SCATTERING OBJECTS
The work carried out in [3] attempts to discover a result for strongly scattering objects along the lines of that found for weakly scattering objects, where the measurement result can in some cases be shown to be the convolution of the scattering potential with the PSF. The authors claim not only to have achieved this, but to have subsequently demonstrated it experimentally [5] by measuring the distortion due to the instrument and applying a correction, allowing the surface of the object to be recovered more accurately. We believe that the conclusions of [3] are erroneous, and that the results shown in [5] demonstrate nothing more than reasonable shift invariance and some noise in their system; the object that they conclude has been measured more accurately is the same one as is used to calibrate the system; thus their calibration of the system just creates a mapping from the measurement results to those that they expect to see, and so it is not a surprise that they subsequently recover this object. As the correction is applied in spatial frequency space, the fact that the object still looks like the one it was calibrated to look like after a shift just shows shift invariance, the slight differences that they obtain can be attributed to a combination of noise and that the shift invariance is not perfect.
One significant problem with the work in [3] is that the surface is not defined uniquely, either in physical space or spatial frequency space; it is just stated that a window function is used to limit this to a measurement of a region of the surface that is of interest. This is problematic for a couple of reasons, one of which is that the result of a measurement is given as [3] Ok 2ΔkH k; (1) whereÕk is the measurement result in spatial frequency space,Δk is the FT of the surface function (corresponding to a cropped section of the total surface, i.e., only the "region of interest" of the air/material interface) andH k is the TF of the instrument that corresponds to the FT of the PSF. IfH k is a fixed property of the instrument and, as in the absence of noise, there is only one result for the measurement of an object Ok, then Eq. (1) requires that the function related to the surface of the object,Δk, has well-defined values for its FT, at least at points whereH k is not equal to zero. If the surface, Δr, were changed, then in general its FT within the support of the TF would change, and the only way in which the measurement result could be remain the same is thatH k must change, which it cannot if it is a fixed property of the instrument. The scattering of light from objects by considering "shadowed" regions and "illuminated" regions has been considered elsewhere for a vectorial case where the transition region from illuminated to shadowed region was found to have a width of the order 2∕kC 1∕3 C m, and where C is a typical radius of curvature of the surface [23,24]. Thus, for a plane wave illumination, the method of considering the illuminated section of the surface as that which contributes most is valid, but this is an approximation. There does not appear to be a clear answer to the question of how to define a surface section that could be FT'd to give a useful result, even for the case of a single illumination angle, and it will be shown later that this section of surface will vary with the illumination angle.

A. Differences as the Surface Section Changes
Before continuing on to the next section, where we will examine the mathematics that was used in [3], we will illustrate that even if the surface function were able to be defined in the manner proposed in [3], then even a slight variation in the section of the surface that we choose to retain leads to significant variation in its FT,Δk. Hence, if a calibration were to be carried out, which consists of calculatingH k from a measurement result, Ok, and the FT of the surface function corresponding to the known object being measured,Δk, such as was proposed in [5], then slight differences in the surfaces used will lead to significantly different results. We stress at this point that we are not saying that Eq. (1) is correct; indeed, we will show later that the surface function of this form is not the correct object function to use; however, as it has been used in prior work, an illustration of its deficiencies seems apt.
There is one final point to make, at the risk of confusing the reader, and this is that while Eq. (1) is not correct, and is at best an approximation of the correct answer, an equation of this type will exist. The only issue is determining the meaning of the function that is multiplied byH k. To see this, the reader just needs to consider the fact that the measurement result can always be FT'd to give a value ofÕk, and that asH k describes the spatial frequencies that the instrument is capable of passing and how they are affected by the instrument itself, then it will have a value independent of the measurement. A relation between the twoÕk F kH k can be found. The meaning of the function, F k, and its relation to the scattering object is not clear, though, and is likely to be highly complex, especially when multiple scattering is present. In the case where the first-order Born approximation is sufficient, then a relation with the scattering potential, and hence the object's refractive index profile, can be found, but this relation is lost when the higher-order terms are needed to describe the scattering. Thus, the fact that there is a function that can be convolved with the PSF in order to produce the measurement result in real space should not be a surprise, as the convolution theorem shows that Or IFTfF kg PSF.
In Fig. 1, a slice of the absolute value of the 3D FT of two different caps of the same infinitesimally thin spherical shell is shown, passing through the symmetry axis k z . These are calculated in the manner described in [6]. Ignoring a multiplicative factor, these objects would take the role of the surface function in Eq. (1). The red outline in the figure marks the limit to the nonzero region of the ideal TF for an instrument with NA 0.5 when the illuminating light is from 500 to 650 nm [25]. Both surfaces correspond to a section of a spherical shell of radius 50 μm, with the section in Fig. 1(a) corresponding to the cap of the spherical shell within a half angle of 30 deg about the symmetry axis, and Fig. 1(b) a half angle of 35 deg. That there are significant differences can be seen by looking at the line of points along the k z axis. It was shown analytically [6] that for these objects, zeros occur in the FT at points along the k z axis separated by 2π∕r 0 1 − cosθ, where r 0 is the radius of the sphere and θ is the half angle within which the cap of the sphere lies. This shows that the location of the zeros depends upon the size of cap that is chosen and, in fact, the phase and amplitude of the FT varies across the region within the TF as the cap size changes.
While the difference in size between the two caps in Fig. 1 was significant, in order to demonstrate the changing FT clearly, even a small change can lead to a significant change in the phase of the elements in the FT. In Fig. 2, the difference in the phase between the spatial frequency elements in the FT of the cap of a sphere of radius 50 μm, within the support of the TF for an instrument of NA 0.5 when the angle subtended by the cap is changed slightly is shown. The sphere has a radius of 50 μm, and is initially reduced to the section that would return light according to geometrical optics, i.e., the cap of a sphere subtended by an angle of 30 deg. This is then compared to the cap of a sphere that subtends an angle of 30.66 deg (corresponding to that imaged according to geometrical optics when the NA is 0.51). The difference in the angle is shown in Fig. 2. This plot shows the magnitude of the difference across a slice through spatial frequency space that passes through the k z axis, about which the results are rotationally symmetric. The phase angle plotted is the magnitude of the minimum possible change in angle, and so is limited to the range 0 to 180 deg. Changes greater than this are mapped into this range, and no difference between a phase lag or advancement is recorded. It can be seen in Fig. 2(a) that the greatest differences in phase occur along the k z axis or very close to it; however, looking at part (b), which just shows the region in part (a) that is away from the k z axis, shows that phase changes of around 12 deg are common. Such errors in a calibration measurement would be associated with the TF, and may lead to significant errors subsequently determined topographies, e.g., by the frequency domain analysis method [11]. Now that it has been shown that the loose definition of the surface used in previous work cannot be sufficient to achieve the stated aims, the next section will examine the assumptions that have been made in [3] when deriving the apparent result and come to a different conclusion to that stating the measurement result is the convolution of the PSF with the surface function.

APPROACH FOR STRONGLY SCATTERING OBJECTS
In [3], the scattering from an object that is strongly scattering, smooth, and slowly varying with respect to the wavelength of the illuminating light is considered. The incident field is a plane wave with unit amplitude, and the scattered field is calculated on a portion of a spherical boundary a long distance away from the object, and corresponds to the portion of the scattered field that would be collected by a given microscope, as is illustrated in Fig. 3.
The scattered field in the region outside of the surface S 2 is given by where the integral is over the surface S 2 which in the following corresponds to the surface of the scattering object,n is a surface normal pointing out from surface S 2 , Gr − r 0 is the free space Greens function, Gr − r 0 expikjr − r 0 j∕jr − r 0 j, and Er 0 is the total field at the surface S 2 , which is equal to the sum of the incident field E i r and the scattered field E s r. A derivation of this result is included in the appendix to allow the reader to clearly see the range of applicability, and the assumptions that are needed. It is at this point in [3] where approximate solutions [26] are introduced for the electric field and its gradient at the object's boundary when the incident field is a plane wave of unit amplitude, E i r exp ik i ·r . These are and which are substituted into Eq. (2). Here R is the Fresnel reflection coefficient of a smooth plane, k i is the wave vector of the incident wave, and r 0 is a position vector corresponding to the values of r at the surface of the object. Equation (3) is just the superposition of the incident and reflected fields from a plane interface, while Eq. (4) is given by differentiating both these fields [26]. Looking at Fig. 3, it can be seen that these expressions will not be valid for the entire boundary, only for the sections of the surface where the wave is incident on the object from the air into the material, and only when the reflection from a plane is a good local representation of the reflection from the surface. In addition, when the incident wave is blocked from falling on a region of the surface by another region of the scattering object, then the field and its gradient will not be given correctly, nor will it be correct if multiple scattering occurs. In [26], Fig. 3. Plane wave E i r is incident upon a strongly scattering, smooth, slowly varying object leading to a scattered field E s r. S 2 is a surface bounding the scattering object, and the arc of a circle represents the portion of a spherical boundary a long distance from the scattering object on which the scattered field is calculated. The distance to the distant boundary must be sufficiently great that the far-field approximation for the Greens function is sufficiently accurate.

Research Article
Vol. 56, No. 10 / April 1 2017 / Applied Optics where these are introduced, they are only used over sections of the surface where this condition would be met, and it is noted that this is an approximate method; indeed, comments are made regarding the truncation of regions of the surface and the use of a "shadowing function" [26].
The authors of [3] substituted Eqs. (3) and (4) into Eq. (2) and used the far-field approximation of the Greens function, Gr b − r expik i r b ∕r b expik irb · r, which implies that the locations of the points r b lie at a distance from the origin that is much greater than both the size of the object and the distance of the object from the origin. This results in being obtained, where k i jk i j and where r b is a position vector of magnitude r b , from the origin to a distant point where the field is being evaluated. The surface that the scattered field is calculated on should not be confused with the bounding surface, S 1 , mentioned in the Appendix A, within which Eq. (2) is valid. As k irb is a wave vector of magnitude k i in the direction given byr b , it will be denoted k s from now on. A limit on the integral is then introduced in [3], reducing it from the full surface to a region on the upper surface made through the use of the function Ar where Ar W r x ; r y δr z − sr x ; r y : Here W r x ; r y is a window function, and δr z − sr x ; r y is a Dirac delta function For E s r b to be the same in Eqs. (5) and (7), the integral over the portion of the surface that is removed must be equal to zero, though that this would be the case is not apparent. If it is not the case, then the E s r b given by Eq. (7) is at best an approximation of the value given by Eq. (5). That said, our experience of light scattering from objects in everyday life implies that this may not be a terrible approximation; after all, if the nonilluminated side of a large strongly scattering object, such as the back surface of a mirror, is changed slightly, we would not expect to see any real change in the light scattered from the front, and indeed, for large objects, ideas from geometrical optics should apply. In fact, as mentioned earlier, the division of the surface of a large scatterer into illuminated and shadow regions has been pursued before [24]. However, it is noted that the transition region is not simple, and only approximate results will be obtained.
If the window function is ignored in Eq. (6), this follows the route laid down in [26] and [27]. However, there is it noted that the termsn s · z in the denominator of the term included corresponds to a projection of the surface area onto a plane when a surface integral is carried out. This surface integral is then converted into a 3D integral to avoid nonlinearities in the simple surface integral form, as was shown in [28]. This implies that these terms are related to the object that is being imaged, and not the TF of the instrument in which they end up [3].
Once an expression for the field on a portion of a distant spherical boundary has been found, the authors of [3] backpropagate a section that corresponds to the region over which a microscope would collect the light in order to see what information remains in the recorded section of the field, and subsequently make a couple of simplifications to the expression obtained. In our view, it is more illuminating to make these simplifications before backpropagating the field.
If we return to Eq. (5), the principle of stationary phase is then applied, as in [3], whereby it is noted that the phase of the exponential term changes in the direction given by k s − k i , and is constant perpendicular to this direction, and so sections of the surface whose normal is in the direction k s − k i should contribute most to the integral. Finally, it is also noted in [3] that in those regions, the term k s k i ·n s is negligible. This leaves At this point, the consequence of the statement that only regions of the surface whose normal are in the direction k s − k i contribute to the integral should be explored further. It can be seen that this condition corresponds to specular reflection. This should not be a surprising result: we have specified that the surfaces are smooth, strongly scattering, and slowly varying, and that we can use the Fresnel equations to describe the scattering from the surface. Taking this to the limit, this would imply that it is only sections of the surface whose normal is in the range where specular reflection would scatter light into the angles over which light is collected that would contribute to the measured field, and therefore it is only these areas where there is any information recorded. An important point to note is that the section of the surface that contributes to the measured field is dependent upon the wave vector of the incident light, specifically the angle at which it is propagating. Within the approximations applied above, we can alter the surface outside of this range of angles in any way that we choose without affecting the measurement result, as long as a) we retain a smooth slowly varying surface, and b) we do not bring any of the surface normals into the range of angles that contribute to the measured field. The section of the surfaces marked in red in Fig. 4 is that which contributes to the measured field on the green portion of the spherical boundaries, while the blue sections can be modified without changing the measurement result. As the blue section of the surface does not contribute to the measured field, and is arbitrary apart from the condition that its surface normal must not fall within the range that would lead to the scattered light contributing to the measured field, then no information about it is recorded by the measurement. Thus, the process of measurement has acted as a filter on the surface, retaining only those sections of the surface that contribute to the measured field. This would be the justification for Eqs. (5) and (7) being equal, but only if Ar filters the surface by angle in the way described. Introducing Ar, calculating the scattered field, and then propagating the field back as was done in [3] misses the point that the measurement of only part of the field has acted as a filtering operation on Ar. It is for this reason that the surface in [3] appears arbitrary.
This introduces a further problem, which is illustrated in Fig. 4(b), which is that the section of the surface that contributes to the measured scattered field is dependent upon the direction from which the light is incident on the scatterer. Instead of A being a function of r, it should be a function of θ i ; ϕ i , the angles that the incident wave vector propagates in, and θ n ; ϕ n the angles that the surface normal is in. Integrating the result obtained for an incident plane wave over the full range of incidence angles then becomes more complex, as instead of the Ar term being constant, it varies over the integration. This invalidates the result gained later in [3] that the scattered field is the convolution of the PSF with a function dependent upon Ar.
If we examine the implications of this, it is easy to see that the expected results from geometric optics are recovered. Using the example of a microscope where the illumination is through the collection optics, as is typical in coherence scanning interferometry, and hence for which the angles that the green cap subtends not only limits the collection angles, but the illumination angles as well, then the range of illumination angles for which a point on the surface contributes to the measured field can be considered. Figure 5 shows the range of illumination angles for which the point marked with an "x" contributes to the measured scattered field. The range of angles that the light is scattered in will cover the same angular range, only propagating out rather than in. For the point at the top of the object, Fig. 5(a), light from all of the illumination angles contribute to the measured field; however, in part (b) the point on the side of the ball contributes only for a limited set of illumination angles, those between the two lines marked with arrows, and Aθ i ; ϕ i ; θ n ; ϕ n 0 for all the rest. As the range of angles that the light will be collected in is the same as the range of illumination angles, it can be seen that the greater the angle of the surface normal with respect to the optical axis of the system, the smaller a contribution it makes to the measurement result, and this now resembles the expected geometrical optics result. This is of course only an approximation of the true result, as numerous approximations have been used to get to this point.

DISCUSSION AND CONCLUSIONS
Previous work suggested that the result of a measurement of a strongly scattering object was given by the convolution of the instrument's PSF with a function that took a value only at a portion of the air/material interface at the surface of the object [3]. In addition, it has been claimed that the calibration of an instrument has been carried out to improve the accuracy of the topography calculated [5]. The results are gained in this prior work only after several approximations have been made, and where the function that describes the surface of the object being measured corresponds to a poorly defined region of the total surface.
In the work presented here, we have carefully examined the route by which the authors obtained this result, and concluded that they did not take into account the filtering effect that the measurement had on the surface function. When this is taken into account, the section of the surface obtained is dependent upon the direction the incident light is propagating in. By considering the effect of this, a result far closer to that corresponding to geometrical optics is obtained. However, it should be noted that this is still only an approximation of the true result, including such things as approximate values for the field and its gradient at the object's surface and the truncation of the full surface integral, and will not have sufficient accuracy in order to allow a microscope to be calibrated.
In addition, it was necessary to address the claims by the previous work that an instrument had been calibrated [5]. We conclude that the paper demonstrates only reasonable shift invariance of their system, as the object that they applied (a) (b) Fig. 5. Illustration of the illumination angles that contribute to the scattered field for the points marked by an "x" in part (a) and (b). The point whose surface normal is aligned with the optical axis of the system, shown by the "x" in part (a), contributes to the measured field for all the illumination angles within the half angle of the lens in a coaxial illumination. In part (b), the point marked with an "x" only contributes for a limited set of illumination angles, those between the two lines marked with arrows.
(a) (b) Fig. 4. Illustration of how the section of the surface that contributes to the measured field varies with illumination angle. Under the approximations mentioned in the text and for the plane wave with wave vector k 1;2 , in part (a) and (b), respectively, it is only the sections of the surface marked in red that contribute to the measured field on the green portion of the distant boundary. The blue section of the surface could be modified without affecting the measured result, as long as the surface normal is not brought into the range corresponding to the red section.
corrections to was the same as the object used to calibrate the instrument. The calibration therefore just creates a mapping from the measurement data to that which they expect to see. We have demonstrated that even if the previous work were correct, the FT of the object function is very sensitive to changes in the size of the region of the surface retained, and this fact, combined with the fact that there are numerous approximations in the prior work, implies that the result obtained by this method would not be sufficiently accurate to calibrate an instrument. In any case, our examination of the mathematics in [3] implies that the surface should not be defined in the manner that was used. As a final note, we hope that the detailed considerations about the method should allow any researchers who wish to build on the work in [3] to quickly understand its strengths, weaknesses, and limitations. We would, however, highly recommend that they look at the more complete work in this area that is in the literature, including work on such things as the Stratton-Chu integral [29].

APPENDIX A
To give clarity to the limits of applicability of the expressions used in the work presented here, a derivation of the expression for the scattered field solely in terms of surface integrals is included here. This is a scalar approximation and very much follows the routes laid down in Chap. 13 of [30] and in [3]. It is worth noting that the result obtained has much in common with the Stratton-Chu integral [29], albeit for a scalar wave instead of the full vectorial solution of Maxwell's equations. This vectorial solution to the problem has been used many times to propagate fields [31][32][33][34] that are due to scattering by objects, such as is being examined in this work.
An expression for the scattered field for points that lie between a surface S 2 , which encloses a scattering object, and a spherical surface with radius R, S 1 , as shown in Fig. 6, is obtained from the scalar wave equation as follows. Starting with the inhomogeneous Helmholtz equation [30] ∇ 2 k 2 Er −4πF rEr; (A1) where E is a single Cartesian component of the electric field, r is a position vector, k 2π∕λ is the wavenumber, and λ is the free space wavelength. The source terms on the right-hand side include the scattering potential F r k 2 n 2 r − 1∕4π, where nr is the refractive index at the point r, and Er is the total electric field, which is the sum of the incident field and the scattered field, Er E i r E s r. The incident field is taken to be a plane wave, which is a solution of the homogeneous Helmholtz equation and a Greens function of the Helmholtz operator will now be introduced, where δ 3 r − r 0 is a 3D Dirac delta function. The Greens function is assumed symmetric with respect to r and r 0 . Subtracting Eq. (A2) from Eq. (A1) gives ∇ 2 k 2 E s r −4πF rEr: (A4) Multiplying Eq. (A4) by Gr − r 0 , and Eq. (A3) by E s r and subtracting gives [30] E s r∇ 2 Gr − r 0 − Gr − r 0 ∇ 2 E s r 4πF rEr − 4πδ 3 r − r 0 : (A5) Interchanging r and r 0 , integrating over a volume V R bound by a sphere of radius R and using Greens theorem to convert the volume integral on the left-hand side into a surface integral and rearranging gives Z Z Z The surface S 1 encloses all the scattering objects, as is shown in Fig. 6. The left-hand side of Eq. (A6) is equal to zero if r lies outside of S 1 , but is equal to E s r if r is within it. As F r 0 is zero outside of the scattering object, the region the volume is integrated over can be reduced to that of the scattering object, V S . If there are multiple scattering objects, the integral can be converted into a set of volume integrals, one over each of the scatterers. If the process above is repeated using Eq. (A1) instead of Eq. (A4) and with the bounding surface this time being S 2 , then the following result is obtained: Z Z Z (A9) which is valid between surfaces S 1 and S 2 . If the radius of S 1 → ∞ with its contribution to E s vanishing, which we will assume to be true, then the result, is obtained.