Enhancing axial localization with wavefront control

Enhancing the ability to resolve axial details is crucial in three-dimensional optical imaging. We provide experimental evidence showcasing the ultimate precision achievable in axial localization using vortex beams. For Laguerre-Gauss (LG) beams, this remarkable limit can be attained with just a single intensity scan. This proof-of-principle demonstrates that microscopy techniques based on LG vortex beams can potentially benefit from the introduced quantum-inspired superresolution protocol.


Introduction
The Abbe-Rayleigh criterion [1,2] defines the minimum resolvable distance between two closely spaced point sources in an image.It is widely recognized, though, that the resulting limit is inadequate for current quantitative imaging [3].In recent years, a collection of methods, which can be gathered under the rough denomination of superresolution microscopy [4][5][6][7][8][9][10], has emerged, surpassing the length scale set by the Abbe-Rayleigh criterion by more than one order of magnitude.
In the broader context of three-dimensional microscopy, accurately determining the position of fluorescent objects along the optical axis is crucial [26].However, the challenge of finding the minimum resolvable separation in the axial direction has only lately been addressed [27,28].The key idea is to resort to the quantum Fisher information [29,30] and the associated Cramér-Rao bound [31][32][33][34] to get a measurement-independent limit.This builds upon the pioneering work of Tsang and coworkers to quantify transverse two-point resolution [35][36][37][38].
Recently, a mode sorter capable of losslessly projecting photons into the radial LG basis set demonstrated quantum-limited estimation of axial separation for incoherent point sources [39].This strategy bears similarities to techniques used in transverse superresolution [40][41][42].Yet the complexity of experimental setups, introducing systematic errors and losses, can undermine the theoretical advantages offered by this optimal scheme.Surprisingly, it has been shown [43] that this ultimate limit can be easily achieved with a single intensity scan when the detector is placed at one of two optimal transverse detection planes.This simplicity and feasibility make it highly valuable for applications requiring extremely precise localization.Additionally, this method has been extended to include vortex-beam illumination [44].The goal of this paper is to experimentally demonstrate the localization limits in this particular scenario.
This paper is organized as follows.In Sec. 2 we revisit the theoretical framework needed to set the ultimate limits in axial resolution and how to achieve such limits.This ideal scenario is complemented in Sec. 3 by taking into consideration detrimental effects, such as finite pixel size or misalignments, which are usually omitted.In Sec. 4 we present the experimental setup used to control the LG beam wavefronts using a spatial light modulator (SLM).In this way, we demonstrate the information gain provided by vortex beams and the experimental validation of the saturation of the quantum Cramér-Rao bound by intensity detection.Finally, our conclusions are summarized in Sec. 5.

Ultimate limits for axial localization
In our quest to unravel the ultimate limits of axial resolution, we approach the challenge as the estimation of the distance traveled by a vortex beam from its beam waist, situated at = 0, to an arbitrary detection plane located at .At this detection plane, we perform an arbitrary measurement.The acquired data are always affected by noise, so they are effectively represented by a stochastic variable denoted by x, enabling us to construct a robust estimator ˆ for the distance .The inference of the parameter is related to the measurement outcomes through some conditional probability density (x| ) that is dictated by the process at hand.
To quantify the maximum information carried by the measured signal, we turn to the concept of Fisher information (FI), defined as In the quantum domain we use a probe state given by the density operator .The propagation encodes the parameter via the transformation = † , with the unitary operator = exp(− ), where is the generator of translations.The measurement is represented by some positive operator-valued measure (POVM) {Π x }, which comprises a set of positive semidefinite selfadjoint operators that resolve the identity [33].By performing this measurement, we obtain a statistical distribution that, according to Born's rule, is given by (x| ) = Tr( Π x ).Afterward, what remains is to obtain the best estimate of given (x| ).
It is natural to ask whether there is an optimal measurement that should be performed on .The quantum Fisher information (QFI) is defined precisely for this purpose: Q ( ) = sup {Π x } F ( ) and depends exclusively on the probe state .By its very same definition, we have Q ( ) ≥ F ( ).
According to the time-honored quantum Cramér-Rao bound (QCRB), the variance of any unbiased estimator ˆ per detection event of the axial distance is bounded by The right hand side of Eq. ( 2) thus represents the ultimate achievable precision regardless of the measurement and depends exclusively on the beam used in the experiment.For definiteness, we take the beam to be in an LG ℓ mode, with complex amplitude given by [45] LG ℓ ( , , where ( , , expression, the following parameters are the radius of wavefront curvature, the beam radius, and the Gouy phase, respectively, at a distance from the beam waist.Here, R = 2 0 /2 is the Rayleigh length and 0 the beam waist radius [47].
The QFI for an arbitrary LG ℓ mode has been recently worked out in Ref. [44]; the result reads which reduces to Q 00 ( ) = 1/ 2 for the fundamental Gaussian mode LG 00 .Note that this QFI is linear in ℓ, and so the axial localization can be significantly improved by using higher LG modes.
Since we are dealing with a single-parameter estimation, the QCRB can always be saturated with a von Neumann measurement projecting the measured signal on the eigenstates of the symmetric logarithmic derivative of the density matrix [29].This implies projecting the signal onto a complete set of optimal modes, which requires a fairly sophisticated and fragile equipment.Therefore, we consider the performance of the possibly inferior, but much more robust scheme of direct intensity detection, as this is the handiest method available in the laboratory.
We assume that the intensity is sampled by a pixelated detector and the signal is dominated by shot noise, which obeys a Poisson distribution [48]: this neglects nonclassical effects, but it is still a suitable model for realistic microscopy.For the time being we ignore the pixel size so any sampling effect can be omitted.We can take ( , | ) = | LG ℓ ( , , )| 2 , which has to be understood as the probability that we observe a given transverse intensity at ( , ), given a known value of .Then, Eq. (1) admits a closed expression; after a lengthy calculation one gets [44]  This function is plotted in Fig. 1 for the beam LG 04 .We appreciate the existence of two well-defined maxima located precisely at opt = ± R .Actually, maxima occur at opt = ± R for any LG ℓ mode.In these two planes, we have maximal wavefront curvature and one can check that Q ℓ ( R ) = F ℓ (± R ); i.e., the quantum limit is saturated.Therefore, in these two detection planes complete information about the axial distance can be extracted with intensityonly detection.

Detrimental effects
In the ideal scenario discussed earlier, certain imperfections and realistic factors were omitted, which can undermine the anticipated metrological advantages.It is crucial to acknowledge and account for these considerations in any practical setup.
A noteworthy observation is that the integrand in Eq. ( 6) can be aptly interpreted as a Fisher information density (FID), we shall denote by ( , | ).This FID reveals the data with highest sensitivity to variations in the coordinate and is plotted in Fig. 2 for the particular case of a LG 04 beam.An intriguing feature is the presence of two highly sensitive regions: the inner and outer tails of the beam, with the latter being particularly informative.We thus conclude that there is no necessity to gather data from an area larger than 5 times the beam radius, as it contributes no valuable information for axial localization.In fact, collecting data beyond this range would only introduce additional noise, further compromising the precision.
Another point to take into account is the finite pixel size.We have to interpret each pixel as an independent measurement channel, so that where is the conditional probability for the -th measurement channel.We numerically examine this effect.Intuitively, one might anticipate a steady decrease in information as the pixel size increases.This is attributed to the integration of the signal across individual pixels, causing the gradients in Eq. ( 6) to blur.However, as revealed in Fig. 3, the situation is not universally characterized by such a monotonic decrease.Remarkably, when considering a LG 04 beam, we observe that using a pixel LG 04 (normalized to the optimal value Q 04 ) as a function of the pixel size (normalized as / ).For the full line the integration starts at the pixel center, whereas for the dotted line the integration starts at the edge of the individual pixel.(Right) Integrated signal for both strategies for the pixel sizes: a) 0.5 / , b) 1.5 / , c) 2.5 / , d) 3.5 / .size 3 times the beam radius actually yields more information than using a pixel size 2 times the beam radius.This counterintuitive phenomenon arises due to specific pixel sizes that integrate the signal from pixels with low FID, while the edge of the adjacent pixel is close or exactly in the area with highest gradient, with the highest FID.Consequently, the complete LG beam is integrated into a single pixel with significantly large pixel sizes, resulting in a loss of information.It is important to emphasize that this behavior is strongly dependent on the sampling mesh grid employed.The blue dotted curve is the Fisher information for the same mode, but with a lateral centroid deviation of 10 mrad.The red dots denote the maxima of both curves, and thus the position of the corresponding optimal detection planes.For the misaligned beam, this optimal plane shifts toward the beam waist (located at = 0).The brown curve is the normalized differential gradient of the Fisher information with and without misalignment.
Another critical issue in the experiment is the misalignment of the mechanical and optical axes.Due to the different propagation angles, a lateral centroid movement is measurable by the detector.The simple model we use here is based on straight light propagation: where Δ is lateral centroid motion at the detection plane and is angle between the optical and mechanical axis.Because Δ = Δ ( ) depends on the estimated parameter, lateral centroid motion becomes an additional information source in the axial localization problem.This is illustrated in Fig. 4, where we have also included the gradient of the FI variation, which turns important to assess the dominant information source: positive values give lateral centroid motion, while negative values are mainly due to diffraction.

Experiment
Our experimental setup is schematized in Fig. 5.A He-Ne laser at wavelength = 633 nm is spatially filtered by an aperture stop (AS) and expanded and collimated by two doublet lenses, 1 and 2 .This collimated beam illuminates the SLM (CRL OPTO XGA3) with 18 m-pixel size.A computer-generated hologram (CGH) is imprinted on the SLM to generate LG vortex beams with a given azimuthal mode index.The pixelated nature of the SLM leads to multiple diffraction orders, which appear in the image plane of the lens 3 .To separate the diffraction orders, we mixed the CGH with a plane wave and filter using a 4 system with an AS.We also effectively relocate the conjugate plane from SLM to a lens 5 , the true imaging lens in the experiment.The 5 has a focal length of 150 mm, and the CGH is 2 mm in diameter.
The measured intensity is proportional to | LG ℓ | 2 .A compact sCMOS camera (pco.edge4.2 LT) with a very high dynamic range (up to 37500 : 1) was used.The detection noise (assumed as Poissonian) is about 100 electrons (with a standard deviation of 8 electrons), so it is negligible.In addition, the saturation level is of the order of 30000 electrons, which provides a fairly good signal-to-noise ratio.
The beam was characterized with a standard Beam Analyzer software.The analysis is based on interpolating the beam profile and measuring the width at 1/ 2 for several (in our case, eight) cross sections.The Rayleigh range is calculated from the measured beam waist.The beam width for LG 00 was compared with the theory at several planes with good match.
We first calibrated the SLM nonlinear response, getting an input-output relation by using several uniform amplitudes.Next, we employed the well-established Gerchberg-Saxton (GS) iterative algorithm [49] to correct the SLM wavefront error.With just one hundred iterations, guided by an initial estimation of the typical astigmatism shape, we successfully restored the wavefront to a satisfactory state.To further refine the system, we employed a modal decomposition in terms of Zernike polynomials [50], minimizing the residual wavefront errors.The results can be observed in Fig. 5.
We observed a deviation in the SLM correction when dealing with azimuthal mode indexes other than ℓ = 1.This discrepancy was attributed to a systematic error arising from changes in the propagation angle on the SLM with varying ℓ values.To mitigate these errors, we developed a customized data processing, ensuring optimal performance and accuracy.
To check the behavior of the QCRB near the Rayleigh range R we suppose small displacements around this plane = R ± , with ≪ R .We manually moved by the sCMOS camera using a linear stage (STANDA 7T167-100Q) with a positioning precision of 5 m.We collected the data stack at five detection planes located at = −200 m, −100 m, 0 m, 100 m, and 200 m from R .Each stack is comprised of 100 intensity frames with a total number of detections per frame about = 7 × 10 5 events in average (with background subtracted).The same procedure was repeated for other LG ℓ beams up to the azimuthal order ℓ = 7, which was the last measurable beam in our setup.In addition, the lateral centroid motion was subtracted from the data.
In the range of measured values, the changes in intensity are small.This suggests to adopt a polynomial basis so that the average intensity Ī ( ) at the th pixel can be expressed as where is the model coefficient matrix.To estimate this matrix we use the five data sets at planes .Let us denote by Ī and the matrices of components Ī = Ī ( ) and = , respectively.We first solve for the model coefficient matrix using the linear inversion estimator: where the superscript (+) indicates the Moore-Penrose pseudoinverse [51,52].This estimator is also known as the ordinary least squares estimator [53] and it is unbiased and consistent.Under the Gauss-Markov assumptions it is also the best linear unbiased estimator [54].
Once Ĉ is known, we can estimate the axial displacement from the set of frames recorded at the Rayleigh plane, denoted as ¯ ( ) , where = 1, . . ., 100.Now, the relation reads This relation can be inverted by the generalized linear squares method [53] and finally, the statistics of displacement estimates ˆ ( ) is evaluated and FI per single detection event estimated by taking In Fig. 6 we represent the resulting FI computed from our experimental data as a function of the azimuthal mode index ℓ.For comparison, we also plot the QFI.We also include a linear fit of the experimental points that confirms the linear growth of the FI with ℓ.This line although quite parallel to the QFI, is a bit below, because of the abnormally low value of the FI for ℓ = 3, due to the experimental difficulties.The rest of the points are close enough to the ultimate limit.Interestingly, we see an apparent violation of the QFI limit for ℓ = 4, but this is surely due to systematic errors.The results in this figure confirm that our protocol is able to determine the axial distance to an accuracy grater than the classical predictions based in the Abbe-Rayleigh approach, which report a depth of focus of the order of the Rayleigh range.If we take, e.g., the data from the LG 07 beam, we get a standard deviation of 0.000432455 m per detection.With = 7 × 10 5 detections registered, this gives 5.1688 m, beating the classical limit (11.67 mm) by three order of magnitude.Similar estimates can be carried out for any LG beam.

Concluding remarks
In summary, we have unveiled the extraordinary potential of pure LG vortex beams for achieving unparalleled precision in axial localization.Through an extensive experimental analysis, we have uncovered the ultimate limits and demonstrated the remarkable advantages of these beams.Our findings reveal an exciting trend: as the azimuthal mode index increases, so does the precision, unlocking a world of possibilities.However, it is crucial to underline the inherent experimental challenges in controlling these vortex-beam wavefronts with SLMs.
Harnessing the power of FI, we have delved into the realm of detrimental effects that have been overlooked in the existing literature.Finite pixel size and the misalignment of optical and mechanical axes, which are often dismissed, have now been brought to the forefront.
Our results offer new insights into the localization problem and open up new avenues for exploration that might be interesting for 3D microscopy.

Fig. 1 .
Fig.1.Fisher information for direct intensity detection with the mode LG 04 as a function of the axial coordinate .At = ± R , the function F 04 reaches its maximal value; these are the optimal planes to place the detector.

Fig. 2 .
Fig.2.Plots of the cross-section of the intensity distribution 04 (red), the FID 04 (green) and the cumulative FI F (blue) for a pure LG 04 beam, as a function of the dimensionless transverse distance / at the optimal detection plane R .Each curve is normalized to unity.

Fig. 3 .
Fig. 3. (Left) Fisher information in direct detection at the optimal plane for the beam LG 04 (normalized to the optimal value Q 04 ) as a function of the pixel size (normalized as / ).For the full line the integration starts at the pixel center, whereas for the dotted line the integration starts at the edge of the individual pixel.(Right) Integrated signal for both strategies for the pixel sizes: a) 0.5 / , b) 1.5 / , c) 2.5 / , d) 3.5 / .

Fig. 4 .
Fig.4.The blue continuous curve represents the Fisher information for intensity detection with the mode LG 04 as a function of the axial coordinate / R , as in Fig.1.The blue dotted curve is the Fisher information for the same mode, but with a lateral centroid deviation of 10 mrad.The red dots denote the maxima of both curves, and thus the position of the corresponding optimal detection planes.For the misaligned beam, this optimal plane shifts toward the beam waist (located at = 0).The brown curve is the normalized differential gradient of the Fisher information with and without misalignment.

Fig. 5 .Fig. 6 .
Fig. 5. (Left) Scheme of our experimental setup.(Right) a) SLM wavefront used for the GS algorithm.The measured corrected beams appear in b) LG 00 , c) LG 01 , and d) LG 04 .The red circles indicate the radii for maximum intensity and corresponding inner and outer radii with intensity drop to 1/ 2 .