\section{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 \cite{rayleigh1879xxxi,goodman2008introduction,born2013principles}. 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 \cite{betzig2006imaging, rust2006sub,hell1994breaking}. 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 \cite{beskrovnyy2005quantum,kolobov2007quantum,kolobov2000quantum}. 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 are needed to enable the superresolution. In addition, this technique, including other approaches that require nonclassical light sources \cite{rozema2014scalable,treps2002surpassing,taylor2014subdiffraction,taylor2013biological,shin2011quantum,schwartz2013superresolution}, cannot be readily applied to incoherent superresolution imaging considered here. While many other methods have been proposed to realize axial super-localization \cite{von2017three}, such as interferometric microscope \cite{backlund2018fundamental,shtengel2009interferometric,schrader19964pi,bewersdorf2006comparison}, point spread function (PSF) engineering \cite{huang2008three, pavani2009three, jia2014isotropic,martinez1995tunable,sales1997fundamental,tamburini2006overcoming}, and multi-plane detection \cite{juette2008three,toprak2007three,dalgarno2010multiplane,abrahamsson2012fast}, these advances can be only used to precisely measure the axial location of 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. \cite{tsang2016quantum} and formulate the estimation of axial separation in the context of quantum metrology \cite{helstrom1969quantum,giovannetti2011advances,holevo2011probabilistic,degen2017quantum,paris2009quantum}. 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\'er-Rao lower bound (CRLB) and characterizes the lower bound of measurement variance for an unbiased estimator \cite{van2004detection,kay1993fundamentals}. To determine the axial location of point sources, the easiest and most commonly used approach, which we refer to as 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 a 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 \cite{yang2016far,paur2016achieving,nair2016far,tham2017beating,tang2016fault, vrehavcek2018optimal,vrehavcek2017multiparameter, tsang2018subdiffraction, tsang2017subdiffractionNJP, zhou2018modern}, which are based on a Hermite-Gaussian mode sorter \cite{zhou2018hermite,paur2016achieving} or mode parity decomposition \cite{nair2016far,tham2017beating}. However, we emphasize that our radial mode sorter comes from very recent advances in spatial mode sorter \cite{zhou2017sorting,zhou2018realization,zhou2018quantum} and cannot be simply realized by mode parity decomposition. Furthermore, while homodyne or heterodyne detection \cite{yang2016far,taylor2014subdiffraction} provides an easier way to implement the spatial mode projective measurement, the shot noise from reference beam fundamentally deteriorates the signal-to-noise ratio, and it has been shown that homodyne or heterodyne detection cannot even outperform direct imaging method when the available photon number is small \cite{yang2017fisher}. 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.

\section{Theory}

\begin{figure}[t]
\includegraphics[width=0.9\linewidth]{setup.pdf}
\caption{Conceptual diagram for (a) direct imaging and (b) sorter-based measurement. A spatial mode sorter can direct different spatial mode components to different locations to perform spatial mode demultiplexing.}
\label{fig:setup}
\end{figure}


The conceptual diagram for direct imaging method and sorter-based measurement is shown in Fig.~\ref{fig:setup}. 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.~\ref{fig:setup}(a). Alternatively, one can detect the optical field in a complete and orthonormal basis set as shown in Fig.~\ref{fig:setup}(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 \cite{tsang2018subdiffraction}, 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 \cite{tham2017beating} and thermal state \cite{tsang2016quantum}. 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 $\ket{\psi}$, where $\braket{r_0}{\psi}=\psi (r_0;z) $ and
\begin{equation}\label{Eq:pupil}
\begin{aligned}
\psi (r_0;z) =\sqrt{2/\pi}\exp(-r_0^2)  \exp(-ikz\text{NA}^2 r_0^2/2),
\end{aligned}
\end{equation}
where NA is the numerical aperture, $z$ is the axial position of the point source, $k=2\pi/\lambda$ is the wavenumber, $\lambda$ is the wavelength, $r_0$ is the normalized radial coordinate in the pupil plane, and $\ket{r_0}$ is the corresponding radial eigenstate. Here we define $r_0=r_p/(f_1 \text{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 \cite{petrov2017measurement}. 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
\begin{equation}
\begin{aligned}\label{Eq:DirectIntensity}
I(r;z) &= \frac{2}{\pi}\frac{1}{w^2(z)}\exp(-\frac{2r^2}{w^2(z)}), \\
w(z) &= \frac{M\lambda}{\pi \text{NA}} \sqrt{ 1+(z\pi \text{NA}^2/\lambda)^2},
\end{aligned}
\end{equation}
where $r$ is the radial coordinate in the image plane, $w(z)$ 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 \cite{tsang2016quantum}, here we assume \textit{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 $\rho =(\ket{\psi_1} \bra{\psi_1}+\ket{\psi_2} \bra{\psi_2})/2$, where $\braket{r_0}{\psi_1}=\psi (r_0;s/2) $ and $\braket{r_0}{\psi_2}=\psi (r_0;-s/2) $. The normalized total intensity at the image plane can be calculated as $I_s(r)\equiv \bra{r}\rho\ket{r}=[I(r;s/2)+I(r;-s/2)]/2$, where $\ket{r}$ 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 only affects the radial profile of pupil function. The radial LG basis in the pupil plane can be denoted as $\ket{\text{LG}_p}$, where $\braket{r_0}{\text{LG}_p}=\text{LG}_p(r_0)$ and
\begin{equation}\label{Eq:LGpupil}
\begin{aligned}
&\text{LG}_p(r_0)=\sqrt{ 2/\pi }  \exp(-r_0^2) L_p( 2r_0^2 ),\\
\end{aligned}
\end{equation}
where $L_p(\cdot)$ is the Laguerre polynomial. While the two-dimensional LG basis involves another azimuthal index $\ell$, this radial subset with $\ell=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.~(\ref{Eq:pupil}). Decomposing the pupil function $\psi (r_0;z)$ to this basis leads to the following radial mode distribution:
\begin{equation}
\begin{aligned}
P(p;z)&=|\braket{\psi}{\text{LG}_p}|^2 = \frac{    4z_R^2 z^{2p}    }{   (4z_R^2+z^2)^{p+1}  },
\end{aligned}
\label{Eq:sorterDistribution}
\end{equation}
where $z_R=\pi w_0^2/\lambda$ and $w_0=\lambda /\pi \text{NA}$ \cite{tsang2016quantum}. For two equally bright sources separated by $s$, the output radial mode distribution becomes $P_s(p) \equiv \bra{\text{LG}_p}\rho\ket{\text{LG}_p}=[P(p;s/2)+P(p;-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 Eq.~(\ref{Eq:DirectIntensity}) and Eq.~(\ref{Eq:sorterDistribution}) 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 the sorter-based measurement by calculating the Fisher information for both techniques. The Fisher information for direct imaging is \cite{tsang2016quantum}
\begin{equation}
\begin{aligned}
\mathcal{J_{\text{direct}}}(s) &= \int_0^{2\pi}d\phi  \int_0^{+\infty}\frac{1}{I_s(r)}\left( \frac{\partial I_s(r)}{\partial s}  \right)^2 rdr \\
&=\frac{4s^2}{(s^2+4z_R^2)^2},
\end{aligned}
\end{equation}
which is independent of the magnification $M$. The Fisher information for the sorter-based measurement is
\begin{equation}
\begin{aligned}
\mathcal{J_{\text{sorter}}}(s) = \sum_{p=0}^{\infty} \frac{1}{P_s\text (p) } \left(  \frac{\partial P_s (p) }{\partial  s}   \right)^2=\frac{4}{s^2+16z_R^2}.
\label{Eq:Csorter}
\end{aligned}
\end{equation}

\begin{figure}[t]
\includegraphics[width=\linewidth]{FisherInfo.pdf}
\caption{(a) The Fisher information (FI) as a function of axial separation for different methods. The sorter and binary sorter can reach the quantum Fisher information for small separation $s$, while the Fisher information of direct imaging drops to zero. (b) The normalized square root of Cram\'er-Rao lower bound (CRLB) for different methods. $N$ is the detected photon number.}
\label{fig:FisherInfo}
\end{figure}

\begin{figure*}[t]
\center
\includegraphics[width=  0.9\textwidth]{Schematic.pdf}
\caption{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.}
\label{fig:Schematic}
\end{figure*}



The quantum Fisher information, i.e. the upper bound of Fisher information of any possible measurements, can be calculated as \cite{Zhixian2018Quantum}
\begin{equation}
\begin{aligned}
\mathcal{K}_{s}=4[\braket{\partial_{s}\psi_1}{\partial_{s}\psi_1}-|\braket{\psi_1}{\partial_{s}\psi_1}|^{2}],
\end{aligned}
\end{equation}
where $\ket{\partial_s \psi_1} = \partial \ket{\psi_1} /\partial s$ and it can be readily shown that $\mathcal{K}_{s}=1/4z_R^2$. We also follow the usual way of using the symmetric logarithmic derivative to calculate the quantum Fisher information and the details are presented in Supplementary Section 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. $\mathcal{J_{\text{sorter}}}(0)= \mathcal{K}_{s} $ [see Eq.~(\ref{Eq:Csorter})], therefore it can be considered to be an optimal measurement for $s$ close to zero. However, in a realistic experiment, a mode sorter can only access a finite-dimensional Hilbert space. Therefore we follow the procedure in \cite{yang2018optimal} to construct other possible optimal measurements that can reach the quantum Fisher information in the limit of $s=0$. In Supplementary Section 2 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
\begin{equation}
\begin{aligned}
P_s^{0}(s)  &= \sum_{p=0}^{\infty} P_s(2p;s)=\frac{1}{2}+\frac{4z_R^2}{8z_R^2+s^2}, \\
P_s^{1}(s)  &= \sum_{p=0}^{\infty} P_s(2p+1;s)=\frac{1}{2}-\frac{4z_R^2}{8z_R^2+s^2}.
\end{aligned}
\end{equation}

Therefore the Fisher information for a binary sorter is
\begin{equation}
\begin{aligned}
\mathcal{J_{\text{binary}}}(s) &= \sum_{q=0}^{1} \frac{1}{P_s^q  (s) } \left(  \frac{\partial P_s^q (s) }{\partial  s}   \right)^2=\\
&=\frac{256z_R^4}{(s^2+8z_R^2)^2(s^2+16z_R^2)}.
\end{aligned}
\end{equation}
The plot of Fisher information for different methods is shown in Fig.~\ref{fig:FisherInfo}(a). It can be readily seen that the Fisher information of direct imaging begins to drop when $s$ is smaller than $2z_R$. In an incoherent imaging microscopy, the axial resolution can be expressed as $\Delta z=2\lambda/\text{NA}^2$ \cite{handbook2006JB}, which can be rewritten as $\Delta z=2\pi z_R$ with our notation. We note that the discrepancy between $2z_R$ and $2\pi 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 \cite{huang2008three} and the result is presented in Supplementary Section 3. 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 is
\begin{equation}
\begin{aligned}
\hat{w}=\sqrt{\frac{2}{N} \sum_{m=1}^{N}r_m^2}, \quad \quad \hat{s}_{\text{direct}}=2z_R \sqrt{\frac{\hat{w}^2}{w_0^2}-1}.
\label{Eq:EstDirect}
\end{aligned}
\end{equation}
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 $\hat{w}$ 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 the algorithm in \cite{smith2010fast} to realize a robust, efficient Gaussian width estimator $\hat{w}$ in our experiment. For binary sorter-based measurement, the maximum likelihood estimator is
\begin{equation}
\begin{aligned}
\hat{Q}=\frac{1}{N}\sum_{q=0}^{1} q\cdot m_q=\frac{m_1}{N}, \quad \hat{s}_{\text{binary}}=2 z_R \sqrt{\frac{2}{1-2\hat{Q}} -2},
\label{Eq:EstBinary}
\end{aligned}
\end{equation}
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 \cite{kay1993fundamentals}
\begin{equation}
\begin{aligned}
\text{Var}(\hat{s}) \geq \frac{(\partial E[\hat{s}]/\partial s)^2}{N\cdot \mathcal{J}(s)},
\label{Eq:CRLB_biased}
\end{aligned}
\end{equation}
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 \cite{tham2017beating} and thermal state \cite{tsang2016quantum}, and the variance that scales as $N^{-1}$ is referred to as the standard quantum limit \cite{giovannetti2004quantum,giovannetti2011advances}. For an unbiased estimator whose expectation is equal to the value of the estimated parameter, i.e. $ E[\hat{s}]=s$, this CRLB reduces to a simpler form as $\text{Var}(\hat{s}) \geq  1/ [N\cdot \mathcal{J}(s)]$, which is just the reciprocal of Fisher information as we plot in Fig.~\ref{fig:FisherInfo}(b).


\section{Experiment}
A schematic for experimental setup is shown in Fig.~\ref{fig:Schematic}. 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 $\mu$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=\pm s/2$ that is generated by SLM separately. A computer-generated hologram is displayed on SLM~1 to generate the desired field at the first diffraction order \cite{mirhosseini2013rapid}. Each time the SLM displays the corresponding hologram to generate either $\psi (r_0;s/2)$ or $\psi (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 \cite{rodenburg2014experimental}. For the Gaussian pupil function we use the parameters of $\text{NA}=0.1$ and $f_1=4$ mm. The calibration data of SLM~1 is presented in Supplementary Section 4.


To construct a binary radial mode sorter we use two polarization-sensitive SLMs (Hamamatsu X10468-02) as shown in the schematic \cite{zhou2017sorting, zhou2018quantum,zhou2018realization}. 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 an arbitrary polarization \cite{zhou2018quantum}, 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-of-principle experiment. To realize a polarization-independent sorter one can use the previously reported interferometric scheme \cite{zhou2017sorting}. 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 $\alpha=\pi/2$ to horizontally polarized light and $\alpha=\pi/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 \cite{zhou2017sorting}. 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 radial mode sorter can be found in \cite{zhou2017sorting, zhou2018quantum,zhou2018realization}. As mentioned earlier, this radial mode sorter cannot be realized by mode parity decomposition based on mirror reflection \cite{nair2016far} or a 4-$f$ system. Moreover, this radial mode sorter is in principle lossless, and the loss of our sorter mainly comes 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 \cite{zhou2018realization}. 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.






\section{Discussion}

The measured separation and the standard deviation as a function of the real separation for different measurement methods are presented in Fig.~\ref{fig:FisherInfoData}. The average detected photon number for each measurement in our experiment is around 2000. As can be seen in Eq.~(\ref{Eq:CRLB_biased}) 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 anti-reflection 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.~(\ref{Eq:EstDirect}) and Eq.~(\ref{Eq:EstBinary}) respectively. One immediate observation from Fig.~\ref{fig:FisherInfoData}(a) is that the measured separation of direct imaging deviates from the real value when the real separation is close to zero. Another observation from Fig.~\ref{fig:FisherInfoData}(b) is that the measured standard deviation does not follow the CRLB and stays finite in the vicinity of $s=0$. Both observations are not due to experimental imperfections as they agree with the Monte Carlo simulation and should be attributed to the bias of the estimator \cite{tham2017beating,kay1993fundamentals}. 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 supplementary section 5 we provide a detailed, analytical calculation of the bias of $\hat{s}_{\text{direct}}$. The expectation value of $\hat{s}_{\text{direct}}$ at $s=0$ can be well approximated as
\begin{equation}
\begin{aligned}
E[\hat{s}_{\text{direct}}]|_{s=0}\approx 0.82N^{-1/4}z_R,
\end{aligned}
\end{equation}
which is $0.123z_R$ for $N=2000$ and very close to the Monte Carlo simulation $0.124z_R$ as shown in Fig.~\ref{fig:FisherInfoData}(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\pi z_R$. A large photon number $N$ can lower the value of bias, which corresponds to the fact that 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 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 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
\begin{figure}[t]
\includegraphics[width=\linewidth]{FisherInfoData.pdf}
\caption{(a) Measured separation and (b) the standard deviation (SD) of $s$ as a function of actual separation for direct imaging method. (c) Measured separation and (d) the SD of $s$ as a function of actual separation for binary sorter-based measurement. The Monte Carlo simulation results and the square root of corresponding CRLB are provided as comparisons. $N$ is the detected photon number.}
\label{fig:FisherInfoData}
\end{figure}
\begin{equation}
\begin{aligned}
\frac{\partial  E[\hat{s}_{\text{direct}}]}{\partial s}\bigg\rvert_{s=0} \approx \frac{0.43N^{1/4}s}{z_R}.
\end{aligned}
\end{equation}
Together with Eq.~(\ref{Eq:CRLB_biased}) it immediately follows that the CRLB becomes $\text{Var}(\hat{s}_{\text{direct}})|_{s=0}\geq 0.74z_R^2/\sqrt{N}$. 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 $\sqrt{\text{Var}(\hat{s}_{\text{direct}})}/(z_R/\sqrt{N}) \geq 5.8$, which is close to the standard deviation in the Monte Carlo simulation result 6.6 as shown in Fig.~\ref{fig:FisherInfoData}(b). It should be noted that Eq.~(\ref{Eq:CRLB_biased}) is an inequality instead of an equality, which explains the discrepancy between 5.8 and 6.6. Furthermore, $ \partial  E[\hat{s}_{\text{direct}}]/\partial s|_{s=0} =0$ implies that the expectation value has a slope of zero when $s=0$ as can be seen in Fig.~\ref{fig:FisherInfoData}(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 $\text{Var}(\hat{s}_{\text{direct}})|_{s=0}\propto N^{-1/2}$, therefore this estimator cannot reach the standard quantum limit when $s$ is small \cite{giovannetti2004quantum,giovannetti2011advances}.



For the sorter-based measurement, it can also be noticed that the standard deviation deviates from the reciprocal of Fisher information and drops to zero when $s$ is small. In the supplement we show that $\partial E[\hat{s}_{\text{binary}}]/\partial s|_{s=0}=0$ and thus $\text{Var}(\hat{s}_{\text{binary}})|_{s=0}\geq 0$, which explains the zero standard deviation that violates the reciprocal of Fisher information as shown in Fig.~\ref{fig:FisherInfoData}(d). It has been pointed out that this so called superefficiency only exists on a set of points with zero measure and the region of superefficiency reduces for more samples \cite{tsang2016quantum, hinkley1979theoretical}. In addition, we have also shown in the supplement that $ E[\hat{s}_{\text{binary}}] |_{s=0}=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 inconsistence 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 associated shot noise cannot be simply eliminated and thus leads to the experimental inconsistence as we describe above. In the supplement we quantify the effect of crosstalk on data processing and the analytical calculation shows that $ E[\hat{s}_{\text{binary}}]|_{s=0}=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 single-pixel detectors, and developing a more robust estimator \cite{smith2010fast}. Despite of these experimental imperfections, it is apparent that the sorter-based measurement can outperform the direct imaging, given the strong bias and higher standard deviation of direct imaging compared to that of the 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}$.



\begin{figure}[t]
\centering
\includegraphics[width= \linewidth]{Centroid.pdf}
\caption{(a) The Fisher information of separation estimation for sorter-based measurement with different centroid positions. The Fisher information for direct imaging with point source pair centroid $z_C=0$ is plotted as a reference. (b) The Fisher information of centroid estimation for astigmatic imaging with different separations.}
\label{fig:Center}
\end{figure}

\begin{figure}[b]
\centering
\includegraphics[width=0.7 \linewidth]{Circ_GaussSorter.pdf}
\caption{The Fisher information of different measurements for an Airy-disk-shaped PSF model.}
\label{fig:Circ}
\end{figure}

To further provide a quantitative description of the improvement of our experiment, we compare the two methods in terms of bias and variance. For the direct imaging, if we want to reduce its bias to $0.049z_R$ that 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, from Fig.~\ref{fig:FisherInfoData} 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 mainly focus 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\rightarrow z/2$, which presents an alternative to the sophisticated interferometric microscopy \cite{backlund2018fundamental,shtengel2009interferometric,schrader19964pi,bewersdorf2006comparison}. In contrast to the interferometric detection scheme which requires nanometer-scale stabilization over a path length on the order of 1 m \cite{von2017three}, 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 Supplementary Section 3 we analyze the effect of misaligned centroid, $i.e.$, centroid of point source pair $z_C \neq 0$. As shown in Fig.~\ref{fig:Center}(a), while the Fisher information drops in the presence of misaligned centroid, the radial mode sorter can provide improved precision for a small separation with $|z_C|<0.3z_R$. In a realistic scenario, an adaptive measurement can be used to estimate both the centroid and separation as discussed in \cite{tsang2016quantum}. However, unlike the case of transverse centroid estimation, the direct imaging does not provide sufficient Fisher information for measuring the axial centroid of point source pair. Here we notice that astigmatic imaging \cite{huang2008three} presents an effective method to overcome this difficulty. The analysis is included in Supplementary Section 3 and the Fisher information of centroid estimation for astigmatic imaging is shown in Fig.~\ref{fig:Center}(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 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 \cite{handbook2006JB, zhang2007gaussian}, for a high-NA imaging system a more accurate PSF model may be needed \cite{gibson1992experimental}. 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 \cite{tsang2018subdiffraction}. Very recently, it has been pointed out theoretically that for the pupil function of a hard-edged aperture, the optimal measurement basis turns to be the Zernike basis \cite{Zhixian2018Quantum, yang2018optimal}, and we discuss other optimal measurements that are easier to implement in Supplementary Section 2. 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 Supplementary Section 6 we calculate the Fisher information of various measurements for an Airy-disk-shaped PSF and the result is shown in Fig.~\ref{fig:Circ}. It can be seen that while the Zernike mode sorter provides the optimal measurement \cite{Zhixian2018Quantum}, 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 only take into account the 0th and 1st order radial modes in calculating the Fisher information of 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 \cite{cheng2013compact, homburg2012gaussian}, 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 \cite{labroille2014efficient} also provide an alternative possible solution for building a Zernike mode sorters. Finally, despite the classical light source used in our experiment, our method can be used for other light sources such as single-photon emitters \cite{tsang2016quantum, tham2017beating}, because the quantum state of photons represents the temporal coherence of light and is generally independent from 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-photons sources \cite{tenne2019super}.



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 Cram\'er-Rao lower bound for an arbitrarily small axial separation. Our method makes three-dimensional superresolution imaging promising and can be potentially useful for enhancing the resolution of optical microscopes.



\section*{Funding}
U.S. Office of Naval Research (No. N000141712443 and N000141512635).

\section*{Acknowledgement}
R. W. B. acknowledges support from Canada Excellence Research Chairs Program and Natural Science and Engineering Research Council of Canada. J.~Y. and A.~N.~J. acknowledge support from National Science Foundation (No. DMR-1506081) and U.S. Army Research Office (No. W911NF-15-1-0496). We thank Mankei Tsang for helpful discussions.
\\
\\
See Supplement 1 for supporting content.

