Towards superresolution surface metrology: Quantum estimation of angular and axial separations

We investigate the localization of two incoherent point sources with arbitrary angular and axial separations in the paraxial approximation. By using quantum metrology techniques, we show that a simultaneous estimation of the two separations is achievable by a single quantum measurement, with a precision saturating the ultimate limit stemming from the quantum Cram\'er-Rao bound. Such a precision is not degraded in the sub-wavelength regime, thus overcoming the traditional limitations of classical direct imaging derived from Rayleigh's criterion. Our results are qualitatively independent of the point spread function of the imaging system, and quantitatively illustrated in detail for the Gaussian instance. This analysis may have relevant applications in three-dimensional surface measurements.

Introduction.-High-resolution imaging is a cornerstone of modern science and engineering, which has enabled revolutionary advances in astronomy, manufacturing, biochemistry, and medical diagnostics. In traditional direct imaging based on classical wave optics, two incoherent point sources with angular separation smaller than the wavelength of the emitted light cannot be resolved due to fundamental diffraction effects [1], a phenomenon recently dubbed "Rayleigh's curse" [2]. Several techniques, including most prominently fluorescence microscopy [3], have been introduced in recent years to overcome this limitation and achieve sub-wavelength imaging [4,5]. Nevertheless, to determine the ultimate limits of optical resolution one needs to resort to a full quantum mechanical description of the imaging process [6]. In this respect, a breakthrough has been reported in a series of works [2,[7][8][9][10][11][12][13][14][15][16][17][18] initiated by Tsang and collaborators [2], who employed techniques from quantum metrology [19][20][21][22] to prove that the achievable error in estimating the angular separation of two incoherent point sources, in the paraxial approximation, is in fact independent of said separation (no matter how small), provided an optimal detection scheme is performed on the image plane. These results, which stem from the fundamental quantum Cramér-Rao bound [19,20] and de facto banish Rayleigh's curse [2], have been corroborated by proofof-principle experiments [23][24][25][26].
The majority of the studies presented so far on quantum superlocalization, however, were limited to the case of point sources aligned on the same object plane, thus neglecting their axial separation. The optical lateral resolution of an imaging system is an important characteristic, but it is not the only figure of merit relevant for the measurement of non-flat surfaces [27]. When probing surface topography, the spacing of the points in an image must be considered, along with the ability to accurately determine the heights of features. In other words, the lateral resolution must be considered in conjunction with the ability of the system to transfer surface amplitudes [28].
To address this issue, here we consider the simultaneous estimation of both angular and axial separations, as well as the 1 2 ҧ ҧ Figure 1. Schematic of the two sources. The four parameters to be estimated are: the angular separation s, the axial separation p, the angular centroid coordinatex, and the axial centroid coordinatez. corresponding centroid coordinates, of two incoherent point sources aligned in general on different object planes. These point sources may model, e.g., two emitters at the edges of a steep section on a rough surface, as indicated by the red dotted outline in Figure 1.
We tackle the problem by resorting to the toolbox of multiparameter quantum metrology, a branch of quantum technology which is attracting increasing interest thanks to its prominent role in fundamental science and applications [19][20][21][22]. We find that Rayleigh's curse does not occur even when the sources have a nonzero axial separation, and both axial and angular distances can be estimated simultaneously and with distance-independent precision by means of a single optimal quantum measurement, meeting the compatibility requirements for saturation of the multiparameter quantum Cramér-Rao bound [29,32]. These results are obtained analytically and are valid for any point spread function of the imaging system obeying the paraxial wave equation. We then specialize to the illustrative case of a Gaussian point spread function, and derive closed formulas for the achievable estimation error and its scaling with the parameters of interest as determined by the quantum Fisher information matrix, showing that in the limit of small angular and axial distances all the parameters, including the centroid coordinates, become statistically independent.
Sources and imaging system model.-We approach the problem of estimating both axial and angular separation of two point sources by following a similar approach to Ref. [2], which is in turn inspired by Rayleigh's work [1]. We assume that the detectable light on the image plane can be described as an incoherent mixture of two quasi-monochromatic scalar paraxial waves, one coming from each source. As shown in Figure 1, our two sources are in general not lying on the same object plane (an 'object plane' is a plane perpendicular to the optical axis z), and they feature an angular separation s and an axial separation p.
Considering thermal sources at optical frequencies, we divide the total emission time into short coherence time intervals τ c , so that within each interval the sources can be assumed weak, i.e., effectively emitting at most one photon. This is a standard approach for modelling incoherent thermal sources [52][53][54][55][56][57][58], and it allows us to describe the quantum state ρ of the optical field on the image plane as a mixture of a zerophoton state ρ 0 and a one-photon state ρ 1 in each time interval (neglecting contributions from higher photon numbers) [59], where ε 1 is the average number of photons impinging on the image plane. In practice, a detectable signal is obtained by measuring the optical field for a time t τ c , so that many coherence time intervals are included, resulting in a non-negligible mean photon number.
We assume in general that the image-plane field amplitude generated by each source takes the form where (x, y) are the image-plane coordinates, (x j , z j ) are the unknown coordinates of the sources j = 1, 2, x j being the coordinate perpendicular to the optical axis and z j the axial distance to the image plane (in this work we assume that the other coordinate y j = 0 is known). The amplitude function ψ(x, y, z) obeys a paraxial wave equation of the form where G is a self-adjoint differential operator featuring only x and y derivatives -for example in free space one would We shall indicate with a (x, y) the field annihilation operator at position (x, y) on the image plane, satisfying the bosonic commutation rule a (x, y) , We can then write the state ρ 1 as the incoherent mixture where the quantum state of the optical field on the image plane corresponding to the emission of one photon by the source r may be expressed as |0 being the field vacuum state. Finally, we may take ψ(x, y, 0) real, which results in some simplifications later on. This can be assumed without loss of generality, as the complex phase of ψ(x, y, 0) may be compensated by a redefinition of a(x, y) that is independent of the source parameters. However, ψ(x, y, z) will have in general a nontrivial phase profile.
Multiparameter estimation and quantum Cramér-Rao bound.-We work under the assumption that the photon statistics of our sources is Poissonian, following a similar approach as in Ref. [2]. We can thus assume that in a single run of the experiment, which lasts for M coherence time intervals, M copies of the state ρ in Eq. (1) are prepared and measured (equivalently, one may consider the input state ρ ⊗M ). On average, this yields Mε photons per run. In order to apply the standard tools of estimation theory, we further assume that ν 1 runs are performed, after which the measurement data are processed to build estimators for the unknown parameters.
In our case, the parameters of interest are the angular and axial relative coordinates and the centroid coordinates of the sources, indicated as s,x, p,z, see Figure 1. We thus write the state ρ as a function of four parameters {λ µ } µ=1,...,4 , where The statistical error (variance) ∆λ 2 µ of any unbiased estimator of the unknown parameter λ µ is lower bounded via the quantum Cramér-Rao bound (qCRb) [19,20] where H is the quantum Fisher information matrix (qFim) of the single-photon state ρ 1 (equivalently, this can be seen as the qFim per coherence time interval per photon). The prefactor on the RHS of Eq. (8) is obtained by exploiting the additivity property qFim(ρ ⊗M ) = M × qFim(ρ), and by approximating that qFim(ρ) ε × qFim(ρ 1 ) at leading order in ε (since the field vacuum state ρ 0 is independent of all source parameters and is always orthogonal to ρ 1 -see also the discussion in the Appendix of Ref. [2]). The resulting linear dependence on the total photon number νMε is characteristic of classical light sources [22,60]. The qCRb suggests that, the higher the qFIm element H µµ , the more precisely (i.e., with lower statistical error) one may be able to estimate the parameter λ µ , by performing a suitable measurement. While for a single parameter the qCRb can always be saturated asymptotically by means of an adaptive procedure [21], this is no longer the case for multiparameter estimation, as the parameters may not always be compatible [32]; we will discuss this issue in detail later in the manuscript.
Results.-We recall that the qFim elements are given by where L µ is the symmetric logarithmic derivative (SLD) for the parameter λ µ , defined implicitly by the equation The following matrix (proportional to the averaged SLD commutators) will also be of interest for our discussion, For the problem under investigation, we have derived general analytical expressions for both matrices H and Γ, as presented in detail in Appendix A. Our derivation relies on the expansion of ρ 1 in the generally non-orthogonal basis followed by standard linear algebraic manipulations. This method results in significant simplifications over previous studies of quantum superlocalization (typically relying on the explicit construction of an orthogonal basis to span the support of ρ 1 and its derivatives, as e.g. in [2]), and may be of independent interest in its own right for the field of multiparameter quantum metrology. Thanks to the representation of |Ψ j given in Eq. (5), it is easy to check that all the scalar products between the above basis vectors only depend on s = x 2 − x 1 and p = z 2 −z 1 , which in turn implies that the qFim is independent of the centroid coordinatesx andz. The corresponding physical interpretation is that the information content of the emitted light is not affected by propagation along the optical axis, or by a redefinition of the image plane origin. Additional simplifications follow from our assumption ψ(x, y, 0) ∈ R, which implies ψ|∂ x ψ = 0. We then find that the qFim is composed of the diagonal elements while the off-diagonal elements are all zero except At the same time, the only nonzero matrix elements of Γ are The following shorthands have been used: where we emphasize that γ = γ(s, p) is the only quantity depending on the source coordinates. A fundamental result can be immediately inferred from Eqs. (13) and below: for any point spread function that satisfies the paraxial wave equation, H ss and H pp are constant. This statement exemplifies how our results provide new insights on the problem of sub-wavelength imaging, while correctly reproducing what is known for p = 0 [2]. We note in particular that Rayleigh's curse does not affect the estimation of the angular separation s nor that of the axial separation p. Taking one step further, we can now investigate how close one can get to the limits imposed by the qCRb in practical experiments. In quantum estimation theory, multiparameter problems embody a nontrivial generalization of the singleparameter case [21,29,31,32]: if an estimation scheme is optimized for a particular parameter, it typically results into an increased error in estimating the others. However, in the best case scenario, such a trade-off does not apply, and one can identify an optimal protocol for the estimation of all the parameters simultaneously. This happens if and only if the parameters are compatible, i.e., they satisfy the following conditions [32]: (i) There is a single probe state yielding the maximal qFim element for each of the parameters; (ii) There is a single measurement which is jointly optimal for extracting information on all the parameters from the output state, ensuring the asymptotic saturability of the qCRb; (iii) The parameters are statistically independent, meaning that the indeterminacy of one of them does not affect the error on estimating the others. We recall also that (ii) holds iff Γ µ,ν = 0 ∀µ ν, while (iii) is equivalent to the condition H µν = 0 ∀µ ν.
In this paper we do not focus on the first condition, since our theory is built around a realistic imaging scenario in which the emission properties of the sources are fixed in advance. Yet, it is worth investigating conditions (ii) and (iii), since they have crucial implications for the actual achievability of the statistical errors given by the qCRb. Remarkably, we find that conditions (ii) and (iii) are always satisfied for the pair of parameters (s, p) -independently of the specifics of the point spread function. In the simplified scenario where (x,z) are estimated independently or known in advance, it is thus possible to construct a physical measurement and estimation strategy for s and p saturating Eq. (8) asymptotically [29,32]. On the other hand, we can see that conditions (ii) and (iii) do not hold in general for the full set of parameters (s, p,x,z). Yet, we shall see in the example below that there is at least one relevant type of point spread function for which conditions (ii) and (iii) are satisfied for all parameters in the limit s → 0, p → 0.
We consider in what follows a Gaussian beam in free space, where z R is a length parameter characterizing the beam, typically assumed of the same order as the wavelength, i.e. ∼ 1/k. Eq. (23) can be obtained, for example, if the fields generated by the two sources are well approximated by Gaussian beams in the vicinity of the image plane [61]. We thus obtain By substituting the above expressions in the qFim elements calculated previously, we find fully analytical closed formulas (as reported in Appendix B) that allow us to perform a comprehensive analysis of the multiparameter estimation problem under investigation. Furthermore, the Gaussian case bears the advantage that it can be easily compared with the existing literature that tackled the estimation of s alone (typically fixing p = 0). To support the solidity of our results, we have indeed checked that, in the limit p → 0, our expressions for H ss and Hxx match the appropriate quantities in Refs. [2,60].
Our results become particularly interesting in the regime ks, kp 1, which is precisely the one of relevance to subwavelength imaging. In the limit we have lim (s,p)→(0,0) lim (s,p)→(0,0) meaning that the four parameters s,x, p,z are approximately statistically independent when the two sources have infinitesimal angular and axial separation. The behaviour of the four diagonal qFim elements H µµ as a function of the separations s and p is illustrated in Fig. 2; the top panel can be compared directly with Fig. 2 of [2]. From the plots and from Eq. (25), we see that the qFIm diagonal elements tend to a nonzero value when s, p → 0. Hence the fundamental lower bound on the total estimation error, ∝ Tr[H −1 ], stays finite even when the two sources are infinitesimally close, instead of diverging as in direct imaging [1,2]. Eq. (26) further suggests that it should be possible to construct a single measurement that is approximately optimal for the estimation of all four parameters when ks, kp 1. The construction of such a measurement will be addressed in future work.
Conclusions.-We determined the ultimate quantum limits to the simultaneous estimation of both angular and axial separations and centroid coordinates of two incoherent point sources on different object planes in the paraxial approximation. Our results indicate that there exists a jointly optimal detection scheme that enables resolving the sources even when arbitrarily close, reasserting that Rayleigh's curse is merely an artefact of classical detection in direct imaging. In practice, a measurement apparatus approaching the optimal precision can be designed by adapting the methods of [15,16,46,47,62], in particular extending the "spatial-mode demultiplexing" or "superlocalization by image inversion interferometry" techniques [2,7] to the axially separated setting considered here.
While some of our findings were illustrated explicitly for Gaussian beams, our framework is general and can be applied to any point spread function that satisfied the paraxial wave equation, thanks to the exact expressions in Eqs. (13)- (20). This leads to qualitatively similar results as those presented here. In particular, the two most important conclusions, namely that the qFim elements for the angular distance s and for the axial distance p are both independent of s and p, and that the joint estimation of s and p fulfils the measurement compatibility condition leading to the saturation of the quantum Cramér-Rao bound in Eq. (8), are in fact valid for any point spread function.
This work constitutes an important application of multiparameter quantum estimation theory to a realistic imaging setting, extending the seminal work of Ref. [2]. Our analysis, combined with the one in [12], yields a quantum enhanced toolbox for full 3D sub-wavelength localization. This paves the way to further experimental demonstrations and innovative metrology solutions in scientific, industrial and biomedical domains, such as sub-nanometre depth mapping in rough surfaces, and dynamical interaction analysis of heterogeneous molecules in a cellular environment [4,5,27,63].
Note added.-Shortly after the initial submission of this work, quantum superresolution of two incoherent sources in three dimensions has been studied independently in Ref. [64], albeit explicit results have been obtained only in the case of a clear circular aperture.
Acknowledgments We start by observing that ρ 1 and all its derivatives -and hence the associated SLDs -are all supported in the subspace spanned by the vectors |Ψ 1 , |Ψ 2 together with We assume that the set {|Ψ j } j is linearly independent provided that x 1 x 2 or z 1 z 2 (this may be always achieved up to appropriate limiting procedures), but the set is not orthonormal in general. Yet, such non-orthogonal basis can still be used to linearly expand any state or operator. The expressions to follow will all depend on the matrix S of scalar products between the basis elements: (S1) Using the representation |Ψ j = exp −iGz j − x j ∂ x |ψ provided in the main text -and exploting that G commutes with all x− and y− derivatives, we find that all the overlaps S i j depend only on the separations s = x 2 − x 1 , p = z 2 − z 1 , and not on the centroid coordinates. For example: where we have also used that ∂ x is anti-Hermitian. Similar simplifications, together with the paraxial wave equation ∂ z j |Ψ j = −iG |Ψ j , allow us to write all the matrix elements of S as per where it is worth noting that only γ = γ(s, p) depends on the source separations, while ∂ x ψ|∂ x ψ , G and G 2 are independent of all source parameters. In the above we exploited the further simplification ψ|∂ x ψ = 0, which follows from the assumption that ψ(x, y, 0) is real (see main text).
In the non-orthogonal basis {|Ψ j } j , the matrix representation of ρ 1 reads Since the above expression may at first sight appear strange, we emphasize that in a non-orthogonal basis, ρ = ρ † does not imply ρ i j = ρ * ji . The derivatives of ρ are instead represented by the matrices We can now employ a symbolic manipulation software (e.g. Mathematica) to solve explicitly the SLD equations L µ ρ + ρL µ = 2∂ µ ρ, which gives us the matrix representation of the SLDs (L x 1 , L x 2 , L z 1 , L z 2 ) in the basis {|Ψ j } j . To find the SLDs corresponding to the variables of interest (s,x, p,z) we then simply apply the rotation Once the SLDs are known, the results presented in the main text follow from the relation H µν + iΓ µν = Tr[ρ 1 L µ L ν ] . (S10)