{'fraction_non_alphanumeric': 0.030239301532676927, 'fraction_numerical': 0.01475321033680655, 'mean_word_length': 4.768057342400294, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum-limited estimation of the axial separation of two incoherent point sources

Spatial mode demultiplexing has recently been proposed as a method to overcome the diffraction limits and can be used to precisely measure the transverse separation between two point sources. Here we show that superresolution can also be achieved along the axial direction through the use of a radial mode sorter, i.e. a device that can decompose the incident field into an orthonormal set of radial modes. We demonstrate axial superresolution in the laboratory through use of a binary radial mode sorter, a device that can sort odd-order and even-order radial modes into two separate output ports. We also show theoretically and experimentally that our method allows for saturating the the quantum Fisher information in the limit of a zero axial separation between two incoherent point sources. Our method presents an important step towards three-dimensional superresolution imaging via spatial mode demultiplexing and can be useful to a variety of applications such as far-field fluorescence microscopy.


INTRODUCTION
Optical microscopy is one of the most important imaging modalities and has been broadly applied in various areas. One crucial metric for an optical microscope is the spatial resolution, which is typically constrained by the diffraction limit, and the Rayleigh criterion is proposed as the resolution limit of an incoherent imaging system [1][2][3]. In recent decades, various methods have been proposed to surpass the diffraction limit. In fluorescence microscopy, a widely used approach is to activate each fluorescence molecule individually, and therefore the overlap between neighboring molecules is avoided, and the localization precision can be improved to tens of nanometers [4][5][6]. This technique usually requires specially prepared samples, and the reconstruction of an image can take a long time due to the sophisticated activation and detection of individual fluorophores. Another superresolution technique is based on decomposing the optical field into the linear prolate spheroidal functions, i.e., the eigenfunctions of aperture in a coherent imaging system [7][8][9]. It is shown that the ultimate limit of resolution of a coherent imaging system is not determined by diffraction but by the signal-to-noise ratio of the measurement. Therefore, a sufficiently large number of photons is needed to enable superresolution. In addition, this technique, including other approaches that require nonclassical light sources [10][11][12][13][14][15], cannot be readily applied to incoherent superresolution imaging considered here. While many other methods have been proposed to realize axial super-localization [16], such as an interferometric microscope [17][18][19][20], point-spread function (PSF) engineering [21][22][23][24][25][26], and multi-plane detection [27][28][29][30], these advances can be used to precisely measure the axial location of only a single point source, and it remains a challenge to determine a small axial separation when two incoherent, simultaneously emitting point sources overlap with each other.
To develop an efficient axial superresolution technique, we follow the procedure in Ref. [31] and formulate the estimation of axial separation in the context of quantum metrology [32][33][34][35][36]. The precision of a measurement method is typically quantified by the Fisher information, and the reciprocal of Fisher information is referred to as the Cramér-Rao lower bound (CRLB) and characterizes the lower bound of measurement variance for an unbiased estimator [37,38]. To determine the axial location of point sources, the easiest and most commonly used approach, which we refer to as the direct imaging method, is to measure the size of the PSF in the image plane and then deduce the axial positions accordingly. However, our calculation in the next section shows that the corresponding Fisher information drops to zero when the axial separation of two incoherent point sources gets close to zero. This result is not surprising because the size of PSF changes slowly when point sources are almost on focus. Nonetheless, a further calculation shows that the quantum Fisher information does not vanish for an arbitrarily small axial separation. The quantum Fisher information is the upper limit of the Fisher information that cannot be exceeded by any possible types of measurement as derived in the quantum metrology theory and can be used to quantify the maximum possible amount of information that can be obtained by a measurement. Given the non-vanishing quantum Fisher information, there should exist a type of measurement that can outperform the direct imaging method and extract the maximum possible amount of information from each photon.
In the following sections, we demonstrate both theoretically and experimentally that the axial superresolution can be achieved at the single-photon level by a radial mode sorter. This radial mode sorter can losslessly project the incident photons into the radial Laguerre-Gaussian basis set. With the same amount of photons, our scheme based on the radial mode sorter can estimate the axial separation with smaller bias and standard deviation. We note that similar strategies have been studied for transverse superresolution [39][40][41][42][43][44][45][46][47][48], which are based on a Hermite-Gaussian mode sorter [40,49] or mode parity decomposition [41,42]. However, we emphasize that our radial mode sorter comes from very recent advances in spatial mode sorters [50][51][52] and cannot be simply realized by mode parity decomposition. Furthermore, while homodyne or heterodyne detection [12,39] provides an easier way to implement the spatial mode projective measurement, the shot noise from the reference beam fundamentally deteriorates the signal-to-noise ratio, and it has been shown that homodyne or heterodyne detection cannot even outperform the direct imaging method when the available photon number is small [53]. Moreover, the mode sorter does not require any active components, such as the local oscillator in heterodyne detection, and thus is more favorable in an experiment.

THEORY
The conceptual diagram for the direct imaging method and sorter-based measurement is shown in Fig. 1. The direct imaging method employs an objective to collect photons and then use a tube lens to form an image of the object as shown in Fig. 1(a). Alternatively, one can detect the optical field in a complete and orthonormal basis set as shown in Fig. 1(b), which can be realized by a spatial mode sorter and is referred to as sorted-based measurement. In the following derivations, we use the Dirac notation to represent the fields and assume a coherent state for each point source. While a semiclassical treatment is sufficient to derive these formalisms inspired by quantum metrology [46], the Dirac notation is convenient to denote the mixed state of the incoherent sources and makes it straightforward to extend the theory to other types of light sources such as single-photon state [42] and thermal state [31]. For a more tractable analysis and experiment, here we assume a Gaussian PSF, and the field distribution at the pupil plane for an on-axis point source is denoted by jψi, where hr 0 jψi ψr 0 ; z and where NA is the numerical aperture, z is the axial position of the point source, k 2π∕λ is the wavenumber, λ is the wavelength, r 0 is the normalized radial coordinate in the pupil plane, and jr 0 i is the corresponding radial eigenstate. Here we define r 0 r p ∕f 1 NA, where r p is the radial coordinate in the pupil plane and f 1 is the objective focal length. This pupil plane field distribution can be viewed as a paraxial, Gaussian approximation to the pupil function of a hard-edged circular aperture [54]. For direct imaging, a tube lens is used to perform a Fourier transform to the pupil function, and the intensity distribution on the image plane becomes where r is the radial coordinate in the image plane, wz denotes the Gaussian beam waist width on the image plane, and M is the magnification of the imaging system. By measuring the beam size, we can estimate the axial position z. Similar to the case of transverse superresolution [31], here we assume a priori knowledge of two on-axis, equally bright incoherent point sources with the centroid located at z 0 plane, and the axial separation between them is s. The density matrix of these two point sources at the pupil plane can be written as ρ jψ 1 ihψ 1 j jψ 2 ihψ 2 j∕2, where hr 0 jψ 1 i ψr 0 ; s∕2 and hr 0 jψ 2 i ψr 0 ; −s∕2. The normalized total intensity at the image plane can be calculated as I s r ≡ hrjρjri I r; s∕2 I r; −s∕2∕2, where jri is the radial eigenstate in the image plane, and the image plane is related to the pupil plane by the Fourier transform. For sorter-based measurement, the incident field is decomposed to an orthonormal basis set, and here we consider the radial Laguerre-Gaussian (LG) basis because we notice that the axial position affects only the radial profile of pupil function. The radial LG basis in the pupil plane can be denoted as jLG p i, where hr 0 jLG p i LG p r 0 and where L p · is the Laguerre polynomial. While the twodimensional LG basis involves another azimuthal index l, this Research Article radial subset with l 0 can still form a complete basis to describe the pupil function because of the rotational symmetry of the pupil function as shown in Eq. (1). Decomposing the pupil function ψr 0 ; z to this basis leads to the following radial mode distribution: where z R πw 2 0 ∕λ and w 0 λ∕πNA [31]. For two equally bright sources separated by s, the output radial mode distribution becomes P s p ≡ hLG p jρjLG p i Pp; s∕2 Pp; −s∕2∕2. It can be noticed that for direct imaging and sorter-based measurement, the two incoherent point sources have the same response because Eqs. (2) and (4) are even functions of z, which suggests that the analysis presented here can also be applied to single-point localization.
We next compare the performance of direct imaging and sorter-based measurement by calculating the Fisher information for both techniques. The Fisher information for direct imaging is [ which is independent of the magnification M . The Fisher information for the sorter-based measurement is The quantum Fisher information, i.e., the upper bound of Fisher information of any possible measurements, can be calculated as [55] where j∂ s ψ 1 i ∂jψ 1 i∕∂s, and it can be readily shown that K s 1∕4z 2 R . We also follow the usual way of using the symmetric logarithmic derivative to calculate the quantum Fisher information, and the details are presented in Supplement 1, which give the same result. The reciprocal of quantum Fisher information is the quantum CRLB, which gives the lower bound of classical CRLB for any possible measurements. We notice that the sorter-based measurement can reach the quantum Fisher information when the separation goes to zero, i.e., J sorter 0 K s [see Eq. (6)]; therefore, it can be considered to be an optimal measurement for s close to zero. However, in a realistic experiment, a mode sorter can access only a finite-dimensional Hilbert space. Therefore we follow the procedure in Ref. [56] to construct other possible optimal measurements that can reach the quantum Fisher information in the limit of s 0. In Supplement 1, we show that a binary radial mode sorter is sufficient to access the quantum Fisher information. A binary sorter has two output ports, and all odd-order radial modes are directed to one output port, while all even-order modes are directed to another output port. Therefore, the photon probability distribution at two output ports is Thus, the Fisher information for a binary sorter is The plot of Fisher information for different methods is shown in Fig. 2(a). It can be readily seen that the Fisher information of direct imaging begins to drop when s is smaller than 2z R . In incoherent imaging microscopy, the axial resolution can be expressed as Δz 2λ∕NA 2 [57], which can be rewritten as Δz 2πz R with our notation. We note that the discrepancy between 2z R and 2πz R comes from our assumption of a Gaussian PSF rather than an Airy disk. However, it can be noticed that the sorter-based measurement stays nonzero and achieve the quantum Fisher information when s approaches zero, which makes it possible to break the diffraction limit. To further illustrate the improvement provided by the radial mode sorter, we calculate the Fisher information of astigmatic imaging [21], and the result is presented in Supplement 1. It is shown that while astigmatism can enhance three-dimensional localization precision of a single point source, it cannot be used directly to resolve the axial separation between two simultaneously emitting point sources without the help of photo-switchable fluorophores.
Having analyzed the performance of each method, now we need to establish the estimator of separation. For direct imaging, it can be verified that the maximum likelihood estimator iŝ where r m is the radial coordinate of m-th photon in the image plane, and N is the total detected photon number. The intuition behind this estimator is to measure the Gaussian width w by detecting the radial coordinates of photons and then use the estimated Gaussian width to calculate the separation s. However, the simple estimator ofŵ does not take into account any experimental imperfections such as detector noise or pixelation, and thus may not be robust in a realistic experiment. Therefore, we apply Research Article the algorithm in Ref. [58] to realize a robust, efficient Gaussian width estimatorŵ in our experiment. For binary sorter-based measurement, the maximum likelihood estimator iŝ where m 0 and m 1 are the photon numbers in the two output ports, and m 0 m 1 N is the total detected photon number. The intuition behind this estimator is to use the photon probability distribution at the output ports of the sorter to estimate the separation. The lower bound of the variance of an estimator for N independent measurements is given by [38] Varŝ ≥ where the right-hand side is referred to as the CRLB, and N is the photon number in the context of our experiment given the Poisson statistics. This formula for CRLB is also applicable to other classical photon states such as single-photon state [42] and thermal state [31], and the variance that scales as N −1 is referred to as the standard quantum limit [33,59]. For an unbiased estimator whose expectation is equal to the value of the estimated parameter, i.e., Eŝ s, this CRLB reduces to a simpler form as Varŝ ≥ 1∕N · J s, which is just the reciprocal of Fisher information, as we plot in Fig. 2(b).

EXPERIMENT
A schematic for the experimental setup is shown in Fig. 3. We use an attenuated laser source to illuminate the spatial light modulator (SLM) to generate the Gaussian pupil function produced by a point source. An acousto-optic modulator (AOM) is driven by a signal generator to produce 3 μs pulses, and the driving signal is also connected to an intensified charge coupled device (ICCD, PI-Max 4 1024i) for synchronization. The average detected photon number in each pulse is around 2000. We use the calibration factor provided by the manufacturer to calculate the photon number in each pixel of the camera. We emulate two incoherent sources by mixing the data for z s∕2, which is generated by SLM separately. A computer-generated hologram is displayed on SLM 1 to generate the desired field at the first diffraction order [60]. Each time, the SLM displays the corresponding hologram to generate either ψr 0 ; s∕2 or ψr 0 ; −s∕2 to simulate a point source located at z s∕2 or z −s∕2, respectively. Since both holograms are never present at the same time, there is no coherence between the two simulated point sources. By using a long exposure time of the camera to incoherently mix the data, we effectively generate two incoherent, simultaneously emitting point sources [61]. For the Gaussian pupil function, we use the parameters of NA 0.1 and f 1 4 mm. The calibration data of SLM 1 are presented in Supplement 1.
To construct a binary radial mode sorter, we use two polarization-sensitive SLMs (Hamamatsu X10468-02) as shown in the schematic [50][51][52]. In our experiment, we use two different areas on a single SLM to act as two SLMs for reduced experimental complexity. Due to the polarization sensitivity of the SLM, this binary mode sorter is designed to work for diagonally polarized light and cannot be directly used for arbitrary polarization [52]; therefore, we use a polarizer before the SLM to filter out undesired polarization. We note that the polarization of photons is not relevant to the theory of superresolution, and thus the use of a polarization-sensitive sorter is permissible for this proof-ofprinciple experiment. To realize a polarization-independent sorter, one can use the previously reported interferometric scheme [50]. A quadratic phase pattern is imprinted on SLM 2 and SLM 3 as the essential ingredient of the sorter. This quadratic phase is identical to the phase of a spherical lens with a focal length of 46.5 cm, and the separation between two SLMs is 65.8 cm. Each SLM performs a fractional Fourier transform of order α π∕2 to horizontally polarized light and α π∕4 to vertically polarized light, respectively. One can check that even-order radial modes remain diagonally polarized and odd-order radial modes become anti-diagonally polarized after passing through both SLMs [50]. Through the use of a half-wave plate (HWP) and a polarizing beamsplitter (PBS), one can efficiently separate odd-and even-order radial modes to distinct output ports. More details about the principle of the radial mode sorter can be found in Refs. [50][51][52]. As mentioned earlier, this radial mode sorter cannot be realized by mode parity decomposition based on mirror reflection [41] or a 4-f system. Moreover, this radial mode sorter is in principle lossless, and the loss of our sorter comes mainly from the limited light utilization efficiency of the SLMs, which can be reduced by using other low-loss devices such as commercially available polarization directed flat lenses [51]. In our experiment, we direct the photons from different output ports to different areas of an ICCD. For direct imaging, we use a 10 cm tube lens to form the image on the ICCD detector plane. For each separation, we repeat the experiment 400 times and calculate the expectation and standard deviation from the collected data based on the maximum likelihood estimators.  Fig. 3. Schematic of the experimental setup. A 633 nm He-Ne laser is attenuated and modulated by an acousto-optic modulator (AOM) to generate weak pulses. A computer-generated hologram is imprinted onto spatial light modulator (SLM 1) to generate the desired pupil function to simulate point sources. Two different methods, the binary sorter-based measurement and the direct imaging method, are used to estimate the separation s. In our experiment, we use a flip mirror to select the measurement method.

DISCUSSION
The measured separation and the standard deviation as a function of the real separation for different measurement methods are presented in Fig. 4. The average detected photon number for each measurement in our experiment is around 2000. As can be seen in Eq. (12), the loss of photons will decrease the detected photon number N and thus increase the variance of measurement and deteriorate the measurement precision. However, since all devices used in our experiment are essentially phase-only elements, the loss can always be reduced to zero by using appropriate antireflection coatings, and in our analysis, we assume a detection efficiency of unity. The Monte Carlo simulation results are provided as comparisons, and they agree well with the experimental data. In the simulation we set the detected photon number to be 2000, and the expectation and standard deviation of both estimators are retrieved by averaging 4000 simulations. We assume a noiseless detector with a sufficiently high spatial resolution in the simulation, and the estimators for direct imaging and binary sorter-based measurement are given by Eq. (10) and Eq. (11), respectively. One immediate observation in Fig. 4(a) is that the measured separation of direct imaging deviates from the real value when the real separation is close to zero. Another observation in Fig. 4(b) is that the measured standard deviation does not follow the CRLB and stays finite in the vicinity of s 0. Neither observation is due to experimental imperfections, as they agree with the Monte Carlo simulation and should be attributed to the bias of the estimator [38,42]. The bias of an estimator is defined as the difference between the estimator's expectation value and the real value of the parameter being estimated. In Supplement 1, we provide a detailed, analytical calculation of the bias ofŝ direct . The expectation value ofŝ direct at s 0 can be well approximated as which is 0.123z R for N 2000 and very close to the Monte Carlo simulation 0.124z R , as shown in Fig. 4(a). It can be noticed that this bias is on the order of z R when N is small, which qualitatively agrees with the conventional axial resolution of 2πz R . A large photon number N can lower the value of bias, which corresponds to the fact that a higher signal-to-noise ratio can enhance the resolution of direct imaging. A simple example is the deconvolution algorithm, which can be used to obtain subdiffraction resolution as long as a sufficiently high signal-to-noise ratio is available. However, the bias of direct imaging scales rather slowly with N as N −1∕4 , and to reduce this bias, a sufficiently large N is needed. While a large photon number is attainable with a bright light source, in a photon-starving experiment such as fluorescence microscopy, it is usually not achievable. The slope of the estimator's expectation is calculated to be Together with Eq. (12), it immediately follows that the CRLB be- . In contrast to the diverging CRLB solely predicted by the reciprocal of Fisher information, the CRLB calculated here takes into account the bias and explains the non-diverging standard deviation as shown in the experiment and Monte Carlo simulation. The scaled CRLB is calculated to be ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Varŝ direct p ∕z R ∕ ffiffiffiffi ffi N p ≥ 5.8, which is close to the standard deviation in the Monte Carlo simulation result 6.6 as shown in Fig. 4(b). It should be noted that Eq. (12) is an inequality instead of an equality, which causes the discrepancy between 5.8 and 6.6. Furthermore, ∂Eŝ direct ∕∂sj s0 0 implies that the expectation value has a slope of zero when s 0, as can be seen in Fig. 4(a). Hence, despite of a finite standard deviation, it is intrinsically unrealistic to use the measured s to recover the real s in the vicinity of s 0 for the direct imaging method, and any attempt to construct an unbiased estimator will lead to a diverging standard deviation. Another observation is that the variance of the estimator for direct imaging scales as Varŝ direct j s0 ∝ N −1∕2 ; therefore, this estimator cannot reach the standard quantum limit when s is small [33,59].
For sorter-based measurement, it can also be noticed that the standard deviation in both simulation and experiment deviates from the reciprocal of Fisher information and drops to zero when s is small. In Supplement 1, we show that ∂Eŝ binary ∕∂sj s0 0 and thus Varŝ binary j s0 ≥ 0, which explains the zero standard deviation that violates the reciprocal of Fisher information as shown in Fig. 4(d). It has been pointed out that this so-called superefficiency exists only on a set of points with zero measure, and the region of superefficiency reduces for more samples [31,62]. In addition, we have also shown in Supplement 1 that Eŝ binary j s0 0, which coincides with the Monte Carlo simulation and suggests that the sorter-based measurement can provide more precise, less biased measurement when s is small. However, we still observe a small, nonzero separation at s 0 in our experiment, and the zero standard deviation is not visible either. We attribute this inconsistency to experimental imperfections, including dark noise of the detector and misalignment of the sorter. At the point of s 0, all photons are supposed to be sorted to the output port of even-order radial modes, and no photons should be detected at the other output port. Nevertheless, when we experimentally characterize our sorter, we observe that 0.28% of detected photons are routed to the wrong output port on average when s 0. In the data processing, we have subtracted this averaged crosstalk before estimating the separation, but the Research Article associated shot noise cannot be simply eliminated and thus leads to the experimental inconsistency, as we describe above. In Supplement 1, we quantify the effect of crosstalk on data processing, and the analytical calculation shows that Eŝ binary j s0 0.043z R , which is very close to the measured value of 0.049z R . There are several ways to further mitigate the crosstalk, such as aligning the sorter more carefully, replacing the ICCD by low-noise singlepixel detectors, and developing a more robust estimator [58]. Despite these experimental imperfections, it is apparent that sorterbased measurement can outperform direct imaging, given the strong bias and higher standard deviation of direct imaging compared to that of sorter-based measurement. In addition, the advantage of sorter-based measurement is supposed to be more obvious with a larger photon number, because the variance of direct imaging scales as N −1∕2 rather than N −1 .
To further provide a quantitative description of the improvement of our experiment, we compare the two methods in terms of bias and variance. For direct imaging, if we want to reduce its bias to 0.049z R , which is obtained in sorter-based measurement, the photon number needs to be on the order of 10 5 to satisfy 0.82N −1∕4 z R 0.049z R . Furthermore, in Fig. 4, it can be seen that the measured standard deviation of direct imaging is approximately twice that of sorter-based measurement. Remember that the variance of direct imaging scales as N −1∕2 ; thus, the standard deviation scales as N −1∕4 , and 16 times more photons are needed to reduce the standard deviation of direct imaging to the level of sorter-based measurement. We note that here we are comparing the experimental data rather than noise-free theoretical predictions. The Monte Carlo simulation shows that the sorter-based measurement has zero bias and zero standard deviation at s 0, and thus the direct imaging needs infinite photons to beat the sorter-based measurement.
In this work, we focus mainly on the superresolution of axial separation for two point sources, but we note that the theory presented above can be directly applied to the localization of axial position z of a single point source as long as the separation s is replaced by s → z∕2, which presents an alternative to sophisticated interferometric microscopy [17][18][19][20]. In contrast to the interferometric detection scheme, which requires nanometer-scale stabilization over a path length on the order of 1 m [16], the common-path radial mode sorter used in our experiment is robust, and no additional stabilization control is needed. We assume that the two point sources are on axis and their center position is known. In Supplement 1, we analyze the effect of misaligned centroid, i.e., centroid of point source pair z C ≠ 0. As shown in Fig. 5(a), while the Fisher information drops in the presence of a misaligned centroid, the radial mode sorter can provide improved precision for a small separation with jz C j < 0.3z R . In a realistic scenario, an adaptive measurement can be used to estimate both the centroid and separation, as discussed in Ref. [31]. However, unlike the case of transverse centroid estimation, the direct imaging does not provide sufficient Fisher information for measuring the axial centroid of a point source pair. Here we notice that astigmatic imaging [21] presents an effective method to overcome this difficulty. The analysis is included in Supplement 1, and the Fisher information of centroid estimation for astigmatic imaging is shown in Fig. 5(b). It can be seen that astigmatic imaging provides appreciable Fisher information over a broad range of centroid and separation. Hence, a hybrid measurement consisting of both the radial mode sorter and astigmatic imaging can be a practical scheme for axial superresolution. Another assumption we make is that a Gaussian PSF is used for more tractable analysis and experiment. While the Gaussian PSF is a widely adopted approximation [57,63], for a high-NA imaging system, a more accurate PSF model may be needed [64]. In this case, one can always establish a complete and orthonormal basis based on the PSF model and construct a sorter accordingly to achieve superresolution accordingly [46]. Very recently, it has been pointed out theoretically that for the pupil function of a hardedged aperture, the optimal measurement basis turns out to be the Zernike basis [55,56], and we discuss other optimal measurements that are easier to implement in Supplement 1. Therefore, based on our result, it can be anticipated that three-dimensional superresolution can be realized as long as a Zernike mode sorter is available. In Supplement 1, we calculate the Fisher information of various measurements for an Airy-disk-shaped PSF, and the result is shown in Fig. 6. It can be seen that while the Zernike mode sorter provides optimal measurement [55], the LG mode sorter as a sub-optimal measurement can still provide nonzero Fisher information at near-zero separation, outperforming the direct imaging measurement. Here we take into account only the zerothand first-order radial modes in calculating the Fisher information of the LG mode sorter and Zernike mode sorter, which should be reasonably achievable in an experiment. Moreover, given the widely used Gaussian-to-tophat laser beam shaper [65,66], it is possible to convert the pupil function of a hard-edged aperture to a Gaussian and then apply the radial mode sorter subsequently. Recent advances in multi-plane light conversion [67] also provide an alternative possible solution for building a Zernike mode sorter. Finally, despite the classical light source used in our experiment, our method can be used for other light sources such as single-photon emitters [31,42], because the quantum state of

Research Article
photons represents the temporal coherence of light and is generally independent of the spatial degree of freedom. Therefore, it is possible to combine the radial mode sorter and intensity correlation measurement to further increase the resolution for single-photon sources [68].
In conclusion, we theoretically and experimentally demonstrate the axial superresolution based on a radial mode sorter. The binary radial mode sorter employed in out experiment can reach the quantum CRLB for arbitrarily small axial separation. Our method makes three-dimensional superresolution imaging promising and can be potentially useful for enhancing the resolution of optical microscopes.