Experimental 3D super-localization with Laguerre–Gaussian modes

Improving three-dimensional (3D) localization precision is of paramount importance for super-resolution imaging. By properly engineering the point spread function (PSF), such as utilizing Laguerre–Gaussian (LG) modes and their superposition, the ultimate limits of 3D localization precision can be enhanced. However, achieving these limits is challenging, as it often involves complicated detection strategies and practical limitations. In this work, we rigorously derive the ultimate 3D localization limits of LG modes and their superposition, specifically rotation modes, in the multi-parameter estimation framework. Our findings reveal that a significant portion of the information required for achieving 3D super-localization of LG modes can be obtained through feasible intensity detection. Moreover, the 3D ultimate precision can be achieved when the azimuthal index l is zero. To provide a proof-of-principle demonstration, we develop an iterative maximum likelihood estimation (MLE) algorithm that converges to the 3D position of a point source, considering the pixelation and detector noise. The experimental implementation exhibits an improvement of up to two-fold in lateral localization precision and up to twenty-fold in axial localization precision when using LG modes compared to Gaussian mode. We also showcase the superior axial localization capability of the rotation mode within the near-focus region, effectively overcoming the limitations encountered by single LG modes. Notably, in the presence of realistic aberration, the algorithm robustly achieves the Cramér-Rao lower bound. Our findings provide valuable insights for evaluating and optimizing the achievable 3D localization precision, which will facilitate the advancements in super-resolution microscopy.


I. INTRODUCTION
Precise three-dimensional (3D) localization is essential for a variety of advanced microscopy techniques, including defect-based sensing [1,2], multiplane detection [3,4], single-particle tracking [5][6][7].Especially, as the position of individual fluorophores can be determined with a precision smaller than the size of the point spread function (PSF), a super-resolution image can be assembled from the estimated positions of the sufficient fluorescent labels [8][9][10][11], to avoid the Abbe-Rayleigh criterion [12,13].The effective achievable resolution is closely related to the precision of individual fluorophore localization.The 3D localization is intimately connected to multi-parameter estimation.Therefore, the advantages of these techniques are better understood in the multi-parameter framework.
Building upon the pioneering work of Tsang and coworkers in quantifying far-field two-point superresolution [14][15][16], the ultimate precision limits of single point source's localization have been extensively investigated [17][18][19].The basic idea is to exploit the quantum Fisher information (QFI) and associated quantum Cramér-Rao bound (QCRB) [20,21].Its theory proves to be effective in establishing under which conditions systems may be preferable to estimate parameters.Com- * lijian.zhang@nju.edu.cnmon imaging strategies rely on complicated experimental setups to achieve the quantum-mechanical bound, such as interferometer arrangement [22,23], mode sorter [24,25], and spatial-mode demultiplexing (SPADE).The multi-parameter cases often involve trade-off relations among the uncertainties on the parameter, since the total information has now to be apportioned [26,27].This makes the multi-parameter optimization problem more involved and more intriguing.
Given the close relationship between QFI and the characteristics of the light field, precision can be significantly enhanced by modifying the system's response, such as through PSF engineering [28][29][30].Recent studies have demonstrated the immense potential of LG modes, as opposed to conventional Gaussian mode, for improving 3D localization precision [31].The superposition of LG modes, specifically the rotation mode characterized by their intensity profiles that rotate on propagation, further enhances the ultimate precision of axial localization.Moreover, the ultimate axial precision of LG modes can be achieved with intensity detection [32,33].The simplicity and feasibility of intensity detection make it extremely valuable for microscopy applications, circumventing potential systematic errors and losses inherent in complex strategies delineated previously [34].However, practical intensity detection introduces pixelated readouts and inherent detection noise, which can compromise the theoretical advantages offered by this simple optimal strategy.Therefore, rigorous mathematical analysis and robust localization algorithms are necessary.
In this work, we rigorously derive the ultimate 3D localization limits of LG modes and rotation modes in the multi-parameter estimation framework.We investigate the accessibility of these limits under ideal intensity detection conditions.Furthermore, we consider the practical limitations imposed by finite pixel size and various detector noise sources.Taking these factors into account, we develop a robust maximum likelihood estimation (MLE) algorithm that iteratively determines the 3D position of a point source.Our algorithm robustly achieves the CRB under low signal-to-noise ratio (SNR) and aberrational conditions.Moreover, it is not limited to symmetric PSFs, but can be extended to accommodate anisotropic PSFs, as well as diverse noise statistics.In this manner, we present comprehensive theoretical and experimental evidence aiming at exhibiting remarkable super-localization capabilities facilitated by LG and rotation modes.

II. THEORETICAL FRAMEWORK
Localization can be regarded as a multi-parameter estimation problem, aiming to determine the 3D coordinates of a point source in the image space, as depicted in Fig. 1 Here, the operators px = −i∂ x and py = −i∂ y represent the lateral displacement as momentum operators, and the axial displace operator is denoted as Ĝz = where k is the wavenumber and ∇ 2 T = ∂ xx + ∂ yy .These operators commute with each other, enabling simultaneous measurement of these unknown parameters [35,36].Through the aforementioned approach, the state in the image space, represented by ρ θ = | Ψ Ψ|, is parameterized with point source 3D coordinates θ = (x e , y e , z e ).The image is a magnified replica of the state |Ψ (x ′ ,y ′ ,z ′ ) in the object space [37].The estimation of 3D coordinates of the object requires a parametric transformation related to the magnification factor of the system.
To proceed further, we consider a shifted normalized LG mode, with a transverse field in the detection plane, given by [38] LG To streamline the derivation, we redefine the coordinate system with the detection plane serving as the origin, denoted as z = 0 and z e represents the distance between the detection plane and the focal plane.
The 3D localization precision is quantified by the covariance matrix Cov(Π, θ) of locally unbiased estimators, which is lower bounded by Cramér-Rao bound (CRB) and QCRB: where N is the number of system copies related to the effective photon counts in each frame.The associated classical Fisher information matrix (CFIm), denoted as F (ρ θ , Π) is defined by where P (X k |θ) is the conditional probability density of observing an outcome X k depending on the underlying 3D position θ of the source, and a specific measurement Π.The associated QFI matrix (QFIm), denoted as Q(ρ θ ), gives the maximum of the CFIm.In the case of pure states, as is the situation we consider, the QFIm is four times the covariance matrix of the generators, which consists solely of diagonal entries.The axial QFI for arbitrary LG modes has been recently worked out in Ref. [33], we extend the results into a 3D scenario, which can be expressed as: (5) The lateral QFI exhibits a linear dependence on p, while the axial QFI shows a quadratic dependence on p, highlighting the significant role of the radial index in effectively enhancing 3D localization precision.These findings are consistent with the numerical results presented in Ref. [31] for low-order LG modes.Notably, the LG 00 mode is the Gaussian mode, serving as a classical benchmark for comparative analysis.
While LG modes exhibit trivial divergence during propagation, more intricate intensity transformations can be achieved by superposing different LG modes.In particular, when employing a set of M constituent LG modes that satisfy the relation [(2p+l) j+1 −(2p+l) j ]/[l j+1 −l j ] ≡ const ≡ V for j = 1, 2, ..., M − 1, the resulting intensity pattern exhibits anisotropy (with non-circular symmetry) and undergoes rotation during the propagation [39,40].As illustrated in Fig. 1(b), the overall rotation angle from the waist to the far field is given by ∆φ(z e = ∞) = V π/2, and ∆φ(z e = −∞) = −V π/2.Notably, half of ∆φ is obtained at the Rayleigh range.These PSFs offer a wider range of applications in superresolution imaging due to their superior localization precision compared to circular symmetric PSFs [11,41,42].To provide a clear physical intuition, we consider the simple example of the superposition of two LG modes (M=2) with l = l ′ and p = p ′ = 0.The 3D QFI for rotation modes can be determined as follows: In contrast to the quadratic precision improvement in the axial localization, the lateral QFI is obtained by summing up the contributions from individual modes.Typically, the QFI is distributed between the phase and intensity variations of the measured beam.Remarkably, by discarding the phase information, the full axial QFI can still be extracted [33].This result prompts us to investigate the efficacy of intensity detection in achieving ultimate 3D localization precision.When ideal detection conditions are assumed, i.e. a detector of infinite area without pixelation and no additional noise sources except for shot noise, the intensity detection projects the quantum state onto the eigenstates of spatial coordinates, represented as Π x,y = |x, y x, y|.In consequence, the probability density is P (X k |θ) = Tr(ρ θ Π x,y ) = | Ψ(xe,ye,ze) | 2 , which corresponds to the (normalized) beam intensity |LG lp | 2 .The summation in Eq. 4 transforms into a twodimensional integral over the spatial domain.For LG modes, after a lengthy calculation, the ideal CFI can be expressed analytically as follows: (7) These results are plotted in Fig. 2. Two detection planes can be found, where full axial QFI and a portion of lateral QFI can be extracted.In the case of certain LG modes with |l| = 0, the lateral CFI can reach the QFI at the beam waist, F x,y (0) = Q x,y , while the axial CFI can reach it at the Rayleigh range, F z (±z R ) = Q z .Intuitively, the precision of axial localization depends on the beam divergence, which determines the rate of beam width variations during propagation.The precision of lateral localization benefits from sharpness of the wave packet.By increasing the radial index p, the sharpness and the beam divergence are enhanced [43], leading to an improvement in the precision of 3D localization.However, as the azimuthal index l is increased, a central dark spot emerges and expands in size.Although this leads to increased divergence, it fails to enhance the sharpness of the beam, thereby hindering improvements in the precision of lateral localization.In the case of rotation modes, we consider the superposition of LG 00 and LG 20 modes as a representative example.Numerical analysis suggests that only a small fraction of the 3D QFI can be extracted with intensity detection.As this specific category lacks radial information, the average lateral CFI is equivalent to that of the LG 00 , see Fig. 2(a).However, the nonstationary rotation behavior significantly enhances axial CFI in the near-focus axial region compared to single LG mode, as shown in Fig. 2(b).These results also underscore the significance of developing quantum-inspired strategies, such as SPADE or mode sorting techniques, to reveal all information about the parameter.Comprehensive derivations of the QFIm and ideal CFIm are included in the Appendix A.

Rotation
LG 01 LG 00 We then incorporate the deteriorating effects in the detection process and present our MLE algorithm.Previous studies have demonstrated that the localization algorithms based on MLE can asymptotically approach the CRB for a few specific scenarios [44][45][46], outperforming nonlinear least squares (NLLS) algorithm [47].
However, discrepancies between the variances of MLE and the precision predicted by the CRB have been observed in the presence of model mismatches and misspecifications [25,46], such as inaccurate noise statistics, PSF mismatch, optical aberrations, and low SNR.These limitations stem from the fact that existing localization MLE algorithms make simplified statistical assumptions for specific scenarios, which inspires us to improve the robustness and generalisability of MLE algorithms by adopting a more refined statistical model.In the case of pixelated intensity detection, the measurements can be described as Π A k = A k |x, y x, y|dxdy.The readout counts in the kth pixel encompass the integrated photon counts N k and contributions from detector noise, which can be expressed as: The photon-electric conversion process in the CCD camera distorts the effective photon counts, resulting in two types of detection noise.The term N b represents signalindependent noise, including background fluorescence, dark current, and readout noise.On the other hand, the variance of N c positively correlates with the signal and arises in the electron amplification process [48].The SNR is determined by the ratio of the average effective photon counts to the noise present at each pixel.These statistical assumptions have been successfully applied in weak measurement scenarios with limited SNR and detector dynamic range [49,50].
Based on the statistical model mentioned above, we describe the MLE algorithm to estimate the parameters with the likelihood function: The estimated results, denoted as θ, maximize the likelihood function.We employ the scoring method to iteratively update the parameter estimates using the inverse of the CFIm and the derivative of the likelihood function: The iterative update scheme is similar to previous approaches in Refs.[46,51], but improves iteration stability and reduces computational complexity [52].
For a system with unknown w 0 , the algorithm simultaneously estimates three unknown parameters: θ = (x e , y e , w(z e )).We assume the nominal axial distance is known, and the estimated beam width is used to determine the system's w 0 and z R .The lateral variance var(x e , ŷe ) can be directly calculated, while axial var(ẑ e ) can be derived using error propagation var(ẑ e ) = var( ŵ(z e ))/[∂ ze ŵ(z e )] 2 .
If the system is fully precalibrated with a known w 0 , the algorithm can estimate z e instead of w(z e ) with slight modification to handle anisotropic PSF, enabling direct 3D position estimation, as demonstrated in the following rotation mode experiment.

III. EXPERIMENT A. Experiment setups
To validate our theoretical framework, we conducted experimental 3D localization using imperfect Gaussian mode, LG modes, and the rotation mode.The experimental setup is schematized in Fig. 3.In these experiments, we use the focus beam waist as a simplified point source realization.
We first assess the performance of the MLE algorithm in the presence of model mismatches and misspecifications by 3D localizing an imperfect Gaussian mode.A He-Ne laser at wavelength λ = 633nm serves as the Gaussian source.A single lens in Module (b) is utilized to form the image, which leads to a beam waist w 0 = 77.48µmand corresponding Rayleigh range z R = 29.8mm.The spatially unfiltered laser and the spherical aberration of the image system lead to the imperfection of Gaussian mode [53].A sequence of images is captured at different axial positions using a scientific ICCD camera (Andor, iStar CCD 05577H) with pixel size 13µm × 13µm.The axial positions range over a span of 60mm with an interval of 3mm.At each position, we acquire 600 intensity images.The pre-calibration noise are characterized by N b ∼ N 515.6, 7.1 2 and N c ∼ N 0, σ 2 c .Here, ln σ 2 c = 1.4 ln (N k ) − 0.7 and N k represents the effective photon counts per pixel.It is worth noting that while Gaussian noise aligns with our CCD response calibration, other statistical distributions can also be accommodated in the algorithm.The effective photon counts per image are approximately N = 1 × 10 4 , obtained by subtracting the detector noise from the total readout.These conditions reflect the typically low SNR encountered in real microscopy scenarios.
We then compare the 3D localization precision of LG modes and the rotation mode, demonstrating their superlocalization capabilities.The desired PSFs are generated using a double-fourier transform optical setup, as depicted in Module (a).Computer-generated holograms (CGH) are imprinted onto the SLM (Meadowlark Optics, P1920-400-800-HDMI), with the desired first-order diffraction selected by a 4f system [54][55][56].The waist radius of CGHs is uniformly set to 500µm.The second lens of the 4-f system is slightly displaced from the focal plane to ensure proper beam focusing.The modulation efficiency of the SLM is adjusted by an HWP.To mitigate the adverse impacts of beam jitter and turbulence, an additional HWP and a BD in Module (c) are employed to create two copies of the output beam, namely the lef t and right beams.These negative effects have less impact considering the previous single-lens imaging systems.By subtracting the estimation results obtained

B. Experiment results
The experimental results related to imperfect Gaussian mode are summarized in Fig. 4. Our algorithm is referred to as Model A. As a comparison, we employ another widely used localization MLE algorithm, referred to as Model B, as described in Ref. [46].The key distinction between the two models lies in Model B's consideration of a Poisson noise source N b while neglecting the contribution of N c .Additionally, Model B treats N and N b as unknown parameters, and offers direct estimation of θ = (x e , y e , w(z e ), N, N b ).As the sub-CFIm of θ = (w(z e ), N, N b ) contains non-diagonal entries, N and N b can often be thought of nuisance parameters [57].
As illustrated in Fig. 4, the ultimate localization precision is given by the QCRB, while CRB of ideal intensity detection reaches QCRB at the focus and the Rayleigh range.The presence of noise diminishes the precision achievable with the ideal CRB.To address this issue, the practical CRB can be obtained using the statistical models, as discussed in detail in the Appendix C. With increasing the distance between the detection plane and the focal plane, the SNR decreases, widening the gap between the practical CRB and the ideal one.By employing a more refined statistical model, Model A accurately predicts the CRB.Conversely, Model B provides a significantly non-tight CRB due to noise misspecification, specifically the Poisson assumption overestimating the variance of the noise.Both algorithms obtain similar precision in the lateral direction, see Fig. 4(a)-(b).However, as the distance increased, substantial discrepancies in the variance of axial localization become apparent, as shown in Fig. 4(c).We infer that these discrepancies may be attributed to nuisance parameters and model mismatch, which often result in reduced algorithm precision [47,58].The inaccurate estimation of w(z e ) occurs due to the biased estimations of N and N b when there exist PSF mismatches caused by aberrations, as LG01 LG11 LG02 LG12 LG20 discussed in the Appendix B. The limitation of Model B in achieving the axial CRB has also been observed in previous works [25,46].These results align with our intuition that a more refined model, taking into account more characteristics of the experimental apparatus, enhances the algorithm's precision and robustness against model misspecification and mismatches.
The experimental results of 3D localizing a set of LG modes are presented in Fig. 5(a)-(c).The images are captured at 10µm intervals within a 100µm range, while the Rayleigh range is approximately 52.71µm.By ensuring modulation accuracy and efficiency, continual improvement in 3D localization precision can be achieved using higher-order modes.Specifically, for the highest-order mode LG 12 generated in our experiment, the variance of the lateral localization is 0.86 × 10 −3 mm 2 and the variance of the axial localization is 0.25 × 10 3 mm 2 .In comparison, the LG 00 mode, which serves as the classical benchmark, achieves a lateral localization precision of 2.03 × 10 −3 mm 2 and an axial localization precision of 7.28 × 10 3 mm 2 under the same detection conditions.The results exhibit an enhancement of up to two-fold in lateral localization precision and up to twenty-fold in axial localization precision when employing LG modes.However, the deleterious effects of pixelation and noise will be exacerbated in higher-order modes.It is imperative to assess these potential limitations in order to achieve the anticipated precision improvement.Additionally, we present the results of the axial localization experiment for the rotation mode in Fig. 5(d), along with the constituent modes LG 00 and LG 20 .Instead of changing the axial distance of the detector plane, a series of holograms are utilized to simulate the propagation of the light field.By performing a pre-calibration of the beam waist, our algorithm enables direct estimation of 3D coordinates.In the vicinity of the focal plane, the CFI of LG modes approaches zero, with their lateral intensity distributions seldom changing during propagation.Moreover, the likelihood function of LG modes exhibits the same values in two axial positions, leading to an ambiguous axial position estimate of ±z.In contrast, rotation modes exhibit a unique and easily detectable rotation angle ∆φ(z e ) as they propagate [33,39].This characteristic eliminates the ambiguity of ±z and greatly improves the precision of axial localization.The experimental results highlight the exceptional precision in the axial localization achieved with the rotation mode in the near-focus region, surpassing that of the LG modes.

IV. CONCLUSION
In summary, our research provides both theoretical and experimental evidence showcasing the exceptional potential of LG and rotation modes for 3D super-localization.To address practical challenges, we develop an iterative MLE algorithm that effectively estimates the 3D positions of point sources with the best possible precision determined by the CRB.By incorporating a refined noise statistic model, our algorithm improves the robustness and generalizability of the localization process, offering significant advantages in scenarios with low SNR and aberrations.
While higher-order or intricate superposition modes demonstrate theoretical advantages in ideal intensity detection, practical experimental imperfections pose challenges in realizing these benefits with PSF engineering and intensity-based strategies.Therefore, when exploring and optimizing the final effective resolution, the practical CRB provided by our algorithm can serve as a reliable benchmark for evaluating precision improvements.Our work builds a bridge between the quantum estimation framework and practical microscopy algorithm, fostering promising advancements in 3D super-resolution microscopy.
In the case of pure states, the QFIm can be expressed as four times the covariance matrix of the generators computed in eigenstates |n + , n − , computed in the eigenstates |n+, n− .The entries of the QFIm for pure LG modes are given by The re-expression of the result in matrix form is given by In the case of rotation modes, we assume an equal superposition of two LG modes with l = l ′ , and p = p ′ = 0: Using the same approach as Eq.A4, the QFIm for rotation modes is given by We proceed to derive the CFIm for ideal intensity detection.In this case, the measurement probability distribution is represented by the normalized intensity distribution of the light field, as Since the detector has infinitesimal pixels, we can obtain the summation form of CFIm calculation through indefinite integration in the spatial domain, which is expressed as: Changing the integration variable t = 2 (x − x e ) 2 + (y − y e ) 2 /w(z e ) 2 and carrying out the derivatives of Laguerre polynomials using the relation ∂ t L l p (t) = −L l+1 p−1 (t), the diagonal elements of CFIm can be obtained in compact expressions: (A9) We first expand the complex integral term as the summation of different Laguerre polynomials: Each integral can be further expanded according to the orthogonal polynomial [61] as The following properties are important for simplifying this summation equation: Complex equations can be simplified and given analytic results as The same routine can be applied to the non-diagonal entries of the matrix, yielding a value of zero.The CFIm of ideal intensity detection for pure LG modes is given by: (A14) We have derived the analytical forms for the QFIm and CFIm as discussed in the main text.However, obtaining the analytical form for the rotational modes proved challenging.To facilitate quantitative analysis, we present numerical results for the rotational modes.Notably, the non-diagonal entries in the x y coordinate matrix are non-zero, indicating mutual influence between the lateral coordinates due to the absence of circular symmetry of the intensity distribution.Consequently, the lateral localization precisions exhibit slight directional differences.

Appendix B: Analyzing Model Mismatch
Model mismatches can arise from various factors, including aberrations in the imaging system or limitations in point spread function (PSF) engineering techniques.The deviation between the observed image and the expected image is commonly quantified using the peak signal-to-noise ratio (PSNR) metric [62,63].Assessing axial aberrations involves examining the beam divergence, which can be evaluated using the M-square parameter [37].
To analyze the lateral mismatch, we employ the Matlab curve fitting tool package for least squares (LS) fitting of the observed image to obtain the expected image.The raw data is averaged and processed with background subtraction and data normalization.The MSE is calculated by comparing the experimental image X of size m × n with the corresponding desired reference image K: The PSNR quality metric is then defined as: where M AX K denotes the maximum value of K.The PSNR results are shown in Fig. 6, revealing high values in the vicinity of the focal plane.LG00(Aberration) LG10 LG01 LG11 LG02 LG12 PSNR [dB] z / z R FIG. 6.An estimation of beam lateral quality, as determined by the peak signal-to-noise ratio for the experimental intensity images.
In the case of generated LG modes, the PSNR diverges into two branches with l = 0 and l = 1, which stems from the fact that images with higher peak values generally yield higher PSNR values.Imperfect Gaussian modes, affected by aberrations, result in lower PSNR values compared to modulated Gaussian modes.Evaluation of axial aberrations is based on the M-square parameter, which represents the ratio of the experimental beam's angular divergence to the theoretical divergence, with the theoretical value given by M 2 lp = 2p + l + 1 [43].The theoretical beam divergence is calculated using the estimated waist radius ŵ0, obtained through LS fitting of the experimentally estimated beam width ŵ(z e ).The performance of the two algorithms differs not only in the axial variance but also in the mean value of the beam width obtained, as depicted in Fig. 7(a).For the imperfect Gaussian mode with l = 0, p = 0, we obtained M-square values of approximately M 2 00 ≈ 1.42 (Model A) and M 2 00 ≈ 1.52 (Model B).These values indicate the presence of significant spherical aberration in our initial single-lens system.In contrast, modes generated through the 4-f system exhibit minimal spherical aberration, demonstrating excellent agreement with the theoretical divergences, as illustrated in Fig. 7(b).LG00 LG10 LG01 LG11 LG02 LG12 (b) Model B Model A FIG. 7. Experimental estimates of beam width for (a) imperfect Gaussian mode and for (b) LG modes at different defocusing distances.The theoretical results (solid lines) are calculated through ŵ0 and the fitting results (dash lines) are obtained using LS fit.

Appendix C: Comparison of Two Localization Algorithms
For pixelated intensity detection, the elements of the CFIm can be obtained by summing the CFIs of each pixel: (C1) In Model A, which falls under the category of mixed noise models involving two types of detector noise N b and N c , the pixel readout probability distribution is given by [49]: where N k represents the effective photon counts in kth pixel, which follows a Poisson distribution P (N k = µ k (x, y)).The detector response, denoted by R(X k |N k ), describes the relationship between the readout counts X k of the kth pixel and the effective photon count N k .The response probability distribution is obtained through the convolution of the probability distributions of N b and N c , which is given by The elements of Model A's CFIm can be written as C4) which involve derivatives of the interest position parameters θ = (x e , y e , w(z e )), and the corresponding precalibrated detector parameters = (N, N b , N c ).The resulting CFIm is a 3x3 matrix given by With a known beam waist w 0 , our program can be configured to directly measure the axial position z instead of the beam radius w.This adjustment can be achieved through the substitution of the partial derivatives associated with the beam radius w within the corresponding matrices, with those pertaining to the axial position z.
In Model B, only involving Poisson noise source N b , the pixel readout probability distribution can be expressed in a simplified form, given by [46]: Utilizing the Stirling approximation (ln n! ≈ n ln n − n for large n), the elements of Model B's CFIM can be simplified as follows: which involve derivatives of five unknown parameters θ = (x e , y e , w(z e ), N, N b ).The resulting CFIm is a 5x5 matrix with a non-diagonal sub-matrix, given by The parameters N and N b can be considered as nuisance parameters, potentially introducing imprecision in waist estimation.In the case of imperfect Gaussian mode localization experiments, our calibration results indicate that the average effective photon counts is N A = 11066.Conversely, direct estimation using Model B yields N B = 9782 and N b ∼ P (516.37), which are biased estimates.The results indicate that the variance of N b in Model B is exaggerated when approximating it as Poisson noise.Consequently, in the presence of aberrations and model mismatch, this bias affects the estimation of the beam width, as illustrated in Fig. 7(a).Pixelation and detection noise have a significant impact on localization precision and the determination of optimal detection plane.In the case of lateral localization, the optimal detection plane coincides with the beam waist, which exhibits the highest SNR.Conversely, for axial localization, the optimal position is at the Rayleigh range, where the SNR is relatively poor.Thus, there exists a trade-off between the SNR and the optimal position, shifting the optimal detection plane towards the waist.The discrepancies between the practical CRB and the ideal CRB exhibit a greater prominence in the axial direction (Fig. 8(b)) as opposed to the lateral direction (Fig. 8(a)), owing to the higher SNR observed at the lateral optimal plane compared to the axial optimal plane.

CRB CRB(Ideal) QCRB Experiment
LG10 LG00 LG01 LG11 LG02 LG12 LG10 LG00 LG01 LG11 LG02 LG12 log[var(z e )] To further investigate the individual effects of pixelation and detector noise, we explore the variation of the ratio of CFI to QFI with increasing waist size in the absence of noise, as illustrated in Fig. 9(a)-(b).It can be observed that the impact of pixelation becomes less significant as the spot size exceeds a certain threshold.Additionally, by employing the waist size utilized in our experimental setup, we analyze the influence of the detector noise in Fig. 9(c)-(d).The ratio of CFI to QFI can be approximated as linearly dependent on the SNR.In contrast to pixelation, detector noise has a more substantial effect on the degradation of localization precision.The limited dynamic range of our detector restricts the achievable SNR in the experiments, resulting in the discrepancies in Fig. 8.These discrepancies can be reduced significantly by improving SNR.
The impact of pixelation and noise is more pronounced The dashed lines mark the waist-to-pixels ratio and the signal-to-noise ratio employed in the experiment.
in higher-order modes.This behavior can be attributed to the narrower central wave packet and the presence of additional side flaps in higher-order modes, rendering them more susceptible to the effects of pixelation.Moreover, the larger patterns associated with higher-order modes also lead to a reduction in the SNR of each pixel.Hence, when dealing with higher-order modes and more sophisticated superposition modes, it becomes even more crucial to conduct a thorough analysis of the actual accuracy based on a refined noise model.

FIG. 1 .
FIG. 1.(a) Schematic illustration of the 3D localization of a point source with an optical microscope-based imaging setup.(b) Schematic illustration of the lateral intensity profiles variation with the rotation mode at different axial positions.

FIG. 2 .
FIG. 2. The (a) lateral CFI and (b) axial CFI for ideal intensity detection as a function of the distance between the detection plane and the focal plane.The average lateral CFI of the rotation mode is equivalent to that of LG00.The horizontal dashed lines indicate the QFI of each mode.Each curve is normalized to the QFI of LG00.

FIG. 3 .FIG. 4 .
FIG. 3. Experimental setups for 3D localization of a point source.The notations used are as follows: (N)PBS, (non) polarized beam splitter; BE, beam expander; SLM, spatial light modulator; AS, aperture stop; HWP, half-wave plate; BD, beam displacer; ICCD, intensified charge-coupled devices.A detailed description of the functions performed by modules (a)-(c) can be found in the text.(d) Measured intensity patterns with cross sections (asterisk) and corresponding fitted profiles (solid curves).
Appendix D: Effect of Pixelation and Noise

FIG. 8 .
FIG. 8. Comparison of the (a) lateral and (b) axial localization precision of different modes in the optimal detection plane.

FIG. 9 .
FIG. 9. (a)  the ratio between lateral CFI and QFI and (b) the ratio between axial CFI and QFI as functions of the waist-topixels ratio.(c) the ratio between lateral CFI and QFI and (b) the ratio between axial CFI and QFI as functions of the SNR.The dashed lines mark the waist-to-pixels ratio and the signal-to-noise ratio employed in the experiment.