Quantum limits of localisation microscopy
Evangelia Bisketzi Dominic Branford Animesh Datta
quant-ph physics.optics








Department of Physics, University of Warwick, Coventry, CV4 7AL, United Kingdom


Localisation microscopy of multiple weak, incoherent point sources with possibly different intensities in one spatial dimension is equivalent to estimating the amplitudes of a classical mixture of coherent states of a simple harmonic oscillator.
This enables us to bound the multi-parameter covariance matrix for an unbiased estimator for the locations in terms of the quantum Fisher information matrix, which we obtained analytically.
In the regime of arbitrarily small separations we find it to be no more than rank two – implying that no more than two independent parameters can be estimated irrespective of the number of point sources.
We use the eigenvalues of the classical and quantum Fisher information matrices to compare the performance of spatial-mode demultiplexing and direct imaging in localisation microscopy with respect to the quantum limits. 


Quantum limits of localisation microscopy
    Animesh Datta
    December 30, 2023
=========================================





§ INTRODUCTION



Precisely locating multiple single emitters is a key challenge in fluorescence microscopy. 
The process of estimating these locations depends on the quality of the image obtained by the microscope.
One of the major limitations to the image quality, known since Abbe and Rayleigh, lies in spatially resolving objects substantially smaller than half the wavelength of the light involved <cit.>.
Known as the Rayleigh limit or diffraction limit, it is a consequence of the diffraction of light due to its wave nature. 

Over the last couple of decades, ways to circumvent the Rayleigh limit in far-field fluorescence microscopy have been invented <cit.>.
Confocal methods such as STED, RESOLFT, and SSIM <cit.> use patterned illumination to spatially modulate the fluorescence pattern of emitters within a diffraction-limited region such that not all of them emit simultaneously, thereby achieving sub-Rayleigh resolution.
Other far-field methods such as PALM, fPALM and STORM <cit.> temporally modulate the fluorescence pattern of emitters with weak laser pulses stochastically such that only a low density of emitters are active within the Rayleigh limit at one time.
Repeating the process many times, images with sub-Rayleigh resolution are reconstructed from the measured positions of individual emitters.
These techniques, with resolution of tens of nanometers, have provided insights into biological processes at the cellular scale that were hitherto unattainable <cit.>.

Though immensely powerful and impressive, none of these methods seek to extract all the information available in the emitted light field.
As in conventional fluorescence microscopy these techniques use `direct imaging'—intensity measurements on the image plane—to extract information from the incident light.
That there is indeed more information in the light field to be extracted was shown by <cit.>. Using methods from classical and quantum estimation theory, it was shown theoretically that two arbitrarily close incoherent point sources may be resolved, and that this may be achieved in practice using a spatial-mode demultiplexing (SPADE) measurement.
In the few years since, theoretical studies have considered different source arrangements or parameters of interest <cit.> in one as well as in two and three spatial dimensions <cit.>.
Other theoretical studies have explored various detection systems that could achieve the ultimate precision in imaging or get close to it <cit.>.
Several experiments have demonstrated some of the principles underlying these detection systems <cit.>.
Advances in this area have been recently reviewed by <cit.>.



Realistic imaging scenarios typically involve more than two point sources or even extended objects.
It has been shown that an extended one-dimensional object much smaller than the Rayleigh limit described only in terms of its centroid and effective radius can be approximated by a two-level quantum system <cit.>.
Theoretical optimality of certain measurement techniques in estimating this effective radius size has also been established in one and two spatial dimensions <cit.>.
Order-of-magnitude bounds on the precision of estimating the normalised moments of extended sources smaller the Rayleigh limit have also been obtained <cit.>.


In this paper, we provide an analytical lower bound on an unbiased estimator's covariance (mean square error) matrix for localisation microscopy – simultaneously estimating the locations of N incoherent, weak point sources of unequal but known intensities in one spatial dimension. The bound is provided by the the quantum Fisher information matrix.
For a Gaussian point spread function (PSF), we first describe the light field on the image plane as a classical mixture of coherent states.
We use this to derive the quantum Fisher information matrix analytically. In the limit of the point sources approaching a single point, we find its rank to be no more than two.
As the inverse of the quantum Fisher information matrix lower bounds the covariance matrix, our result implies that no more than two independent parameters can be estimated in localisation microscopy in the limit of arbitrarily small separations.
In this limit, we provide a mathematical explanation for our observation in terms of an approximation of the light field involving only the first two Hermite-Gauss modes.
Finally, we compare performance of conventional direct imaging and the recently proposed SPADE <cit.> in localisation microscopy with the quantum bounds we obtain. 
In the limit of the point sources approaching a single point, we find the classical Fisher information matrices for both these detection systems to be rank one. Furthermore, in the sub-Rayleigh limit, SPADE does not attain the quantum limit for localisation microscopy. For the subset of parameters where scalings may be optimal, we find SPADE to be short of the quantum limit in absolute precision.

This paper is organised as follows: In Section <ref> we provide a quantum mechanical description of localisation microscopy.
The appropriate framework to study the quantum limits of the localisation problem is quantum estimation theory, the toolbox of which is described in Section <ref>,   leading to the definition of the quantum Fisher information matrix (QFIM).
In Section <ref> we provide an analytic expression of the QFIM for localisation microscopy, our main technical result. We then draw conclusions about its rank and its implications for localisation microscopy. 
We end in Section <ref> with further insights and discussions about the sinc PSF and the potential of detection systems attaining the quantum limits of localisation microscopy.








§ QUANTUM DESCRIPTION OF LOCALISATION MICROSCOPY



We consider localisation microscopy – the problem of estimating the locations of N incoherent point sources or emitters located in a one-dimensional spatial configuration as in Fig. <ref>.
As we assume them to be weak, such that on average no photons arrive on the image place within a coherence time with probability (1-ϵ), where ϵ≪ 1 and one photon arrives with probability ϵ.
We also assume the optical field on the image plane to be quasi-monochromatic and paraxial <cit.>.
The quantum state of this optical field is then

    ρ_opt≈ (1-ϵ)ρ_vac+ϵρ,

where we have neglected terms of second and higher orders in ϵ and ρ_vac=|vac⟩⟨vac| is the vacuum state and ρ is the one-photon state.
 
 The one-photon density matrix on the object plane is an incoherent mixture of position eigenstates ρ= ∑_i=1^N w_i |χ_i⟩⟨χ_i|, where w_i are the relative intensities with ∑_i=1^N w_i=1.
An imaging system maps ĉ_x^†, the creation operator producing one photon in the position x on the object plane, to the corresponding image plane operator ĉ_i^† <cit.> 

    ĉ_i^†= ∫ dx Ψ_(x-χ_i) ĉ_x^†,

where χ_i is the position on the source on the object plane and ψ_(x) is the PSF.
On the image plane this becomes

    ρ=∑_i=1^N w_i|ψ_i⟩⟨ψ_i|,

where

    |ψ_i⟩=∫ dx Ψ_(x-χ_i) |x⟩,

as follows from Eq. (<ref>).


An ideal imaging system with Ψ_(x) = δ(x) is free of any Rayleigh limit as it transmits all spatial frequencies from the object to the image plane.
In practice, a Gaussian PSF
 
    ψ_(x)= 1/(2 πσ^2)^1/4 e^-x^2/4σ^2,

with σ= λ/(2 πNA), where NA is the numerical aperture of the imaging system is a good approximation for quasimonochromatic paraxial light <cit.> and also allows us to obtain analytical results.
For such a PSF, the state of Eq. (<ref>) has an intensity distribution of the form illustrated in Fig. <ref>.
For a Gaussian PSF, the |ψ_i⟩ can be expanded in the Hermite-Gauss (HG) basis as (See Appendix <ref>)

    |ψ_i⟩ =∑_k=0^∞α_i^k/√(k!) e^-α_i^2/2  |ϕ_k⟩≡|α_i⟩,

where |ϕ_k⟩ are the HG modes
[Unlike the conventional quantum optical coherent states which reside in the phase space of the electromagnetic field, our coherent states reside in physical space on the image plane.
This mathematical form was also identified by <cit.> but only used for numerical calculations.
]
This has the same mathematical form as the coherent states, produced by the displacement operator 𝒟(α_i) = e^α_ia^†+α_i^*a <cit.> acting on the ground state of the harmonic oscillator with α_i=χ_i/2σ∈ℝ the dimensionless positions of the sources.
Thus the one-photon state on the image plane is
 
    ρ≡ρ_α = ∑_i=1^N w_i|α_i⟩⟨α_i|,

 a classical mixture of coherent states in the HG basis.

The above is a quantum optical rendition of localisation microscopy—a classical optics problem.
It enables us to harness the mathematical formalism associated with coherent states and provides a basis that spans the space of the quantum state as well as its derivative.
The latter is an essential ingredient of deriving the quantum Fisher information matrix analytically in Section <ref>.
We also hope that this description will provide insights into the quantum limits to localisation microscopy in the presence of shot noise and assist in designing detection systems that attain these quantum limits. 

 




§ QUANTUM ESTIMATION THEORY






Localisation has long been treated as an estimation problem with the unknown locations of the sources χ≡{χ_i}, i=1,…,N being the parameters to be estimated <cit.>.
In our formulation, the limits to the localisation of the point sources are the same as estimating the amplitudes α≡{α_i}, i=1,…,N of the coherent states in Eq. (<ref>). 
Let these estimates be α̃≡{α̃_i}.
The precision of our estimate is then given by the covariance (or mean square error) matrix defined as

    Cov[α]  = ∑_z p(z|α) (α - α̃)^T(α - α̃),

where p(z|α) is the probability distribution of the collected data labelled by, for instance, the pixel z on the image plane.
Cov[α] is a positive symmetric matrix whose i-th diagonal element denotes the variance of an estimator of α_i given the data collected.
The (i,j)-th off-diagonal element denote the covariance in the estimation of α_i and α_j. 

Given the data collected, the maximum amount of information that can be extracted from it to obtain the most precise estimate of the locations is given by the Cramér-Rao bound (CRB) <cit.>. 
For an unbiased estimator, this bound is given by

    Cov[α] ≥1/Mϵ[𝒞(ρ_α, Π_z)]^-1,

 where M is the number of coherence times over which the data is collected, making Mϵ the total photon count.
This inequality is saturable but generally only in the asymptotical limit of many repetitions <cit.>. 
𝒞(ρ_α, Π_z) is the classical Fisher information matrix (CFIM) whose elements are given by <cit.> 

    [𝒞(ρ_α, Π_z)]_μν =
    ∑_z 1/p(z|α )∂ p(z|α)/∂α_μ∂ p(z|α)/∂α_ν.

The probability distribution p(z|α) = Tr( ρ_αΠ_z) results from detecting the light on the image plane using a specific detection system Π_z.
Fluorescence microscopy typically employs intensity detectors Π_z = {n_z}, n = 0, 1, ⋯, at each pixel x, known as direct imaging.
It is then evident that the CFIM depends on the detection system used, and not surprising that it determines the amount of information that can be extracted from the light field at the image plane. 

To identify the quantum limit on the precision of localisation microscopy, the CFIM must be maximised over all possible physically allowed detection systems.
This set is given by positive operator-valued measures (POVMs) <cit.> and the maximisation is bounded as <cit.>.

    max _{Π_z} 𝒞(ρ_α, Π_z) ≤𝒬(ρ_α),

by the quantum Fisher information matrix (QFIM). 
Its matrix elements are given by

    [𝒬(ρ_α)]_μν=Tr[ ρ_αL^μL^ν+L^νL^μ/2],

with L_μ being the symmetric logarithmic derivative (SLD) corresponding to the parameter α_μ.
The SLD is determined by the Lyapunov equation 

    2∂ρ_α/∂α_μ = ( ρ_αL^μ +L^μρ_α).


The quantum limit to localisation microscopy is thus given by

    Cov[α] ≥1/Mϵ[𝒞(ρ_α, Π_z)]^-1≥1/Mϵ[𝒬(ρ_α)]^-1.

The QFIM depends only on the light field on the image plane and determines the maximum amount of information that can be extracted from it using detection systems allowed by quantum mechanics.
Deriving an analytical expression for 𝒬(ρ_α) for state in Eq. (<ref>) is our main result, which we present in the next section.


A practical issue following the identification of the quantum limit is its attainability.
For cases where a single parameter is unknown then a measurement can be found to satisfy the equality of Eq. (<ref>), which involves projecting onto the eigenstates of the SLD <cit.>. However this strategy does not generalise to multiple parameters, as in localisation microscopy, in general.

For multi-parameter estimation the attainability is tantamount to saturating the second inequality in Eq. (<ref>). A necessary condition for the saturability of any scalar form of Eq. (<ref>) is the satisfaction of weak commutativity <cit.>
 
    Tr(ρ_α [L^μ,L^ν])= 2Tr(Im( ρ_α L^μL^ν)) = 0.

Moreover, through the quantum theory of asymptotic normality <cit.>, this condition becomes sufficient with the application of collective measurements over multiple copies of ρ_α  <cit.>.

Any scalar function of the covariances can be bounded by the inverse of QFIM with the lower bound following from the spectral decomposition of QFIM. To that end, calculating the eigenvalues of the QFIM and their scaling is of importance for the multi-parameter estimation.
For localisation microscopy, ρ_α as in Eq. (<ref>) as well as its derivative are real matrices. Thereby, L^μ,L^ν are also real and the above condition is always satisfied
[Since the density matrix and its derivatives are real-valued in the orthonormal {|ϕ_k⟩} basis, Eq. (<ref>) is a system of equations with real coefficients.
Hence L^ν must be real as well, and so ρ_α L^μL^ν is real-valued. We thank Ben Wang for bringing this to our attention].
The quantum limit for localisation microscopy is therefore attainable, at least in principle, although collective measurements over multiple copies <cit.> of the light field on the image plane may be required.


Alternative parameterisations of the system—where the new parameters α' are functions of the old parameters α—can be dealt with by a transformation of the QFIM.
Given the transformation matrix B with elements B_ij=∂α_i/ ∂α'_j, the QFIM of the transformed parameters is <cit.>
 
    𝒬'=B  𝒬  B^T,

provided the transformation is non-singular. This can be used to recast our results in terms of, for instance, the moments of the point source distribution. 



 




§ RESULTS
 

We now present our main result – the analytical expression of the QFIM for localisation microscopy.
This expression allows us to conclude that the QFIM is a rank two matrix as α_i → 0.
Eq. (<ref>) then implies that the eigenvalues of Cov[α] remains finite for no more than two independent parameters.
Thus, no more than two independent parameters can be estimated from the entire set α as α→ 0. 

We lack a fully satisfactory physical explanation for this restriction on the number of estimable parameters, but provide an explanation involving only the first two Hermite-Gauss modes for α_i ≪ 1.  




 §.§ Analytical expression of QFIM
 


The state in Eq. (<ref>) can be expressed in the basis of {|α_i⟩, â^†|α_i⟩} as

    ρ_α   = A [ D_w   0;   0   0 ] A^†≡ A ρ_A A^†

where

    A = [    |α_1⟩    |α_2⟩        ⋯    |α_N⟩ a^†|α_1⟩        ⋯ a^†|α_N⟩ ]

and D_w = ( w_1 , w_2 , ⋯ , w_N ) denotes a diagonal matrix.
Although the basis used in Eq. (<ref>) is non-orthogonal this representation can still be used to evaluate the QFIM <cit.>.
The coherent states  {|α_i⟩} are linearly independent and span the support of the state in Eq. (<ref>).
The support of the derivative is spanned by {|α_i⟩} and {â^†|α_i⟩},  which are also linearly independent.

The Grammian matrix

    Υ = A^† A,

whose elements consist of the scalar products between the basis vectors ⟨α_j | α_k|⟩, ⟨α_j | a^† | α_k|⟩, ⟨α_j | a | α_k|⟩, and ⟨α_j | aa^† | α_k|⟩
is in block form,

    Υ =	[  Υ_αα Υ_α d; Υ_d α Υ_d d ],

where 

    (	Υ_αα)_ij   = ⟨α_i | α_j|=⟩ e^ - (α_i-α_j)^2/2 , 
    ( 	Υ_α d )_ij= Υ_d α^†   = ⟨α_i | a^† | α_j|=⟩α_i e^ -(α_i-α_j)^2/2=D_αΥ_αα, 
    ( Υ_d d)_ij   = ⟨α_i | aa^† | α_j|=⟩ (α_i α_j + 1) e^ -(α_i-α_j)^2/2 
       =D_αΥ_αα D_α + Υ_αα,

and  D_α = ( α_1 , α_2 , ⋯ , α_n ).


Since ∂_α|α⟩ = (a^† - α) |α⟩ for real α, 
the derivative of the quantum state is

    ∂_j ρ_α = A  w_j[ -2α_j E_j       E_j;       E_j         0 ]  A^T≡ A (∂_jρ)_A A^†,
 
where ∂_j denotes the derivative with respect to α_j and (E_j)_kl = δ_jkδ_jl.
Similarly, the SLD L_A^i can be written in the generic form 

    L^j = A L_A^j A^T = A [  L_αα^j L_α d^j;  L_dα^j  L_dd^j ] A^T

where L_αα^j corresponds to the elements ⟨α_i| L^j |α_j⟩, L_α d^j to ⟨α_j| L^j â^†|α_i⟩ etc.
The Lyapunov equation Eq. (<ref>) can be now rewritten as

    2∂_jρ_A = ρ_AΥ L_A^j+L_A^jΥρ_A,

and the QFIM elements from Eq. (<ref>) as

    [𝒬(ρ_α)]_jk= ∂_jρ  L^k=∂_jρ_AΥ  L^k_AΥ.


Using the Tracy-Singh block kronecker product  and the block column "vecb" operator <cit.> defined as

    L_A^j =[  L_αα^j; L_d α^j; L_α d^j;  L_dd^j ],

where X=vec(X) is the column vectorisation of a matrix and X its transpose.
Eq. (<ref>) can be blockwise vectorised to

    2 ∂_jρ_A =( 𝕀 (ρ_AΥ ) + (ρ_AΥ ) 𝕀) L_A^j

with 𝕀 being the identity matrix.
Using the matrix identity <cit.>
    A^T B C D^T=A^T^T(D  B) C,

the QFIM elements from Eq. (<ref>) can be re-expressed as

    [ 𝒬(ρ_α)]_ij   = ∂_iρ_A^T (ΥΥ) L_A^j
       = w_i[ -2 α_iE_i       E_i       E_i         0 ][  Γ_αα^j; Γ_d α^j; Γ_α d^j;  Γ_dd^j ],

where we have defined 

    (ΥΥ) L^j_A 
    = Γ^j 
    = [  Γ_αα^j; Γ_d α^j; Γ_α d^j; Γ_dd^j, ]

which is the outstanding quantity to be determined.

We now recast Eq.(<ref>) and (<ref>) as

    2 ∂_iρ_A    =( Υ^-1ρ_A + ρ_AΥ^-1) (ΥΥ) L_A^i
       =( Υ^-1ρ_A + ρ_AΥ^-1) Γ^i.

Putting it all together, we obtain

    [ -4 w_iα_iE_i;     2 w_iE_i;     2 w_iE_i;            0 ] = [       0;   𝔸   0;       0; 0 0 0 0 ] [  Γ_αα^j;  Γ_dα^j; Γ_α d^j;  Γ_dd^j ],

 where

    𝔸=[ D_w⊗υ_αα+υ_αα⊗ D_w          D_w⊗υ_α d         υ_α d⊗ D_w;           D_w⊗υ_dα           D_w⊗υ_dd                  0;          υ_dα⊗ D_w                  0          υ_dd⊗ D_w ],

and {υ_αα, υ_α d, υ_dα, υ_dd} defines the inverse of Υ via

    Υ^-1 = 
    [  υ_αα υ_α d;  υ_dα  υ_dd ].

Note that the inverse Υ^-1 always exists since Υ is the Grammian matrix of linearly independent vectors.
The elements of Υ^-1 can be found using the formula of blockwise inversion (See Appendix <ref>).


Noticing that Γ_dd^j does not contribute in Eq. (<ref>), Eq. (<ref>) can be reduced to

    [ -4 w_iα_iE_i;     2 w_iE_i;     2 w_iE_i ] 
    = 𝔸[  Γ_αα^j;  Γ_dα^j; Γ_α d^j ]

where
𝔸 is a 3N^2× 3N^2 invertible matrix unless α_i =α_j for some i,j, which is a singular case for which the rank of the density matrix reduces. 
Hence the unique solution to Eq. (<ref>) is 

    [  Γ_αα^j;  Γ_dα^j; Γ_α d^j ]= 𝔸^-1[ -4 w_jα_jE_j;     2 w_jE_j;     2 w_jE_j ],

where the block matrices that compose the 𝔸^-1 can be found by using the formulas for blockwise inversion (See Appendix <ref>).

Substituting Eq. (<ref>) into Eq. (<ref>) gives us the QFIM elements

    [𝒬(ρ_α)]_ij   = 2 w_i w_jE_i[
    𝕀⊗Υ_dαΥ_αα^-1 + Υ_dαΥ_αα^-1⊗𝕀 - 2α_i 𝕀⊗𝕀] S^-1[
    𝕀⊗Υ_αα^-1Υ_α d + Υ_αα^-1Υ_α d⊗𝕀 - 2α_j 𝕀⊗𝕀] E_j
       16mu +
    4 w_i δ_ij[ 1+α_i^2 - (Υ_ααD_αΥ_αα^-1D_αΥ_αα)_ij]

where S^-1 = ( Υ_αα^-1⊗ D_w + D_w⊗Υ_αα^-1)^-1 is an N^2× N^2 matrix and Υ_αα is the inverse of the submatrix of Υ which exists, as it is the Grammian matrix of linearly independent vectors {|α_i⟩}. 
Eq. (<ref>) is an analytic expression for the QFIM elements for localisation microscopy and our main result.



Fig. <ref> shows the elements of the QFIM for the localisation microscopy of three point sources.
We choose them to be equidistant, that is, (α_1,α_2,α_3)=(x,2x,3x) and w_1 = w_2 = w_3 = 1/3 for illustration purposes.
Note the non-zero off-diagonal elements evidencing correlations in the precision around and below the Rayleigh limit of x ∼ 1. 

While the diagonal elements are all non-vanishing, more crucially as x → 0 the diagonal and off-diagonal elements combine to make the QFIM singular.
This is revealed by a closer analysis of the QFIM matrix as in Fig. <ref> which shows that only two of its eigenvalues remain non-zero as the sources approach each other. 
This is in spite of all the diagonals elements of the QFIM remaining non-zero even as x → 0, as Fig. <ref> shows.

This behaviour of only two non-zero eigenvalues also holds for other values of N.
 We have explicitly checked this for N=4,… ,10 as well as when the sources are not equally spaced.
In Fig. <ref> in Appendix <ref> we plot the eigenvalues of the QFIM for N=4,5 as further examples. 
In the case of different relative intensities the results are the same except of the limiting case of one extremely bright source w_j≫ 1, w_i≠ j≪ 1, where the rank of the QFIM is approximately one (Fig. <ref> in Appendix <ref>).


Since the QFIM has rank two as x → 0, its inverse is ill-defined except on a two-dimensional subspace.
This implies that the N × N covariance matrix for localisation microscopy, as per Eq. (<ref>), will also be unbounded except on a two-dimensional subspace.
Thus, no more than two independent parameters can be estimated in localisation microscopy as the point sources approach each other.



In other words, the rank-deficient nature of the QFIM shows that a form of the Rayleigh limit resurfaces for any N>2.
This had been suggested by previous works based on order-of-magnitude bounds for the diagonal elements on the CFIM <cit.> or uppers bounds on the diagonal elements of the QFIM <cit.>.
Our analytical expression for the full QFIM—its diagonal and off-diagonal elements for any N—shows that this rank two behaviour is truly quantum mechanical in origin.
Furthermore, knowing the full QFIM matrix allows us to uncover the nature in which N-2 of the eigenvalues approach zero. We return to the behaviour in which this rank deficiency or Rayleigh limit emerges in Sec. <ref>.






 §.§ Why rank two?


We now provide an explanation for the rank deficiency of the QFIM in the regime of small separations which can be seen as the re-emergence of the Rayleigh limit.
To that end, we expresses the state in Eq. (<ref>) in terms of the real-valued displacement operator 𝒟(α_i)= e^α_iĉ_i^†-α_iĉ_i as 

    ρ=∑_i√(w_i) 𝒟(α_i) |0⟩⟨0|√(w_i) 𝒟^†(α_i).

In the limit of very small separations (α_i ≪ 1), the displacements are approximately

    𝒟(α_i)=𝕀+α_i( â^†-â) + 𝒪(α_i^2)  ,

where 𝕀 is the identity operator and the displacement α_i is real.
Up to the second order in α_i, the normalised quantum state of the light field on the image place is then

    ρ_α^(2)=
     [ 1-𝒞_2   𝒞_1;   𝒞_1   𝒞_2;       ],

where 𝒞_i are the first two moments

    𝒞_1=∑_i=1^N w_iα_i,   𝒞_2=∑_i=1^N w_iα_i^2.


Eq. (<ref>) describes the state of two-level quantum system—the two levels being the first two HG modes.
A similar approximation which described the state relative to a PSF centred at a fixed reference point was used in Ref. <cit.> to estimate the centroid and the effective radius of a distribution of incoherent point sources. 
We now consider the more general problem of estimating the location of N point sources.

The QFIM for α (See Appendix <ref>) is

    𝒬(ρ_α^(2)) ≡𝒬=1/𝒜( 𝐈  α) 
    ℳ[ 𝐈^T; α^T ],

with

    ℳ = [ M_11 M_12; M_21 M_22 ],

where  M_11= (𝒞_2-1)𝒞_2,
M_12=M_21=𝒞_1( 1-2  𝒞_2),
M_22=4 𝒞_1^2 -1,
𝒜=( 𝒞_2-1)𝒞_2+𝒞_1^2,
and  I=( 1   1 … 1 )^T.

The QFIM 𝒬 is an N × N matrix, which is a product of three matrices of dimensions N × 2, 2 ×2 and 2 × N.
Since (AB) ≤min{ (A),  (B) }, and the matrix ℳ has rank 2, the QFIM 𝒬 has rank no more than two.
Although a two-level quantum system has the potential of estimating three real parameters, localisation microscopy in this limit can estimate only two as the two-level system possesses a real density matrix
[As the localisation parameters α are real, Tr( ρ_α^(2) σ_y)=2  Im( 𝒞_1) =0, where σ_y is the Pauli Y matrix.]
This is another way of arguing that as the point sources get closer, the light field on the image plane has enough information to estimate only two parameters.
A physical reason for this observation would be highly desirable.




§ DISCUSSION


Our analytical expression for the QFIM for localisation microscopy has enabled us to show that as point sources get closer, no more than two independent parameters can be estimated.
A rank-deficient QFIM occurs when the quantum state does not contain enough information to permit the estimation some of the parameters or combinations thereof.
The parameters that can be estimated correspond to the non-zero eigenvalues of the QFIM.
Without additional knowledge of the source distribution this restricts us to estimating functions of the first two moments f(𝒞_1,𝒞_2) only deep in the sub-Rayleigh limit. As Eq. (<ref>) shows, when all {α_j } are unknown as in localisation microscopy, there is vanishing information about any single α_i itself.
This is in contrast to the scalar QFI [𝒬(ρ_α)]_ii for α_i which is non zero, but assumes that all the other {α_j} are known.
The manner in which the eigenvalues of the QFIM tend to zero is of interest in the search for optimal detection systems for localisation microscopy.
Numerical fitting in Fig. <ref> shows the vanishing eigenvalues of the QFIM approach zero polynomially.
The degree of the polynomial is given by d=2⌊μ -1/2⌋, where μ is the order the eigenvalue when arranged in descending order and ⌊·⌋ is the floor function.
These scalings are now extracted from the elements of the full QFIM of the localisation parameters α– rather than from bounds on estimating the various moments independently as in previous works <cit.>.
















Unlike the latter, we can now compare the absolute performance of detection systems for localisation microscopy relative to its quantum limit. Indeed, while Fig. <ref> shows the 2n-th eigenvalue of the QFIM closely parallel to the n-th eigenvalue of the CFIM for SPADE <cit.>, there is a large gap in the absolute terms. This could be due to the sub-optimality of SPADE for estimating the ⌊ N/2⌋ parameters it is sensitive to[Conventional SPADE is not sensitive to all the parameters needed to describe the sources' distribution, only its even moments <cit.>.]. Similar scalings were observed with detection using superpositions of the conventional SPADE basis <cit.> that are sensitive to the other half of the moments. 
For reference over a range of separations,  Fig. <ref> in Appendix <ref> shows the eigenvalues of the CFIM for SPADE as well as direct imaging. Note that for both, the CFIM tends towards a rank one matrix. 



Finally, although our analytical result is derived with a Gaussian PSF, we expect the rank deficiency of the QFIM to be present in a more general family of PSFs.
To that end, Fig. <ref> shows the numerically obtained eigenvalues of the QFIM for three equidistant point sources of equal intensities under a sinc PSF (See Appendix <ref>) defined as

    ψ_PSF(x)=1/√(σ)( π x/σ)

This PSF is the exact form for diffraction through a sharp one-dimensional slit which in its principal peak is well-approximated as Gaussian.

An approximation involving the first two spherical Bessel modes as in Sec. <ref> can be performed for a sinc PSF as well, leading to similar insights.
A proof of this rank deficiency for arbitrary PSFs and a physical explanation remains an open question.


To conclude, we have obtained several insights into the quantum limits of localisation microscopy via an analytical expression for the QFIM. In particular, the behaviour of the eigenvalues of the QFIM deep in the sub-Rayleigh limit revealed that only two parameters are eventually estimable. It also enabled us to compare the performance of known detection systems relative to the quantum limit in absolute terms, a question left open in the literature <cit.>. The gap identified by us should motivate the search for detection systems, ideally on a single copy of the light field on the image plane, seeking to reduce or eliminate it. 


§ ACKNOWLEDGEMENTS

We thank Francesco Albarelli, Jamie Friel, Mankei Tsang, Liang Jiang, Lijian Zhang, Alex Retzker for enlightening discussions.
This study has been supported by the UK EPSRC (EP/K04057X/2),
the UK National Quantum Technologies Programme (EP/M01326X/1, EP/M013243/1) and the University of Warwick Global Partnership Fund.



*equationsection

§ EXPRESSING THE DENSITY MATRIX IN THE HG BASIS

 
The density matrix is written in terms of the kets |ψ_i⟩, which are expressed in the position space as in Eq. (<ref>).
We assume a normalised Gaussian point spread function (PSF) of the form
 
    ψ_(x)=1/(2 πσ^2)^1/4 e^-x^2/4σ^2,

 and so
 
    |ψ_i⟩ = ∫ dx ψ_(x-χ_i) |x⟩.

The kets |ψ_i⟩ can be expressed in terms of the complete Hermite-Gauss modes as
 
    |ψ_i⟩ = ∑_q=0^∞⟨ϕ_q|ψ_i||%s⟩⟩ϕ_q,

where |ϕ_q⟩ are the Hermite-Gauss modes, which can be expressed in the position space as <cit.>
    |ϕ_q⟩=1/(2πσ^2)^1/41/√(2^qq!)∫ dx H_q( x/√(2)σ) e^- x^2/4 σ^2|x⟩,

where H_q(x) are the Hermite polynomials.
The coefficients of the expansion Eq. (<ref>) are 
 
    ⟨ϕ_q|ψ_i|   ⟩=1/√(2πσ^2)1/√(2^qq!)∫ dx dx' H_q( x/√(2)σ) e^- x^2/4 σ e^-(x'-χ_i)^2/4σ^2⟨x|x'|
    ⟩   =e^--χ_i^2/8σ^2/√(2πσ^2 2^qq!)∫ dxH_q( x/√(2)σ) e^- (x/√(2)σ -χ_i/√(2)σ )^2
       = ( χ_i/2 σ)^qe^-1/2(χ_i/2 σ)^2/√(q!)

Setting χ_i/2 σ= α_i we get
 
    |α_i⟩≡|ψ_i⟩=∑_q=0^∞α_i^q/√(q!) e^-α_i^2/2|ϕ_q⟩

 which has the same mathematical form as the coherent states with {|ϕ_q⟩} forming the Fock basis <cit.>.
 
 
The state in Eq. (<ref>) can be also written in terms of the displacement operators 𝒟(α_i)=e^α_i (a^†-a), with α_i = χ_i/2 σ∈ℝ
    ρ_α=∑_i√(w_i)𝒟(α_i) |0⟩⟨0|√(w_i)𝒟^†(α_i)

where 𝒟(α) is the displacement operator.

The derivative of each coherent state with respect to its real amplitude α is given by

    ∂|α⟩/∂α    = ∂ D(α)/∂α|0⟩= (â^† -α)|α⟩, 
    ∂⟨α|/∂α   = ∂ D^†(α)/∂α⟨0|= (â-α)|α⟩,

which yields the formula Eq. (<ref>).




§ ANALYTIC RESULTS FOR N SOURCES

The Tracy-Singh product <cit.> defined for matrices A and B subdivided into blocks A_ij and B_kl
is A  B where the (i,j)-th block of A  B is A_ij B whose (k,l)-th block is in turn A_ij⊗ B_kl.
That is if A,B are block matrices with 
    A = [ A_11 A_12; A_21 A_22 ],
     and 
    B = [ B_11 B_12; B_21 B_22 ],

then the Tracy-Singh product is

    A ⊙ B = [ [c|c]
    A_11⊙ B       A_12⊙ B;       A_21⊙ B       A_22⊙ B ] 
    =
    [ [cc|cc]
    A_11⊗ B_11         A_11⊗ B_12         A_12⊗ B_11         A_12⊗ B_12;         A_11⊗ B_21         A_11⊗ B_22         A_12⊗ B_21         A_12⊗ B_22;         A_21⊗ B_11         A_21⊗ B_12         A_22⊗ B_11         A_22⊗ B_12;         A_21⊗ B_21         A_21⊗ B_22         A_22⊗ B_21         A_22⊗ B_22 ]


Using the above definition, the matrix of Eq. (<ref>) is found to be

    ( Υ^-1ρ_A + ρ_AΥ^-1)=
    	[ D_w⊗υ_αα+υ_αα⊗ D_w          D_w⊗υ_α d         υ_α d⊗ D_w                  0;           D_w⊗υ_dα           D_w⊗υ_dd                  0                  0;          υ_dα⊗ D_w                  0          υ_dd⊗ D_w                  0;                  0                  0                  0                  0 ].

where the elements of Υ^-1 can be found using the formula of blockwise inversion:

    Υ^-1=[  υ_αα υ_α d;  υ_dα  υ_dd ]

with

    υ_αα   = Υ_αα^-1+Υ_αα^-1D_αΥ_αα( Υ_dd -Υ_ααD_αΥ_αα^-1D_αΥ_αα)^-1Υ_ααD_αΥ^-1
    υ_α d   =-Υ_αα^-1D_αΥ_αα( Υ_dd -Υ_ααD_αΥ_αα^-1D_αΥ_αα)^-1
    υ_dα   =-( Υ_dd -Υ_ααD_αΥ_αα^-1D_αΥ_αα)^-1Υ_ααD_αΥ_αα ^-1
    υ_dd   =( Υ_dd -Υ_ααD_αΥ_αα^-1D_αΥ_αα)^-1

The inverse of the block matrix Υ_αα exists, because it is the Gramian matrix of the linear independent vectors |α_i⟩. 

For the QFIM elements we need to evaluate the inverse of the top left 3N^2 × 3N^2 part of the matrix of Eq. (<ref>) which we denote 𝔸.
In order to obtain the inverse of 𝔸, we need to further partition 𝔸 as

    𝔸=[ ε ϑ; φ ϖ ]

with

    ε   = [ D_w⊗υ_αα+υ_αα⊗ D_w ]   ϑ   = [  D_w⊗υ_α d υ_α d⊗ D_w ]
    φ   = [  D_w⊗υ_dα; υ_dα⊗ D_w ]   ϖ   = [  D_w⊗υ_dd         0;         0 υ_dd⊗ D_w ]

The inverse of ϖ is

    ϖ^-1 = 
    [  D_w^-1⊗υ^-1_dd               0;               0 υ^-1_dd⊗ D_w^-1;                 ]

The elements of 𝔸^-1 will be given by the formulas

    (𝔸^-1)_11   = (ε - ϑϖ^-1φ)^-1 = S^-1
    
    (𝔸^-1)_12   = - S^-1ϑϖ^-1
    
    (𝔸^-1)_21   =  - ϖ^-1φ S^-1
    
    (𝔸^-1)_22   = ϖ^-1+ ϖ^-1φ S^-1ϑϖ^-1

After calculations and by substituting the Υ^-1 elements from Eq. (<ref>), we derive the explicit form of 𝔸^-1 elements:

    (𝔸^-1)_11   = S^-1 = ( Υ_αα^-1⊗ D_w +D_w⊗Υ_αα^-1)^-1
     
    (𝔸^-1)_12   = S^-1[ 𝕀⊗ (Υ_αα^-1Υ_α d) (Υ_αα^-1Υ_α d) ⊗𝕀 ]
    
    (𝔸^-1)_21   = [ 𝕀⊗ (Υ_dαΥ_αα^-1); (Υ_dαΥ_αα^-1) ⊗𝕀 ] S^-1
    
    (𝔸^-1)_22   = [  D_w^-1⊗υ^-1_dd               0;               0 υ^-1_dd⊗ D_w^-1;                 ] 
    + [ 𝕀⊗ (Υ_dαΥ_αα^-1); (Υ_dαΥ_αα^-1) ⊗𝕀 ] S^-1[ 𝕀⊗ (Υ_αα^-1Υ_α d) (Υ_αα^-1Υ_α d) ⊗𝕀 ]


The QFIM elements are then obtained from Eq. (<ref>) and (<ref>)
    𝒬_ij   = 2 w_i w_j[ -2 α_iE_i       E_i       E_i ]𝔸^-1[ -2 α_jE_j;       E_j;       E_j ]
       = 2 w_i w_jE_i[ -2α_i 𝕀⊗𝕀       𝕀⊗𝕀       𝕀⊗𝕀 ]𝔸^-1[ -2α_j 𝕀⊗𝕀;       𝕀⊗𝕀;       𝕀⊗𝕀 ]E_j
       = 2 w_i w_jE_i[
    𝕀⊗Υ_dαΥ_αα^-1 + Υ_dαΥ_αα^-1⊗𝕀 - 2α_i 𝕀⊗𝕀] 
     S^-1[
    𝕀⊗Υ_αα^-1Υ_α d + Υ_αα^-1Υ_α d⊗𝕀 - 2α_j 𝕀⊗𝕀] 
    E_j
       16mu +
    2 w_i w_j
    E_i E_w^-1⊗υ^-1_dd + υ^-1_dd⊗ E_w^-1E_j
       = 2 w_i w_jE_i[
    𝕀⊗Υ_dαΥ_αα^-1 + Υ_dαΥ_αα^-1⊗𝕀 - 2α_i 𝕀⊗𝕀] 
     S^-1[
    𝕀⊗Υ_αα^-1Υ_α d + Υ_αα^-1Υ_α d⊗𝕀 - 2α_j 𝕀⊗𝕀] 
    E_j
       16mu +
    4 w_i δ_ij
    (υ^-1_dd)_ij
       = 2 w_i w_jE_i[
    𝕀⊗Υ_dαΥ_αα^-1 + Υ_dαΥ_αα^-1⊗𝕀 - 2α_i 𝕀⊗𝕀] S^-1[
    𝕀⊗Υ_αα^-1Υ_α d + Υ_αα^-1Υ_α d⊗𝕀 - 2α_j 𝕀⊗𝕀] E_j
       16mu +
    4 w_i δ_ij[ 1+α_i^2 - (Υ_ααD_αΥ_αα^-1D_αΥ_αα)_ij]


Finally, to complement the discussion in the main text,we  present some further examples of the QFIM eigenvalues for N=4,5 spurces and in Fig. <ref> we present the eigenvalues of the QFIM for 3 sources in the case of unequal weights (relative intensities) Fig.(<ref>). 




§ ANALYTIC RESULTS FOR X<<Σ



The state in the sub-diffraction regime is given by Eq. (<ref>). The derivative can be calculated immediately from this formula and it is

    ∂/∂α_iρ =-2α_i|0⟩⟨0|+ [ |0⟩⟨1|+|1⟩⟨0|] + 2α_i|1⟩⟨1|=[ -2α_i     1;     1  2α_i ]

By solving the SLD equation ∂_a_iρ_α = (ρ_α L_i+L_iρ_α), we can determine the SLDs in the {|0⟩,|1⟩} basis:

    L_i=2/(𝒞_2 -1)𝒞_2+𝒞_1^2[        𝒞_2𝒞_1 +(𝒞_2-𝒞_1^2)α_i  (𝒞_2 -1)𝒞_2+(𝒞_1-2𝒞_2𝒞_1)α_i;  (𝒞_2 -1)𝒞_2+(𝒞_1-2𝒞_2𝒞_1)α_i 𝒞_1-𝒞_1𝒞_2 +(2𝒞_1^2+𝒞_2-1)α_i ]

Knowing the SLDs, we can obtain the QFIM of Eq. (<ref>).

As already mentioned in the main text, the rank of the QFIM only depends on the matrix 

    [ M_11 M_12; M_21 M_22 ]

of Eq. (<ref>), with the elements of this matrix given by Eq. (<ref>). 
The eigenvalues μ_1, μ_2 of the matrix Eq. (<ref>) are

    μ_1   = 1/2(𝒞_2 ^2-√(((𝒞_2 -1) 𝒞_2 +4 𝒞_1 ^2-1)^2+4 ((𝒞_2 -1) 𝒞_2 +𝒞_1 ^2))-𝒞_2 +4 𝒞_1 ^2-1),
    μ_2   = 1/2(𝒞_2 ^2+√(((𝒞_2 -1) 𝒞_2 +4 𝒞_1 ^2-1)^2+4 ((𝒞_2 -1) 𝒞_2 +𝒞_1 ^2))-𝒞_2 +4 𝒞_1 ^2-1)

The condition for the eigenvalues to be zero is

    (0=3 𝒞_1 ^2-1-√((3 𝒞_1 ^2-1)^2)𝒞_2 =1/2(1-√(1-4 𝒞_1 ^2)))(0=3 𝒞_1 ^2-1-√((3 𝒞_1 ^2-1)^2)𝒞_2 =1/2(1+√(1-4 𝒞_1 ^2)))

The first part 0=3 𝒞_1 ^2-1-√((3 𝒞_1 ^2-1)^2) is always true, as it reduces to the identity (3 𝒞_1 ^2-1)^2=(3 𝒞_1 ^2-1)^2. For the second part we have

    𝒞_2 =1/2(1±√(1-4 𝒞_1 ^2)) ⇔( 2𝒞_2 -1 )^2=1-4𝒞_1^2⇔𝒞_2^2-𝒞_2 +𝒞_1^2 =0

Substituting 𝒞_2 and 𝒞_1 we get

    ( ∑_i=1^Nα_i^2)^2 -∑_i=1^Nα_i^2 +(∑_i=1^Nα_i) ^2 =0 ⇔( ∑_i=1^Nα_i^2)^2 -∑_i=1^Nα_i^2+∑_i=1^Nα_i^2  +2∑_i,j=1,i≠ j^Nα_iα_j =0 
       ( ∑_i=1^Nα_i^2)^2 +2∑_i,j=1,i≠ j^Nα_iα_j =0

Since α_i are strictly positive, except one that can be zero, this sum of positive terms cannot be equal to zero. Therefore, this statement is always false.
Thus, the Eq. (<ref>) becomes (1  0)  (1  0)=0, which means that the two eigenvalues can never be zero and the QFIM will be rank 2. 




§ CALCULATION OF THE QFI FOR THE SINC PSF

 
The expansion of the Sinc function on the HG modes is not ideal for numerical calculations.
Instead we use the spherical Bessel function of the 1^st kind and express the states onto those modes in which we then truncate.
If the PSF is a  function, the |ψ_i⟩ are

    |ψ_i⟩=1/√(σ)∫_-∞^∞( π (x-X_i)/σ) |x⟩

We can use the identity <cit.>
    ( π (x-x')/σ) =∑_q=0^∞ (2q+1) J_q( π x/σ) J_q( π x'/σ),

where J_q(x) is the spherical Bessel function of the 1^st kind. The spherical Bessel function are orthogonal in all ℝ
    ∫_-∞^∞ dx   J_q(x) J_p(x) = π/2q+1δ_qp,

therefore we can define the orthonormal basis

    |j_q⟩= √(2q+1/σ)∫_-∞^∞ dx   J_q( π x/σ) |x⟩

The set of the spherical Bessel functions is a basis in ℝ, but is not complete since it is not a resolution of identity as we can see from Eq. (<ref>).
Hence, we can expand the sinc function on the bessel function basis, using the identity Eq. (<ref>):

    |ψ_i⟩   =1/√(σ)∫_-∞^∞∑_q=0^∞ (2q+1) J_q( π x/σ) J_q( π X_i/σ) |x⟩
       =∑_q=0^∞√(2q+1)  J_q( π X_i/σ)  |j_q⟩

Using the identity for the Bessel functions

    ∂ J_q(x)/∂ x=J_q-1(x)-q+1/2J_q(x)

we can also have an expression for the derivative of |ψ_i⟩
    ∂|ψ_i⟩/∂ X_i=π/σ( J_q-1( π X_i/σ)-q+1/2J_q( π X_i/σ) )

We see that both the state ρ and its derivatives are completely expressed within the basis |j_q⟩.
This means that we can use the definition of the SLD (Eq. <ref>) and express the SLD in the same basis.

    2 ∂ρ/∂α_μ= ρ L^μ + L^μρ

In this way the fact that the specific basis is not complete does not affect our calculations.

For the numerical calculations we have to truncate our state in the appropriate amount of modes. From Figs. <ref> and <ref>, we can see that our conclusions do not change with the use of a non-Gaussian PSF. 


§ EIGENVALUES OF THE CFIM FOR SPADE AND DIRECT IMAGING


Finally, we present the eigenvalues of the CFIM for SPADE and direct imaging fir a large range of separations.