Axial superlocalization with vortex beams

Discerning the nanometer-scale details of living cells, tissues, and materials is of utmost importance for many modern research efforts. A trail toward to this holy grail was opened with the appearance of a set of methods, dubbed with the generic term of superresolution microscopy [1, 2], capable of bypassing the diffraction limit [3–5]: a barrier that was traditionally thought to be impenetrable.

A number of these techniques can also reveal three-dimensional (3D) structural details: relevant examples include stimulated-emission-depletion microscopy [6], PSF engineering [8–12], photoactivated-localization microscopy [7], and multiplane detection [13–15], to cite only but a few. All of them rely on a very accurate localization of point sources; they differ in how point objects are excited and how the corresponding emitted photons are collected.

For 3D imaging, the emitter is fluorescently labeled and determining its axial position is an indispensable ingredient. This problem has been throughly examined and some impressive results have been demonstrated so far [16]. However, the fundamental depth precision achievable by any such engineering method has been considered only lately [17–19]. The rationale behind is to make a systematic use of the quantum Fisher information (QFI) [20] and the associated quantum Cramér–Rao bound (QCRB) to get a measurement-independent limit [21, 22]. This is much along the lines of the work of Tsang and co-workers to quantify the transverse two-point resolution [23–27], which has led to the dispelling of the Rayleigh curse [28–31].

In a recent work [32], the ultimate precision in axial localization using Gaussian beams has been established. This limit can be attained with just one intensity scan, as long as the detection plane is placed at one optimal position.

In this paper, we generalize those results and derive quantum limits for axial localization with Laguerre–Gauss (LG) beams, which carry quantized orbital angular momentum [33]. Here, the beam waist acts as a realization of the light emitted by a point source after e.g. mode conversion. Another relevant situation is the reflection of the beam from a surface in surface topology measurements, etc. By linearly superposing different LG modes, one can realize beams with amplitude, phase, and intensity patterns that simply rotate, under free space propagation, maintaining the transverse shape. These rotating structures lie at the core of a variety of sensing techniques [34–37].

We demonstrate that a meager part of the full (quantum) information is available in intensity scans and only a small fraction of this can be attributed to the rotation. This clearly confirms the potential of modal expansions inspired by quantum information protocols [38], which allow for reaching the QCRB and thus can be considered as the optimal measurement. Our results make 3D superresolution imaging more feasible and potentially useful for improving the resolution of optical microscopes.

The problem we address here is to estimate the distance traveled by a vortex beam from the beam waist to an arbitrary detection plane. We thus consider the beam waist as an object whose axial distance is to be determined. In what follows, we shall represent the fields using Dirac notation, for it makes it straightforward to expand the theory to other types of light states.

We take the beam to be represented by the pure state |Ψ(0)⟩, where z = 0 denotes the position of the object plane. The axial displacement is thus characterized by a unitary operation

$$\begin{equation}\vert {\Psi}\left(z\right)\rangle ={\mathrm{e}}^{\mathrm{i}Gz}\vert {\Psi}\left(0\right)\rangle ,\end{equation} \tag{ 2.1 }$$

where the Hermitian operator G is the generator. To pinpoint the specific form of G, it is appropriate to use the transverse-position representation Ψ(x, y; z) = ⟨ x, y|Ψ(z)⟩. It then follows directly from equation (2.1) that

$$\begin{equation}\mathrm{i}{\partial }_{z}{\Psi}\left(x,y;z\right)=-G{\Psi}\left(x,y;z\right).\end{equation} \tag{ 2.2 }$$

On the other hand, vortex beams are solutions of the paraxial Helmholtz equation [39, 40]

$$\begin{equation}2\mathrm{i}k{\partial }_{z}{\Psi}\left(x,y,z\right)={\nabla }_{T}^{2}{\Psi}\left(x,y,z\right),\end{equation} \tag{ 2.3 }$$

where k is the wavenumber and ${\nabla }_{T}^{2}={\partial }_{xx}+{\partial }_{yy}$ is the transverse Laplacian. A direct comparison leads us to

$$\begin{equation}G{\mapsto}-\frac{1}{2k}{\nabla }_{T}^{2}.\end{equation} \tag{ 2.4 }$$

The detection plane is placed at z, wherein we perform an arbitrary measurement. Given the formal analogies between spatial modes in wave optics and pure states in quantum theory, and also the mathematical similarities in describing evolution and detection of such objects, the amount of information about the axial distance z carried by the measured signal is quantified by the QFI. For pure states, as it is our case, the QFI reduces to [21]

$$\begin{equation}\mathcal{Q}\left(z\right)=4\enspace \mathrm{Var}\left(G\right).\end{equation} \tag{ 2.5 }$$

Except for the factor 4, the QFI is the variance of the generator G computed in the initial state |Ψ(0)⟩. According to the time-honored QCRB, the variance of any unbiased estimator $\hat{z}$ of the axial distance z satisfies

$$\begin{equation}\mathrm{Var}\left(\hat{z}\right){\geqslant}\frac{1}{\mathcal{Q}\left(z\right)},\end{equation} \tag{ 2.6 }$$

whose saturation provides the ultimate precision in axial distance estimation.

To proceed further, we take the structure of the transverse field to correspond with LG modes

$$\begin{align}\hfill {\mathrm{L}\mathrm{G}}_{pl}\left(r,\phi ,z\right)& =\langle r,\phi ,z\vert p,l\rangle =\sqrt{\frac{2p!}{\pi \left(p+\vert l\vert \right)!}}\frac{1}{w\left(z\right)}{\left[\frac{\sqrt{2}r}{w\left(z\right)}\right]}^{\vert l\vert }\hfill \\ \hfill & \quad {\times}{L}_{p}^{\vert l\vert }\left(\frac{2{r}^{2}}{w{\left(z\right)}^{2}}\right){\text{e}}^{-\frac{{r}^{2}}{w{\left(z\right)}^{2}}}\enspace \mathrm{exp}\left(\mathrm{i}\left[\frac{k{r}^{2}}{2R\left(z\right)}-l\phi -{\psi }_{pl}\left(z\right)\right]\right),\hfill \end{align} \tag{ 2.7 }$$

where (r, ϕ, z) are cylindrical coordinates, ${L}_{p}^{\vert l\vert }\left(\cdot \right)$ is the generalized Laguerre polynomial, l ∈ {0, ±1, ±2, ... } is the azimuthal mode index and p ∈{0, 1, 2, ... } is the radial index. The parameters R(z), w(z), and ψ_pl (z) are

$$\begin{align}\hfill R\left(z\right)& =z\left[1+{\left({z}_{\text{R}}/z\right)}^{2}\right],\hfill \\ \hfill {w}^{2}\left(z\right)& ={w}_{0}^{2}\left[1+{\left(z/{z}_{\text{R}}\right)}^{2}\right],\hfill \\ \hfill {\psi }_{pl}\left(z\right)& =\left(2p+\vert l\vert +1\right)\mathrm{arctan}\left(z/{z}_{\text{R}}\right),\hfill \end{align} \tag{ 2.8 }$$

and represent the radius of curvature of the wave front, the beam radius, and the Gouy phase [41], respectively, at an axial distance z from the beam waist located at z = 0. Here, ${z}_{\mathrm{R}}=k{w}_{0}^{2}/2$ is the Rayleigh length and w₀ the beam waist radius [42].

To facilitate the derivation of the QFI corresponding to axial displacements of vortex beams it is advantageous to use an established correspondence between eigenstates of a two-dimensional harmonic oscillator and paraxial beams [43]. For the case of LG beams, one defines

$$\begin{equation}{a}_{{\pm}}=\frac{1}{\sqrt{2}}\left({a}_{\xi }\mp \text{i}{a}_{\eta }\right),\end{equation} \tag{ 3.1 }$$

where a_ξ and a_η are dimensionless bosonic operators for each independent amplitude of the oscillator and obey the standard commutation relations $\left[{a}_{s},{a}_{{s}^{\prime }}^{{\dagger}}\right]={\delta }_{s{s}^{\prime }}$ (s, s' ∈ {ξ, η}). Similar relations are obeyed by a_±.

The eigenstates of the harmonic oscillator |n₊, n₋⟩, generated by the action of ${a}_{+}^{{\dagger}}$ and ${a}_{-}^{{\dagger}}$ on the vacuum, are nothing but LG modes with azimuthal and radial indices given by

$$\begin{equation}l={n}_{+}-{n}_{-},\hspace{25.0pt}p=\mathrm{min}\left({n}_{+},{n}_{-}\right).\end{equation} \tag{ 3.2 }$$

Applying the usual definition of the momentum operator ${p}_{s}=\frac{1}{\sqrt{2}\mathrm{i}}\left({a}_{s}-{a}_{s}^{{\dagger}}\right)$, we have

$$\begin{equation}{\tilde {\nabla }}_{T}^{2}={p}_{\xi }^{2}+{p}_{\eta }^{2}=\left({p}_{\xi }+\text{i}{p}_{\eta }\right)\left({p}_{\xi }-\text{i}{p}_{\eta }\right)=\left({a}_{+}-{a}_{-}^{{\dagger}}\right)\left({a}_{+}^{{\dagger}}-{a}_{-}\right).\end{equation} \tag{ 3.3 }$$

If we recall (2.4), and take into account that, in the proper units, $\left(x,y\right){\mapsto}\left(\sqrt{2}\xi /{w}_{0},\sqrt{2}\eta /{w}_{0}\left. \right\}\right)$, we get ${\nabla }_{T}^{2}{\mapsto}2{\tilde {\nabla }}_{T}^{2}/{w}_{0}^{2}$. In this way, the QFI of a pure LG mode $\left\vert {n}_{+},{n}_{-}\right\rangle$ reads:

$$\begin{equation}\mathcal{Q}\left(z\right)=\frac{1}{{z}_{\mathrm{R}}^{2}}\left(2{n}_{+}{n}_{-}+{n}_{+}+{n}_{-}+1\right)=\frac{1}{{z}_{\mathrm{R}}^{2}}\left[2p\left(p+\vert l\vert \right)+2p+\vert l\vert +1\right].\end{equation} \tag{ 3.4 }$$

For the particular case of the Gaussian mode LG₀₀ we get

$$\begin{equation}\mathcal{Q}\left(z\right)=\frac{1}{{z}_{\mathrm{R}}^{2}},\end{equation} \tag{ 3.5 }$$

that is, the quantum bound (per single detection) is precisely the Rayleigh range [44]. Note that the QFI is linear in |l|, which means that axial localization can be improved by using LG beams with large OAM.

Other sets of Gaussian transverse modes can be characterized using the so-called Hermite–Laguerre sphere [45, 46]. These modes are represented by a point, of spherical coordinates (θ, ϕ), on that sphere, and they are generated by the rotated operators

$$\begin{align}\hfill {a}_{1}\left(\theta ,\phi \right)& ={a}_{+}{\text{e}}^{-\text{i}\phi /2}\enspace \mathrm{cos}\left(\frac{\theta }{2}\right)+{a}_{-}{\text{e}}^{\text{i}\phi /2}\enspace \mathrm{sin}\left(\frac{\theta }{2}\right),\hfill \\ \hfill {a}_{2}\left(\theta ,\phi \right)& =-{a}_{+}{\text{e}}^{-\text{i}\phi /2}\enspace \mathrm{sin}\left(\frac{\theta }{2}\right)+{a}_{-}{\text{e}}^{\text{i}\phi /2}\enspace \mathrm{cos}\left(\frac{\theta }{2}\right).\hfill \end{align} \tag{ 3.6 }$$

In particular, θ = π/2 gives rise to Hermite–Gauss (HG) modes, and θ = 0, π to LG modes analyzed above. Combining equations (3.3) and (3.6), we get a direct generalization of (3.4); viz,

$$\begin{equation}\mathcal{Q}=\frac{1}{{z}_{\mathrm{R}}^{2}}\left[4+{n}_{1}+{n}_{2}\left(3+{n}_{2}\right)+{n}_{1}\left(3+4{n}_{2}\right)+\left({n}_{1}-{n}_{1}^{2}+{n}_{2}+4{n}_{1}{n}_{2}-{n}_{2}^{2}\right)\mathrm{cos}\left(2\theta \right)\right],\end{equation} \tag{ 3.7 }$$

where n₁ and n₂ denote the eigenvalues of the corresponding number operators. Notice that the QFI (3.7) is independent of ϕ and is optimized by (3.4). This proves that LG modes are better than their HG counterparts, as axial localization is concerned.

Apart from a trivial divergence, transversal intensity profiles of pure vortex beams do not change on propagation. More complex intensity transformations can be realized by superposing two or more vortex beams. In particular, rotating structures are of interest in microscopy [34–37]. Rotation of highly symmetric spots is easily detected and the corresponding defocusing parameter can be estimated from the measured rotation angle. This is roughly sketched in figure 1.

Zoom In Zoom Out Reset image size

To simplify the details as much as possible, we take the simple example of the superposition of two LG modes with different azimuthal numbers l ≠ l', and p = p' = 0

$$\begin{equation}\vert {{\Psi}}_{l{l}^{\prime }}\rangle =\frac{1}{\sqrt{2}}\left(\vert {\mathrm{L}\mathrm{G}}_{0l}\rangle +\vert {\mathrm{L}\mathrm{G}}_{0{l}^{\prime }}\rangle \right)=\frac{1}{\sqrt{2}}\left(\left\vert {n}_{+},0\right\rangle +\left\vert {n}_{+}^{\prime },0\right\rangle \right).\end{equation} \tag{ 3.8 }$$

In this case, the QFI about the axial position of the source reads

$$\begin{equation}\mathcal{Q}=\frac{1}{{z}_{\mathrm{R}}^{2}}\left[4+2\left(\vert l\vert +\vert {l}^{\prime }\vert \right)+{\left(\vert l\vert -\vert {l}^{\prime }\vert \right)}^{2}\right].\end{equation} \tag{ 3.9 }$$

Considering OAM as a resource for axial localization and using the maximum available OAM for one of the components of the superposition, l = ±|l_max|, the QFI is maximized when the second component is in the fundamental mode l' = 0, whereupon ${\mathcal{Q}}_{\mathrm{max}}=\left[4+2\vert {l}_{\mathrm{max}}\vert +\vert {l}_{\mathrm{max}}{\vert }^{2}\right]/{z}_{\mathrm{R}}^{2}$,and the quantum bound becomes quadratic in |l|. At first sight, the better performance of vortex superpositions (3.9) over pure vortex beams (3.4) is somehow related to beam rotation. This might suggest a proxy for the measured axial distance traveled from the waist to the detection plane. However, as we show in the next section, things are not that simple.

In single-parameter estimation, the QFI can always be accessed and the corresponding QCRB saturated with a von Neumann measurement projecting the measured signal on the eigenstates of the symmetric logarithmic derivative of the density matrix [21]. In our context, the practical implementation of such measurements requires a spatial mode demultiplexer/sorter that performs simultaneous projection of the signal onto a complete orthonormal set of modes [47, 48]. Mode separation is usually achieved with a sequence of spatial light modulators implementing a suitable unitary transformation. Unavoidable systematic errors and losses introduced by such complicated experimental setups may ruin any theoretical advantage offered by optimal strategies.

In consequence, we consider the performance of the possibly inferior, but much more robust intensity detection, because it is the simplest method at hand for the experimentalist. As the information (3.4) and (3.7) about the axial distance is carried by both intensity and phase of the measured beam, we might aptly ask how much information is sacrificed by completely ignoring the phase.

As usual, due to the noise, the detection can be considered as a random process. In consequence, the (normalized) beam intensity,

$$\begin{equation}p\left(r,\phi \vert z\right)=\vert {\Psi}\left(r,\phi ;z\right){\vert }^{2},\end{equation} \tag{ 4.1 }$$

can be seen as the probability density of a detection conditional on the axial distance z. We take the detection as dominated by shot noise, which obeys a Poisson distribution [49]: although this neglects nonclassical effects, it is still a suitable model for realistic microscopy.

The classical Fisher information about z, per single detection, thus reads [50]

$$\begin{equation}\mathcal{F}={\int }_{0}^{\infty }{\int }_{0}^{2\pi }\frac{{\left[{\partial }_{z}p\left(r,\phi \vert z\right)\right]}^{2}}{p\left(r,\phi \vert z\right)}\enspace r\mathrm{d}r\enspace \mathrm{d}\phi ,\end{equation} \tag{ 4.2 }$$

and it is a suitable tool to quantify the information content about axial displacements accessible from the detected transversal intensity profile. For simplicity, we take the pixel size negligibly small, so that any sampling effect can be ignored. We also define the radial and azimuthal Fisher informations

$$\begin{align}\hfill {\mathcal{F}}_{r}& ={\int }_{0}^{\infty }\frac{{\left[{\partial }_{z}{\int }_{0}^{2\pi }p\left(r,\phi \vert z\right)\mathrm{d}\phi \right]}^{2}}{{\int }_{0}^{2\pi }p\left(r,\phi \vert z\right)\mathrm{d}\phi }\enspace r\enspace \mathrm{d}r,\hfill \\ \hfill {\mathcal{F}}_{\phi }& ={\int }_{0}^{2\pi }\frac{{\left[{\partial }_{z}{\int }_{0}^{\infty }p\left(r,\phi \vert z\right)r\enspace \mathrm{d}r\right]}^{2}}{{\int }_{0}^{\infty }p\left(r,\phi \vert z\right)\mathit{\text{r}}\enspace \mathrm{d}\mathit{\text{r}}}\enspace \mathrm{d}\phi ,\hfill \end{align} \tag{ 4.3 }$$

respectively. They quantify the information of the radial and the azimuthal intensity components due to the z-dependence. For example, beam rotation only contributes to ${\mathcal{F}}_{\phi }$, whereas beam divergence contributes to ${\mathcal{F}}_{r}$.

For pure |LG_pl ⟩ modes, Ψ(r, ϕ; z) = LG_pl (r, ϕ, z), the classical Fisher information can be obtained in a closed form. First notice the integrand of (4.2) is independent of ϕ. Carrying out z-derivatives of Laguerre polynomials using the relation ${\partial }_{z}{L}_{p}^{l}\left(z\right)=-{L}_{p-1}^{l+1}\left(z\right)$, and changing the integration variable r yields

$$\begin{equation}\mathcal{F}=\frac{4p!}{\left(\vert l\vert +p\right)!}{\left[\frac{{\partial }_{z}w\left(z\right)}{w\left(z\right)}\right]}^{2}\underset{0}{\overset{\infty }{\int }}{\text{e}}^{-t}{t}^{\vert l\vert }{\left[2t{L}_{p-1}^{\vert l\vert +1}\left(t\right)+\left(t-\vert l\vert -1\right){L}_{p}^{\vert l\vert }\left(t\right)\right]}^{2}\enspace \mathrm{d}t.\end{equation} \tag{ 4.4 }$$

We are thus left with evaluating six different integrals involving products of Laguerre polynomials. Compact expressions for each case can be obtained from the generalized orthogonality relation [51]

$$\begin{equation}{\int }_{0}^{\infty }{\mathrm{e}}^{-t}{t}^{\mu }{L}_{p}^{l}\left(t\right){L}_{{p}^{\prime }}^{{l}^{\prime }}\left(t\right)\mathrm{d}t={\left(-1\right)}^{p+{p}^{\prime }}{\Gamma}\left(\mu +1\right)\sum\limits _{k=0}^{\mathrm{min}\left(p,{p}^{\prime }\right)}\left(\genfrac{}{}{0pt}{}{\mu -l}{p-k}\right)\left(\genfrac{}{}{0pt}{}{\mu -{l}^{\prime }}{{p}^{\prime }-k}\right)\left(\genfrac{}{}{0pt}{}{\mu +k}{k}\right).\end{equation} \tag{ 4.5 }$$

Straightforward simplifications leave us with the final result

$$\begin{equation}\mathcal{F}=4\left[2p\left(p+\vert l\vert \right)+2p+\vert l\vert +1\right]{\left[\frac{{\partial }_{z}w\left(z\right)}{w\left(z\right)}\right]}^{2}=\frac{2p\left(p+\vert l\vert \right)+2p+\vert l\vert +1}{R{\left(z\right)}^{2}/4}.\end{equation} \tag{ 4.6 }$$

For the planes of maximal wavefront curvature z_opt = ±z_R, where $R{\left({z}_{\mathrm{R}}\right)}^{2}=4{z}_{\mathrm{R}}^{2}$, the quantum limit is saturated with intensity sensitive detection

$$\begin{equation}{\mathcal{F}}_{\text{opt}}=\mathcal{Q}.\end{equation} \tag{ 4.7 }$$

Thus, for any pure LG beam and any beam waist location, two detection planes can be found, where complete information about the axial distance can be extracted with intensity-only detection. Hence, full potential of high-order vortex beams for axial metrology can be exploited with direct detection techniques. This generalizes the results obtained for Gaussian beams [32].

However, numerical analysis suggests that this result does not hold for superpositions of LG modes: a single transversal intensity scan is no longer optimal for any detector position z. In figure 2 we plot the classical Fisher information (normalized to the optimal QFI) for different vortex superpositions. First, notice that the optimal detector position is no longer at z = ±z_R, as it was for pure LG beams. Second, the larger the angular momentum l carried by the signal, the smaller portion of the total QFI can be extracted with intensity measurements.

Zoom In Zoom Out Reset image size

This is further illustrated in figure 3. Locating intensity detectors at the optimal detection planes for every value of l, the classical Fisher information grows sublinearly with l, so that progressively smaller portion of quadratic QFI is available from intensity data and the gap between the quantum bound and the performance of the best estimator from intensity data worsens with l.

Zoom In Zoom Out Reset image size

For example, for l = 2, only about 17% of the QFI is available from intensity scan at the optimal detector position, as we can see in figure 4. Of this, only a small fraction is due to changes in the azimuthal intensity profile, as arises, e.g., in rotations, except for very small propagation distances, where azimuthal profile accounts for up to 10% of QFI. This ratio quickly approaches zero with increasing l.

Zoom In Zoom Out Reset image size

The information content of rotating beam intensity distributions is surprisingly low. As bad as it sounds, this is not a negative result, but encouraging news for quantum metrology. It highlights inadequacy of simple detection techniques in this particular metrological scenario, where intensity detection fails to reveal all potentially accessible information about the parameter of interest. The hidden Information has to be accessed with advanced detection techniques. Here a generic tool is projecting/sorting the signal into optimal set of spatial modes derivable from the quantum detection theory. Hence, we uncover a huge potential of axial superlocalization based on spatial-mode projections applied to higher order vortex beam superpositions [44].

We have established the ultimate quantum limits for axial localization using vortex beams. For pure LG beams, this limit is attained with an intensity scan with the detector located at one of two optimal planes. For superpositions of LG beams, in particular of those with intensity profiles rotating on defocusing, this is no longer true. This means that microscopy methods based on rotating vortex beams may benefit from replacing traditional intensity scans with advanced mode-sorting techniques.

No new data were created or analysed in this study.

We thank Robert W Boyd and Aaron Z Goldberg for helpful discussions. We acknowledge financial support of the project ApresSF, under the QuantERA programme, which has received funding from the European Union's Horizon 2020 research and innovation programme, the program H2020 (project StormyTune), Grant Agency of the Czech Republic (Grant No. 18-04291S), the Palacký University (Grant No. IGA_PrF_2020_004), and the Spanish Ministerio de Ciencia e Innovacion (Grant PGC2018-099183-B-I00).