## Abstract

Fluorescence phase-shifting interferometry (FPSI) is an optical technique that coherently combines the phase-shifted 4π steradian emission wavefronts of a single fluorescent emitter to obtain multiple interferograms from which the emitter axial displacement can be retrieved with high precision. Here, we study the axial displacement sensitivity in 4-step FPSI within the framework of maximum-likelihood (ML) phase estimation. Using Monte-Carlo simulations, we show that regardless of the method used to preprocess the measured interferograms, the variance of the ML estimate of the axial displacement approaches the Cramér-Rao lower bound and is closely limited from above by the variance of the classical 4-step phase shifting estimator. The difference between these lower and upper bounds depends on the interferogram visibility and signal-to-noise-ratio (SNR), with a percentage change of up to 29% that yields an absolute change of a few sub-nanometers to some nanometers for SNRs larger than ~5. Our results suggest that for these levels of SNR, the use of the computationally simpler classical 4-step phase shifting estimation can be adequate to accurately determine axial displacements in FPSI, as we also experimentally verified. FPSI interferograms with lower SNR but good visibility can benefit from the use of the ML estimator, provided that the spatiotemporal phase stability of the FPSI system is high.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Precise localization of fluorescent emitters in three dimensions is an important task in the biophysical sciences, and which has enabled key applications, such as super-resolution imaging and tracking [1–6]. While localizing single emitters with high precision in the lateral dimension is an established process [7–9], precisely assessing the position of the emitters along the axial direction requires dedicated measurement and analysis methods. Amongst these methods are widefield imaging techniques including, for example, multiple z-slice de-convolution imaging [10], evanescent-field intensity decay analysis [11], defocused and multi-focal imaging [3, 12–14], axially asymmetric point-spread-function engineering [1, 4, 5], and fluorescence interference contrast microscopy [15]. All these techniques can localize multiple, optically resolved, fluorescent particles to within <~100 nm and over a depth-of-field of ~0.1–20 μm, depending on the photon yield and excitation beam profile used. Another class of three-dimensional (3D) particle localization approaches are focus locking and volumetrically/circularly laser scanning that lock the position of a single particle to follow its motion with ~5–150–nm axial localization precision over a large axial distance of ~5−70 μm, depending on the working distance of the objective lens and the axial travel range of the microscope stage used [2, 16–20].

Recently, a new widefield imaging technique that can precisely localize multiple, optically resolved, fluorescent emitters along the optic axis using fluorescence phase-shifting interferometry (FPSI) has been theoretically investigated and experimentally demonstrated with ~1–10–nm axial localization precision in thin (<~250 nm) and thick (~1 μm) samples [21–25]. In FPSI, multiple (at least three) phase shifts *δ _{n}*,

*n*= 1,2,…,

*N*are first introduced between the 4π steradian emission wavefronts of a fluorescent emitter positioned at the point (

*x*

_{0},

*y*

_{0},

*z*

_{0}) in space, as illustrated in Fig. 1(a). This can be accomplished using, for example, a phase–shifting mirror in a triangular two–opposing objectives interferometer [see inset of Fig. 1(a)]. The wavefronts are then combined in the interferometer and subsequently imaged on a two–dimensional sensor array to produce phase-shifted interferograms, as depicted in Fig. 1(b) for

*N*= 4.

The lateral position of the emitter (*x*_{0}, *y*_{0}) can be assessed with high precision from the dc intensity of the interferograms using established processing algorithms, such as least-squares (LS) and maximum-likelihood (ML) fits [7, 8]. The axial displacement of the emitter *z*_{0} can be evaluated from the interferograms by measuring the emitter wavefront phase difference Δ*ϕ* (modulo 2π) with classical *N*–step algorithms followed by linearly mapping of Δ*ϕ* to *z*_{0} [21]. The use of the classical algorithm in FPSI enabled to retrieve axial displacement of the emitters in *<~λ*/2– and ~2*λ*–thick samples, depending whether procedures that remove the ambiguity in the axial displacement calculation (due to phase wrapping) were used [22–25]. One such procedure applies inner and outer spatial weighting masks on the measured phase-shifted interferograms to retrieve Δ*ϕ* from both the inner and outer sections of the interferograms, enabling the axial localization of an emitter over an extended depth–of–field (DOF) of ~2*λ* [23, 24]. When retrieving Δ*ϕ* only from the inner section of the interferograms, a limited DOF of *~λ*/2 can be achieved. This preprocessing procedure was used on the phase-shifted interferograms simulated or measured throughout the work reported here.

The classical *N*-step algorithm is known to yield an unbiased and efficient ML estimate of the axial displacement under a Gaussian noise model with the same noise variance across all pixels and in all interferograms [26]. However, the noise model in FPSI is approximately Poisson with a parameter determined by the signal strength and the readout/background noise in each pixel [22, 24]. Consequently, the classical *N*-step algorithm is suboptimal when estimating the axial displacement from FPSI measurements.

The main aim of this work is to put on a firm ground the adoption of classical PSI approaches in FPSI nanoscopy by performing a fundamental and comprehensive study of the estimation of the axial displacement in FPSI using estimation theory tools. In particular, we investigate the ML estimator for the axial displacement in 4-step FPSI. Using high numerical aperture (NA) vectorial imaging theory and Monte-Carlo simulations, we compare the performances of the ML estimator and the classical 4-step algorithm (in terms of their standard deviations) over *~λ*/2 and ~2*λ* depth-of-fields, and their dependence on the number of photons detected, readout/background noise level, interferogram visibility, and the spatial weighting masks used to preprocess the interferograms. We show that regardless of the preprocessing schemes employed the relative difference in the standard deviations of the ML and classical estimators can be as high as ~29%, resulting in maximum absolute differences of a few sub-nanometers to some nanometers for practical interferogram signal-to-noise ratios (>~5). Our results indicate that, in practice, the use of the computationally simpler classical 4-step phase shifting estimation is adequate to accurately determine axial displacements in FPSI, as also experimentally verified. FPSI interferograms with lower SNR but good visibility (>~0.7) can benefit from the use of the ML estimator, provided that the spatiotemporal phase stability of the FPSI instrument is high.

This paper is organized as follows: in section 2, we outline the theoretical framework of ML and classical 4-step estimations in FPSI; in section 3, we present the simulation and experimental results and discussion, and conclusions are drawn in section 4.

## 2. Maximumlikelihood (ML) and classical estimation for the axial displacement of a fluorescent emitter in 4-step fluorescence phase-shifting interferometry (FPSI)

We first present the statistical model for the noise in FPSI interferograms, providing the basis for the mathematical derivation of the ML estimator and Cramér-Rao lower bound of the axial displacement of a fluorescent emitter in FPSI. The noise model includes contributions of both shot and readout/background noise. Next, we detail the derivation of the log-likelihood equations given four phase-shifted measured interferograms for binary or Gaussian weighting preprocessing schemes of the interferograms. These equations were then maximized numerically to obtain the ML estimates. We also describe the derivation of the Cramér-Rao lower bound, which expresses a lower bound on the variance of estimators of the emitter's axial displacement. Finally, we outline the classical 4-step phase shifting algorithm and its variance for the estimation of the axial displacement of an emitter. These estimation theory tools allowed us to comprehensively study and compare the performance of the ML and classical estimators in 4-step FPSI, as presented in section 3.

#### 2.1 Noise model in FPSI

In FPSI, two additive statistically independent noise sources—the Poisson shot noise and the Gaussian readout noise—corrupt the intensity of individual pixels ** ρ** = (

*x*,

*y*) = (

*ρ*cos

*φ*,

*ρ*sin

*φ*) in the

*n−*th phase-shifted interferogram of an optically resolved fluorescent emitter positioned at (

*ρ*_{0},

*z*

_{0}) and imaged onto plane

*z*,

*I*(

_{n}**,**

*ρ**z*|

*ρ*_{0},

*z*

_{0}). Mathematically,

*I*(

_{n}**,z|**

*ρ*

*ρ*_{0},

*z*

_{0}) is given by

*λ*) is a Poisson distribution with parameter

*λ*, and N(

*μ*,

*σ*

^{2}) is a Gaussian distribution with mean

*μ*and variance

*σ*

^{2}. Also,

*Ī*represents the spatial distribution of expected photons in

_{n}^{shot}*I*

_{n}when

*n*(

_{readout}**) = 0, and is expressed as**

*ρ**I'*and

*I”*being the dc intensity and fringe modulation of the measured interferogram, respectively, Δ

*ϕ*the imaged wavefront phase difference of the fluorescent emitter (from which

*z*

_{0}is retrieved), and

*δ*the

_{n}*n−*th phase shift between the 4π steradian emitter emission wavefronts, which equals to

*n*π/2,

*n*= 0,1,2,3 in 4-step FPSI. Also,

*ρ*_{0}

*= (*

^{s}*Mx*

_{0},

*My*

_{0}),

*z*

_{0}

*=*

^{s}*M*

^{2}

*z*

_{0}with

*M*being the optical magnification. Note that

*I'*can include a background noise term. Approximating the Gaussian readout noise distribution by a Poisson distribution [27–29],

*I*– (μ

_{n}*− σ*

_{ro}^{2}

*) is distributed as*

_{ro}#### 2.2 ML estimator for the axial displacement of a fluorescent emitter in 4-step FPSI

Neglecting edge and pixelation effects of the image sensor arrays [22, 24], the measurement data set for the ML estimation of *z*_{0} in FPSI constitutes, in general, *N* weighted phase-shifted interferogram intensities *I _{n}^{weighted}*(

*z*|

*z*

_{0}) given by [22–24]

*w*(

**) represents a spatial weighting mask.**

*ρ*For a binary spatial mask—explicitly, *w*(** ρ**) = 1,

**∈**

*ρ**A*and

*w*(

**) = 0,**

*ρ***∉**

*ρ**A*where

*A*is a region around an optically resolved fluorescent emitter in the image plane—and as the intensity of individual pixels in a phase-shifted interferogram is statistically independent, the use of Eq. (3) allows to write the probability mass function of

*I*(

_{n}^{weighted}*z*|

*z*

_{0}) – (μ̅

*– σ*

_{ro}^{̅2}

*) as*

_{ro}*δ*=

_{n}*n*π/2,

*n*= 0,1,2,3 and known μ̅

*, σ*

_{ro}^{̅2}

*given the parameter vector*

_{ro}**= [**

*θ**Ī'*,

*Ī”*, $\overline{\Delta \varphi}$] is the product of each probability mass function of

*I*(

_{n}^{weighted}*z*|

*z*

_{0}) – (μ̅

*– σ*

_{ro}^{̅2}

*), and reads as*

_{ro}

*I**(*

^{weighted}*z*|

*z*

_{0}) = [

*I*(

_{0}^{weighted}*z*|

*z*

_{0}),

*I*(

_{1}^{weighted}*z*|

*z*

_{0}),

*I*(

_{2}^{weighted}*z*|

*z*

_{0}),

*I*(

_{3}^{weighted}*z*|

*z*

_{0})] and

**= [**

*k**k*

_{0},

*k*

_{1},

*k*

_{2},

*k*

_{3}]. Given the measured interferograms

*I**and for known μ̅*

^{weighted}*, σ*

_{ro}^{̅2}

*, the ML estimator of*

_{ro}*z*

_{0},

*ẑ*

_{0}

^{ML}, is the value of

*z*

_{0}that maximizes the logarithm of the likelihood function given by Eq. (8), and can be computed numerically by solving the following log-likelihood optimization equation set for

*θ**ẑ*

_{0}

^{ML}[21].

In the case of an arbitrary mask *w*(** ρ**), the random variable

*I*(

_{n}^{weighted}*z*|

*z*

_{0}) – μ̅

*represents a weighted sum of independent random variables, and, by virtue of the central limit theorem, its probability density function can be well approximated by a Gaussian distribution of mean*

_{ro}*δ*=

_{n}*n*π/2,

*n*= 0,1,2,3 and known μ̅

*, σ*

_{ro}^{̅2}

*, σ*

_{ro}^{̅2}

*given the parameter vector*

_{In}^{shot}**= [**

*θ**Ī'*,

*Ī”*, $\overline{\Delta \varphi}$] is then expressed as the product of each Gaussian probability density function of

*I*(

_{n}^{weighted}*z*|

*z*

_{0}) – μ̅

_{ro}

*I**(*

^{weighted}*z*|

*z*

_{0}) = [

*I*(

_{0}^{weighted}*z*|

*z*

_{0}),

*I*(

_{1}^{weighted}*z*|

*z*

_{0}),

*I*(

_{2}^{weighted}*z*|

*z*

_{0}),

*I*(

_{3}^{weighted}*z*|

*z*

_{0})] and

**= [**

*i**i*

_{0},

*i*

_{1},

*i*

_{2},

*i*

_{3}]. Given the measured interferograms

*I**and for known μ̅*

^{weighted}*, σ*

_{ro}^{̅2}

*, σ*

_{ro}^{̅2}

*, we obtain*

_{In}^{shot}*ẑ*

_{0}

^{ML}by first solving numerically the log-likelihood optimization equation set for

*θ**ẑ*

_{0}

^{ML}[21]. We note that the log likelihood optimization equations [Eq. (9) and (14)] were obtained by differentiating the logarithm of Eq. (8) or Eq. (13) with respect to the parameter vector

**and setting the result equal to zero. A Nelder-Mead simplex algorithm was then used to solve iteratively Eqs. (9) and (14) for**

*θ***in Matlab.**

*θ*#### 2.3 Cramér-Rao lower bound for the axial displacement of a fluorescent emitter in 4-step FPSI

For a binary spatial mask, the maximum precision σ^{2}_{ẑ0}^{ML} with which the axial displacement of a fluorescent emitter *z*_{0} can be estimated is computed by mapping the Cramér-Rao lower bound for ${\widehat{\overline{\Delta \varphi}}}^{\text{ML}}$

*ẑ*

_{0}

^{ML}. We note that Eq. (15) assumes that the quality of the ML estimate of $\overline{\Delta \varphi}$ is essentially not degraded when the parameters

*Ī'*,

*Ī”*are unknown. By substituting the probability mass function of Eq. (8) in Eq. (15) and performing the ensemble average operator, σ

^{2}

_{Δ}

_{ϕ}^{ML}(and σ

^{2}

_{ẑ0}

^{ML}) can be obtained as

*Ī”*, $\overline{\Delta \varphi}$,

*Ī*are given by Eq. (6), and σ

_{n}^{shot}^{̅2}

*by Eq. (7), respectively.*

_{ro}In a similar way, σ^{2}_{ẑ0}^{ML} can also be derived for an arbitrary spatial weighting mask by substituting the probability density function of Eq. (13) in Eq. (15) and applying the ensemble average operator, resulting in Eq. (16) with *Ī _{n}^{shot}* and σ

^{̅2}

*replaced by σ*

_{ro}^{̅2}

*and σ*

_{In}^{shot}^{̅2}

*of Eq. (11), and*

_{ro}*Ī”*, $\overline{\Delta \varphi}$ by those defined in Eq. (10).

#### 2.4 Classical N-step phase shifting algorithm and its variance for the estimation of the axial displacement of a fluorescent emitter in 4-step FPSI

The classical 4-step phase shifting algorithm for evaluating $\overline{\Delta \varphi}$ for an arbitrary spatial weighting mask can be expressed as [24]

*z*

_{0}, we map ${\widehat{\overline{\Delta \varphi}}}^{\text{C}}$to

*ẑ*

_{0}

^{C}with the same linear transformation used to map ${\widehat{\overline{\Delta \varphi}}}^{\text{ML}}$ to

*ẑ*

_{0}

^{ML}in section 2.2. We note that the variance of ${\widehat{\overline{\Delta \varphi}}}^{\text{C}}$ can be approximated by linearization of Eq. (17) as

*Ī*(which equals to

_{n}^{weighted}*Ī*) and

_{n}^{shot}*Ī”*are defined in Eq. (10), and ${\sigma}_{{I}_{n}^{weighted}}^{2}$ equals to

*Ī*+ σ

_{n}^{shot}^{̅2}

*for a binary spatial mask [Eqs. (6) and (7)] and to σ*

_{ro}^{̅2}

*+ σ*

_{In}^{shot}^{̅2}

*for an arbitrary mask [Eq. (11)].*

_{ro}## 3. Results and discussion

In this section, we present two scenarios for the ML and classical estimations of the axial displacement of a fluorescent emitter *z*_{0} measured by 4-step FPSI. The first scenario involves the computation of *z*_{0}−estimations over *~λ*/2 depth-of-field, whereas the second scenario employs phase unwrapping to calculate *z*_{0}−estimations over an extended depth-of-field of ~2*λ*. We detail the simulation comparisons of the performance of these *z*_{0}−estimations as a function of the number of photons detected, readout/background noise level, interferogram visibility, and the spatial weighting masks used to preprocess the interferograms. The simulation data for these comparisons comprised noisy 4-step FPSI measurements of a stationary fluorescent emitter generated by high−NA vectorial imaging and Monte-Carlo simulations as described in our previous work [24]. We also present experimental data acquired by a temporal fluorescence phase-shifting interferometry setup. Our findings suggest that for practical SNRs (>~5) the use of the computationally simpler classical 4-step phase shifting estimation is adequate to accurately determine axial displacements in FPSI. FPSI interferograms with lower SNR but good visibility (>~0.7) can benefit from the use of the ML estimator, provided that the spatiotemporal phase stability of the FPSI instrument is high.

#### 3.1 ML and classical z_{0}−estimations from 4-step FPSI simulations over a depth-of-field (DOF) of ~λ/2

_{0}

Figure 2(a) shows the bias, *E*{*ẑ*_{0}^{ML}}−*z*_{0}, of the ML estimator of the axial displacement *z*_{0} of a fluorescent emitter using simulation data of 4-step FPSI measurements with binary (top panel) and Gaussian (bottom panel) spatial masks at a signal-to-noise (SNR) of 10.6 and an effective visibility (EV) of 0.6 for the binary mask and 0.67 for the Gaussian mask. These figures-of-merit are defined in the legend of the figure. As observed from this figure, *ẑ*_{0}^{ML} is effectively an unbiased estimator for both masks. The classical estimator of 4-step FPSI *ẑ*_{0}^{C} is also unbiased for both masks (data not shown; see [30] for a proof).

Figures 2(b) and 2(c) presents the corresponding sample standard deviations of the ML and classical estimators for the binary [Fig. 2(b)] and Gaussian [Fig. 2(c)] spatial masks. Also shown for each mask type are the Cramér-Rao lower bound (CRLB) and the analytical standard deviation of the classical estimator calculated by Eqs. (16) and (18), respectively. We see from these graphs that for both the binary and Gaussian spatial masks (i) the standard deviation of the unbiased ML estimator of *z*_{0} was oscillatory with the axial displacement of the emitter, (ii) the modulation depth of these oscillations, given by the expression 1–{1–0.5⋅EV^{2}}^{½} derived from Eq. (16), was calculated to be ~10% for the EV values used in these simulations, (iii) the standard deviation of the ML estimator of *z*_{0} was lower than that of the classical estimator for all values of *z*_{0}, except for *z*_{0} ≈0 or ± 90 nm. At these axial displacement values, Eqs. (16) and (18) become equal to 0.5⋅SNR^{−2} as $\overline{\Delta \varphi}=0\text{or}\pm \pi $, (iv) the standard deviations of the ML and classical estimators increased with the axial displacement of the emitter (due to the decrease in the number of photons in the central area of the imaged interferograms as defocus increases), (v) the unbiased ML estimator of *z*_{0} is effectively an efficient estimator (because its variance appeared to achieve the CRLB) and is therefore the minimum variance unbiased estimator. We note that the differences in the standard deviations of the ML and classical estimators between the binary and Gaussian masks were small, implying that the simpler binary mask could be a proper choice for retrieving *z*_{0}. As a final remark, we point out that the square of the mean squared error of the above estimators, *E*{[*ẑ*_{0}^{ML/C}−*z*_{0}]^{2}}^{½}, which represents the axial localization precision achieved by the estimators, equals to their standard deviations, as these estimators are unbiased.

To study the dependence of the performances of the ML and classical estimators of *z*_{0} on the SNR and EV of the FPSI interferograms, we plotted their standard deviations at increasing SNR and EV levels by varying the readout/background noise, as presented in Figs. 3(a) and 3(b) for the binary and Gaussian spatial masks, respectively. As observed, for both mask types, the standard deviations of both estimators improved with lower readout/background noise and their percentage change (or modulation depth) increased. We note that the standard deviations of both estimators slightly differed between the binary and Gaussian masks, with higher (lower) values for the binary mask at the lowest (highest) SNR. Also, the modulation depth of the CRLB computed using the binary mask appears to be slightly smaller than that calculated with the Gaussian mask for all SNRs. These differences stem from small differences in the SNR and EV values obtained by these two masks. Similar results were observed when varying the number of photons detected instead of the readout/background noise level (data not shown). Next, we drew in Figs. 3(c) and 3(d) the standard deviations of the ML and classical estimators of *z*_{0} for the binary and Gaussian spatial masks, respectively, at increasing SNRs but similar EVs. This was achieved by changing both the number of photons detected and readout/background noise levels such that EVs were maintained similar between the various SNRs. As in Figs. 3(a) and 3(b), for both mask types, the standard deviations of both estimators improved with higher SNR, but here their modulation depth remained similar across the different SNR levels. These results show that regardless of the spatial weighting mask used to preprocess the interferograms the SNR regulates the standard deviation of the classical estimator of *z*_{0} as well as the peak values of the standard deviation of the ML estimate of *z*_{0}, whereas the EV determines the modulation depth (or percentage change) of the standard deviation of the ML estimator over *z*_{0}.

The percentage change, expressed using Eq. (16) as 1–{1–0.5⋅EV^{2}}^{½} with possible EV values ranging from 0 to 1, can obtain values as high as ~29%, depending on the visibility of the acquired interferograms. This maximum percentage change represents the best possible improvement in axial localization accuracy that the ML estimation can offer over the classical estimation of *z*_{0} in FPSI. For practical interferogram SNRs (>5), this maximum percentage change translates into an absolute difference smaller than 3 nm. To be able to measure this level of difference, the spatiotemporal phase stability of the FPSI instrument needs to be high; otherwise, the percentage change reduces to lower levels, resulting in undetectable sub-nm absolute differences between the ML and classical estimators at high interferogram SNRs, as shown by our experimental results in section 3.3.

#### 3.2 ML and classical z_{0}−estimations from 4-step FPSI simulations over a DOF of ~2λ

_{0}

To extend the applicability of the ML and classical estimates of *z*_{0} to a DOF larger than *λ*/2, we used a pair of weighted inner and outer imaged wavefront phase differences ${\overline{\Delta \varphi}}^{i/o}$, which can unambiguously encode the axial position of a fluorescent emitter over a DOF of ~2*λ* [23, 24]. The ML and classical estimators of ${\overline{\Delta \varphi}}^{i/o}$as well as their CRLBs and variances were calculated using the same procedures outlined in sections 2.2–2.4, but with a set of (possibly overlapping) four–weighted–inner and four–weighted–outer phase-shifted interferograms. Mapping of these estimators to *ẑ*_{0} and their CRLBs and variances to σ^{2}_{ẑ0} was obtained via the use of ${\overline{\Delta \varphi}}^{i}$versus${\overline{\Delta \varphi}}^{o}$lookup curves of a static fluorescent emitter located at (0, 0, *z*_{0}), as detailed in [23, 24]. Figure 4a presents the bias *E*{*ẑ*_{0}^{ML}}−*z*_{0} of the ML estimator of *z*_{0} over an extended DOF using simulation data of 4-step FPSI measurements with binary (top panel) and Gaussian (bottom panel) spatial masks at SNR of ~10 and EV of 0.67. As observed, *ẑ*_{0}^{ML} is effectively an unbiased estimator for both masks. The classical estimator of 4-step FPSI *ẑ*_{0}^{C} is also unbiased for both masks (data not shown). Figures 4(b) and 4(c) shows the corresponding sample standard deviations of the ML and classical estimators for the binary [Fig. 4(b)] and Gaussian [Fig. 4(c)] masks. Also depicted for each mask type are the CRLB and the analytical standard deviation of the classical estimator.

Similarly to the results in Figs. 2(b) and 2(c), the CRLB and the standard deviation of the ML estimator also exhibited here an oscillatory pattern that, however, faded at the edges of the DOF due to the lower EV as defocus increases. Furthermore, as also seen in Figs. 2(b) and 2(c), the standard deviations of the ML and classical estimators deteriorated with increasing axial displacement of the emitter due to the decrease in SNR as defocus increases. Note that the first-order analytical approximation for the standard deviation of the classical estimator (in gray dashed line) is not a tight approximation at the edges of the extended DOF, and higher-orders terms are likely necessary to improve it. Also, it is interesting to note that while the ML estimator of *z*_{0} was still bounded from below over the entire extended DOF by the CRLB and from above by the standard deviation of the classical estimator, the variance of the ML estimator of *z*_{0} started to deviate from the CRLB at larger axial displacements (|*z*_{0}|>~*λ*/3). The degree of this deviation increased with the spatial overlap of the inner and outer masks as seen in Figs. 4(b) and 4(c) (Fig. 4(b), no overlap between the inner and outer binary masks; Fig. 4(c), existing overlap between the inner and outer Gaussian masks). These results show that the ML estimator may not be efficient over the entire extended DOF, so that another unbiased non−ML estimator with a variance lower than that of the ML estimator may exist. It is noteworthy that the differences in Figs. 4(b) and 4(c) in the standard deviations of the ML and classical estimators between the binary and Gaussian masks stem from the smaller size of the inner binary mask that guaranteed no overlap between the inner and outer binary masks, but resulted in a slightly lower SNR, and hence in a slightly poorer axial localization precision than that achieved by the Gaussian masks. Finally, we point out that the performances of the ML and classical estimators of *z*_{0} as a function of the SNR and EV over the extended DOF were similar to those described in section 3.1 for the ~*λ*/2–DOF (data not shown).

#### 3.3 Experimental comparison between the ML and classical z_{0}−estimations in 4-step FPSI

_{0}

The simulation results in sections 3.1–3.2 indicate that the main improvement in the axial localization precision of the ML estimator is within the axial range of |*z*_{0}|<~*λ*/4 to levels of up to ~29% compared to the axial localization precision of the classical 4-step FPSI estimator. To examine the ML estimator with experimental 4-step FPSI data, we built a triangle self-referencing interferometer with two-opposing high-NA microscope objectives (63 × /1.15, Zeiss) and a phase-shifting mirror (MCL) in one of the interferometer's arms. The implementation of this setup is simpler than that of existing FPSI systems [23, 25]. By translating the phase-shifting mirror in steps corresponding to π/2, the system recorded four phase-shifted interferograms of 170−nm fluorescent dark red beads (excitation: 660 nm, emission: 680 nm; Invitrogen) on a sCMOS camera with an exposure time of 4 ms at ~100 frames−s^{−1}. The beads were excited with a wide beam from a 633−nm helium neon laser, and were mounted in between the two objectives in a #1 coverslip sandwich placed onto a nanopositioning stage (Newport), which was used to adjust their axial displacements *z*_{0}, as illustrated schematically in Fig. 1(a). An emission filter centered at ~677 nm with a bandwidth of ~30 nm (Semrock) was attached to the camera to eliminate spurious laser excitation light.

Figures 5(a) and 5(b) shows the mean, minimum and maximum values of the standard deviations (localization precisions) of the ML and classical estimators of the axial displacement of fluorescent beads over |*z*_{0}|<~*λ*/4 using binary and Gaussian spatial weighting masks on the measured interferograms, respectively. Also depicted in Figs. 5(a) and 5(b) are the mean (in gray solid line), and the minimum and maximum values (in gray dashed lines) of the analytical standard deviations of the classical estimator calculated using Eq. (18), with parameters evaluated from the measured interferograms and camera readout noise.

We can see from these figures that regardless of the spatial weighting mask used the performances of the ML and classical estimators were comparable over the entire range of axial displacements. This result stems from the high SNRs (>100) and relatively low EVs (~0.55) of the measured interferograms that led to undetectable sub-nanometer modulation depth levels of the standard deviation of the ML estimator. Such levels were undetectable due to the large dispersion in the values of the standard deviation of both estimators at each axial displacement. This dispersion is probably associated with the variations in the number of photons detected and readout/background noise levels, and with the limited spatiotemporal phase stability of our system. We also observe from Figs. 5(a) and 5(b) that the standard deviations of both estimators were well within the range of standard deviations predicted from Eq. (18), indicating the usefulness of this equation for approximating the axial localization precision. Similar results were obtained when using the ML and classical estimators of *z*_{0} over an extended DOF of ~1.25*λ*, as presented in Figs. 5(c) and 5(d). We note that acquisition of experimental data from a larger DOF was possible, but at the expense of a lower statistical precision in the estimations of *z*_{0} due to photobleaching of the beads.

Finally, we point out that in order that the ML estimator derived here for 4-step FPSI be advantageous over the classical PSI estimate, FPSI interferograms with SNRs lower than ~5 and with good visibility (close to 1 for achieving the maximum improvement in the axial localization precision) should be obtained. Such measurement conditions would yield an absolute improvement of several nanometers in the axial localization precision of the ML estimator compared to the classical estimator. Experimentally, this would require the development of an FPSI system with a high spatiotemporal phase stability and a careful control of the SNR of the fluorescent emitter.

## 4. Conclusions

The overall aim of this paper was to put on a firm ground the adoption of classical PSI approaches in FPSI nanoscopy by performing a fundamental and comprehensive study of the estimation of the axial displacement in FPSI using estimation theory tools. Specifically, we derived the ML estimator of the axial displacement of fluorescent emitters from 4-step FPSI measurements plagued by shot noise and readout/background noise. The ML estimator was developed for two cases: when the axial displacement is within a DOF of ~*λ*/2 and when it is within an extended DOF of ~*2λ*. Using Monte-Carlo simulations of FPSI, we have shown that regardless of the scheme adopted to preprocess the measured interferograms the ML estimator is characteristically oscillatory, unbiased (i.e., its mean is equal to the actual value of the axial displacement) over the DOF, and efficient (i.e., the variance on the error of the estimate is given by the CRLB) over part of the extended DOF.

By comparing the performance of the ML estimate and the classical 4-step classical phase shifting estimator, and regardless of the spatial weighting masks used to preprocess the FPSI interferograms, we found that the relative difference in the standard deviations of these two estimators can be as high as ~29% (depending on the EV of the measured interferograms), yielding maximum absolute differences of a few subβ−nanometers to some nanometers for practical interferogram SNRs (>5). These results suggest that in practice the use of the computationally simpler classical 4-step phase shifting estimation together with simple binary weighting interferogram preprocessing schemes is adequate to accurately determine axial displacements in FPSI, as we also validated experimentally. FPSI interferograms with lower SNR but good visibility (>~0.7) can benefit from the use of the ML estimator, provided that the spatiotemporal phase stability of the FPSI system is high.

As a final remark, the ML estimation framework developed here for FPSI can be extended to other wide-aperture interferometry methods used for quantitative phase measurements, such as quantitative phase imaging and optical coherence microscopy [29, 30].

## Funding

Israel Science Foundation (ISF) (1599/12); FP7 People: Marie-Curie Actions (293932).

## Acknowledgments

We thank Elad Arbel and Itay Remer from Agilent for the initial design and construction of the experimental setup.

## References and links

**1. **H. P. Kao and A. S. Verkman, “Tracking of single fluorescent particles in three dimensions: Use of cylindrical optics to encode particle position,” Biophys. J. **67**(3), 1291–1300 (1994). [CrossRef] [PubMed]

**2. **V. Levi, Q. Ruan, and E. Gratton, “3-D particle tracking in a two-photon microscope: Application to the study of molecular dynamics in cells,” Biophys. J. **88**(4), 2919–2928 (2005). [CrossRef] [PubMed]

**3. **E. Toprak, H. Balci, B. H. Blehm, and P. R. Selvin, “Three-dimensional particle tracking via bifocal imaging,” Nano Lett. **7**(7), 2043–2045 (2007). [CrossRef] [PubMed]

**4. **M. A. Thompson, J. M. Casolari, M. Badieirostami, P. O. Brown, and W. E. Moerner, “Three-dimensional tracking of single mRNA particles in Saccharomyces cerevisiae using a double-helix point spread function,” Proc. Natl. Acad. Sci. U.S.A. **107**(42), 17864–17871 (2010). [CrossRef] [PubMed]

**5. **Y. Shechtman, L. E. Weiss, A. S. Backer, S. J. Sahl, and W. E. Moerner, “Precise three-dimensional scan-free multiple-particle tracking over large axial ranges with tetrapod point spread functions,” Nano Lett. **15**(6), 4194–4199 (2015). [CrossRef] [PubMed]

**6. **A. von Diezmann, Y. Shechtman, and W. E. Moerner, “Three-dimensional localization of single molecules for super-resolution imaging and single-particle tracking,” Chem. Rev. **117**(11), 7244–7275 (2017). [CrossRef] [PubMed]

**7. **A. Yildiz, J. N. Forkey, S. A. McKinney, T. Ha, Y. E. Goldman, and P. R. Selvin, “Myosin V walks hand-over-hand: single fluorophore imaging with 1.5-nm localization,” Science **300**(5628), 2061–2065 (2003). [CrossRef] [PubMed]

**8. **L. S. Churchman, Z. Ökten, R. S. Rock, J. F. Dawson, and J. A. Spudich, “Single molecule high-resolution colocalization of Cy3 and Cy5 attached to macromolecules measures intramolecular distances through time,” Proc. Natl. Acad. Sci. U.S.A. **102**(5), 1419–1423 (2005). [CrossRef] [PubMed]

**9. **L. Holtzer and T. Schmidt, “The tracking of individual molecules in cells and tissues,” in *Single Particle Tracking and Single Molecule Energy Transfer*, C Bräuchle, DC Lamb, J Michaelis eds. (Wiley-VCH, Weinheim, 2010).

**10. **D. Li, J. Xiong, A. Qu, and T. Xu, “Three-dimensional tracking of single secretory granules in live PC12 cells,” Biophys. J. **87**(3), 1991–2001 (2004). [CrossRef] [PubMed]

**11. **R. M. Dickson, D. J. Norris, Y. L. Tzeng, and W. E. Moerner, “Three-dimensional imaging of single molecules solvated in pores of poly(acrylamide) gels,” Science **274**(5289), 966–968 (1996). [CrossRef] [PubMed]

**12. **M. Speidel, A. Jonás, and E. L. Florin, “Three-dimensional tracking of fluorescent nanoparticles with subnanometer precision by use of off-focus imaging,” Opt. Lett. **28**(2), 69–71 (2003). [CrossRef] [PubMed]

**13. **S. Ram, P. Prabhat, J. Chao, E. S. Ward, and R. J. Ober, “High accuracy 3D quantum dot tracking with multifocal plane microscopy for the study of fast intracellular dynamics in live cells,” Biophys. J. **95**(12), 6025–6043 (2008). [CrossRef] [PubMed]

**14. **S. Ram, D. Kim, R. J. Ober, and E. S. Ward, “3D single molecule tracking with multifocal plane microscopy reveals rapid intercellular transferrin transport at epithelial cell barriers,” Biophys. J. **103**(7), 1594–1603 (2012). [CrossRef] [PubMed]

**15. **B. Nitzsche, F. Ruhnow, and S. Diez, “Quantum-dot-assisted characterization of microtubule rotations during cargo transport,” Nat. Nanotechnol. **3**(9), 552–556 (2008). [CrossRef] [PubMed]

**16. **T. Ragan, H. Huang, P. So, and E. Gratton, “3D particle tracking on a two-photon microscope,” J. Fluoresc. **16**(3), 325–336 (2006). [CrossRef] [PubMed]

**17. **H. Cang, C. S. Xu, D. Montiel, and H. Yang, “Guiding a confocal microscope by single fluorescent nanoparticles,” Opt. Lett. **32**(18), 2729–2731 (2007). [CrossRef] [PubMed]

**18. **W. Supatto, S. E. Fraser, and J. Vermot, “An all-optical approach for probing microscopic flows in living embryos,” Biophys. J. **95**(4), L29–L31 (2008). [CrossRef] [PubMed]

**19. **M. F. Juette and J. Bewersdorf, “Three-dimensional tracking of single fluorescent particles with submillisecond temporal resolution,” Nano Lett. **10**(11), 4657–4663 (2010). [CrossRef] [PubMed]

**20. **V. Levi and E. Gratton, “Three-dimensional particle tracking in a laser scanning fluorescence microscope,” in *Single Particle Tracking and Single Molecule Energy Transfer*, C. Bräuchle, D. C. Lamb, and J. Michaelis, eds. (Wiley-VCH, 2010).

**21. **G. Shtengel, J. A. Galbraith, C. G. Galbraith, J. Lippincott-Schwartz, J. M. Gillette, S. Manley, R. Sougrat, C. M. Waterman, P. Kanchanawong, M. W. Davidson, R. D. Fetter, and H. F. Hess, “Interferometric fluorescent super-resolution microscopy resolves 3D cellular ultrastructure,” Proc. Natl. Acad. Sci. U.S.A. **106**(9), 3125–3130 (2009). [CrossRef] [PubMed]

**22. **C. von Middendorff, A. Egner, C. Geisler, S. W. Hell, and A. Schönle, “Isotropic 3D Nanoscopy based on single emitter switching,” Opt. Express **16**(25), 20774–20788 (2008). [CrossRef] [PubMed]

**23. **D. Aquino, A. Schönle, C. Geisler, C. V. Middendorff, C. A. Wurm, Y. Okamura, T. Lang, S. W. Hell, and A. Egner, “Two-color nanoscopy of three-dimensional volumes by 4Pi detection of stochastically switched fluorophores,” Nat. Methods **8**(4), 353–359 (2011). [CrossRef] [PubMed]

**24. **E. Arbel, A. Praiz, and A. Bilenca, “Fluorescence phase-shifting interferometry for axial single particle tracking: a numerical simulation study,” Opt. Express **22**(16), 19641–19652 (2014). [CrossRef] [PubMed]

**25. **F. Huang, G. Sirinakis, E. S. Allgeyer, L. K. Schroeder, W. C. Duim, E. B. Kromann, T. Phan, F. E. Rivera-Molina, J. R. Myers, I. Irnov, M. Lessard, Y. Zhang, M. A. Handel, C. Jacobs-Wagner, C. P. Lusk, J. E. Rothman, D. Toomre, M. J. Booth, and J. Bewersdorf, “Ultra-high resolution 3D imaging of whole cells,” Cell **166**(4), 1028–1040 (2016). [CrossRef] [PubMed]

**26. **E. W. Rogala and H. H. Barrett, “Phase-shifting interferometry and maximum-likelihood estimation theory,” Appl. Opt. **36**(34), 8871–8876 (1997). [CrossRef] [PubMed]

**27. **F. Aguet, D. Van De Ville, and M. Unser, “A maximum-likelihood formalism for sub-resolution axial localization of fluorescent nanoparticles,” Opt. Express **13**(26), 10503–10522 (2005). [CrossRef] [PubMed]

**28. **F. Huang, T. M. Hartwich, F. E. Rivera-Molina, Y. Lin, W. C. Duim, J. J. Long, P. D. Uchil, J. R. Myers, M. A. Baird, W. Mothes, M. W. Davidson, D. Toomre, and J. Bewersdorf, “Video-rate nanoscopy using sCMOS camera-specific single-molecule localization algorithms,” Nat. Methods **10**(7), 653–658 (2013). [CrossRef] [PubMed]

**29. **S. Chen, C. Li, and Y. Zhu, “Sensitivity evaluation of quantitative phase imaging: a study of wavelength shifting interferometry,” Opt. Lett. **42**(6), 1088–1091 (2017). [CrossRef] [PubMed]

**30. **A. Dubois, “Phase-map measurements by interferometry with sinusoidal phase modulation and four integrating buckets,” J. Opt. Soc. Am. A **18**(8), 1972–1979 (2001). [CrossRef] [PubMed]