Far-field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit
FFar-field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit
Ranjith Nair ∗ and Mankei Tsang
1, 2 Department of Electrical and Computer Engineering,National University of Singapore, 4 Engineering Drive 3, Singapore 117583 Department of Physics, National University of Singapore, 2 Science Drive 3, Singapore 117551 (Dated: August 21, 2018)We obtain the ultimate quantum limit for estimating the transverse separation of two thermalpoint sources using a given imaging system with limited spatial bandwidth. We show via the quan-tum Cram´er-Rao bound that, contrary to the Rayleigh limit in conventional direct imaging, quantummechanics does not mandate any loss of precision in estimating even deep sub-Rayleigh separations.We propose two coherent measurement techniques, easily implementable using current linear-opticstechnology, that approach the quantum limit over an arbitrarily large range of separations. Ourbound is valid for arbitrary source strengths, all regions of the electromagnetic spectrum, and forany imaging system with an inversion-symmetric point-spread function. The measurement schemescan be applied to microscopy, optical sensing, and astrometry at all wavelengths.
PACS numbers: 42.30.-d, 42.50.-p, 06.20.-f
The Rayleigh criterion for resolving two incoherent op-tical point sources [1] is the most widely used benchmarkfor the resolving power of an imaging system. Accordingto it, the sources can be resolved by direct imaging onlyif they are separated by at least the diffraction-limitedspot size of the point-spread function of the imaging sys-tem. While the criterion is heuristic and does not takeinto account the intensity of the sources or the measure-ment shot noise, recent work [2–5] has made it rigorousby taking as resolution measure the classical Cram´er-Raolower bound (CRB) of estimation theory [6] on the meansquared error (MSE) of any unbiased estimate of theseparation of the sources using spatially-resolved image-plane photon counting. These works showed that if thedetected average photon number per mode N s ≪
1, theMSE of any unbiased estimator based on direct imag-ing diverges as the source separation decreases to zeroover an interval comparable to the Rayleigh limit. Thisphenomenon, dubbed Rayleigh’s curse in [7], stems fromthe indistinguishability between the photons coming fromthe two sources and imposes a fundamental limitationof direct imaging in resolving sources much closer thanthe spot size, even when the measured photon numberis taken into account. Recent developments in far-fieldmicroscopy [8] sidestep Rayleigh’s curse by preventingmultiple sources from emitting simultaneously, but con-trol over the emission properties of sources is unavailablein target sensing or astronomical imaging.While the development of novel quantum states of lightand measurement techniques has given rise to the vastfield of quantum imaging [9], fundamental quantum lim-its in resolving two incoherent sources have been largelyneglected since the early days of quantum estimation the-ory [10, 11]. Recently, the coherent [12] and incoherent[7] two-source resolution problems were revisited usingthe quantum Cram´er-Rao bound (QCRB) [11, 13] thataccounts for all (unbiased) measurement techniques al- lowed by quantum mechanics. Under a weak-source as-sumption similar to that in [2–5], it was found in [7] thatthe QCRB showed no dependence on the separation ofthe sources. Linear optics-based measurements that ap-proach the bound were also proposed [7, 14]. Subsequentdemonstrations of superresolution [15–18] have substan-tiated the feasibility of these proposals. Nevertheless,the classical treatments [2–5] and the quantum treatment[7] neglect multi-photon coincidences and bunching, phe-nomena that figure prominently in quantum optics [19].While such an approximation leads to correct conclusionsfor weak sources, e.g., at optical frequencies [20], it isproblematic for intense sources, e.g., in the microwave tofar-infrared regimes, for high-temperature astronomicalsources, and for optical demonstrations using pseudother-mal light generated from laser sources [21]. As such, aquantum-optically rigorous derivation of the resolutionlimit is as yet unavailable.In this paper, we solve these problems and obtain theQCRB for estimating the separation of two thermal pointsources of arbitrary strength using rigorous quantum op-tics and estimation theory, and show that resolution isnot fundamentally compromised at sub-Rayleigh sepa-rations. We then propose two schemes that approachthe QCRB. The finite spatial-mode demultiplexing (fin-SPADE) scheme performs photon counting in a finitenumber of suitably chosen transverse spatial modes ofthe field. The interferometric pixelated superlocalizationby image inversion interferometry (pix-SLIVER) schemeuses pixelated detector arrays in the two interferometeroutputs. The two schemes approach the QCRB overgreater ranges of the separation as the number of ac-cessed modes (fin-SPADE) or the number of pixels (pix-SLIVER) is increased.
Source and system model:
Consider two thermalpoint sources being imaged under paraxial conditions bya spatially-invariant unit-magnification imaging system a r X i v : . [ qu a n t - ph ] N ov x y x y z A + A - O I d/ - d/ d/ - d/ Imager
FIG. 1. A spatially-invariant imaging system: Point sourcesat (± d / , ) of the object plane O have images centered at (± d / , ) of the image plane I but spread out by the PSF ofthe system. (Fig. 1) – such an assumption entails no essential loss ofgenerality [22]. We assume that the system’s amplitudepoint-spread function (PSF) ψ ( ρ ) (∫ I d ρ ∣ ψ ( ρ )∣ = ) isinversion-symmetric, i.e., ψ (− ρ ) = ψ ( ρ ) , where ρ = ( x, y ) is the transverse coordinate in the image plane I . Mostimaging systems, e.g., those with circular or rectangularentrance pupils, satisfy this assumption [22].Two incoherent thermal point sources, each of effectivestrength (average photon number) N s [23], are describedby a pair of dimensionless amplitudes A = ( A + , A − ) ∈ C with the probability density [19, 24]: P N s ( A ) = ( πN s ) − exp [− (∣ A + ∣ + ∣ A − ∣ ) / N s ] . (1)In order to focus on the essential physics of the problem,we assume that the centroid (midpoint) of the sources isimaged at the optical axis and that the line joining thesources is aligned with the x -axis, so that images of thesources are centered at d ± = (± d / , ) respectively in theimage plane. Estimating the centroid of two incoherentsources by direct imaging is subject to much less strin-gent bounds than the separation [2, 7, 11] and may bedone using a portion of the available signal [7]. We alsoassume that a single quasimonochromatic temporal mode ξ ( t ) (∫ T d t ∣ ξ ( t )∣ = ) is excited over the observation in-terval [ , T ] . Extensions to multiple temporal modes canbe made using standard techniques [11].Conditioned on the value of A , the electromagneticfield in the image plane, described by the positive-frequency field operator ˆ E ( ρ , t ) [25], is in a pure coher-ent state ∣ ψ A,d ⟩ that is an eigenstate of ˆ E ( ρ , t ) with theeigenfunction ψ A,d ( ρ , t ) given by:ˆ E ( ρ , t ) ∣ ψ A,d ⟩ = ψ A,d ( ρ , t ) ∣ ψ A,d ⟩ ; (2) ψ A,d ( ρ , t ) = [ A + ψ ( ρ − d + ) + A − ψ ( ρ − d − )] ξ ( t ) , (3)where we have used the spatial invariance of the imagingsystem to write (3). The unconditional quantum state ρ d then has the P -representation [19]: ρ d = ∫ C d A + d A − P N s ( A ) ∣ ψ A,d ⟩ ⟨ ψ A,d ∣ . (4) Source separation d / < N o r m a li ze d F i s h e r i n f o r m a ti on QFI N s = 0.1QFI N s = 1QFI N s = 5DI N s = 0.1DI N s = 1DI N s = 5 FIG. 2. (Color online) The QFI of Eq. (7) (solid lines), thelower bound of Eq. (12) (dash-dotted lines) on spatially- andnumber-resolved direct imaging (DI) for the Gaussian PSF(10). The plots are normalized to the respective maximumvalues N s / σ of the QFI and are independent of the PSFhalf-width σ . Fundamental quantum bound:
The quantum Fisherinformation (QFI) K d of the state family { ρ d } determinesthe quantum Cram´er-Rao bound (QCRB) E [ ˇ d − d ] ≥ K − d (5)on the MSE of any estimator ˇ d of the separation derivedfrom an unbiased measurement POVM [11, 13, 26]. Ourderivation of K d proceeds by calculating the quantum fi-delity F ( ρ d , ρ d ) = Tr √√ ρ d ρ d √ ρ d between the (non-commuting) states (4) for two neighboring separations d and d and employing the relation K d = × lim d ,d → d − F ( ρ d , ρ d )( d − d ) (6)between the fidelity and the QFI [26, 27]. The details ofthe derivation are given in the Appendix, with the result: K d = − β ( ) N s − γ ( d ) [ ( + N s ) N s ( + N s ) − δ ( d ) N s ] , (7)where δ ( d ) = ∫ I d ρ ψ ∗ ( ρ ) ψ ( ρ − ( d, )) (8)is the overlap function of the PSF for translations in the x − direction, γ ( d ) = ∂δ ( d )/ ∂d , and β ( d ) = ∂γ ( d )/ ∂d [28].In particular, − β ( ) = ∫ I d ρ ∣ ∂ψ ( ρ ) ∂x ∣ ≡ ( ∆ k x ) , (9)the mean-squared spatial bandwidth of the PSF in the x -direction, and is independent of orientation for circular-symmetric PSFs. The first term in (7) – identical to theresult in [7] – is independent of d and dominates in the N s ≪ N s , this value is attainedin the large- d limit ( γ ( d ) → d → ∞ ) but also for x y z I d/ - d/ N Q N N Multimode fiber ˆ E ( ρ ,t ) FIG. 3. Fin-SPADE: The image-plane field is coupled intoa multimode fiber and separated into its components in theHermite-Gaussian TEM q modes of order 0 ≤ q ≤ Q by evanes-cent coupling to single-mode fibers supporting those modes.Detectors record the photon number in each mode. d =
0, so that Rayleigh’s curse is evaded. The QFI suf-fers a dip at intermediate values whose relative depthincreases with increasing N s . This is the net effect ofcorrecting the overestimation of the single-photon prob-ability and neglect of multi-photon events in the weak-source model of [7]. The QFI (7) and a lower bound onthe FI of spatially-resolved direct detection (see follow-ing) are shown in Fig. 2 for a system with the circularGaussian PSF ψ G ( ρ ) = ( πσ ) − / exp [−∣ ρ ∣ / ( σ )] , (10)for which − β ( ) = /( σ ) .Theoretical results guarantee the existence of multi-step POVMs that attain the QFI [29], but we now givetwo linear-optics schemes that closely approach it. Fin-SPADE:
For a system with the Gaussian PSF(10), consider the separation of the image-plane fieldˆ E ( ρ , t ) into its components in the TEM q Hermite-Gaussian (HG) basis [30] { ψ q ( ρ )} q with ψ G ( ρ ) ≡ ψ ( ρ ) , followed by number-resolved but not necessar-ily time-resolved photon counting over [ , T ] in eachof the modes with order 0 ≤ q ≤ Q . The coupling tothe TEM q modes can be accomplished (Fig. 3) in thesame way as SPADE [7]. Mathematically, fin-SPADEimplements a simultaneous measurement of the opera-tors { ˆ N q = ˆ a † q ˆ a q } Qq = withˆ a q = ∫ T d t ∫ I d ρ ˆ E ( ρ , t ) ψ ∗ q ( ρ ) ξ ∗ ( t ) , (11)resulting in a ( Q + ) -vector N = ( N , . . . , N Q ) T of thenumber of counts in each mode.The statistical correlations among the HG modes inthe state (4) make a direct calculation of the FI J d [ N ] of fin-SPADE difficult. We turn instead to a general lowerbound on the FI J x [ Y ] on an arbitrary parameter x ofany vector observation Y = ( Y , . . . , Y M ) T ∈ R M depend-ing on x . For µ = (⟨ Y ⟩ x , . . . , ⟨ Y M ⟩ x ) T the mean vector Source separation d / < N o r m a li ze d F i s h e r i n f o r m a ti on QFIDI Q=0Q=1Q=2 Q=3 Q=5 Q=9Q=15
FIG. 4. Fin-SPADE performance: The QFI (solid), the lowerbound (12) on the FI of fin-SPADE (dashed) for various Q ,and of direct imaging (DI) (dashed-dotted). The GaussianPSF (10) is assumed and N s = . N s / σ of the QFI andare independent of σ . The DI bound assumes a detector ofwidth 17 σ with P d =
50 pixels at 100% fill factor and is stableto increase in P d . Number-resolving unity-quantum-efficiencydetectors are assumed for all the measurement schemes. and C = ⟨( Y − µ )( Y − µ ) T ⟩ x the covariance matrix of Y evaluated at x , we have [31]: J x [ Y ] ≥ ˙ µ T C − ˙ µ , (12)where ˙ µ = ∂ µ / ∂x . Formally similar expressions have ap-peared in the quantum estimation literature [32].The mean and covariance of N in the state ρ d for thefin-SPADE measurement can be calculated using semi-classical photodetection theory [33] as detailed in the Ap-pendix. The resulting bound (12) is plotted in Fig. 4 fora representative value of N s = . d ≳ σ – in this regime, interference between the sourcesis minimal and the QCRB follows that for localizing asingle source [7, 11]. We see that measuring the first 6HG modes already achieves the quantum bound (7) overthe range d = − σ and that increasing Q widens theregion of saturation of the quantum bound. Pix-SLIVER:
Consider a PSF that is reflection-symmetric about the y − axis, i.e., ψ (− x, y ) = ψ ( x, y ) ,but otherwise arbitrary. Fig. 5 shows a schematic ofpix-SLIVER. Using an extra reflection in one arm of abalanced Mach-Zehnder interferometer, we separate theimage-plane field into its symmetric ( s ) and antisymmet-ric ( a ) components with respect to inversion of the image-plane field operator in the x -axis. The output field oper-ators areˆ E ( s ( a )) ( x, y, t ) = [ ˆ E ( x, y, t ) ± ˆ E (− x, y, t )] / + [ ˆ E v ( x, y, t ) ∓ ˆ E v (− x, y, t )] / , (13) W DA (sym) x y z I d/ - d/ BP P PR W DA (asym)B: 50:50 BeamsplitterP: Plane mirrorR: Phase retarderDA: Detector arrayB ˆ E v ( ρ ,t ) ˆ E ( ρ ,t ) FIG. 5. Pix-SLIVER: The image-plane field is separatedinto its symmetric and antisymmetric components (13) usinga balanced Mach-Zehnder interferometer with an extra reflec-tion in one arm before detecting the two outputs using iden-tical detector arrays of width W pixelated in the x -direction. where ˆ E v ( ρ , t ) is the (vacuum-state) field operator en-tering the empty port of the first beam splitter in Fig. 5.The two outputs are detected using two detector arrayspixelated along the x -direction. Each array consists of P pixels of equal x -width. To show that super-resolution ispossible without number-resolving detectors, we assume on-off detection in each pixel. For a pixel p ∈ { , . . . , P } in the α ∈ { s, a } array, such a measurement correspondsto measuring the operator ˆ K ( α ) p = f ( ˆ N ( α ) p ) , whereˆ N ( α ) p = ∫ T d t ∫ A ( α ) p d ρ ˆ E ( α ) † ( ρ , t ) ˆ E ( α ) ( ρ , t ) (14)is the total photon number operator measured over thepixel area A ( α ) p of array α , and f ( x ) = x = K = ( K ( s ) , . . . , K ( s ) P , K ( a ) , . . . , K ( a ) P ) are calculatedin the Appendix. For the Gaussian PSF (10), the lowerbound on the FI J d [ K ] of pix-SLIVER is plotted inFig. 6 for various values of P , showing how the QFIcan be approached more and more closely over the entirerange of separation values by increasing P . Discussion:
The sensitivity of our schemes at sub-Rayleigh separations can be intuitively understood as fol-lows. Information on d is encoded in the energy distribu-tion in any basis of spatial modes on I , each of which is ina thermal state. The FI of any one mode scales roughly as ∼ [ N ′ ( d )/ N ( d )] [14], for N ( d ) the mean energy in themode and N ′ ( d ) = ∂N ( d )/ ∂d , and is large if N ( d ) ∼ mode, most of the FI is contributed by theTEM mode (Fig. 4). Direct imaging is a poor way toestimate the energy in the latter, since the much largerenergy in the TEM mode acts like background noise.Similarly, in pix-SLIVER, the antisymmetric component(comprising the odd modes in any basis of modes withdefinite parity about the centroid) carries the most infor-mation at sub-Rayleigh separations (Fig. 6). Source separation d / < N o r m a li ze d F i s h e r i n f o r m a ti on QFIDI P=1 P=3P=5P=9P=15P=40symasym
FIG. 6. Pix-SLIVER performance: The QFI (solid), thelower bound (12) on the FI of pix-SLIVER using on-off de-tection with various P values(dashed lines), the lower bound(12) for DI (dash-dotted line) with number-resolved detection,and the contributions of the symmetric (sym) and antisym-metric (asym) field components to (12) for P =
40 (dottedlines). The Gaussian PSF (10) is assumed and N s = . N s / σ of the QFI and are independent of σ . The lower bounds as-sume detector array(s) of width 17 σ and 100% fill factor. TheDI bound assumes an array with P d =
50 pixels and is stableto increase in P d . While the QCRB can be approached by the maximum-likelihood (ML) estimator in the limit of a large num-ber of measurements [6], suboptimal estimators can alsoevade Rayleigh’s curse [15–18]. That small values of P achieve a substantial fraction of the QFI in pix-SLIVERis in line with work on detecting beam displacementsusing pixelated detectors [34]. The optical componentsused in pix-SLIVER have counterparts in other regionsof the electromagnetic spectrum, leading to potential ap-plications from the microwave to the gamma-ray regions[35]. Generalizations to 2D-separation estimation [36]and variants of pix-SLIVER using image inversion de-vices [37], can be envisaged. Recently developed tech-niques [38] may help to generalize our quantum limit tomultiple parameters and to unequal source strengths. Acknowledgements:
This work is supported by theSingapore National Research Foundation under NRFGrant No. NRF-NRFF2011-07 and the Singapore Min-istry of Education Academic Research Fund Tier 1Project R-263-000-C06-112.
Author contributions:
R.N. developed the sourcemodel, calculated the QFI, and invented pix-SLIVER.M.T. and R.N. bounded the FI of fin-SPADE, and R.N.applied (12) to all the detection schemes.
Note added:
During this work, we became aware ofan alternative derivation by Lupo and Pirandola [39] ofa more general quantum bound applicable to arbitraryquantum states, including our bound Eq. (7) for thermalsources as a special case. Our proposal of concrete mea-surement schemes and their near-optimality for a broadrange of source separations, however, are unique resultshere. ∗ Corresponding author: [email protected][1] Lord Rayleigh, The London, Edinburgh, and DublinPhilosophical Magazine and Journal of Science , 261(1879); M. Born and E. Wolf, Principles of Optics:Electromagnetic Theory of Propagation, Interference andDiffraction of Light (Cambridge University, 1999).[2] E. Bettens, D. Van Dyck, A. Den Dekker, J. Sijbers, andA. Van den Bos, Ultramicroscopy , 37 (1999).[3] S. V. Aert, A. den Dekker, D. V. Dyck, and A. van denBos, Journal of Structural Biology , 21 (2002).[4] S. Ram, E. S. Ward, and R. J. Ober, Proceedings ofthe National Academy of Sciences of the United Statesof America , 4457 (2006).[5] J. Chao, E. S. Ward, and R. J. Ober, J. Opt. Soc. Am.A , B36 (2016).[6] H. Cram´er, Mathematical Methods of Statistics (PMS-9) , Vol. 9 (Princeton university press, 2016); C. R. Rao,Bull. Calcutta Math. Soc. , 81 (1945); H. L. Van Trees, Detection, Estimation, and Modulation Theory : Part I ,1st ed. (Wiley-Interscience 1st Ed, 2001).[7] M. Tsang, R. Nair, and X.-M. Lu, Phys. Rev. X ,031033 (2016).[8] S. Weisenburger and V. Sandoghdar, ContemporaryPhysics , 123 (2015).[9] L. A. Lugiato, A. Gatti, and E. Brambilla, Journal of Op-tics B: Quantum and Semiclassical Optics , S176 (2002);Y. Shih, IEEE Journal of Selected Topics in QuantumElectronics , 1016 (2007); M. I. Kolobov (ed.), Quan-tum imaging (Springer Science & Business Media, 2007).[10] C. W. Helstrom, IEEE Transactions on Information The-ory , 389 (1973).[11] C. W. Helstrom, Quantum Detection and EstimationTheory (Academic Press, 1976).[12] M. Tsang, Optica , 646 (2015).[13] A. S. Holevo, Probabilistic and Statistical Aspects ofQuantum Theory (Edizioni della Normale, Pisa, Italy,2011).[14] R. Nair and M. Tsang, Optics Express , 3684 (2016).[15] Z. S. Tang, K. Durak, and A. Ling, Opt. Express ,22004 (2016).[16] F. Yang, A. Tashchilina, E. S. Moiseev, C. Simon, andA. I. Lvovsky, Optica , 1148 (2016).[17] W. K. Tham, H. Ferretti, and A. M. Steinberg, (2016),arXiv:1606.02666 [physics.optics].[18] M. Pa´ur, B. Stoklasa, Z. Hradil, L. L. S´anchez-Soto, andJ. Rehacek, Optica , 1144 (2016).[19] L. Mandel and E. Wolf, Optical Coherence and QuantumOptics (Cambridge University, 1995).[20] J. Zmuidzinas, J. Opt. Soc. Am. A , 218 (2003);M. Tsang, Phys. Rev. Lett. , 270402 (2011).[21] L. E. Estes, L. M. Narducci, and R. A. Tuft, J. Opt.Soc. Am. , 1301 (1971).[22] J. W. Goodman, Introduction to Fourier Optics , 3rd ed.(Roberts and Company Publishers, 2005), Sec. 5.3.[23] To be precise, N s is the average number of photons fromeach source reaching the image plane , allowing us to writeEqs. (2)-(3). [24] J. W. Goodman, Statistical Optics (John Wiley & Sons,1985).[25] We assume a single polarization and that the quasi-monochromatic scalar field operator ˆ E ( ρ , t ) has been castin units of √ photons ⋅ m − ⋅ s − .[26] M. Hayashi, Quantum Information (Springer, 2006).[27] S. L. Braunstein and C. M. Caves, Physical Review Let-ters , 3439 (1994).[28] Inversion symmetry of the PSF entails that δ ( d ) = δ (− d ) = δ ∗ ( d ) [28]. Note that δ ( d ) ≤ δ ( ) =
1, so that γ ( ) = Asymptotic Theory of Quantum Statis-tical Inference: Selected Papers (World Scientific, 2005);A. Fujiwara, Journal of Physics A: Mathematical andGeneral , 12489 (2006).[30] A. Yariv, Quantum Electronics (Wiley, New York, 1989).[31] M. Stein, A. Mezghani, and J. A. Nossek, IEEE SignalProcessing Letters , 796 (2014).[32] M. Hotta and M. Ozawa, Phys. Rev. A , 022327 (2004);W. Zhong, X.-M. Lu, X. X. Jing, and X. Wang, Journalof Physics A: Mathematical and Theoretical , 385304(2014).[33] J. H. Shapiro, IEEE Journal of Selected Topics in Quan-tum Electronics , 1547 (2009).[34] G. C. Knee and E. M. Gauger, Phys. Rev. X , 011032(2014).[35] R. Adam, et al. , (2015), arXiv:1502.01582 (Planck Col-laboration); K. C. Patel and S. R. Spath, Proc. SPIE , 112 (2004); M. C. Weisskopf, The Chandra X-Ray Observatory Great Science with a Great Observa-tory , Tech. Rep. MSFC-E-DAA-TN27074 (NASA Tech.Rep., 2015); T. C. Weekes,
Very high energy gamma-rayastronomy (CRC Press, 2003).[36] S. Z. Ang, R. Nair, and M. Tsang, (2016),arXiv:1606.00603 [quant-ph].[37] K. Wicker and R. Heintzmann, Optics Express , 12206(2007); K. Wicker, S. Sindbert, and R. Heintzmann,Optics Express , 15491 (2009).[38] P. Marian and T. A. Marian, Phys. Rev. A , 022340(2012); L. Banchi, S. L. Braunstein, and S. Pirandola,Phys. Rev. Lett. , 260501 (2015).[39] C. Lupo and S. Pirandola, (2016), arXiv:1604.07367[quant-ph] (to appear in Phys. Rev. Lett.)[40] C. Gerry and P. Knight, Introductory Quantum Optics (Cambridge University Press, 2005).
Fundamental quantum limit on transverse resolution
In this Section, we give the details of the derivation ofthe QCRB for separation estimation.As in Eq. (4) of the main text, the quantum state ofthe electromagnetic field in the image plane is given bythe coherent-state decomposition ρ d = ∫ C d A + d A − P N s ( A ) ∣ ψ A,d ⟩ ⟨ ψ A,d ∣ . (15)Here P N s ( A ) = ( πN s ) exp (− ∣ A + ∣ + ∣ A − ∣ N s ) , (16)is the probability density of the source field amplitudes A = ( A + , A − ) and the conditional state ∣ ψ A,d ⟩ is an eigen-vector of the image-plane field operator ˆ E ( ρ , t ) witheigenfunction ψ A,d ( ρ , t ) = [ A + ψ ( ρ − d / ) + A − ψ ( ρ + d / )] ξ ( t ) , (17)where d = ( d, ) . This eigenfunction is simply the semi-classical complex field amplitude that results from the su-perposition of the images of the two sources conditionedon the amplitude vector A .In order to evaluate the fidelity F ( ρ d , ρ d ) in Eq. (6)of the main text, we need to first choose transverse spatialmodes in which to express the quantum states ρ d and ρ d . Transverse spatial modes
For an arbitrary vector a = ( a x , a y ) in the image plane,consider the overlap function δ ( a ) ∶= ∫ I d ρ ψ ∗ ( ρ ) ψ ( ρ − a ) . (18)The Cauchy-Schwarz inequality implies that ∣ δ ( a )∣ ≤ δ ( ) =
1. We have δ ∗ ( a ) = ∫ I d ρ ψ ∗ ( ρ − a ) ψ ( ρ ) (19) = ∫ I d ρ ψ ∗ ( ρ ) ψ ( ρ + a ) (20) = δ (− a ) . (21)For an inversion-symmetric PSF, we can say more.Changing variables to σ = − ρ with d σ = d ρ , we have δ ∗ ( a ) = ∫ I d σ ψ ∗ (− σ ) ψ (− σ + a ) (22) = ∫ I d σ ψ ∗ ( σ ) ψ ( σ − a ) (23) ≡ δ ( a ) , (24)where we have used inversion-symmetry ψ (− ρ ) = ψ ( ρ ) of the PSF in the last step. For such PSFs, the overlapfunction is thus real-valued for all a ∈ I . We make theinversion-symmetry assumption throughout this paper.Since we are considering only the estimation of the x -component of the separation between the sources, weslightly abuse the above notation to define the overlapfor a scalar argument as δ ( d ) ∶= δ (( d, )) . (25)We then have δ ( d ) = δ ∗ ( d ) = δ (− d ) ≤ d of the x -separation. Consider two different values d and d of the separa-tion. For d = ( d , ) , the functions χ ( ρ ) = ψ ( ρ − d / ) + ψ ( ρ + d / )√ N χ ( ρ ) = ψ ( ρ − d / ) − ψ ( ρ + d / )√ N (27)with normalization constants given by N = + δ ( d ) , (28) N = − δ ( d ) , (29)are orthonormal over the image plane I . The functions(27) will be two of our mode functions. In like manner,for d = ( d , ) , the functions ̃ χ ( ρ ) = ψ ( ρ − d / ) + ψ ( ρ + d / )√ N (30) ̃ χ ( ρ ) = ψ ( ρ − d / ) − ψ ( ρ + d / )√ N (31)are orthonormal over the image plane with the normal-ization constants N = + δ ( d ) , (32) N = − δ ( d ) . (33)Using (26), we can readily verify that ̃ χ is orthogonalto χ and ̃ χ is orthogonal to χ . However ̃ χ is not ingeneral orthogonal to χ and neither is ̃ χ orthogonal to χ . In order to obtain an orthonormal set of transversespatial modes, the Gram-Schmidt process can be used todefine χ ( ρ ) = ̃ χ ( ρ ) − µ s χ ( ρ )√ − µ s , (34) χ ( ρ ) = ̃ χ ( ρ ) − µ a χ ( ρ )√ − µ a , (35)with µ s = ∫ I d ρ χ ∗ ( ρ )̃ χ ( ρ ) = δ [( d − d )/ ] + δ [( d + d )/ ]√N N , (36) µ a = ∫ I d ρ χ ∗ ( ρ )̃ χ ( ρ ) = δ [( d − d )/ ] − δ [( d + d )/ ]√N N . (37)The set { χ , χ , χ , χ } is an orthonormal set of trans-verse spatial modes that span the same space as { ψ ( ρ ± d / ) , ψ ( ρ ± d / )} . Note that inversion symmetry of thePSF implies that χ and χ are symmetric with respectto inversion about ρ = χ and χ are antisym-metric under inversion. Density operators ρ d and ρ d Equation (16) implies that the incoherent thermal sourceamplitudes A are circular-complex Gaussian randomvariables satisfying the relations: E [ A µ ] = E [ A µ A ν ] = E [ A ∗ µ A µ ] = N s (40) E [ A ∗+ A − ] = µ, ν ∈ {+ , −} ranging over the two sources. Define thesum and difference amplitudes S = A + + A − , (42) D = A + − A − (43)which satisfy the relations E [ S ] = E [ D ] = E [ S ] = E [ D ] = E [ SD ] = E [ S ∗ S ] = E [ D ∗ D ] = N s E [ S ∗ D ] = statistically independent circular-complexGaussian random variables. Clearly, specifying the pair ( S, D ) is equivalent to specifying A = ( A + , A − ) . Therandom variables ∣ A + ∣ and ∣ A − ∣ are independent andare both distributed exponentially with mean N s [24].Analogously, the random variables ∣ S ∣ and ∣ D ∣ are alsoindependent and are both distributed exponentially withmean 2 N s . Consider the coherent-state decomposition (15) for ρ d .Conditioned on the source amplitudes, the eigenfunction(17) can be rewritten as ψ A,d ( ρ , t ) = ⎛⎝ S √ N χ ( ρ ) + D √ N χ ( ρ ) ⎞⎠ ξ ( t ) , (45)in terms of the spatial modes defined in the previoussubsection. Since S and D are i.i.d. circular-Gaussianvariables, we may write, given the P -representation (15)[19, 33]:- ρ d = ρ th (N N s ) ⊗ ∣ ⟩ ⟨ ∣ ⊗ ρ th (N N s ) ⊗ ∣ ⟩ ⟨ ∣ , (46)where ρ th ( N ) = ∞ ∑ n = N n ( N + ) n + ∣ n ⟩ ⟨ n ∣ (47) = πN ∫ C d α exp (− ∣ α ∣ N ) ∣ α ⟩ ⟨ α ∣ (48)is the single-mode thermal state of N average pho-tons (written above in its number-state and coherent-state decompositions). The four spatiotemporalmodes in the above representation are respectively χ ( ρ ) ξ ( t ) , χ ( ρ ) ξ ( t ) , χ ( ρ ) ξ ( t ) , and χ ( ρ ) ξ ( t ) , and wehave omitted including an infinity of other spatiotempo-ral modes which are in the vacuum state for all values ofthe separation and are not useful for estimating it.Consider now the coherent-state decomposition (15)for ρ d . Conditioned on the source amplitudes, the eigen-function (17) can be rewritten as ψ A,d ( ρ , t ) = ⎛⎝ S √ N ̃ χ ( ρ ) + D √ N ̃ χ ( ρ )⎞⎠ ξ ( t ) , (49) = ⎧⎪⎪⎨⎪⎪⎩ S √ N [ µ s χ ( ρ ) + √ − µ s χ ( ρ )] + D √ N [ µ a χ ( ρ ) + √ − µ a χ ( ρ )]⎫⎪⎪⎬⎪⎪⎭ ξ ( t ) (50)The unconditional density operator ρ d can then be writ- ten in the same set of modes used for writing (46), asfollows:- ρ d = { U s [ ρ th (N N s ) ⊗ ∣ ⟩ ⟨ ∣] U † s } ⊗ { U a [ ρ th (N N s ) ⊗ ∣ ⟩ ⟨ ∣] U † a } , (51)where U s is the two-mode beam-splitter unitary (see, e.g., ref. [40]) whose action on coherent states is U s ( ∣ α ⟩ ∣ β ⟩ ) ↦ ∣ µ s α − √ − µ s β ⟩ ∣ µ s β + √ − µ s α ⟩ (52)and on the number state-vacuum product is U s ( ∣ n ⟩ ∣ ⟩ ) ↦ n ∑ k = √( nk ) µ ks ( − µ s ) n − k ∣ k ⟩ ∣ n − k ⟩ . (53)Similarly, U a is the two-mode beamsplitter unitary whoseaction on coherent states is U a ( ∣ α ⟩ ∣ β ⟩ ) ↦ ∣ µ a α + √ − µ a β ⟩ ∣ µ a β − √ − µ a α ⟩ (54)and on the number state-vacuum product is U a ( ∣ n ⟩ ∣ ⟩ ) ↦ n ∑ k = √( nk ) µ ka ( − µ a ) n − k ∣ k ⟩ ∣ n − k ⟩ . (55) State fidelity
The quantum fidelity between ρ d and ρ d is given by F ( ρ d , ρ d ) = Tr √√ ρ d ρ d √ ρ d . (56)Since both density operators (46) and (51) factorize into aproduct of density operators on the symmetric (spannedby the modes χ ( ρ ) ξ ( t ) and χ ( ρ ) ξ ( t ) ) and the anti-symmetric modes (spanned by the modes χ ( ρ ) ξ ( t ) and χ ( ρ ) ξ ( t ) ), we can mutliply the fidelities for each pair.Considering the symmetric modes first, let r ∶= N N s + N N s , (57) r ∶= N N s + N N s , (58)so that the symmetric components of the density opera-tors under each hypothesis are ρ ( sym ) d = ( − r ) ∞ ∑ n = r n ∣ n ⟩ ⟨ n ∣ ⊗ ∣ ⟩ ⟨ ∣ , (59) ρ ( sym ) d = ( − r ) ∞ ∑ n = r n U s ( ∣ n ⟩ ⟨ n ∣ ⊗ ∣ ⟩ ⟨ ∣) U † s . (60)Then √ ρ ( sym ) d ρ ( sym ) d √ ρ ( sym ) d (61) = ( − r )( − r ) ∞ ∑ n,n ′ ,n ′′ = r n + n ′′ r n ′ ∣ n ⟩ ⟨ n ∣ U s ∣ n ′ ⟩ ⟨ n ′ ∣ U † s ∣ n ′′ ⟩ ⟨ n ′′ ∣ , (62) = ( − r )( − r ) ∞ ∑ n,n ′ ,n ′′ = r n + n ′′ r n ′ ∣ n ⟩ µ n ′ s δ n n ′ µ ∗ n ′ s δ n n ′′ ⟨ n ′′ ∣ (63) = ( − r )( − r ) ∞ ∑ n = r n r n ∣ µ s ∣ n ∣ n ⟩ ⟨ n ∣ , (64)where we have used Eq. (53) to evaluate the matrix elements in Eq. (62). Consequently, F ( ρ ( sym ) d , ρ ( sym ) d ) = Tr √√ ρ ( sym ) d ρ ( sym ) d √ ρ ( sym ) d (65) = ( − r ) / ( − r ) / − ∣ µ s ∣ √ r r (66) = [√( + N N s ) ( + N N s ) − ∣ µ s ∣ √N N N s ] − . (67)In similar fashion, we find F ( ρ ( asym ) d , ρ ( asym ) d ) = [√( + N N s ) ( + N N s ) − ∣ µ a ∣ √N N N s ] − , (68)resulting in the expression F ( ρ d , ρ d ) = [√( + N s [ + δ ( d )]) ( + N s [ + δ ( d )]) − N s ∣ δ [( d − d )/ ] + δ [( d + d )/ ] ∣] − × [√( + N s [ − δ ( d )]) ( + N s [ − δ ( d )]) − N s ∣ δ [( d − d )/ ] − δ [( d + d )/ ] ∣] − (69)for the overall fidelity. Quantum Cram´er-Rao bound
Let d = d and d = d + ∆ d . The quantum Fisher infor-mation (QFI) K d on d is given by [26, 27] K d = × lim ∆ d → − F ( ρ d , ρ d + ∆ d )( ∆ d ) = − ∂ F ( ρ d , ρ d ) ∂d ∣ d = d . (70)Since the symmetric and antisymmetric modes are intensor-product states, K d is the sum of the QFIs K sym d and K asym d from the respective subsystems [26]. Defining γ ( d ) = δ ′ ( d ) ,β ( d ) = γ ′ ( d ) , (71)the QFI from the symmetric modes is found after somealgebra to be: K sym d = [ β ( d ) − β ( )] N s − N s γ ( d ) + N s [ + δ ( d )] . (72)Similarly, the QFI from the antisymmetric modes isfound to be K asym d = −[ β ( d ) + β ( )] N s − N s γ ( d ) + N s [ − δ ( d )] , (73)giving a total QFI K d = K sym d + K asym d (74) = − β ( ) N s − γ ( d ) [ ( + N s ) N s ( + N s ) − N s δ ( d ) ] , (75)which is Eq. (7) of the main text. Here β ( ) = − ∫ I d ρ ∣ ∂ψ ( ρ ) ∂x ∣ ≡ −( ∆ k x ) , (76)For circularly symmetric PSFs, this quantity is inde-pendent of the direction of the x -axis and is the mean-squared spatial bandwidth of the PSF. Fisher Information lower bounds for concretemeasurements
In this Section, we give the derivation of the lowerbound on the Fisher information for direct imaging, fin-SPADE, and pix-SLIVER. Consider a vector random variable Y = ( Y , . . . , Y M ) T ∈ R M whose probability density P Y ∣ X ( y ∣ x ) depends on anunknown parameter x . The classical Fisher information(FI) J x [ Y ] of Y on x [6] is typically difficult to computeunless the components of Y are statistically independent.However, a general lower bound J x [ Y ] ≥ ˙ µ T C − ˙ µ (77)was recently derived in [31]. Here µ =(⟨ Y ⟩ x , . . . , ⟨ Y M ⟩ x ) T is the mean observation vec-tor, C = ⟨( Y − µ )( Y − µ ) T ⟩ x is the covariance matrixof Y , and ˙ µ = ∂ µ / ∂x . All the above quantities arefunctions of x . The bound (77) is very convenient as itdepends only on the first two moments of the observationvector, which are easier to compute. In contrast, the FI J x [ Y ] depends on the full joint probability density of Y (conditioned on x ).We compute this lower bound for various measure-ments below. Since all the measurements involve atmost linear-optical processing prior to photodetection,the classicality (in the sense of having a non-negative P -representation [19, 33]) of the incoming state ρ d is pre-served. It is well known that, for such states, the quan-tum theory of photodetection gives the same quantita-tive statistics as the semiclassical theory of photodetec-tion [19, 33]. Let the input field E ( ρ , t ) be subjected toarbitrary linear-optics processing and the resulting out-put field E det ( ρ , t ) impinge on an ideal continuum pho-todetector surface. Semiclassical photodetection theorydictates that, conditioned on the source amplitudes A ,the incident field generates a space-time Poisson randomprocess at the photodetector output with the rate func-tion (or intensity) ∣ E det ( ρ , t )∣ . Unconditional statisticscan then be obtained by averaging over the source dis-tribution using (16). We will follow this approach in thesequel. Lower bound on direct imaging
Consider first the case of direct detection in the im-age plane with a pixelated detector array centered atthe origin and of width W in the x -direction. For sim-plicity, we assume it to be infinite in the y -direction,but pixelated in the x -direction with P d pixels of width W / P d . We assume ideal unity-quantum-efficiency andnoiseless number-resolved photon counting in each pixel.0Let p ∈ { , . . . , P d } be the pixel index and let pixel p be defined by the region A p = {( x, y ) ∶ l p ≤ x ≤ r p , −∞ ≤ y ≤ ∞} (78)of the image plane. The observation consists of the vector N = ( N , . . . , N P d ) T of measured counts in each pixel.Conditioned on A , the intensity function I A ( ρ , t ) inthe image plane is, using (17), I A ( ρ , t ) = ∣ ψ A,d ( ρ , t )∣ (79) = {∣ A + ∣ ∣ ψ ( ρ − d / )∣ + ∣ A − ∣ ∣ ψ ( ρ − d / )∣ + [ A ∗+ A − ψ ∗ ( ρ − d / ) ψ ( ρ + d / )]} ∣ ξ ( t )∣ . (80)The conditional photocounts N p ∣ A on the detectors p ∈{ , . . . , P d } integrated over the observation interval [ , T ] are then independent Poisson random variables with themeans µ p ∣ A = ∫ T d t ∫ A p d ρ I A ( ρ , t ) . (81)We now suppose the PSF has the Gaussian form ψ G ( ρ ) = ( πσ ) / exp (− ∣ ρ ∣ σ ) , (82)although the treatment is readily generalized to arbitraryPSFs. We obtain µ p ∣ A = ∣ A + ∣ α p + ( A ∗+ A − ) β p + ∣ A − ∣ γ p , (83)where α p = Q ( l p + d / σ ) − Q ( r p + d / σ ) ,β p = ( − d σ ) [ Q ( l p σ ) − Q ( r p σ )] , (84) γ p = Q ( l p − d / σ ) − Q ( r p − d / σ ) , and Q ( x ) = √ π ∫ ∞ x d t exp ( − t ) (85)is the Q-function.The mean photocount µ p = E [ N p ] = E A [ µ p ∣ A ] is then µ p = N s ( α p + γ p ) , (86)where we have used eqs. (38)-(41). We then have˙ µ p = ∂µ p ∂d = N s √ πσ ⎧⎪⎪⎨⎪⎪⎩ exp ⎡⎢⎢⎢⎣ − ( l p − d / ) σ ⎤⎥⎥⎥⎦ − exp ⎡⎢⎢⎢⎣ − ( r p − d / ) σ ⎤⎥⎥⎥⎦ + exp ⎡⎢⎢⎢⎣ − ( r p + d / ) σ ⎤⎥⎥⎥⎦ − exp ⎡⎢⎢⎢⎣ − ( l p + d / ) σ ⎤⎥⎥⎥⎦⎫⎪⎪⎬⎪⎪⎭ . (87)The ( p, p ′ ) -th element of the covariance matrix of N equals E [ N p N p ′ ] − µ p µ p ′ . Now E [ N p N p ′ ] = E A [ µ p ∣ A µ p ′ ∣ A ] (88) = { E A [ µ p ∣ A ] E A [ µ p ′ ∣ A ] if p ≠ p ′ E A [ µ p ∣ A ] if p = p ′ . (89)1Straightforward computations using the relations (38)-(41) and (83) give the matrix elements C p p ′ = { N s ( α p + β p + γ p ) + N s ( α p + γ p ) if p = p ′ ,N s ( α p α p ′ + β p β p ′ + γ p γ p ′ ) if p ≠ p ′ . (90)In obtaiing the elements of the covariance matrix, wehave also used the fact that E [∣ A + ∣ ] = E [∣ A − ∣ ] = N s ,which follows from the exponential statistics of ∣ A + ∣ and ∣ A − ∣ (see Sec. ). Using eqs. (87) and (90), the lowerbound (77) can be evaluated numerically for any givensystem parameters – see Figs. 2, 4, and 6 of the maintext. The limit of continuum image-plane photodetectionis achieved for P d → ∞ , but it was observed that the theFI lower bound did not change discernibly for P d ≳ P d =
50 was used in plotting the direct imaging curvesin Figs. 2, 4, and 6 of the main text.
Lower bound on fin-SPADE performance
Suppose the PSF has the Gaussian form ψ G ( ρ ) = ( πσ ) / exp (− ∣ ρ ∣ σ ) . (91)As discussed in the main text, the fin-SPADE measure-ment measures the photon number in each Hermite-Gaussian mode TEM q (with profile ψ q ( ρ ) ) of theimage-plane field for 0 ≤ q ≤ Q over the interval [ , T ] .This results in a ( Q + ) -vector N = ( N , . . . , N Q ) T of thenumber of counts in each mode. The moments of N canbe found using the semiclassical photodetection theoryas follows.Conditioned on A , the amplitude B q ∣ A in the q -th chan-nel can be written (cf. Eq. (11) of the main text):- B q ∣ A = ∫ T d t ∫ I d ρ ψ A,d ( ρ , t ) ψ ∗ q ( ρ ) ξ ∗ ( t ) . (92)As shown in [7], the integrals may be associated to theprobability amplitudes of a coherent state in the Fockbasis so that B q ∣ A = κ q / exp (− κ / )√ q ! R q , (93)where R q = { S (if q even) D (if q odd) , (94)and κ = d σ . (95) Conditioned on A , the photocounts N q ∣ A in each q -channel are independent Poisson random variables withthe means µ q ∣ A = ∣ B q ∣ A ∣ = κ q exp (− κ ) q ! ∣ R q ∣ (96) ≡ f q ∣ R q ∣ , (97)where f q is the Poisson probability of mean κ . For theunconditional mean, we have µ q ∶= ⟨ N q ⟩ = E A [ µ q ∣ A ] (98) = E A [ f q ∣ R q ∣ ] (99) = N s f q , (100)since ∣ S ∣ and ∣ D ∣ are i.i.d. random variables distributedexponentially with mean 2 N s . We also need ∂µ q ∂d = N s d σ κ q − [ q − κ ] exp (− κ ) q ! (101) = N s d σ ( f q − − f q ) , (102)where we define f − = q = q ′ , we have E [ N q ] = E A [ E [ N q ∣ A ]] (103) = E A [ f q ∣ R q ∣ + f q ∣ R q ∣ ] (104) = N s f q + N s f q , (105)where we have used the fact that N q ∣ A is Poisson-distributed. If q ≠ q ′ but q − q ′ is even, R q = R q ′ , sowe get E [ N q N q ′ ] = E A [ E [ N q ∣ A N q ′ ∣ A ]] (106) = E A [ µ q ∣ A µ q ′ ∣ A ] (107) = E A [ f q f q ′ ∣ R q ∣ ] (108) = N s f q f q ′ . (109)If q ≠ q ′ and q − q ′ is odd, E A [∣ R q ∣ ∣ R q ′ ∣ ] = E A [∣ R q ∣ ] E A [∣ R q ′ ∣ ] , so that E [ N q N q ′ ] = E A [ E [ N q ∣ A N q ′ ∣ A ]] (110) = E A [ µ q ∣ A µ q ′ ∣ A ] (111) = E A [ f q f q ′ ∣ R q ∣ ∣ R q ′ ∣ ] (112) = N s f q f q ′ . (113)Thus, the covariance matrix C of N has the qq ′ − th entry C qq ′ = ⎧⎪⎪⎪⎨⎪⎪⎪⎩ N s f q + N s f q if q = q ′ , N s f q f q ′ if q ≠ q ′ and q − q ′ is even,0 if q ≠ q ′ and q − q ′ is odd.(114)From Eqs. (101) and (114), the lower bound (77) can benumerically evaluated, as displayed in Fig. 4 of the maintext.2 Lower bound on pix-SLIVER performance
Consider the pix-SLIVER setup of Fig. 5 of the maintext with identical detector arrays in the symmetric (s)and antisymmetric (a) output ports. The overall dimen-sions of the arrays are as in Sec. , except that we consider P pixels in each array. For a conservative comparison, wetake P < P d . In addition, we also assume on-off (Geigermode) detection in each pixel, so that each component ofthe observation K = ( K ( s ) , . . . , K ( s ) P , K ( a ) , . . . , K ( a ) P ) is 0(if the corresponding pixel did not fire) or 1 (if it did). Incontrast, we allowed number-resolved detection in directimaging (see Sec. ).We now assume that the PSF is symmetric relative to reflection about the y − axis, i.e., ψ (− x, y ) = ψ ( x, y ) forall x and y – circular symmetry of the PSF is clearly asufficient condition for this to hold. Conditioned on A ,the semiclassical field amplitude in the two interferometeroutputs is given by (cf. Eq. (13) of the main text):- E ( s ( a )) A ( x, y, t ) = [ ψ A,d ( x, y, t ) ± ψ A,d (− x, y, t )] / . (115)Since the field ˆ E v ( ρ , t ) is in vacuum, the open input portof the first beam splitter does not contribute to the fieldamplitude. We can rewrite the above as E ( s ) A ( x, y, t ) = S [ ψ ( x + d / , y, t ) + ψ ( x − d / , y, t )] , (116) E ( a ) A ( x, y, t ) = D [ ψ ( x − d / , y, t ) − ψ ( x + d / , y, t )] . (117)where we have used the reflection symmetry of the PSF.The resulting conditional intensity patterns on the two detectors are I ( s ) A ( x, y, t ) = ∣ S ∣ [∣ ψ ( x − d / , y, t )∣ + ∣ ψ ( x + d / , y, t )∣ ]+ ∣ S ∣ [ ψ ∗ ( x − d / , y, t ) ψ ( x + d / , y, t )] , (118) I ( a ) A ( x, y, t ) = ∣ D ∣ [∣ ψ ( x − d / , y, t )∣ + ∣ ψ ( x + d / , y, t )∣ ]− ∣ D ∣ [ ψ ∗ ( x − d / , y, t ) ψ ( x + d / , y, t )] . (119)The integrated intensity I ( α ) p ∣ A on pixel p ∈ { , . . . , P } ofthe α ∈ { s, a } detector array over the observation interval [ , T ] is then I ( α ) p ∣ A = ∫ T d t ∫ A p d ρ I ( α ) A ( x, y, t ) . (120)Specializing to the Gaussian PSF (91), these integralsevaluate to I ( s ) p ∣ A = ∣ S ∣ [ α p + γ p + β p ] ≡ ∣ S ∣ f ( s ) p , (121) I ( a ) p ∣ A = ∣ D ∣ [ α p + γ p − β p ] ≡ ∣ D ∣ f ( a ) p , (122)where α p , γ p , and β p are defined in Eq. (84) and theabove equations serve to define the quantities { f ( α ) p } .Conditioned on A , the probability of a detector clickin the ( α, p ) -th pixel is simply the probability that oneor more photons impinge on the pixel: E [ K αp ∣ A ] ≡ µ ( α ) p ∣ A = − exp (− I ( α ) p ∣ A ) . (123)Consequently, µ ( α ) p ≡ E [ K αp ] (124) = E A [ K ( α ) p ∣ A ] (125) = − E A [ exp (− I ( α ) p ∣ A )] (126) = f ( α ) p N s + f ( α ) p N s , (127)where we have used the fact that ∣ S ∣ and ∣ D ∣ are ex-ponentially distributed with mean 2 N s to evaluate theexpectation over A . It follows that˙ µ ( α ) p = f ( α ) p N s ( + f ( α ) p N s ) , (128)for3˙ f ( s ( a )) p = √ πσ ⎧⎪⎪⎨⎪⎪⎩ exp ⎡⎢⎢⎢⎣ − ( l p − d / ) σ ⎤⎥⎥⎥⎦ − exp ⎡⎢⎢⎢⎣ − ( r p − d / ) σ ⎤⎥⎥⎥⎦ + exp ⎡⎢⎢⎢⎣ − ( r p + d / ) σ ⎤⎥⎥⎥⎦ − exp ⎡⎢⎢⎢⎣ − ( l p + d / ) σ ⎤⎥⎥⎥⎦⎫⎪⎪⎬⎪⎪⎭∓ ( d √ πσ ) exp ( − d σ ) [ Q ( l p σ ) − Q ( r p σ )] . (129)For the second moments E [ K ( α ) p K ( α ′ ) p ′ ] , three cases arise. If p = p ′ and α = α ′ , E [ K ( α ) p K ( α ′ ) p ′ ] = E [ K ( α ) p ] (130) = E A [ µ ( α ) p ∣ A ] (131) = µ ( α ) p . (132)If α ≠ α ′ (so that the pixels are in different detector arrays), the independence of S and D ensures that K ( α ) p and K ( α ′ ) p ′ are independent also so that E [ K ( α ) p K ( α ′ ) p ′ ] = µ ( α ) p µ ( α ′ ) p ′ . (133)Finally, if α = α ′ but p ≠ p ′ , E [ K ( α ) p K ( α ′ ) p ′ ] = E A [ E [ K ( α ) p ∣ A K ( α ) p ′ ∣ A ]] (134) = E A [ µ ( α ) p ∣ A µ ( α ) p ′ ∣ A ] (135) = E A [( − exp (− I ( α ) p ∣ A )) ( − exp (− I ( α ) p ′ ∣ A ))] (136) = − + f ( α ) p N s − + f ( α ) p ′ N s + + ( f ( α ) p + f ( α ) p ′ ) N s , (137)where we again use the exponential distribution of ∣ S ∣ and ∣ D ∣ to evaluate the expectation over A . Fromthese second moments, means (127), and (128), the lowerbound (77) can be numerically evaluated, with the resultsdisplayed in Fig. 6 of the main text. Note that Eq. (133) implies that the covariance matrix C is a direct sum ofmatrices for the symmetric and antisymmetric outputs,so that the lower bound (77) is also the sum of corre-sponding terms – these are shown separately in Fig. 6 ofthe main text for the case of P ==