The point spread function in interferometric scattering microscopy (iSCAT). I. Aberrations in defocusing and axial localization

Interferometric scattering (iSCAT) microscopy is an emerging label-free technique optimized for the sensitive detection of nano-matter. Previous iSCAT studies have approximated the point spread function in iSCAT by a Gaussian intensity distribution. However, recent efforts to track the mobility of nanoparticles in challenging speckle environments and over extended axial ranges has necessitated a quantitative description of the interferometric point spread function (iPSF). We present a robust vectorial diffraction model for the iPSF in tandem with experimental measurements and rigorous FDTD simulations. We examine the iPSF under various imaging scenarios to understand how aberrations due to the experimental configuration encode information about the nanoparticle. We show that the lateral shape of the iPSF can be used to achieve nanometric three-dimensional localization over an extended axial range on the order of 10$\,\mu$m either by means of a fit to an analytical model or calibration-free unsupervised machine learning. Our results have immediate implications for three-dimensional single particle tracking in complex scattering media.


Introduction
Interferometric scattering (iSCAT) microscopy is a powerful label-free technique well suited for the sensitive detection and tracking of nanoscale matter such as individual viruses, colloidal nanoparticles, proteins and single molecules [1][2][3][4][5]. While iSCAT is devised for detection of weak scattering inherent to nanoparticles, it belongs to the broader family of interferometric microscopies which include holography, interference reflection microscopy, phase contrast microscopy and quantitative phase imaging [2,3].
In interferometric microscopies, the signal of interest is accompanied by an imaging background which oftentimes is that of a random speckle-like pattern. Such backgrounds present a challenge to remove, especially if they are dynamic. However, with appropriate subtraction or mitigation, ultra-weak signals can be identified down to the shot noise limit. With special sample preparation conditions even single molecules have been detected with iSCAT [6][7][8]. The high detection sensitivity of interferometric microscopy permits localization of single particles to high precision in all four spatio-temporal dimensions, and hence the technique has found application in tracking the fast nanoscale dynamics of various entities [5,[9][10][11][12][13].
Positional information which is encoded within the interferometric point spread function (iPSF) has been exploited through different methods. In digital holographic microscopy, for example, the positions of large beads are oftentimes extracted through numerical reconstruction of the hologram or through fitting based on Mie theory [14,15]. In the absence of a model or a fit, experimental point spread functions can be interpreted against a calibration data set, as is common within optical and magnetic tweezing microscopies [16][17][18]. It is only very recently that imaging above and below the focal plane within in-line digital holography has been modeled for an index-matched sample by explicit inclusion of the effect of the lens. Doing so opens the possibility to extend the depth of field to which one can image colloidal particles over a range of many microns [19].
In this work we set out to investigate the iPSF through both theory and experiment. We present a robust vectorial diffraction model to treat the imaging of a nanoparticle arbitrarily positioned about the focus in a wide-field reflection iSCAT microscope. Our model describes the image formation in the presence of aberrations due to the refractive index mismatch between the sample and coverslip, which is common in biological experiments. Furthermore, we discuss the origin of the remarkably long range and fine resolution of iSCAT particle tracking in the axial direction and explain how the iPSF behaves differently from conventional intensity-based PSFs in the diffraction-limited focal region. We also demonstrate how an understanding of the iPSF and application of the model can enable tracking the diffusion of lipids on a large artificial cell under physiological conditions. Finally, we present the use of unsupervised machine learning to extract the axial information encoded in the lateral iPSF profile.

The interference equation
Central to iSCAT microscopy is the principle that the electric field scattered from the nanoscopic scatterer of interest (E sca ) is superposed on the detector with some portion of a reference beam (E ref ), rendering the detected intensity as where φ dif = φ ref − φ sca is the differential phase accumulated between the phases of the scattered and reference fields at the point of detection. The interferometric nature of the cross-term (last term in Eq. (1)) permits detection of extremely weak scattering signals on top of a large background with sensitivity down to the shot noise limit. The most popular imaging modality for iSCAT is illustrated in Fig. 1(a), and is composed of a common-path wide-field illumination scheme with camera-based detection. In this mode, the reference beam consists of the reflection of the illumination at the sample interface. The most sensitive imaging is achieved with an oil-immersion objective of high numerical aperture. As it becomes clear below, it is important to realize that these components are typically designed for imaging objects directly placed on the coverslip.
The iPSFs as shown in Fig. 1(b) exhibit a characteristic series of rings of increasing size and oscillating contrast caused by the mismatch between the spherical scattered and the near-planar reference wavefronts. These features are important as they encode a wealth of information about the material of the scatterer and its three-dimensional position to nanometer precision. Although we recently benefited from this rich information for tracking transmembrane proteins [13], to date, there has not been a complete description of the iSCAT imaging system. In particular, the important role of omnipresent aberrations on iPSFs has been missing. Such an understanding is critical for accurate and quantitative interpretation of the information contained in iSCAT images.

How geometric aberrations contribute to the iPSF and the extended depth of focus
The complex iPSF at the plane of detection results from the differing wavefront curvatures of the reflected and scattered fields as well as the phase difference accumulated between them. Typically, the reference field is considered to have a planar wavefront generated by focusing the illuminating light beam in the back focal plane of the microscope objective (see Fig. 1(a)). The scattered light has a quasi-spherical wavefront emanating from the nanoparticle. Under normal conditions the scattered light is additionally subjected to multiple phase-altering aberrations that a nanosca�erer lens objec�ve beam splitter The sample is typically a sub-wavelength scatterer such as a gold nanoparticle, which is taken to be freely positionable with respect to the coverslip and the objective focus. The inset shows that the scattered light with approximately spherical wavefront mixes with planar reflections off the coverslip. (b) Exemplar images on the camera showing the point spread function of a nano-scatterer at different positions. Each image is 2µm × 2µm and normalized to the extremum of its own contrast. arise due to the substrate interface where a refractive index mismatch occurs. The most critical issue in modeling the iPSF is to accurately describe the aberrations introduced by the stratified layers of differing refractive indices which affect both the illuminating and scattered fields. This is not only important for a quantitative description of iSCAT images, but as we shall see, it also provides very valuable axial information.
We start with a qualitative discussion on the origin of the critical aberrations which affect the phase and the amplitude of the scattered field. Fig. 2(a) illustrates a green ray emanating from a point source such as a nanoparticle scatterer. In the absence of an interface, these rays would seamlessly propagate (dashed red line). On the imaging side, the collected unaberrated rays are focused into a tight diffraction-limited spot which would possess on axis a sharp amplitude profile. In the presence of interfaces, e.g., that of a coverslip, the green ray emitted from the source becomes distorted through refraction ( Fig. 2(a)). When tracing such rays on the image side, one finds an aberrated extended focus region because the degree of refraction at interfaces depends on the angle of the wave vectors.
The traveling phases of the illumination and scattering fields lead to a rapidly varying phase shift, which amounts to a round-trip phase of π for a displacement of λ m /4 where λ m is the reduced illumination wavelength in the medium (i.e., in the order of 100 nm). This phase shift has been previously exploited to achieve an exquisite axial sensitivity in iSCAT [9,12,20]. In addition, because the wavefront is composed of a distribution of spatially-confined wave vectors, the phase fronts undergo a cumulative π shift along the direction of propagation, a phenomenon known as the Gouy phase shift [21,22]. For the tight focus of the unaberrated beam, the Gouy phase shift occurs rapidly in proportion to the extent of the focal volume. For the aberrated rays, however, the Gouy phase accumulates much more slowly across the focal depth. In the next section, we present a vectorial diffraction model to describe the imaging of a scattering dipole with explicit inclusion of these aberrations for a wide-field reflection interferometric microscope.  Here the '0' labels the nominal design values. Note that the perpendicular polarization component of the light is maintained while the rays are refracted by the interfaces opposite to the parallel polarization vector that changes direction for every passing through an interface. In many experimental measurements, the system under study deviates from the intended design scenario, for example by having the nanoparticle situated off the coverslip (c).

A vectorial diffraction model for the iPSF
The vectorial description of diffraction has proven to be a powerful and accurate framework to analytically compute the light emitted by a dipole through stratified media [23,24]. This forms the conceptual basis for our model. We first begin by recasting the interference equation (Eq. (1)) to a vectorial format to accommodate an illumination field with an arbitrary polarization. We write E inc = E inc + E inc ⊥ , where and ⊥ denote two orthogonal polarization states in the plane normal to the propagation direction. The detected intensity is then given by where E denotes a vectorial field, and the argument γ is a vector containing all the physical parameters in the coherent imaging system, which we will shortly discuss and for brevity omit from writing explicitly. Experimentally, one often views the iPSF by firstly having subtracted the imaging background from I det and then normalizing again by this background. Hence, we write In the following, we consider each field contribution in the interference equation.
The reflected reference field: E ref The reference field arises from a reflection of the incident illumination at the interface between the coverslip and the sample and is, in principle, polarization sensitive. The strength and phase of this field are given by the reflectance R and complex amplitude reflection coefficientr obtained from Fresnel coefficients If multiple interfacial layers exist, one can readily cascade the Fresnel coefficients.
The scattered field: E sca The amount of light scattered by a nanoparticle is determined by the strength of the induced dipole polarization encapsulated by the polarizability [25], Here, the particle semiaxes are a 1 ≤ a 2 ≤ a 3 , L i denotes the depolarization factor along each axis, and 1 and m are the wavelength-dependant permittivities of the nanoparticle and the medium, respectively. For a spherical nanoparticle of radius a, L i =1/3, yielding the familiar expression for polarizability, and the material-dependent scattering phase The scattering cross-section for a spherical particle is then given by where k is the wave vector of the illuminating field. It follows that the scattering amplitude for each polarization component is described by where T is the transmittance of the interface. We have introduced η = 1/π arcsin[min(NA/n m , 1)] as a collection efficiency factor to take into account the angular extent to which the scattered light is collected by the objective. To ensure that the collection angle does not exceed the limit of π/2, the minimum function sets the maximum limit of the sin function to be 1. Here, we have introduced E sca andÊ sca as the complex amplitude and a unit vector for the scattered field E sca , respectively. The former is a complex number accounting for the absolute scattering amplitude and material-specific scattering phase of the nanoparticle. The latter term describes the three-dimensional spatial distribution of the vectorial scattered field.
We take as our foundation the vectorial diffraction framework for a radiating dipole emitting through numerous stratified layers [26]. This treatment accounts for the refractive aberrations introduced in Fig. 2(b) and (c). We begin with a complete description of an unaberrated or 'ideal' scenario, depicted in Fig. 2(b) and then apply a correction factor to the phase to account for the effect of aberrations, as formulated by Ref. [27]. This correction factor is in the form of an optical path difference (OPD). It should be emphasized that most modern high-NA objectives are designed to perform aberration-free imaging of a point source seated directly upon the coverslip at the focus of the lens. Hence, we describe this scenario as unaberrated or the nominal design scenario.
Following Ref. [26], the vectorial electric field of a radiating dipole situated in the sample medium at the position x p = (x p , y p , z p ) with a moment directed along the unit vectorê p = (sin θ p cos φ p , sin θ p sin φ p , cos θ p ) is given by: where θ p and φ p are the zenith and azimuthal angles of the dipole (the optical axis is taken as z-axis shown in Fig. 2), and E m is the magnitude of the field. The superscript m onê m ⊥ stands for the sample medium. The unit vectorsê m andê ⊥ are the p-and s-polarized components of E m in the sample medium. We note that since the polarization direction of the s-polarized component e ⊥ is preserved in the optical system (see Fig. 2(b)), we do not assign a superscript to it.
After propagating through the layers in the object space and the objective lens, the dipole field impinging on the aperture is given by ⊥ are the Fresnel transmission coefficients for p-polarized and s-polarized light from layer j to layer j + 1. This takes into account a glass slab for the coverslip and an immersion oil layer with thickness and refractive index pairs (z p , n m ), (t g , n g ) and (t i , n i ) respectively, but it may be generalized to include an arbitrary number of strata. In our case, superscript indices 1 and 2 label the interfaces between water-glass and glass-oil, respectively. The unit vectors of the polarization directions in spherical coordinates arê A ray in the immersion layer described by a wave vector k i has the zenith and azimuthal angles of θ i and φ, respectively, yielding The field impinging on the aperture E a (Eq. (10)) is collected up to the angular extent of NA, and is subsequently imaged on the detector using the Richards-Wolf integral [26,28,29] (expressed in coordinates of the object space), Here, A 0 is a scaling factor that incorporates E m . The OPD correction term (Λ) is given by [26] where the vector τ = (n i , n i0 , n g , n g0 , n m , t i0 , t g , t g0 ) describes the optical parameters of the setup with '0' denoting the nominal design parameters. The field E d is a function of the observation , while z focus denotes the focal position with respect to the origin of the coordinate system, which lies at the interface between the coverslip and the sample layer.
To apply this framework laid out for a radiating dipole field to that of a scattering dipole, we make the following modifications. Firstly, we point out that under realistic experimental conditions (i.e., non-design conditions), the incident light traverses an extra optical path given by [(n m z p + n g t g + n i t i ) − (n g0 t g0 + n i0 t i0 )] before impinging on the particle. These additional phase terms can be readily included within the OPD to account for the excitation phase. We then define Λ sca as the aberration term for a scattering dipole as Thus, to obtain the electric field of the scatterer at the detector plane, we consider the electric field E d for a radiating dipole (Eq. (13)), replace the OPD term Λ (Eq. (14)) with that of a scattering dipole source Λ sca (Eq. (15)), and implement an induced dipole moment along the polarization of the incident lightê inc in place of the dipole momentê p [23,26,30].

The extraordinary extended depth of focus in iSCAT and the uniqueness of the defocused iPSF
Leading on from the qualitative introduction to the role of aberrations introduced earlier in Section 2.2, here we present an explicit examination of the complex scattered light and the interesting features of it such as structure, asymmetry and uniqueness. These aberration-induced contributions richly tailor the iPSF and yield localization of a nanoparticle over an extended axial depth of focus, a feature that is inherent to iSCAT microscopy. We begin by presenting Fig. 3 where we consider the example of a gold nanoparticle (GNP) held above the coverslip and where the focal plane is swept throughout the sample space, depicted in Fig. 3(a). In Fig. 3(b-d) we present the meridional iPSF profile, the amplitude and differential phase components respectively for the scattered light for the experimental arrangement shown in Fig. 3(a). In Fig. 3(e) and (f) we show the amplitude and the differential phase of the unaberrated scattered field, i.e. the case where the GNP lies in a medium with the refractive index similar to that of the coverslip, hence in absence of any interface. A notable distinctive feature of the iPSF is the strong asymmetry in the distribution of both amplitude and phase about the position of the GNP. In our formulation, we consider the difference between the phase of the scattered light to that of the reflected light. Therefore, one directly sees the Gouy phase in the phase plots of Fig. 3. Examining the phase evolution reveals a slow variation when the focus is below the particle position, but this phase then rapidly accelerates close to the true position of the GNP where thereafter one observes a dramatic continuum of ripples in a honeycomb-like pattern that extends onward for many micrometers. One interesting observation in the amplitude distribution is that the maximum iSCAT signal in the image space Fig. 3. The effect of aberrations on the iPSF. (a) The geometry of imaging in which the particle is at a fixed height above the coverslip (z GNP = 3 µm) and the focal plane is swept starting from the coverslip (z focus = 0) up to 7 µm above the particle. (b) The meridional profile of the iPSFs as a function of the focal position z focus . With an ever-increasing defocus, additional side lobes begin to appear around the main peak which continue to spread out radially. These features keep the profiles of the iPSFs unique over a long axial range.
The blue dashed lines act as a guide to the eye for one such feature. (c,d) The aberrated amplitude and differential phase of the scattered light.
(e,f) The corresponding amplitude and differential phase of the scattered light in an ideal interface-free configuration.
occurs for a particle location shortly below the focus position. Moreover, as shown in Fig. 3(c), the aberrations also affect the absolute strength of the scattered field, which is found to be roughly 40% lower than its value in the absence of interfaces (see Fig. 3(e)). The observed asymmetry proves to be a great asset for removing the inherent ambiguity that is caused by the periodic traveling phase along the axial direction in interferometric signals [9,12,20]. We note that coverslip-induced aberrations stem from the sample geometry and hence are common not only to iSCAT but also to other microscopy modalities. However, interferometric microscopies uniquely benefit from these aberrations because one detects the scattered signal through its field amplitude rather than intensity, and hence as the former decays more slowly across space, there persists a larger detectable signal over a longer axial range. In the case of intensity-based microscopies with a limited depth of focus (typically of about 0.5 µm at high NA), one needs to employ additional phase-based engineering means to extend the axial range [31,32].
In particle tracking experiments, the focus of the microscope objective is often set at a fixed height within the sample depth while the particle can travel freely in all three dimensions. We now investigate the iPSF profiles for this configuration theoretically by considering different axial locations of a particle and focus positions, shown in Fig. 4(a). In brief, through defocusing, the corresponding iPSFs accumulate a slowly varying Gouy phase shift in the scattered field over many micrometers, leading to profiles typified by the red-boxed iPSFs in Fig. 4(a). In comparison, when the particle is free to change position the main phase parameter at play is now the accumulated traveling phase between the incident and scattered light with respect to the reference. This phase difference rapidly accumulates over tens of nanometers rendering rapidly changing iPSFs -depicted in blue in Fig. 4(a). In Fig. 4(b), we further highlight the remarkable sensitivity of the iPSF to the axial position of the nanoparticle. Owing to aforementioned aberrations, the radial profiles of the iPSFs show strong distinct structure with axial steps of as small as 10 nm. This feature was recently exploited to track a gold nanoparticle on the plasma membrane of a live cell [13].
To take a broader view of the parameter space of Fig. 4(a), in Fig. 4(c) we plot the peak c z focus (μm)  centroid contrast as a continuous function of focal position and particle height. In doing so, we can identify three regimes labeled Ω 1 , Ω 2 and Ω 3 that show distinct trends in their contrast dependence. Regime Ω 1 describes the case where the focus is below the particle and displays a weak contrast due to the flat variation in amplitude as well as the accumulation in the Gouy phase being exceptionally slow. Regime Ω 2 displays the strongest contrast as the focus lies close to, and approximately one micrometer above, the physical height of the particle as the aberrated wavefronts begin tightening to a focus which consequently begins to ramp up the accumulated Gouy phase. Finally, there is regime Ω 3 which is similar to the imaging scenario of regime Ω 1 inasmuch as that a large defocus is present, but where now the defocus lies in the space beyond the particle on the non-objective side. As shown in Fig. 3, this region is richly structured by numerous Gouy phase shift ripples and a slow decay in amplitude that extends outwards for many micrometers far beyond the axial location of the particle. To summarize, in interferometric detection scheme, one boosts the useful information encoded by the aberrations in the amplitude and phase of the scattered light by mixing it with a strong planar reference light. In this way, iSCAT yields two advantages. Firstly, the recorded information is doubly imprinted by both the phase and amplitude of the field which is not possible in the intensity-based measurements, and secondly, iSCAT signals win in sensitivity through the structural features of iPSFs.

Verification of iPSF model against experiment and simulation
We evaluated the calculated iPSF against those obtained from an experimental wide-field reflection iSCAT microscope as depicted in Fig. 1. This microscope uses an oil-immersion objective (100x, NA = 1.4) and a sample prepared on standard glass coverslip of 170 µm thickness. Gold nanoparticles were chosen as model sub-wavelength scattering particles. We used a laser at λ = 525 nm for illumination and a camera as detector. Numerical simulations replicating the aforementioned experimental arrangement were performed using finite difference time-domain (FDTD) calculations with the commercially available package 3D Electromagnetic Simulator (Lumerical Inc.).

A particle moving in the axial direction
An important application of iSCAT microscopy is in three-dimensional tracking of nanoscopic objects, e.g. in fluidic channels [9] or biological cells [13]. Thus, it is important to acquire a quantitative mastery of the iPSF as a function of the axial position of the particle in the sample. As an example of a realistic laboratory application, we tracked a single GNP (diameter 40 nm) diffusing on a model lipid membrane of a giant unilamellar vesicle (GUV), which is often considered as an artificial cell system [33]. Here, we investigated a GUV with a diameter of 10 µm upon which a GNP was bound and free to diffuse on its spherical surface. As shown in Fig. 5(a), the GUV was held by a micropipette above the coverslip. Further details regarding the procedure to perform single particle tracking on a GUV can be found in Ref. [5,34].
The 3D tracking of the GNP is achieved by fitting the model to the experimental iPSF extracted from a raw iSCAT video. Here, the wide-field image of the GUV delivers direct information on its spherical shape and size, which is used as an independent means to verify the height of the GNP above the coverslip [5]. The nonlinear iterative fitting of the model can be initiated with a coarse initial guess of the lateral position of the particle. This guess can be obtained, for example, using a localization algorithm based on radial symmetry [35]. Similarly, for the initial coarse guess of the axial position of the particle one can either use information from the sample geometry, as is possible here with a GUV, or alternatively one can roughly match the recorded iPSF with that of a model-generated template, such as the series presented in Fig. 4.  sufficiently accurate to describe the sensitive and delicate iPSF formation encountered under realistic experimental conditions.
Lateral and axial localization of the GNP in each video frame yields 3D trajectories spanning a range of many micrometers in each of the dimensions, as displayed at the bottom of Fig. 5(a). Fig. 5(c) presents a histogram of the axial and lateral localization errors, which amount to only a few nanometers. The bimodal distribution in the lateral localization error (blue, Fig. 5(c)) result from the radial asymmetry in the experimental iPSFs (see Fig. 5(b)). This asymmetry disproportionately hampers the fitting of the model function to the iPSFs which rely upon the secondary lobe information, that is, iPSFs whose central peak is naturally weak. We note, however, that the asymmetry in the experimental iPSFs may originate from other sources of aberration, and if desired could be described within our analytical model description via inclusion of the appropriate Zernike coefficients. Fig. 5(d) plots the maximum of the central lobe in the iPSF as a function of the axial position of the GNP. As predicted by our model discussed in Fig. 3, one finds a shift between the axial position of the objective focal plane and the height at which the maximum peak contrast is observed. Furthermore, the focus region has an extended character and the Gouy phase is accumulated in a smooth fashion when the imaging plane is located below the particle. However, the Gouy phase has discrete distortions when the image plane lies above the particle and this in turn leads to the observed irregularities in the oscillation of contrast profile. These data iterate the asymmetry about the imaging plane caused by aberrations and the extended depth of field in iSCAT microscopy. For a visual comparison, we mark the shallower depth of focus achieved with intensity-based imaging from high NA objectives with the gray band in Fig. 5(d).
We note that in determining the ground truth for height assignment in the experimental trajectory, we assume the GUV to be a perfect sphere. The GUV, however, does not possess a uniform radius of curvature. This leads to slight disagreement with the height assignment found via the central contrast of the modeled iPSFs for GNP positions in larger defocusing regions.
To avoid the complications of realistic systems, we also devised a textbook scenario, where a nanoparticle could be positioned only in the axial direction under full control. Here, we placed a GNP (diameter 50 nm) upon the 50 µm tip of a rounded quartz tip fabricated by melting a quartz rod (diameter 1 mm) using an oxyacetylene flame to produce a rounded end (see Fig. 6(a)). A dilute solution of GNPs was drop cast upon the plasma-cleaned quartz tip which resulted in several GNPs being deposited across the spherical tip surface at a low density. The quartz tip was independently positionable with respect to the coverslip by a piezoelectric actuator while both the coverslip and objective focus remained static throughout (see Fig. 6(a)). To suppress unwanted scattering from the quartz tip, the sample medium was index matched to quartz (n med = 1.461) by use of a 91.5% (w/v) glycerol solution. We remark that a faint background signal from incomplete suppression of the tip scattering remains, but it does not affect extraction of the iPSF. A thin layer (25 nm) of TiO 2 was deposited upon the coverslip via atomic deposition to compensate the reduction in reflectivity from the coverslip interface. We note that the added TiO 2 layer is thin enough to not perturb the wavefronts. The focus of the microscope objective here was positioned at a height 1.5 µm above the coverslip. A stack of measured, modeled and simulated iPSFs is shown in Fig. 6(b) for GNPs placed at various axial positions. To examine the agreement among these three case studies, we plot the central contrast of each iPSF over the full axial range in Fig. 6(c). Again, the degree to which all three sets agree with one another reflects the accuracy of the model despite the high sensitivity of the aberrations to the details of the experimental arrangement such as thicknesses and indices of the individual layers. We emphasize that even for the iPSFs with a weak peak contrast, which may often be regarded as giving no detectable signal, the side lobes are highly structured and provide a readily usable alternative. We also draw attention to the comparison of the axial contrast trend here to that of Fig. 3, with a greater symmetry in the distribution of contrast above and below focal plane (z focus = 1.5 µm). This is due to the higher refractive index of glycerol, which is more closely matched to that of the glass coverslip than that of the water, hence weakening the effect of aberrations. This effect is also evident in the better match between the objective focal plane and the height at which the maximum iPSF contrast takes place. We point out, however, that an extended focus is still maintained.
It is also meaningful to compare the model against the simulation. While both agree compellingly with the experimental results, we note numerous advantages to a model framework over simulation. The model offers a physical picture of the formation of the iPSFs, including all the layers of varying refractive indices in a real experimental configuration and regardless the thickness of these interfacial layers, it samples the wave vectors at the exit pupil and analytically calculates the associated aberrations. The computational time for the model is negligible (milliseconds) whereas a full electromagnetic simulation requires times on the order of hours to days. Simulations are particularly helpful if the geometry of the system is complex [36]. However, in the simulation of the electromagnetic fields over a range of several micrometers, one must sample the field close to the object and then numerically propagate it further owing to the expensive computational resources and time this would otherwise require. Indeed, simulations might be susceptible to sampling artifacts and convergence problems for challenging sample geometries. In our case, the limited accuracy of our simulations can be seen in a somewhat poorer agreement with the experiment than when compared to the results of the model (see Fig. 6(c)).

Focusing through a nanoparticle held at a fixed height
Another situation of practical interest concerns the case where the focus of the microscope objective is tuned through the sample. To study this situation in a model system, we examined a GNP placed at the water-glass interface as well as a GNP embedded within gelatin at 4 µm above the coverslip. Furthermore, to investigate the effect of the material-dependent phase shift upon scattering, we also compared these results with those of a dielectric nanoparticle at the water-glass interface (see Fig. 7).
Panel (a) of Fig. 7 considers the elemental case of a dielectric particle upon the coverslip, e.g. as in applications of protein detection [8]. Here, it is important to understand the iPSF in order to adjust the focal plane during measurements or to maximize the detected contrast. We find that when the focal plane coincides with the position of the dielectric particle sitting on the coverslip, the contrast is maximally destructive, testifying to the Gouy phase shift, wherein the focused point source of light is out of phase with a plane wave of initially the same starting phase by π/2. Here, it should be borne in mind that since the imaginary part of the dielectric function for the particle is negligible, the phase of the scattered field is the same as that of the illumination. We also observe, accordingly, that the observed contrast has a strong symmetry about the position of the particle.
Panel (b) considers a GNP, which introduces an additional material and wavelength-dependent phase shift, φ sca on the incident light. This shift renders a GNP with a more weakly negative contrast in the focus as compared to the dielectric particle in (a), requiring defocusing to maximize the contrast. In panel (c) we consider the case of a GNP displaced off the coverslip, thus invoking the need for the introduction of additional phase-shifting aberrations. Indeed, as observed in Fig. 3, we find that the envelope of the amplitude is asymmetric, and the contrast experiences many cycles over the axial range, while decaying much more slowly with progressive defocusing above the position of the particle.
Upon closer examination, the axial trends of the iPSFs presented in panels (a) and (b) of Fig. 7 also reveal some asymmetry about the focus position. These are due to slight deviations of the coverslip and immersion oil properties (thickness, refractive index) with respect to the design parameters (see details in the caption of Fig. 7). The asymmetry observed in panel (c), which  ), but for a static GNP situated off the coverslip. The height of the particle was fitted to be z GNP =4 µm, in water. The medium here was a dilute mixture of gelatin in water with n m = 1.33. Each image is 2 µm × 2 µm and normalized to the extremum of its own contrast.
is due to the aberrations caused by defocusing through the water-glass interface, however, is much more pronounced. We wish to emphasize that in a truly aberration-free microscope, the Gouy phase shift and the amplitude distribution are symmetric about the focus and hence lead to ambiguity in height assignment [37,38].

Feature-based height-assignment of iPSFs through machine learning
Conventional methods such as iterative algorithms widely used in interferometric imaging techniques, e.g. in digital holography, achieve accurate results at high processing speeds for strong signals from large particles (typically micron-sized). However, they tend to face difficulties tackling noisy distorted signals and have, thus, employed machine learning tools for various image analysis tasks [39,40]. Considering that efficient and accurate 3D tracking of nano-matter has enormous potential in modern biology, biophysics and soft condensed-matter physics, we now discuss methods to characterize measured iPSFs from small nanoparticles. The iPSF vectorial model developed in the previous section can now be used as a fit function for experimental measurements using a nonlinear optimizer. This would allow one to deduce the position of the nanoparticle in all three spatial dimensions. However in the presence of a background and a high level of noise in the iPSF signal, axial tracking becomes challenging. Machine learning techniques perform reliably in such conditions. Moreover, as opposed to nonlinear fitting algorithms, such learning-based tools do not require initial guesses of the axial position of the particle. We recall that the challenge in tracking particles axially in earlier iSCAT studies was the short accessible axial range of about 100 nm limited by the periodicity of the contrast in the central lobe of the iPSF: when attempting to assign a height to the centroid contrasts (see the region outlined with the dashed line in Fig. 8(a)), one finds multiple solutions for the axial position of the GNP, z GNP . In Fig. 4(b), however, we showed that the ring features of the iPSF provide an exquisite sensitivity to axial displacements. We now show that this sensitivity can be exploited in an unsupervised machine learning scheme to resolve the motion of a GNP. In Fig. 8(a), we present an example of a GNP that travels axially a known distance of 300 nm along the optical axis. In Fig. 8(b), the iPSFs associated with the three distinct color-coded regions of Fig. 8(a) illustrate how their unique ring structure resolves this ambiguity. We now investigate the application of unsupervised machine learning tools for nanoparticle tracking. The proof of principle work shown here can be extended to even longer axial ranges using more complex deep learning algorithms. We therefore quantify the natural branching in the iPSFs to facilitate benchmarking for the development of follow-up learning algorithms.  The first step of the procedure is to group similar looking iPSFs. This can be done by defining a so-called distance function as a quantitative measure to compare different iPSFs. This can be done by computing the pixel-wise difference between two iPSFs under study and then summing up the squares of the outcome for all pixels. We then employ a K-means-based unsupervised clustering algorithm [41] to sort all iPSFs within the region of interest (the region outlined with the dashed line in Fig. 8(a)) by mutual similarity into an optimal number of groups -referred to as 'clusters'. In Fig. 8(a), we color-code three such clusters. Here we used the default K-means++ algorithm [42] from the Statistics and Machine Learning Toolbox in MATLAB. For each cluster, the iPSFs are averaged to obtain their 'cluster center'. The group of iPSFs are then fitted with the model and averaged within each cluster to obtain their mean fit iPSF. In Fig. 8(b), we plot the three cluster centers next to the mean fitted iPSFs for comparison. The clustered iPSF representatives show excellent agreement with the outcome of the complete physical model. We note that because the iPSFs are radially symmetric, we can extract the mean radial profile of the iPSFs over N × N-pixels to reduce the complexity of this computation from N 2 down to N/2 and, thus, significantly boost the computational performance with a minimum memory footprint. After a successful clustering we can go from 1D profiles back to the full 2D iPSF shapes while using the labels and groupings that we obtained for the 1D case. Interestingly, the results presented in Fig. 8(b) show that we can nevertheless go beyond the structural rings of the iPSF to even account for slight imperfections in the imaging setup.
Next, we benchmark the clustering accuracy through a 'silhouette value' as a measure for the consistency within clusters [43]. The silhouette value of one data point after being clustered is defined as the difference between its smallest mean inter-cluster dissimilarity and its mean intra-cluster dissimilarity. The former quantity is a measure of how dissimilar (i.e., far) this point is compared to the other points in other clusters than its own cluster. The latter quantity denotes the average dissimilarity of this point compared to all other data points within the same cluster. The silhouette coefficient is normalized to ± 1. When an iPSF is strongly matched to its own assigned cluster, but poorly matched to other clusters, the silhouette coefficient tends to +1. If a silhouette value is larger than 0.5, one considers the data to tend to cluster. In Fig. 8(c), we present the silhouette coefficients obtained when sorting the experimental data presented in panel (a) with an increasing number of cluster groups. The curve starts at a value about 0.7 and then shows a sudden drop when more than three cluster groups are considered. This suggests the optimal number of iPSF clusters to be three. In this case, the SNR of the experimental iPSFs were about 2.
To examine the robustness of this approach to the SNR, we synthesized iPSFs for the same axial range and added different noise levels. We performed three rounds of modeling with SNR levels of 1,2 and 4 and followed it with the silhouette analysis shown in Fig. 8(d). The dashed-dotted curve shows the SNR case of 1. In this case, the silhouette value starts with a value below 0.5 and falls almost linearly suggesting that there is no meaningful grouping in the iPSFs as all the discriminatory features between them are buried in noise. The solid line presents the case for a SNR of 2. The curve starts with a silhouette value of about 0.6 and descends for more than three clusters, similar to what we observed for the experimental data in panel (c). The case of a higher SNR of 4, yields the silhouette curve shown with the dashed line. The curve suggests that the discriminatory features of iPSFs at this SNR are strong enough to show a meaningful grouping even for more than three clusters. We note that for higher SNR levels more iPSF side lobes lie above the noise level and provide additional discriminatory sources of information.
The iPSFs extracted through clustering can be used together with a simple cosine function to build a template for the iPSFs in a given axial range. The iPSF centroid value of the green cluster is solely used as the maximum value of the cosine term cos(2k m z GNP ) to obtain the dashed gray line shown in Fig. 8(e). The blue and red iPSFs sit opposite to each other on the guiding curve due to their dissimilarity in the clustering process. Using the 2D shapes of the clustered iPSFs and their relative axial distance, one can interpolate and extrapolate iPSFs for every single nanometer of this axial range. The centroid contrasts of the estimated iPSFs plotted with the solid black line overlaps strongly with the guiding dashed gray line. The 2D estimated iPSFs can be fitted independently using the iPSF model to obtain the axial position of the particle. In Fig. 8(f), we plot the results from the unsupervised machine learning method against the height of the particle using the model, revealing an excellent agreement.

Discussion and outlook
Optical microscopy has played a very important role in the development of many disciplines of nanoscience and nanotechnology over the past three decades and continues to be a dynamic field of research [44]. One of the most valuable modern methods in optical microscopy is the realization that one can localize single molecules and nanoparticles with much better precision than the diffraction limit, thus, giving access to the investigation of nanoscopic trajectories. While most of the efforts in single particle tracking have been limited to lateral movements on surfaces, there is an acute need for three-dimensional trajectories, as methods become powerful enough to study more complex geometries such as live biological cells and tissues. The intrinsically confined nature of the axial focus profile of the fluorescence intensity and the low SNR in single-molecule fluorescence microscopy render such endeavors very challenging. One of the approaches to combat this problem is based on PSF engineering, where one sculpts the PSF to provide axial information over a range of 2 µm, reaching a localization precision of about 10-20 nm [31,45].
By examining a rigorous model, numerical simulations and comparison with controlled measurements, we have shown how the use of iSCAT microscopy not only provides a high SNR and temporal resolution [2], but it also provides a uniquely powerful avenue for recording three-dimensional trajectories with nanometer spatial precision. Here, the aberrated scattered field manifests as an extended and asymmetric field distribution along the optical axis. This distorted confinement also extends the distance over which the Gouy phase evolves. Together, these effects produce an iPSF consisting of an ample amount of light, allowing an unprecedented 3D tracking of a nanoparticle over the range of up to 10 µm with nanometric precision. A similar result was also recently reported by formulating the effect of the lens in imaging an index-matched sample using an in-line digital holographic microscope [19]. These by far outperforms conventional dark-field and PSF-engineered fluorescence microscopies.
The iPSF model presented here also opens new strategies for removal of the unwanted imaging background which accompanies the signal of interest. Previous works have exploited features unique to the PSF to separate the weak signal of nanoparticles from the speckle background of glass coverslips [46][47][48] and for the dynamic backgrounds such as those encountered in living biological systems such as cells [13]. Deconvolution with a deterministic iPSF or spectroscopic discrimination of the scattered field through a wavelength-dependent polarizability [49] will further improve the separation of probe and background. Furthermore, while the axial extent of the iPSF is already large, we envisage improvements based on PSF engineering [50], permitting deep tissue imaging with interferometric sensitivity and precision. We speculate that additional aberrations originating from a highly scattering media, such as biological tissue, might also encode a wealth of information about the scatterer and its environment which could be retrieved over an extended axial range necessary for deep tissue imaging.
The model framework presented here can be extended to describe the position and orientation of non-symmetric particles such as nano-ellipsoids. It can also model the iPSF for a nano-scatterer illuminated by non-planar wavefronts and is, therefore, readily applicable to other common microscopy modalities such as confocal and bright-field transmission schemes. Furthermore, phase plates used to boost the detection contrast of particularly weak scatterers can be treated [51,52].
Another promising avenue for advanced image analysis aided by the iPSF model is the extraction of important physical quantities through machine learning. To benefit from fully automatic and robust deep learning techniques, one needs reliable physical models to generate ground truth iPSFs to train such networks. Our work provides a model which can generate data for use in training of supervised machine learning methods or in validating the outcomes of unsupervised machine learning.
In conclusion, we have used experimental observation and numerical simulation to confirm that the vectorial diffraction model developed in this article accurately describes the operation of the popular wide-field reflection iSCAT microscope. We have shown that the interferometric point spread function, assisted by geometry-induced aberrations in the scattering wavefront, contains a wealth of information about the scattering source, allowing nanometrically accurate localization over an extended axial range in all three dimensions. We emphasize, however, that the radiation pattern of a scatterer is modified when situated at a sub-wavelength distance from the sample-substrate interface [53]. This effect, which also modifies the iPSF profile, has been neglected in our current work, but it will be the topic of a sequel publication.