Optimal observables and estimators for practical superresolution imaging
Giacomo Sorelli, Manuel Gessner, Mattia Walschaers, Nicolas Treps
OOptimal observables for practical super-resolution imaging
Giacomo Sorelli, Manuel Gessner, Mattia Walschaers, and Nicolas Treps
Laboratoire Kastler Brossel, Sorbonne Universit´e, ENS-Universit´e PSL,CNRS, Coll`ege de France, 4 Place Jussieu, F-75252 Paris, France (Dated: February 11, 2021)Recent works identified resolution limits for the distance between incoherent point sources. How-ever, it is often unclear how to choose suitable observables and estimators to reach these limits inpractical situations. Here, we show how estimators saturating the Cram´er-Rao bound for the dis-tance between two thermal point sources can be constructed using an optimally designed observablein the presence of practical imperfections, such as misalignment, crosstalk and detector noise.
Introduction –
Our capability of resolving small detailsin microscopy, and remote objects in astronomy is de-termined by the achievable precision of optical imagingdevices. Historical resolution limits – as those of Abbe[1] or Rayleigh [2] – address the effect of diffraction, butthey can be overcome when the signal-to-noise ratio ishigh enough [3]. In fact, the last decade provided us witha variety of super-resolution techniques to beat these lim-its by fluorescence microscopy [4–6], homodyne measure-ments [7–9] or intensity measurements in an appropriatebasis [10, 11].A paradigmatic imaging problem is the estimation ofthe separation of two incoherent point sources. By usingtools from quantum metrology [12–17], this can be op-timally solved by spatial-mode demultiplexing [18]. Ex-tensions of these results to thermal sources [19, 20], totwo-dimensional imaging [21] and (for faint sources) tomore general scenarios [22–25] are also available.Several experiments [26–29] implemented a binary ver-sion of this demultiplexing technique, that distinguishesbetween the fundamental and the first excited modes.Modern light-shaping techniques, such as multi-plane-light-conversion [30] or wave-front shaping [31], have re-cently enabled experiments demultiplexing into multiplemodes [32]. However, in all experiments performed sofar, the distance estimation was performed by measur-ing one mode at the time. To push these experimentstowards their ultimate resolution limits, it is crucial todetermine a practical estimation strategy that optimallycombines the information contained into multiple demul-tiplexed modes.In this letter, we identify such an estimation strategy.In particular, we show how an estimator for the sepa-ration between two, arbitrarily bright, thermal sourcescan be constructed from an optimized linear combinationof demultiplexed intensity measurements. Our estima-tor takes into account misalignment [18], crosstalk [33]and detector noise [34, 35], and therefore it is directlyrelevant for practical applications. Even in presence ofthe aforementioned imperfections, for faint sources, thisestimator is efficient, i.e. it saturates the Cram´er-Raobound. Moreover, we demonstrate that our strategy is,in the noiseless case, also optimal for arbitrarily brightsources, if sufficiently many modes are measured.
Source model –
We estimate the transverse separation d between two point sources located at positions ± r ,with r = ( d cos θ, d sin θ ) /
2. After propagation througha diffraction-limited imaging system (possibly with finitetransmissivity κ ) the sources are described by the modes u ( r ± r ), with u ( r ) the point spread function (PSF)of the imaging system. The modes u ( r ± r ) are nonorthogonal, therefore they cannot be used to properlyrepresent the quantum state of the sources in the imageplane. This issue is solved by introducing the orthogonalmodes u ± ( r ) = ( u ( r + r ) ± u ( r − r )) / (cid:112) ± δ ), with δ = Re (cid:90) d r u ∗ ( r + r ) u ( r − r ) , (1)the overlap between the images of the two sources. Thefield operators ˆ b ± associated with the modes u ± ( r ) canbe related to the operators ˆ s ± = (ˆ s ± ˆ s ) / √
2, with ˆ s / the field operators associated with the modes generatedby the sources, according to [20]ˆ b ± = √ κ ± ˆ s ± + (cid:112) − κ ± ˆ v ± , (2)where ˆ v ± are field operators associated with auxiliarymodes, that we can assume to be in the vacuum state,and κ ± = κ (1 ± δ ). Equation (2) allows to propagatethe quantum state ˆ ρ of the sources into its image ˆ ρ ( d, θ )after transmission through the imaging system [36].We assume the sources to be in the state ˆ ρ = ˆ ρ s ( N ) ⊗ ˆ ρ s ( N ), with ˆ ρ a ( N ) a thermal state with mean photonnumber N in the mode associated with the field operatorˆ a . Thermal sources with unequal brightnesses are dis-cussed in [36]. According to Eq. (2), ˆ ρ is mapped intoˆ ρ ( θ, d ) = ˆ ρ b + ( N + ) ⊗ ˆ ρ b − ( N − ) in the image plane, withthe information on the parameter d encoded both in theshape of the modes u ± ( r ) and the mean photon numbers N ± = N κ ± [36]. Construction of optimal observable and estimator –
Toestimate the distance d between the two thermal sources,we use the method of moments [17]. Following this ap-proach an estimator of d is obtained from the samplemean x µ = (cid:80) µi =1 x i /µ of a series of µ independent mea-surement results x i of a given observable ˆ X . In partic-ular, the separation estimator is given by the parameter a r X i v : . [ qu a n t - ph ] F e b d d s u u u ... u QQ v v v ... v QQ N N N ... N QQ + N D + N D + N D + N DQQ m m m ... m QQ ¯ x µ d (cid:104) ˆ X (cid:105) ∆ X ˜ d ± ∆ d Data acquisition Parameter estimationFIG. 1. Schematic representation of the estimation proce-dure. In the data acquisition block, the image of the twosources enters into an (eventually misaligned) demultiplex-ing device that performs a mode decomposition affected bycrosstalk. The intensity of each mode is then measured withnoisy detectors. Parameter estimation is performed (in post-processing) linearly combining, the measured intensities withoptimal coefficients, and comparing the result with a calibra-tion curve. value ˜ d for which the sample mean x µ equals the expecta-tion value (cid:104) ˆ X (cid:105) [37] of the measurement operator ˆ X . Thedependence of (cid:104) ˆ X (cid:105) on the separation d could be knowneither from theory or from a previous calibration exper-iment (see the parameter estimation block of Fig. 1).For sufficiently large values of µ , according to the cen-tral limit theorem, x µ is normally distributed with meanvalue (cid:104) ˆ X (cid:105) and variance (∆ ˆ X ) /µ . Consequently, the es-timation error is given by the error propagation formula(∆ d ) = (∆ ˆ X ) /µ (cid:16) ∂ (cid:104) ˆ X (cid:105) ∂d (cid:17) .Let us now assume that we can measure the intensityof K spatial modes v k ( r ) with associated field operatorsˆ a k . In other words, we have access to the measurement operators ˆ N = ( ˆ N , . . . , ˆ N K ), with ˆ N k = ˆ a † k ˆ a k , and wecan measure arbitrary linear combinations ˆ X ˜ m = ˜ m · ˆ N ,with ˜ m the measurement-coefficients vector. In thiscase, it is possible to perform an analytical optimiza-tion over all possible linear combinations of the accessibleoperators, i.e. to calculate the measurement sensitivity M [ d, θ, ˆ N ] = max ˜ m (cid:16) ∂ (cid:104) ˆ X ˜ m (cid:105) ∂d (cid:17) / (∆ ˆ X ˜ m ) . This optimiza-tion yields [38] M [ d, θ, ˆ N ] = D [ d, θ, ˆ N ] T Γ[ d, θ, ˆ N ] − D [ d, θ, ˆ N ] , (3)with Γ[ d, θ, ˆ N ] the covariance matrix with elementsΓ k,l [ d, θ, ˆ N ] = (cid:104) ˆ N k ˆ N l (cid:105)− N k N l , and D [ d, θ, ˆ N ] the deriva-tives vector, with components D k [ d, θ, ˆ N ] = ∂N k ∂d , wherewe denoted the mean photon number in the measurementmodes as N k = (cid:104) ˆ N k (cid:105) . The optimum given by Eq. (3) isobtained for ˜ m = m [38], with m = η Γ − [ d, θ, ˆ N ] D [ d, θ, ˆ N ] , (4)and η a normalization constant. The measurement sensi-tivity M [ d, θ, ˆ N ] obeys the following chain of inequalities M [ d, θ, ˆ N ] ≤ F [ d, θ, ˆ X m ] ≤ F Q [ d, θ ]. Here, F [ d, θ, ˆ X ] de-notes the Fisher information that bounds the sensitivityof the estimation of d from measurements of ˆ X accordingto the Cram´er-Rao lower bound (∆ d ) ≥ ( µ F [ d, θ, ˆ X ]) − [12–17]. Finally, F Q [ d, θ ] = max ˆ X F [ d, θ, ˆ X ] is the quan-tum Fisher information which determines the ultimatemetrological sensitivity [13].To calculate the quantities (3) and (4), we extend themodes u ± ( r ) to a complete orthonormal basis, and weexpand the field operators ˆ a k in this basis. By meansof this expansion, the mean photon number N k and thecorrelations (cid:104) ˆ N k ˆ N l (cid:105) can be fully expressed in terms of theexpectation values (cid:104) ˆ b †± ˆ b ± (cid:105) = N ± , (cid:104) ˆ b †± ˆ b ± ˆ b †± ˆ b ± (cid:105) = 2 N ± and (cid:104) ˆ b †± ˆ b ± ˆ b †∓ ˆ b ∓ (cid:105) = N + N − [36]. Finally, we obtain (cid:16) Γ[ d, θ, ˆ N ] (cid:17) k,l = ( N κ ) ( | f + ,k | | f + ,l | + | f − ,k | | f − ,l | + 2 Re (cid:2) f + ,k f ∗ + ,l f − ,k f ∗− ,l (cid:3) ) + δ k,l N k , (5)( D [ d, θ, ˆ N ]) k = ∂N k ∂d = 2 N κ Re (cid:18) f ∗ + ,k ∂f + ,k ∂d + f ∗− ,k ∂f − ,k ∂d (cid:19) , (6)with the mean photon number N k = N κ ( | f + ,k | + | f − ,k | ), and f ± ,k = (cid:82) d r v ∗ k ( r ) u ( r ∓ r ). Optimality of the estimator –
We now study the opti-mality of this strategy by comparing the measurementsensitivity M [ d, θ, ˆ N ] (3) with the Fisher information.To this goal, and for the rest of the letter, we consider aGaussian PSF, u ( r ) = (cid:112) / ( πw ) exp (cid:0) −| r | /w (cid:1) . More-over, let us assume ideal intensity measurements in theHermite Gaussian (HG) basis v k ( r ) = u k =( n,m ) ( r ) with u ( r ) = u ( r ). In this scenario, the covariance ma-trix (5) can be inverted analytically and it is possibleto derive exact expressions for the measurement sensitiv-ity M [ d, θ, ˆ N ] and coefficients m for an arbitrary num-ber K of measured modes [36]. In particular, if wemeasure intensity in the full HG basis, our estimationstrategy saturates the quantum Cram´er-Rao bound cal-culated in [19, 20], i.e. lim K →∞ M [ d, θ, ˆ N ] = F Q [ d, θ ][36]. When the number of received photons is low FIG. 2. Measurement sensitivity M [ d, θ, ˆ N ] for intensitymeasurements into HG modes u nm ( r ) with n, m ≤ Q =1 , K = ( Q + 1) ) including (a) misalignment ( d s / w =0 . , θ s = π/ (cid:104)| c ij | (cid:105) = 0 . σ k = 0 . ∀ k ), and (d) all three imperfections com-bined. In panel (b) and (d) solid lines and bands represent themean and one standard deviation computed over 500 crosstalkmatrices. Dashed lines show the results for ideal measure-ments, while the green solid line describes direct imaging re-sults [36]. For all plots, we assumed Nκ = 1 . θ = π/ ( N κ (cid:28) F Q [ d, θ ] = 2 N κ/w [18], accordingly small distances canbe resolved as well as large ones. Moreover, for N κ (cid:28) M [ d, θ, ˆ N ] ≈ (cid:88) k =1 N k (cid:18) ∂N k ∂d (cid:19) , (7)that coincides with the Fisher information for demulti-plexing [18], i.e. in the low brightness regime our esti-mator is efficient for arbitrary K . On the other hand,increasing the sources brightness such that N κ (cid:38) d , a reduction whichis more pronounced the larger N κ [19, 20]. For ther-mal sources of arbitrary brightness, an expression for theFisher information for finite K is unknown. Noise sources –
In the following, we describe howexperimental imperfections (whose impact on intensitymeasurements is illustrated in the data acquisition blockof Fig. 1) can be included in our model.First, the assumption that the demultiplexing modebasis v k ( r ) is perfectly centred with respect to the cen-troid of the two sources is often not true in prac-tice. In particular, an imperfect knowledge of the cen-troid position leads to a misalignment of the sourcesand the measurement basis. A two-dimensional shift r s = ( d s cos θ s , d s sin θ s ) of the sources with respect tothe ideal HG measurement basis can be readily includedin our model by substituting f ± ,k = β k ( ± r − r s ) = (cid:82) d r u ∗ k ( r ) u ( r ∓ r + r s ) in Eqs. (5) and (6). Figure 2 (a) shows that alignment must be done with a precisionof the order of the sources’ separation. When this is thecase, the impact of misalignment on the measurementsensitivity M [ d, θ, ˆ N ] can be ignored.Let us now consider the impact of imperfections in themode decomposition. In particular, we model crosstalkbetween the detection modes as a unitary matrix c kl that maps the ideal HG mode basis u l ( r ) into the ac-tual measurement basis v k ( r ) = (cid:80) l c kl u l ( r ). In prac-tically relevant scenarios, crosstalk is generally weak,namely the off diagonal elements of the coupling ma-trix are much smaller than the diagonal ones [32, 33].The overlap functions to be used in Eqs. (5) and (6)are now given by f ± ,k = (cid:80) l c kl β l ( ± r − r s ). To assessthe impact of weak crosstalk on our estimator, following[33], we numerically generate random K × K crosstalkmatrices c ij resulting into an average crosstalk probabil-ity (cid:104)| c ij | (cid:105) = (cid:104) (cid:80) Ki (cid:54) = j =1 | c ij | /K ( K − (cid:105) , where (cid:104)·(cid:105) repre-sents here an ensemble average [39]. Figure 2 (b) showsthe measurement sensitivity M [ d, θ, ˆ N ] averaged over 500random crosstalk matrices.The last noise source we consider is electronic noiseat the detection stage, i.e. dark counts. We modelthis effect, by adding to the quantum mechanical pho-ton number operators in the measurement modes ˆ N k a classical random variable ξ k which is thermally dis-tributed with mean value N ( D ) k = 2 N κσ k . The meanphoton number N (cid:48) k in the detection modes as well asthe covariance matrix Γ (cid:48) [ d, θ, ˆ N ] are now calculated notonly taking quantum mechanical expectation values, butalso classical averages over the probability distributionof ξ k . We therefore obtain N (cid:48) k = N k + N ( D ) k , and since N ( D ) k does not depend on d , the derivative vector (6) isunaffected by dark counts. On the other hand, the co-variance matrix (5) acquires an additional diagonal termΓ (cid:48) [ d, θ, ˆ N ] = Γ[ d, θ, ˆ N ] + δ kl N ( D ) k (2 N ( D ) k + 1). The influ-ence of dark counts on the sensitivity M [ d, θ, ˆ N ] is shownin Fig. 2 (c), assuming σ k to be the same for all modes.The combined effect of misalignment, crosstalk anddark counts on the measurement sensitivity M [ d, θ, ˆ N ]is shown in Fig. 2 (d) for experimentally relevant imper-fections [40]. Considering all noise sources together, wehave that misalignment and crosstalk affect the overlapfunctions f k, ± , while dark counts affect only the diago-nal of the covariance matrix. Accordingly, Eq. (7), whichcoincides with the Fisher information in the N κ (cid:28) N k with N (cid:48) k . As aconsequence, in the low brightness regime, our estimatorremains efficient even in the presence of noise.Independently of the number of received photons N κ ,Fig. 2 shows that all noise sources cause M [ d, θ, ˆ N ] tovanish for d →
0, and therefore make it harder to resolvesmall distances. However, even when all imperfectionsare combined, demultiplexing allows to outperform directimaging for small separations. The minimal distance at
FIG. 3. Smallest sources’ separation d at which ideal directimaging outperforms demultiplexing (into HG modes u nm ( r )with n, m ≤
2) as a function of the (a) misalignment, (b)crosstalk and (c) dark counts strengths for a fixed finite valueof the other two imperfections (same as in Fig. 2), and dif-ferent brightnesses Nκ = 1 . which ideal direct imaging outperforms imperfect demul-tiplexing into HG modes u nm ( r ) with n, m ≤ Optimal observable –
Let us finally discuss the observ-able that practically achieve the sensitivity bounds dis-cussed above. To this goal, in Fig. 4, we present the co-efficients m ij , Eq. (4), of the optimal linear combinationof intensity measurements in the HG modes u ij ( r ) with i, j ≤ d . Comparingthe top three panels in Fig. 4, we see that, for small sep-arations and all noise sources, the coefficients are weaklydependent on d . Accordingly, in the, arguably, most in-teresting range of small separations d , the observable thatmakes the best use of the available measurements doesnot change with the real value of the parameter. Furtherillustrations of this behaviour are provided in [36].Let us now have a look at the amplitudes of the variouscoefficients. First, the fundamental mode u containsno information on d for small separations, accordingly m = 0. In the ideal case (blue bars in Fig. 4), for every k ≤ Q all coefficients m i,k − i (with i ≤ k ) are degenerate,and their amplitude increase with k . In fact, in the ab-sence of noise, the optimal observable amplifies the smallsignals in the higher order modes to extract the most in-formation on the parameter out of them. Different noisesources modify this behaviour. Misalignment (red barsin Fig. 4) is influential for d (cid:46) d s and tends to increasehigher-order coefficients. On the contrary, both crosstalkand dark counts (green and orange bars in Fig. 4) addnoise to the higher order modes. Accordingly, the co-efficients for these modes get strongly depleted, and forsmall separation the ultimate sensitivity can be achieved FIG. 4. Optimal oefficients m ij for measurements in theHG basis u ij ( r ) with i, j ≤
1. The modes’ intensity dis-tributions are plotted below the corresponding coefficients.Different noise sources are considered: (blue) none, (red)misalignment ( d s / w = 0 . θ s = π/ (cid:104)| c ij | (cid:105) = 0 . σ k = 0 . ∀ k ).Green bars for crosstalk are averaged over 500 crosstalk ma-trices. All plots correspond to Nκ = 1 . θ = π/ by only measuring u ( r ) and u ( r ). For larger sepa-rations, even in presence of noise, the optimal observ-able gets significant contributions also from higher ordermodes (e.g. orange bars in the bottom panel of Fig. 4). Conclusions –
We presented a procedure to estimatethe separation between two thermal sources which makesan optimal use of a demultiplexed image. In the limitingcase of ideal intensity measurements into infinitely manyHG modes, our approach reaches the quantum Cram´er-Rao bound for arbitrary sources’ brightness and sepa-rations. For faint sources, the sensitivity M [ d, θ, ˆ N ] ofour method saturates the Fisher information even in thenoisy scenario. In other words, in the N κ (cid:28) d , the weak dependence of the optimalobservable on d allows to identify an estimation strategyvalid over a wide range of separations.Our results say that a good calibration is crucial toachieve optimality. In fact, the optimal coefficients inEq. (4) can be calculated either by measuring the covari-ance matrix and the derivatives vector with the help oftwo test sources, or from theory after characterizing thedifferent noise sources. In the latter case, a precise mea-sure of the crosstalk matrix of a demultiplexer, and thedark count level of detectors are not problematic. On theother hand, misalignment errors are mainly due to an im-perfect knowledge of the sources’ centroid. These errorsmay be limited by scanning the multiplexer position.Finally, we point out that with a fixed apparatus differ-ent measurement coefficients can be chosen at the estima-tion stage. Therefore, even in presence of a dynamicallychanging parameter, our approach allows to select at anytime the optimal observable in post processing. Acknowledgement –
GS acknowledges financial supportof ONERA - the French aerospace lab. MG acknowl-edges funding by the LabEx ENS-ICFP:ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*. This work was par-tially funded by French ANR under COSMIC project(ANR-19-ASTR-0020-01). This work received fundingfrom the European Union’s Horizon 2020 research and in-novation programme under grant agreement No 899587.This work was supported by the European Union’s Hori-zon 2020 research and innovation programme under theQuantERA programme through the project ApresSF. [1] E. Abbe, Beitr¨age zur Theorie des Mikroskops und dermikroskopischen Wahrnehmung, Archiv f¨ur Mikroskopis-che Anatomie , 413 (1873).[2] F. Lord Rayleigh, Investigations in optics, with specialreference to the spectroscope, The London, Edinburgh,and Dublin Philosophical Magazine and Journal of Sci-ence , 261 (1879).[3] J. W. Goodman, Statistical optics (John Wiley & Sons,2015).[4] S. W. Hell and J. Wichmann, Breaking the diffrac-tion resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy, Opt. Lett. , 780 (1994).[5] T. A. Klar, S. Jakobs, M. Dyba, A. Egner, and S. W.Hell, Fluorescence microscopy with diffraction resolutionbarrier broken by stimulated emission, Proceedings of theNational Academy of Sciences , 8206 (2000).[6] E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lind-wasser, S. Olenych, J. S. Bonifacino, M. W. Davidson,J. Lippincott-Schwartz, and H. F. Hess, Imaging intra-cellular fluorescent proteins at nanometer resolution, Sci-ence , 1642 (2006).[7] M. T. L. Hsu, V. Delaubert, P. K. Lam, and W. P.Bowen, Optimal optical measurement of small displace-ments, Journal of Optics B: Quantum and SemiclassicalOptics , 495 (2004).[8] V. Delaubert, N. Treps, C. C. Harb, P. K. Lam, and H.-A. Bachor, Quantum measurements of spatial conjugatevariables: displacement and tilt of a gaussian beam, Opt.Lett. , 1537 (2006).[9] O. Pinel, J. Fade, D. Braun, P. Jian, N. Treps, andC. Fabre, Ultimate sensitivity of precision measurementswith intense gaussian quantum light: A multimodal ap-proach, Phys. Rev. A , 010101 (2012).[10] C. Helstrom, Resolution of point sources of light as an-alyzed by quantum detection theory, IEEE Transactionson Information Theory , 389 (1973).[11] M. Tsang, Resolving starlight: a quantum perspective,Contemporary Physics , 279 (2019).[12] C. W. Helstrom and C. W. Helstrom, Quantum detectionand estimation theory , Vol. 3 (Academic press New York,1976).[13] S. L. Braunstein and C. M. Caves, Statistical distanceand the geometry of quantum states, Phys. Rev. Lett. , 3439 (1994).[14] A. S. Holevo, Probabilistic and statistical aspects of quan-tum theory , Vol. 1 (Springer Science & Business Media, 2011).[15] M. G. A. Paris, Quantum estimation for quantum tech-nology, International Journal of Quantum Information , 125 (2009).[16] V. Giovannetti, S. Lloyd, and L. Maccone, Advances inquantum metrology, Nature Photonics , 222 (2011).[17] L. Pezz`e and A. Smerzi, Quantum theory of phase es-timation, in Proceedings of the International School ofPhysics ”Enrico Fermi” , Course 188, Varenna, edited byG. M. Tino and M. A. Kasevich (IOS Press, Amsterdam,2014) pp. 691 – 741.[18] M. Tsang, R. Nair, and X.-M. Lu, Quantum theory ofsuperresolution for two incoherent optical point sources,Phys. Rev. X , 031033 (2016).[19] R. Nair and M. Tsang, Far-field superresolution of ther-mal electromagnetic sources at the quantum limit, Phys.Rev. Lett. , 190801 (2016).[20] C. Lupo and S. Pirandola, Ultimate precision bound ofquantum and subwavelength imaging, Phys. Rev. Lett. , 190802 (2016).[21] S. Z. Ang, R. Nair, and M. Tsang, Quantum limit for two-dimensional resolution of two incoherent optical pointsources, Phys. Rev. A , 063847 (2017).[22] J. ˇRehaˇcek, Z. Hradil, B. Stoklasa, M. Pa´ur, J. Grover,A. Krzic, and L. L. S´anchez-Soto, Multiparameter quan-tum metrology of incoherent point sources: Towards re-alistic superresolution, Phys. Rev. A , 062107 (2017).[23] M. P. Backlund, Y. Shechtman, and R. L. Walsworth,Fundamental precision bounds for three-dimensional op-tical localization microscopy with poisson statistics,Phys. Rev. Lett. , 023904 (2018).[24] Z. Yu and S. Prasad, Quantum limited superresolutionof an incoherent source pair in three dimensions, Phys.Rev. Lett. , 180504 (2018).[25] C. Napoli, S. Piano, R. Leach, G. Adesso, andT. Tufarelli, Towards superresolution surface metrology:Quantum estimation of angular and axial separations,Phys. Rev. Lett. , 140505 (2019).[26] M. Pa´ur, B. Stoklasa, Z. Hradil, L. L. S´anchez-Soto, andJ. Rehacek, Achieving the ultimate optical resolution,Optica , 1144 (2016).[27] Z. S. Tang, K. Durak, and A. Ling, Fault-tolerant andfinite-error localization for point emitters within thediffraction limit, Opt. Express , 22004 (2016).[28] F. Yang, A. Tashchilina, E. S. Moiseev, C. Simon, andA. I. Lvovsky, Far-field linear optical superresolution viaheterodyne detection in a higher-order local oscillatormode, Optica , 1148 (2016).[29] W.-K. Tham, H. Ferretti, and A. M. Steinberg, Beat-ing rayleigh’s curse by imaging using phase information,Phys. Rev. Lett. , 070801 (2017).[30] J.-F. Morizur, L. Nicholls, P. Jian, S. Armstrong,N. Treps, B. Hage, M. Hsu, W. Bowen, J. Janousek, andH.-A. Bachor, Programmable unitary spatial mode ma-nipulation, J. Opt. Soc. Am. A , 2524 (2010).[31] S. Rotter and S. Gigan, Light fields in complex media:Mesoscopic scattering meets wave control, Rev. Mod.Phys. , 015005 (2017).[32] P. Boucher, C. Fabre, G. Labroille, and N. Treps, Spa-tial optical mode demultiplexing as a practical tool foroptimal transverse distance estimation, Optica , 1621(2020).[33] M. Gessner, C. Fabre, and N. Treps, Superresolution lim-its from measurement crosstalk, Phys. Rev. Lett. , , 1941015 (2020).[35] C. Lupo, Subwavelength quantum imaging with noisy de-tectors, Phys. Rev. A , 022323 (2020).[36] G. Sorelli, M. Gessner, M Walschaers, N. Treps,Manuscript in preparation.[37] All expectation values are intended with respect to thestate of the sources in the image plane ˆ ρ ( θ, d ), e.g. (cid:104) ˆ X (cid:105) =Tr (cid:104) ˆ X ˆ ρ ( θ, d ) (cid:105) .[38] M. Gessner, A. Smerzi, and L. Pezz`e, Metrological non-linear squeezing parameter, Phys. Rev. Lett. , 090503(2019).[39] Here, we assumed K − dimensional cross-talk matrices, however by considering matrices of size D > K thismodel allows to account for losses into modes that arenot measured [33].[40] We assumed an alignment of the demultiplexer up to 2%of the PSF diameter ( d s / w = 0 . (cid:104)| c ij | (cid:105) =0 . σ k = 0 . Nκ = 1 . M [ d, θ, ˆ N ] for brightsources ( Nκ (cid:38)(cid:38)