Localising two sub-diffraction emitters in 3D using quantum correlation microscopy

The localisation of fluorophores is an important aspect of the determination of the biological function of cellular systems. Quantum correlation microscopy (QCM) is a promising technique for providing diffraction unlimited emitter localisation that can be used with either confocal or widefield modalities. However, so far, QCM has not been applied to three dimensional localisation problems. Here we show that quantum correlation microscopy provides diffraction-unlimited three-dimensional localisation for two emitters within a single diffraction-limited spot. By introducing a two-stage maximum likelihood estimator, our modelling shows that localisation precision scales as $1/\sqrt{t}$ where $t$ is the total detection time. Diffraction unlimited localisation is achieved using both intensity and photon correlation from Hanbury Brown and Twiss measurements at as few as four measurement locations.


I. INTRODUCTION
The quest for super-resolution imaging and localisation of fluorophores is one of the most important challenges for optical microscopy [1][2][3][4][5], a convenient imaging modality for both functional and live-cell imaging.However, the length scale for visible light restricts imaging to of order 500 nm because of the diffraction limit.This length scale is large compared with many biological regions of interest, and this limitation has spurred research into techniques that go beyond the diffraction limit [6][7][8], and are ideally diffraction unlimited [9,10].
As well as improving resolution, capturing the three-dimensional structure of the cellular and sub-cellular environment is also vital to understanding its function.Several techniques exist for three-dimensional imaging, some of which incorporate super-resolution techniques in whole or in part.Examples of three-dimensional imaging techniques are z-stack [11][12][13], light sheet microscopy [14][15][16] and total internal reflection fluorescence microscopy [17][18][19].
Both light sheet and total internal reflection fluorescence microscopy require complex optical configurations/components.
Our goal has been to study the role that quantum effects play on resolution.In particular, we have been studying quantum correlation microscopy (QCM) [20][21][22].This approach uses both intensity and photon correlation from Hanbury Brown and Twiss measurements (HBT), to improve localisation and/or imaging, and has been shown to be useful in both widefield [23,24] and confocal [6,25] imaging modalities.Our approach has more in common with z-stack in that we consider a confocal arrangement using discrete focal points in threedimensional space.
We first develop the analytics of three-dimensional localisation using quantum correlation microscopy, based on which we develop a maximum likelihood estimator (MLE) using an approximated likelihood function.We confirm our results using Monte-Carlo (MC) simulations for the diffraction unlimited two-emitter localisation process.Our methodology allows the efficient usage of information from each photon to minimize the detection time and total intensity, and thereby to mitigate possible phototoxicity observed in super-resolution approaches [26].Our method is particularly suited for high photostability single-photon emitters such as nitrogen-vacancy (NV) centres in diamond as biomarkers [27].This paper is organised as follows.First we describe the theory of localisation of two emitters in three dimensions, using both intensity and HBT information.We then develop a MLE estimator to realize the 3D localisation algorithm Using the MLE, we run MC simulations for the 3D localisation of two randomly generated emitters, including different detection configurations.Finally we exhibit the scaling law of resolution, demonstrating the diffraction-unlimited potential of the approach.

II. SYSTEM
We consider two single photon emitters of unequal brightness, positioned in three dimensions within one diffraction limit.Fig. 1 shows the schematic of 3D HBT measurement.The emitter positions are at For simplicity, we assume that the emitters are excited via an external light source so that the two emitters are equally illuminated.Collection is performed by a single confocal microscope objective that is capable of being translated in the x, y, and z directions, which leads to translation of the focal point.We define the focal points as ξ j ∈ R 3 , where j ranges of over all of the focal points used for a particular localisation scheme.This configuration means that the probability of exciting each emitter is equal, though the probability of collecting photons from the emitters varies depending on the precise location of each emitter within the microscope point spread function.Photons from the two emitters are directed via a 50/50 beamsplitter to two single photon counting modules, e.g.avalanche photodiodes (APD), as shown in Fig. 1.The measurements from the APDs are read out both independently, to obtain intensity, and in coincidence, to obtain the HBT signal.Assuming each single-photon-emitter emits only one photon during its emission lifetime, correlation events can only occur when each emitter fires one photon and these are received by the two APDs simultaneously.Both the intensity and second-order correlation function of the two-emitter system at zero time delay g 2 (τ = 0) are used for the two-emitter localisation.
We define the intrinsic brightness of emitter i to be P 0,i ; this is the maximum probability of detecting a photon from the emitter when it is located at the centre of the detection point spread function (PSF), i.e. when x i = ξ j for one of the ξ j .The probability of detecting a photon from emitter i is defined by the microscope point spread function, which we treat as Gaussian where r i = x 2 i + y 2 i is the transversal distance from emitter i to the detection optical axis, z is the distance propagated from the plane where the wavefront is assumed flat.Note that this definition of the detection probability also takes into account any additional loss processes in the system.
The beam waist w(z) is the radius of the laser beam at which the irradiance is 1/e 2 (13.5%) of the centre of the beam at z: where w 0 is the Gaussian beam waist in the focal plane and the Rayleigh range, z R is the value of z where the cross-sectional area of the beam is doubled: The coincidence rate is expressed as a function g (2) (τ ), where τ is the time delay between photon detections at each of the APDs.For our purposes we are solely interested in the second order correlation at zero time delay τ = 0 (i.e.simultaneous detection at each APD) as this provides the most quantum information about the relative brightness of the emitters: where α = P 2 /P 1 is the brightness ratio of two emitters (0 < α < 1).In what follows we use g (2) as simplified notation for g (2) (τ = 0).
Each measurement includes an intensity and g (2) .These are both functions of emitter coordinates and relative brightness.Since we are inferring seven independent variables x 1 , y 1 , z 1 , x 2 , y 2 , z 2 , α in total, at least eight equations, that is, four measurement locations are required to solve the problem, i.e. the minimum number of ξ j is j = 1...4.To reduce the possibility that the 3D searching descends into a local minimum, we introduce a scanning function to localise the maximum intensity area, with which we run an intensity scan in 3D and select the top two maximum intensity locations as initial guesses for the optimizer.

III. MAXIMUM LIKELIHOOD LOCALISATION
We model the detected photon counts from emitter i via a Poisson distribution where t is the detection time i.e. the photon acquisition time, and Pois (P i t) is the Poisson distribution with rate P i t.
The total normalised intensity from the two emitters can then be expressed as: The 2nd order correlation is: where c 12 is the detected photon counts from both emitters simultaneously (coincidence), c 11 and c 22 are the uncorrelated events generated by the same emitter 1 or 2, In this paper, the maximum likelihood estimator (MLE) is used to estimate the positions of emitters via measurements I and g (2) .Specifically, let P 0,1 = 1, P 0,2 = α ∈ (0, 1), f be the probability density functions (PDF) of I j and g j , respectively, given the locations of emitters x 1 and x 2 , and the focusing positions, i.e. the maximum of PSFs of the j-th For simplicity of notation, g (2) and g (2) j is written as g and g j , respectively, when they are used as subscripts, e.g.f g (2) = f g and Then, the joint log-likelihood function given the measurements is and, as a result, the MLE is 1:m ).(10) The PDFs of I and g (2) are not available in closed form, so the analytic solution to Eq.( 10) is intractable.However these PDFs approach a Gaussian distribution with increasing t (see Section ?? for more detail).As a result, the likelihood function is approximated by where µ I j (•), σ 2 I j (•), µ g j (•) and σ 2 g j (•) are the mean and variance of the intensity I and the g (2) , as defined in Eq.(??).
Finding the MLE using Eq. ( 11) is still non-trivial, as the log-likelihood involves highly non-linear functions of the parameters.Lagarias et al. used the Nelder-Mead simplex algorithm (NMSA) to solve an equivalent equation to Eq. ( 11) [28].In our case, however, the nonlinearity of the likelihood function results in multiple local minima.As a result, the optimization results given by NMSA are highly dependent on the selection of initial guesses, from which the algorithm typically converges to a local (non-global) minimum.To mitigate this problem, we first roughly estimate the {x 1 , x 2 , α} using a (1-st order) method of moments estimator (MME) [29], i.e.
as the seed for the NMSA.
Since the objective function in the MME described in Eq. ( 12) is simpler than the likelihood function in terms of nonlinearity, the application of MME in this case is less sensitive to the initial guess.The final estimate is found by applying NMSA to find the MLE using the estimate given by the MME as the initial guesses.

IV. EFFECTIVE PSF CHARACTERISATION
In an optical system we have a point spread function (PSF), which is typically defined by the physical properties of the optical arrangement and therefore to determine parameters like the uncertainty in position.Although in our study there is no PSF in this definition, there are probabilities associated with it and that makes it possible to develop an effective PSF that plays the same role on a statistical basis, i.e.where are the emitters most likely to be and how this probability varies as a function of space.However in any given simulation of our study, all we obtain is a single estimate of the emitter locations and their relative brightness.
We use MC simulations to localise two emitters close with respect to the diffraction limit.
After the MC simulation we use equal-cluster k-means to group the inferred positions into two clusters and determine their centroids [30].These centroids of the two k -means clusters are the estimated emitter locations.We define the effective PSF size by determining the radial distance between each cluster centroid and the [1 − 1/ √ e] th of MC inferred locations, ρ 39.5% .If this distribution were exactly Gaussian, without any noise fluctuations, then ρ 39.5% would correspond to the standard deviation of the PSF w 0 , with 2ρ 39.5% being the width of the PSF.We then determine the volume that includes the closest ρ 39.5% of the cluster centroids: The effective width of the point spread function weff is the average of two emitter estimations w eff,i : Fig. 2(a) shows the volumes enclosed the ρ 39.5% of each emitter at acquisition t = 10 4 .To make the visualization clear we project the 3D image into the xoy plane as contours, shown as Fig. 2(b).The contours are the projection of the closest ρ 39.5% inferred locations defined by Eq. ( 13), and the blue and red ones corresponding to emitter 1 and 2 respectively.It shows the evolution of weff with photon acquisition time t, where t is in units of expected time between detections, assuming only one photon is fired from each emitter within one detection time.As t increases more photons, that is, information, would be measured.The where n is the bin value.The histograms show the probability of the model giving the inferred locations with a specific weff .The mode of each histogram is shown as the white solid line overlaying on the histogram map in Fig. 3(b).From the histograms we can see that as the acquisition time increases, the peak i.e. the mode of the histogram gradually shifts to smaller weff , indicating that the resolution of the model increases with the acquisition time t.
There is considerable freedom to choose the focal points, ξ j , and each configuration will give rise to different resolution, with different results for different emitter system.For this reason, we explore several 3D measurement configurations.Motivated by the idea that the measurements could be realized in practice via two sets of two-dimensional HBT trilateration, we present two six-detector configurations [Fig.show the histograms and fitting results for the mode of the distribution for the various measurement configurations in 3D space.

VI. SCALING LAW OF RESOLUTION
To explore the resolution weff scaling with the acquisition time t, we extract the modes of the histograms.The inset Fig. 6 shows the total counts of the histogram as the function of t (here we take the detection configurations with four, six, and nine detections as examples).
We investigate the histogram with the weff limitation [0, 1]w 0 , when a count falls into this region the model successfully achieves localisation resolution under the diffraction limit.
The plateau shows that the majority weff of the 600 random cases are localised with the sub-diffraction resolution.In the short time regime, the counts are low because the low correlation events are insufficient to realize sub-diffraction localisation.And when t > 10 6 , we suspect that limited by the computation power the MC statistics (M = 1000 in this study) is not enough to realise the localisation.Therefore, we study the scaling law within the plateau region in Fig. 6 inset, roughly 10 4 to 10 6 .We then run the linear fit of the modes weff and t in log-log scale as Eq.( 17), which gives the scaling law shown in Eq. (18).The dots are the modes of the histograms in Figs. 4 and 5, and the lines are linear fitting.
The inset shows the time window we choose to fit, within which the model could achieve sub-diffraction localisation for the majority of 600 cases.Fig. 6 shows the modes of histograms weff as scatter dots and the fitting is shown as the solid lines.The scaling law fitting results are shown in Table I.According to [31] the superresolution enhancement goes with a 1/ N photon scaling, where N photon is the total number of photons collected in the entire image spot during the observation time window.Because N photon is proportional to the acquisition time, here we are expecting the resolution scaling law 1/ √ t, so that the slope in Eq. ( 17) b ≈ −0.5.Results in Fig. 6 show that with fewer detections i.e. minimum four, the scaling is worse than 1/ √ t, and increasing the detection numbers could improve the scaling therefore better localisation.The possible causes of the scaling better than 1/ √ t come from both the fitting of the histogram modes and the linear fitting, and the binning in histograms would also change the value of b slightly, however it does not contradict that the theoretical precision scaling with the number of photons collected.Our results show diffraction unlimited localisation with increasing the detection time.We develop an approximated likelihood function as the analysis tool, and investigate several practical detection configurations by applying the estimator to many random ground truths, and the modes in statistics show the majority of the MC simulations generate 1/ (t) scaling in localisation precision.
With the combination of intensity and second-order correlation information, our approach shows that it is achievable to localise sub-diffraction single-photon-emitters in three dimensions.We demonstrate the localisation of only two emitters, but the same approach could possibly be extended to more and ultimately track multiple emitters in three dimensions.
The configurations that we have considered are relatively natural, and require only minor modification of a more conventional z-stack style confocal scan, albeit with relatively few focal points as the minimum number of focal points is four, in two focal planes.We also notice that the precision could be improved by increasing the measurements or adopting more efficient detection configurations, which could possibly be optimised via the investigation of Fisher information matrix in 3D.

FIG. 1 :
FIG. 1: (a) Schematic of HBT system of HBT localisation system.(b) Two emitters (stars) within one 3D Gaussian PSF, which corresponds to the photon collection probability.Excitation is assumed to be uniform.The green plane indicates the objective focal plane, and the red surfaces are the equal intensity surfaces: the outside one is the intensity I = 1/e 2 n i P 0,i to show the 3D Gaussian beam waist, the inside one is the intensity I = 0.5 n i P 0,i to show the 3D focal spot.Through the confocal microscopy photons from the two emitters are directed via a 50/50 beamsplitter into two APDs, and both intensity and HBT signal are obtained.
4(c)], or z-scan microscopy [Fig.5(c)].Figs. 4 and 5 show different measurement configurations identified by the positions of the maximum of the microscope confocal point spread function as it is translated through space.These maxima are indicated by the blue dots in the schematic.From Figs. 4 and 5 we gradually increase the number of detections.Right columns of Figs. 4 and 5

FIG. 4 :FIG. 5 :
FIG. 4: Three of the detection configurations analysed, with corresponding w eff histograms.(a) z-augmented trilateration configuration, with three focal points in the xoy plane as a two-dimensional trilateration and two along z axis.(c) Orthogonal trilateration configuration with two perpendicular trilaterations.(e) Orthogonal trilateration plus: seven focal points with two parallel trilaterations plus one at the origin.(b) (d) and (f) are the corresponding weff histogram maps of the right column configurations, with 600 random ground truths and 1000 MC simulations to localise each ground truth at each time point.

weff = 10 a * t b ( 18 FIG. 6 :
FIG.6: Scaling of the mode of the histogram of weff as a function of total acquisition time.

TABLE I :
The scaling law fitting results of different detection configurations.