Fundamental limit of resolving two point sources limited by an arbitrary point spread function

TL;DR

This paper extends the analysis of the fundamental limits of resolving two point sources beyond Gaussian PSFs to arbitrary PSFs, revealing optimal measurement modes and challenging traditional resolution limits.

Contribution

It generalizes the Fisher Information analysis to arbitrary PSFs and identifies the optimal measurement modes for superresolution imaging.

Findings

Optimal modes are sinc-Bessel for hard apertures.

Resolution limit is not constrained by Rayleigh criterion.

Exact minimum mean squared error calculations support the analysis.

Abstract

Estimating the angular separation between two incoherently radiating monochromatic point sources is a canonical toy problem to quantify spatial resolution in imaging. In recent work, Tsang {\em et al.} showed, using a Fisher Information analysis, that Rayleigh's resolution limit is just an artifact of the conventional wisdom of intensity measurement in the image plane. They showed that the optimal sensitivity of estimating the angle is only a function of the total photons collected during the camera's integration time but entirely independent of the angular separation itself no matter how small it is, and found the information-optimal mode basis, intensity detection in which achieves the aforesaid performance. We extend the above analysis, which was done for a Gaussian point spread function (PSF) to a hard-aperture pupil proving the information optimality of image-plane sinc-Bessel modes, and generalize the result further to an arbitrary PSF. We obtain new counterintuitive insights on energy vs. information content in spatial modes, and extend the Fisher Information analysis to exact calculations of minimum mean squared error, both for Gaussian and hard aperture pupils.

Figures

Equations

$$m_f(\theta )$$

$$+\frac{1}{2}{\left|{\int }_{-\infty }^{+\infty }\overline{f(x)}A(x-\theta )dx\right|}^2,$$

$$m_q(\theta )$$

$$\forall x\in \mathbb{R},\sum _{q=0}^{\infty }(1+2q)j_q(x{)}^2=1\Rightarrow \sum _{q=0}^{\infty }{m}_q(x)=1,$$

$$\mathcal{I}(\theta )=\frac{f^{\prime }(\theta {)}^2}{f(\theta )}.$$

$${\mathcal{I}}_q(\theta )$$

$$=\frac{4{\pi }^{2}N}{{\sigma }^2}(1+2q){\left(\frac{q\sigma }{\pi \theta }{j}_q\left(\frac{\pi \theta }{\sigma }\right)-j_{q+1}\left(\frac{\pi \theta }{\sigma }\right)\right)}^2.$$

$$\forall x\in I,\sum _{q=0}^{\infty }$$

$$+(1+2q)j_{q+1}(x{)}^2=\frac{1}{3}.$$

$$\forall \theta \in I,\mathcal{I}(\theta )=\sum _{q=0}^{\infty }{\mathcal{I}}_q(\theta )=\frac{4{\pi }^{2}N}{3{\sigma }^2}.$$

$$m_q(\theta )$$

$${\mathcal{I}}_q(\theta )$$

$$\mathcal{I}(\theta )=\sum _{q=0}^{\infty }{\mathcal{I}}_q(\theta )$$

$${\mathcal{I}}_{\text{Direct}}={\int }_{-\infty }^{\infty }\frac{I^{\prime }(x,\theta {)}^2}{I(x,\theta )}dx,\text{mbox}with{I}^{\prime }(x,\theta )=\frac{\partial I(x,\theta )}{\partial \theta },$$

$$\theta \in \mathbb{R},{\mathcal{I}}_S(\theta )=\frac{m_S^{\prime }(\theta {)}^2}{m_S(\theta )}\le \sum _{q\in S}{\mathcal{I}}_q(\theta )=\sum _{q\in S}\frac{m_q^{\prime }(\theta {)}^2}{m_q(\theta )}.$$

$$\exists c>0,\text{s.t.}c{m}_g(\theta )\underset{\theta \to 0^{+}}{\sim }{m}_g^{\prime }(\theta {)}^2.$$

$${\mathcal{I}}_{\text{0-BinSPADE}}(\theta )$$

$${\mathcal{I}}_{\text{1-BinSPADE}}(\theta )$$

$${\mathcal{I}}_{\text{0-BinSPADE}}(\theta )$$

$${\mathcal{I}}_{\text{1-BinSPADE}}(\theta )$$

$${\mathrm{MSE}}_{\hat{\theta }}(\theta )\ge \frac{(1-B_{\hat{\theta }}^{\prime }(\theta ){)}^2}{\mathcal{I}(\theta )}+B_{\hat{\theta }}(\theta {)}^2,$$

$$\hat{\theta }(y_0,\dots ,y_Q)=\mathrm{argmax}\left\{\prod _{q=0}^Q\mathcal{P}(y_q|\theta )\right\}.$$

$$\hat{\theta }(y_q)=m_q^{-1}\left(y_q/N\right).$$

$$\hat{\theta }(y_q,y_{q,r})=m_q^{-1}\left(y_q/(y_q+y_{q,r})\right).$$

$${\Gamma }_A(x^{\prime })={\int }_{-\infty }^{+\infty }\overline{A(x)}A(x+x^{\prime })dx.$$

$$m_{A,0}(\theta )={\left|{\Gamma }_A\left(\frac{\theta }{\sigma }\right)\right|}^2,m_{A,0,r}(\theta )$$

$${\mathcal{I}}_{\text{0-BinSPADE}}(\theta )=\frac{4N}{{\sigma }^2}\frac{\Re {\overline{{\Gamma }_A^{(1)}\left(\frac{\theta }{\sigma }\right)}{\Gamma }_A\left(\frac{\theta }{\sigma }\right)}^2}{{\left|{\Gamma }_A\left(\frac{\theta }{\sigma }\right)\right|}^2\left(1-{\left|{\Gamma }_A\left(\frac{\theta }{\sigma }\right)\right|}^2\right)},$$

$${\Gamma }_A(x)\underset{\theta \to 0}{=}1+i\beta x-\frac{\alpha }{2}{x}^2+\mathcal{O}(x^3),$$

$${\mathcal{I}}_{\text{0-BinSPADE}}(\theta )\underset{\theta \to 0}{\to }\frac{4N}{{\sigma }^2}(\alpha -{\beta }^2),$$

$${\int }_{-\infty }^{+\infty }\overline{A^{(q)}(x)}A(x+x^{\prime })dx=(-1{)}^q{\Gamma }_A^{(q)}(x^{\prime })$$

$$\frac{1}{\sqrt{\sigma }}{M}_q\left(\frac{x}{\sigma }\right)=\sum _{k=0}^q(-1{)}^k\frac{{\omega }_{k,q}}{{\sigma }^{q+\frac{1}{2}}}{A}^{(q)}\left(\frac{x}{\sigma }\right),$$

$$m_{A,q}(\theta )=\frac{1}{2}{\left|\sum _{k=0}^q\frac{\overline{{\omega }_{k,q}}}{{\sigma }^k}{\Gamma }_A^{(k)}\left(\frac{\theta }{\sigma }\right)\right|}^2+\frac{1}{2}{\left|\sum _{k=0}^q\frac{\overline{{\omega }_{k,q}}}{{\sigma }^k}{\Gamma }_A^{(k)}\left(-\frac{\theta }{\sigma }\right)\right|}^2$$

$$m_{A,1}(\theta )$$

Full text

I Introduction and Background

Consider estimating the angular separation $2\theta$ between two incoherently-radiating $\lambda$-wavelength quasi monochromatic point sources in the far field that are symmetrically disposed about the line of sight. The aperture of the camera has diameter $D$, and during the integration time the total mean photon number collected is denoted $N$. A conventional camera uses a lens in the plane of the aperture pupil to focus the image in an image plane, and detects the image-plane intensity pattern using a detector pixel array. The field amplitude in the image plane is an aperture-blurred version of the true object profile, i.e., a scaled version of convolution of the object-plane field (two independent delta functions for the above problem) with the amplitude spread function (ASF) $A(x)$ of the camera’s aperture. It is well known that no matter what $\theta$ is, the minimum mean squared error (MMSE) of estimating $\theta$ can be made arbitrarily small by using a long enough exposure, i.e., by taking $N\to\infty$. Rayleigh showed that, for a conventional camera that measures the intensity in the image plane, even if that intensity measurement is done with infinitely many infinitesimally-tiny shot-noise-limited detector pixels, when $\theta$ decreases below $\sim\lambda/D$ [1], the mean squared error (MSE) of estimating $\theta$ drastically degrades (increases) for the integration time (hence $N$) held fixed. In a recent breakthrough result, Tsang et al. showed that, assuming a Gaussian ASF, intensity measurement in the infinite Hermite-Gauss (HG) basis in the image-plane coordinates attains a Fisher Information ${\cal{I}}_{\rm HG}(\theta)$ that is independent of $\theta$ no matter how small is $\theta$, and equals the high-$\theta$ MSE attained by conventional image-plane direct detection [5]. They also showed that the quantum Fisher Information (QFI) ${\cal{I}}_{\rm Q}(\theta)$ for estimating $\theta$—which makes no assumptions on how the optical field collected by the aperture gets pre-processed and detected—equals ${\cal{I}}_{\rm HG}(\theta)$, establishing that a linear spatial mode sorting prior to detection, which separates mutually-orthogonal HG modes and detects each with individual detector pixels, is an optimal detector for this problem. This showed that Rayleigh’s criterion is an artifact of the conventional philosophy of intensity measurement in the image plane. There is rich information content in the phase of the image to extract, which optimally using a shot-noise-limited intensity measurement, one must use a non-trivial spatial-mode transformation to the aperture field prior to detection to manipulate the post-detection shot noise so as to maximize the information content about $\theta$ in the noisy detection outcomes.

This result opens up a variety of interesting questions, some important ones being: (a) what is the information-optimal mode basis for the two-point-source problem with a hard-aperture pupil (sinc ASF) and for other general ASFs, (b) what is the right minimal set of modes that carry almost all the relevant information about a passive imaging problem, (c) what is the actual advantage in MSE (Fisher Information provides a lower bound on the MMSE, via the Cramer Rao lower bound (CRLB), which is not always achievable), and (d) how does this theory generalize to more complex imaging problems, and to broadband light.

In this paper, we focus on the two-point-source problem described above. Our contributions are summarized below:

1. With a rectangular hard aperture, i.e. sinc ASF, we show that measuring in the sinc-Bessel (SB) mode basis achieves the QFI and the Fisher Information is independent of $\theta$, analogous to measuring in the HG mode basis with a Gaussian aperture.

2. We evaluate the exact MMSE of estimating $\theta$ with the optimal mode basis and compare with the CRLB.

3. We illustrate a counterintuitive distribution of energy vs. information in the individual modes of the optimal mode basis, which provide insights on efficient measurement design. In particular we find that if one were to extract and detect a single mode and its orthogonal complement (e.g., the binary SPADE of [5]), the $1^{\rm st}$ mode is optimal rather than the $0^{\rm th}$ mode, in the deep sub-Rayleigh limit. We discuss information loss due to leaky mode separation.

4. We provide the optimal binary SPADE measurement for a general ASF which attains the QFI in the low $\theta$ limit, and new insights into constructing information optimal mode bases.

II Optimal modes for hard aperture

Let $A(x)$ be the (generally complex-valued) energy-normalized Amplitude Spread Function (ASF) of the aperture, i.e., $\int_{-\infty}^{+\infty}|A(x)|^{2}\ dx=1$. The image plane field is an incoherent sum of two symmetrically-shifted copies of the ASF at $\pm\theta$, each of which is perfectly self coherent. If one projects the image plane field onto the complex-valued spatial mode $f(x)$, $\int_{-\infty}^{+\infty}|f(x)|^{2}\ dx=1$, the fraction of the intensity in the image-plane field that appears in the $f(x)$ mode is given by:

[TABLE]

where $\overline{f(x)}$ denotes complex conjugate. We will call $m_{f}(\theta)$ the measurement function for mode $f$. In [5], the authors showed that if $A(x)$ is Gaussian, then projecting the image-plane field simultaneously onto the infinite Hermite Gauss (HG) orthonormal mode basis functions $f_{q}(x)$, $q=0,1,\ldots$, attains a vector of measurements whose Fisher Information content on $\theta$ is quantum optimal, and independent of $\theta$. We will develop a similar strategy but for the practically relevant case of a space-limited (hard) aperture. This produces a cardinal sine ($\operatorname{sinc}$) ASF, i.e., $A(x)={\operatorname{sinc}}(x)$.

In order to match the energy distribution of the PSF, we will choose an orthonormal basis for which the ASF is the first basis function. This is the case of the Spherical Bessel Functions of the First Kind, $M_{q}(x)=\sqrt{1+2q}\,j_{q}(\pi x),q\in\mathbb{N}$, which are all either even or odd. Note that $j_{0}(\pi x)=\operatorname{sinc}(x)$. It is simple to verify that $\int M_{q}(x)\operatorname{sinc}(x^{\prime}-x)\ dx=M_{q}(x^{\prime})$. We will introduce a spatial scale factor $\sigma$ of dimensions of length in the image plane coordinate ($x$) to capture the actual ‘length’ of the ASF in the image plane. It will depend upon the diameter of the aperture and the focal length of the imaging system. The fraction of the total energy collected in the $q^{\rm th}$ mode is given by:

[TABLE]

From [8], equation 1.10.50 :

[TABLE]

shows that the sinc Bessel (SB) modes capture all the energy in the image-plane field. With $N$ being the total mean photon number collected over the camera’s integration time, the number of photons in the $q^{\rm th}$ SB mode is $N\,m_{q}(\theta)$.

Each of separated SB modes is detected using a shot-noise-limited detector. The total number of orthogonal temporal modes in the collected field $M\approx T(\Delta\nu)$, where $T$ is the integration time and $\Delta\nu$ is the bandwidth of the light around its center frequency. The number of photons per temporal mode $N_{0}\ll 1$ at optical frequencies. $N=MN_{0}$. It is simple to show that with $N_{0}\ll 1$ and $M\gg 1$, and with the photon statistics in the individual modes being distributed with the thermal (geometric, or Bose Einstein) distribution, that the total number of photons in the $q^{\rm th}$ spatial mode is Poisson distributed with mean $N\,m_{q}(\theta)$. In order to calculate the Fisher Information of $\theta$ in the $q^{\rm th}$ SB mode, we will rely on the following result:

Lemma 1** (Fisher Information For A Poisson Corrupted Process)**

Let $f$ be a $\mathcal{C}_{1}$ function with values in $\mathbb{R}^{+}$ mapping the variable of interest $\theta$ to a measurement $Y\sim\mathcal{P}\left(y|f(\theta)\right)$ where $\mathcal{P}$ is the Poisson distribution. Then the Fisher Information associated to the process can be written as:

[TABLE]

This expression is similar to that obtained for a Gaussian-corrupted measurement process, $\mathcal{I}(\theta)=f^{\prime}(\theta)^{2}/\eta^{2}$, where $f^{\prime}(\theta)$ is the sensitivity and $\eta^{2}$ is the noise variance. The Fisher Information in the $q^{\rm th}$ SB mode evaluates to:

[TABLE]

Lemma 2** (Series of Spherical Bessel Function Of The First Kind)**

We have the following result on a series of Spherical Bessel Function of the First Kind on any finite interval $I$ of $\mathbb{R}$ and containing [math] (see Appendix for proof):

[TABLE]

The measurement outputs on any orthogonal mode set are statistically independent Poisson random variables. So, the total Fisher Information for the vector-parametrized estimator is equal to the sum of the individual Fisher Informations from each mode. Using Lemma 2, we deduce that measuring all the SB modes leads to a $\theta$-independent Fisher Information, i.e.,

[TABLE]

In [5], it was shown that $\mathcal{I}(\theta)$ in (8) is the QFI for estimating $\theta$ with a hard aperture. Hence, we now have a proof that a specific SB mode sorting prior to direct detection achieves the QFI, in the sense that the classical Fisher Information of the SB-mode measurement exactly matches the QFI.

III Energy vs. Information Content in Modes

With the Gaussian ASF $A(x)=(2\pi\sigma^{2})^{-\frac{1}{4}}\exp(-x^{2}/4\sigma^{2})$, the measurement function and Fisher information in the $q^{\rm th}$ image-plane HG mode, $q\in\mathbb{N}$, are respectively given by [5]:

[TABLE]

The total Fisher information from measuring all the (infinitely many) HG modes equals the QFI bound for any $\theta$, i.e. [5],

[TABLE]

The Fisher Information attained by infinite-spatial-resolution image-plane direct detection is given by,

[TABLE]

where $I(x,\theta)=|A(x,\theta)|^{2}$ is the normalized spatial distribution of energy in the image plane, also the probability density function of measuring a photon at spatial position $x$, conditioned on $\theta$. $\mathcal{I}_{\text{Direct}}$ approaches the QFI in (11) for $\theta\to\infty$, but goes to zero as $\theta\to 0$. So, the information advantage of the HG mode measurement over image-plane direct detection is maximum at small $\theta$ (sub Rayleigh regime) [5].

Comparing (9) and (10) with  (2) and (6), we see several analogous trends. First, measuring all the modes (HG or SB, respectively) in either case captures all the image-plane energy for any $\theta$. Hence, any spatial mode orthogonal to the span of the respective mode sets would neither have any energy nor any information content. Further, in both cases the $q=0$ mode captures all of the energy at $\theta=0$. Finally, as noted in [6], only the $q=1$ HG mode contributes to the total Fisher information ${\mathcal{I}}(\theta)$ at low $\theta$ for the Gaussian aperture case. With a hard aperture, the $q=1$ SB mode has that exact same property.

In what follows, we will consider information contributions from individual modes and in sets of modes (without resolving modes in the set). We first consider the following lemma:

Lemma 3** (Fisher Information Inequality On Aggregated Measurements)**

Consider measurement functions $\left\{m_{q}(\theta)\right\}$ corresponding to an orthonormal family of modes. Let us say we make an aggregated measurement where we project the image plane field on to a collection of modes $S$, i.e., an effective measurement function $m_{S}(\theta)=\sum_{q\in S}m_{q}(\theta)$. This measurement cannot give us more information than the sum of the information in the individual modes in the set, i.e.,

[TABLE]

We now want to find a single mode $g$ whose information content does not go to zero (i.e., goes to a constant $c$) as $\theta\to 0$. To see the requirement on $g$, let us consider the following:

Proposition 4** (Insensitivity Property)**

Given any properly normalized and continuous mode $g(x)$ over $\mathbb{R}^{+}$ and any ASF $A(x)$, the first order derivative $m_{g}^{\prime}(\theta)$ of the associated measurement function $m_{g}(\theta)$ goes to [math] as $\theta\to 0^{+}$.

Hence, in order for the information content in mode $g$, $I_{g}(\theta)=m_{g}^{\prime}(\theta)^{2}/m_{g}(\theta)$ to go to a constant as $\theta\to 0^{+}$, the following equivalence relation must be satisfied:

[TABLE]

For the above to hold, it is necessary (but not sufficient in general) for $m_{g}\to 0$ as $\theta\to 0^{+}$. Thus, $m_{g}(\theta)$ should neither have sensitivity nor should it capture any energy at $\theta=0$. Yet, if (14) is satisfied, measuring $g$ will produce non zero information for $\theta\to 0^{+}$. Note that the insensitivity property applies to $I^{\prime}(x,\theta)$ in (12) for the direct measurement as well: it equals [math] regardless of the ASF. So an image-plane direct measurement provides no information about $\theta$ when $\theta\to 0^{+}$.

A binary SPADE (Bin-SPADE) receiver measures a single mode $g(x)$ and its orthogonal component (i.e., the leftover energy) leading to a simple implementation [5]. For both Gaussian and hard apertures, BinSPADE receivers constructed for $g(x)$ being the respective $q=0$ mode (the ASF mode) and one with the $q=1$ mode ($q=1$ HG or SB mode, respectively) attain a non-zero information at $\theta\to 0^{+}$, which equals the respective QFI limit. We will refer to these two receivers as 0-BinSPADE and 1-BinSPADE, respectively. The Fisher Information for these two measurements for the Gaussian ASF are given by:

[TABLE]

whereas, for sinc ASF (and SB modes), the respective Fisher Informations are given by:

[TABLE]

Even though both BinSPADE receivers attain the QFI for $\theta\to 0$, the former significantly outperforms the latter for higher $\theta$ (see Fig. 1). However 1-BinSPADE is much more robust to imperfect implementation. Let us say $\epsilon>0$ is a leakage parameter such that the BinSPADE receiver measures $(1-\epsilon)m_{g}(\theta)$ and its orthogonal complement. Then, absolutely any $\epsilon>0$ results in $\mathcal{I}_{\text{0-BinSPADE}}(\theta)$ to collapse to [math] at $\theta=0$. On the contrary, the performance of 1-BinSPADE at small $\theta$ remains fairly stable and leakage tolerant, with the loss in Fisher information being proportional to the energy loss in the detected mode, thus continuing to satisfy the condition stated above and hence retaining a non-zero information at $\theta\to 0^{+}$.

IV Mean Squared Error (MSE) Analysis

Let us recall the following inequality between the Mean Square Error ($\operatorname{\text{MSE}}_{\widehat{\theta}}$) of an estimator $\widehat{\theta}$, its bias $B_{\widehat{\theta}}$, and the Fisher Information ${\mathcal{I}}(\theta)$:

[TABLE]

where $B_{\widehat{\theta}}^{\prime}(\theta)$ is the first derivative of the bias function with respect to $\theta$. For comparing the MSE attained by mode-sorting receivers and image-plane detection, we will use the Normalized Root Mean Square Error, $\operatorname{\text{NRMSE}}_{\widehat{\theta}}=\sqrt{\operatorname{\text{MSE}}_{\widehat{\theta}}(\theta)}/\theta$ and the adapted Cramer-Rao bound defined as $1/\theta\sqrt{\mathcal{I}(\theta)}$ [3, 4].

For a set of independent and identically distributed (i.i.d.) measurements $\left\{y_{q}\right\},q=0,\ldots,Q$ drawn from the conditional distribution $\mathcal{P}(y_{q}|\theta)$ (as will be the case when a set of orthogonal image-plane spatial modes are detected simultaneously), the Maximum Likelihood Estimator (MLE) is given by:

[TABLE]

Assuming the prior knowledge of the total photon number $N$ collected during the integration time, the MLE for a measurement that just measures the $q(x)$ mode, is given by:

[TABLE]

For the $q$-BinSPADE receiver (i.e., one that measures mode $q$ and its orthogonal complement), the total photon number can be inferred from the sum count over the two measurements:

[TABLE]

When $N$ is high, $y_{q}+y_{q,r}$ becomes a good estimator for $N$ and the two previous expressions become similar.

In Figs. 2 and 3, we show the performance of MLE for the 0 and 1-BinSPADE receivers for the Gaussian and hard apertures, respectively, for three choices of the total number of detected photons. The MLE performs very close to the CRB at high SNR and generally follows the same trend across all $\theta$. As noted by Tsang et al. in [5], the estimator can become super efficient when it is biased, which happens at very low and large $\theta$. In both graphs, we also plot the numerically-evaluated CRB for image-plane direct measurement using (12). The MSE attained by the BinSPADEs are seen to outperform image-plane detection by a clear margin for $\theta$ much smaller than the Rayleigh separation ($\theta/\sigma=1$). Despite the superior Fisher information of the BinSPADE receivers, the NRMSE highlights the unavoidable precision decay at low $\theta$. Both the BinSPADEs and direct-measurement receivers can reach arbitrary precision in estimating $\theta$ provided they can collect a large enough photon flux (from increased exposure time). But the BinSPADEs remain more efficient than the direct measurement strategy in the sub Rayleigh regime.

V Efficient Binary SPADE For Arbitrary ASF

Consider a general (complex-valued) ASF, $A(x/\sigma)/\sqrt{\sigma}$, $\int_{-\infty}^{\infty}|A(x)|^{2}\ dx=1$, in $\mathcal{C}^{\infty}$, and its autocorrelation function,

[TABLE]

Using (1) the measurement functions associated with measuring the $A(x)$ mode and its orthogonal complement are:

[TABLE]

The Fisher Information for this generalized 0-BinSPADE receiver evaluates to:

[TABLE]

where $\Gamma_{A}^{(1)}$ is the first derivative of $\Gamma_{A}$. As for the two cases (Gaussian and hard aperture) studied before, this generalized 0-BinSPADE receiver collects all the image-plane energy in the [math] (ASF) mode for $\theta\to 0^{+}$. This follows from the peak property of the autocorrelation function and the normalization of the ASF, which gives us $m_{A,0}(0)=1$. It is interesting to note that $\mathcal{I}_{\text{0-BinSPADE}}(\theta)$ is independent of any phase present in the aperture function and is solely based on its intensity profile. Thus, most optical aberrations such as defocusing, spherical aberration or coma, among others, if taken into account into the projection mode, do not degrade the information, unlike the image-plane direct measurement.

Assuming $\Gamma_{A}(x)$ admits a second order expansion near $\theta=0$, one can write a Taylor series as follows:

[TABLE]

with $\alpha\geq 0$, and $\beta\in\mathbb{R}$. One can further show that:

[TABLE]

which is a $\theta$-independent constant as before. Note that a real ASF benefits from the fact that $\beta=0$. Also, from the Wiener-Khinchin theorem, increasing $\alpha$ is equivalent to greater spatial variations in the ASF profile, i.e., sharper or numerous edges. It is simple to verify that (27) reduces to the corresponding formulas for the QFI for the Gaussian ($\alpha=1/4$, $\beta=0$) and hard rectangular apertures ($\alpha=\pi^{2}/3$, $\beta=0$), attained by the respective [math]-BinSPADE receivers as discussed above.

Let us now construct an orthonormal basis for a general ASF by using derivatives of $A(x)$. Let us first note the following identity relating the $q^{\text{th}}$ derivative of $A(x)$ to that of $\Gamma_{A}$:

[TABLE]

Clearly, the functions $(A^{(q)})_{q}$ need not be orthogonal, and hence cannot be used for parallel measurements. We do a Gram-Schmidt orthogonalization by selecting weights $\omega_{k,q}\in\mathbb{C}$ such that each measurement mode $M_{q}(x)$ can be written as a linear combination of the $(A^{(q)})_{q}$ functions,

[TABLE]

such that $\left\{M_{q}\right\}$ forms an orthonormal basis (after removing some possibly identically null functions). It is easy to verify that $\omega_{0,0}=1$ (due to the energy normalization property of $A(x)$) and hence the first mode ($M_{0}$) simply equals the ASF. The measurement functions for a simultaneous measurement of the $M_{q}$ modes can be expressed as:

[TABLE]

We consider a generalized $1$-BinSPADE receiver associated with the $M_{1}$ mode. The corresponding measurement function is given by:

[TABLE]

and the associated Fisher Information is given by:

[TABLE]

Assuming $\Gamma_{A}(x)$ can be expanded as in (26), we find that the generalized $1$-BinSPADE receiver attains the same Fisher Information as that of the generalized [math]-BinSPADE receiver: $4N/\sigma^{2}(\alpha-\beta^{2})$, at $\theta\to 0^{+}$. As discussed before in the context of Gaussian and hard rectangular apertures, a similar robustness advantage (for 1 over 0-BinSPADE) exists in the general case. Whether or not $\sum_{q=0}^{\infty}{\mathcal{I}}_{q}(\theta)=4N/\sigma^{2}(\alpha-\beta^{2}),\forall\theta$ holds for a general complex-valued $A(x)$ with our construction of the measurement modes, is left open. In future work, it will also be interesting to investigate a systematic generalization of the optimal modes for a pre-detection mode-sorting based receiver for more complex imaging problems, and to prove that the quantum Fisher Information limit for any incoherent-light imaging problem can be attained by an appropriate receiver that applies a pre-detection linear mode transformation.

Appendix

Proof:

Let $f$ be a function $\mathbb{R}\rightarrow\mathbb{R}^{+}$ modeling the output of a process on a variable $\theta$ that is corrupted by Poisson noise:

[TABLE]

such that $\ln(p_{f}(y|\theta))$ is twice differentiable with respect to $\theta$. We note $f^{\prime}$ and $f^{\prime\prime}$, respectively its first and second derivatives. Then we can write the corresponding Fisher Information as:

[TABLE]

∎

Proof:

If we consider the output of multiple independent measurement functions $(f_{0}(\theta),\ldots,f_{Q}(\theta))$, all twice differentiable with respect to $\theta$, and their respective outputs $y_{0},\ldots,y_{Q}$, we have for the Fisher Information:

[TABLE]

∎

Proof:

We consider the two following series of functions over a finite interval $I$ of $\mathbb{R}$ containing zero, for some fixed positive integer $b$ :

[TABLE]

We have the following loose upper bound on the Spherical Bessel Functions, from their series definition:

[TABLE]

We can pick $X$ such that, $\forall x\in I,\ |x|\leq X$ and we can test that, for the first series (40) we have the convergence bound:

[TABLE]

As the underlying series is positive, monotonic and converges (by the ratio test), we can choose a $P$ large enough to reduce the remainder, and subsequently the bound, to any small $\epsilon$ we wish. In other words:

[TABLE]

We can thus conclude that the function series is uniformly Cauchy and converges uniformly over $I$.

For the second series (41) we have:

[TABLE]

where the same conclusion applies.

Finally, we consider the new set of series of functions, related respectively to $A_{b,Q}$ and $B_{b,Q}$:

[TABLE]

We note that their point-wise convergence can easily be established at $x=0$, as $\forall q>0,\ j_{q}(0)=0,\ j_{0}(x)=1$ and, considering the previous results, we can apply the Differentiation Theorem to obtain their respective uniform convergence as well as the relations:

[TABLE]

∎

Proof:

We will use the following two results from [8] (Equations 1.10.50 and 1.10.52), where $\operatorname{Si}(x)$ denotes the Sine Integral.

[TABLE]

By combining them, we obtain:

[TABLE]

Thanks to the Differentiation results of the previous lemma (51), we can write from the derivative of (53):

[TABLE]

We then find the following as the solution of a first order linear differential equation involving (55) and (57) as well as the property (52) and the constraint that the series is equal to 0 at $x=0$:

[TABLE]

With this result, we can compute the derivative of (55) to express the following series and its respective shifted version:

[TABLE]

The next series is also found as the solution of another differential equation involving (61) and (62), with the same value constraint at $x=0$ than previously:

[TABLE]

And after deriving (61) we obtain:

[TABLE]

Finally, we can combine the results (57) to (64) to get:

[TABLE]

∎

Proof:

Considering any continuous and differentiable ASF $A(x)$ as well as mode $g(x)$, the measurement function $m_{g}$ for two point sources separated by $2\theta$ can be written as:

[TABLE]

We have for the limit of its first derivative:

[TABLE]

We can proceed similarly for the Fisher Information in the case of direct detection. We recall the expression:

[TABLE]

In the case of two point sources, the normalized intensity profile can be written with the previous ASF as : $I(x)=(|A(x+\theta)|^{2}+|A(x-\theta)|^{2})/2$. One can write for its first derivative:

[TABLE]

∎

Proof:

Let $m_{S}$ be an aggregated measurement over a collection $S$ of measurement functions from orthogonal modes, i.e. $m_{S}(\theta)=\sum_{q\in S}m_{q}(\theta)$, all corrupted by Poisson noise. At a location $\theta$ where $\forall q,\ m_{q}(\theta)>0$, we have the following measurement function and Fisher Information:

[TABLE]

We can simplify the notations for the current location into sets $(m_{q})$ and $(m^{\prime}_{q})$ and observe that:

[TABLE]

Here, if $m^{\prime}_{q}$ or $m^{\prime}_{r}$ are equal to zero or their product is negative, and with the previous positivity constraint, the sum is clearly positive. Otherwise, one can write:

[TABLE]

As we have, $m_{r}m_{q}^{\prime}/m_{q}m_{r}^{\prime}>0$ and $\forall x>0,\ x+1/x\geq 2$. Thus, (72) is always positive and we can conclude with the inequality:

[TABLE]

∎

Remark 5** (Generalized projection modes for a BinSPADE receiver)**

One can expand the generalized mode expression (29) in the case of the two first modes for a normalized ASF $A(x)$ into :

[TABLE]

where $A^{(1)}(x)$ is the first derivative of the ASF, $\Gamma_{A}$ is the autocorrelation of $A$ and $\Gamma_{A}^{(1)}$, $\Gamma_{A}^{(2)}$ are respectively its first and second derivatives. Here, one can note, as the autocorrelation is hermitian, that $\Gamma_{A}^{(1)}(0)$ is purely imaginary. Thus in the case of a purely real ASF it is equal to zero and the first mode can be constructed with only $A^{(1)}(x)$.

Remark 6** (Fisher Information of a leaky BinSPADE receiver)**

Considering a leakage factor $0<\epsilon<1$, we define the efficiency of a BinSPADE receiver as $\rho=1-\epsilon$. We have the following expressions of the Fisher Information for 0 and 1-BinSPADE respectively; first, for a Gaussian aperture:

[TABLE]

For a rectangular aperture:

[TABLE]

Finally, for a generic ASF $A(x)$, properly normalized and in $\mathcal{C}^{2}$:

[TABLE]

[TABLE]

One can notice from the common normalization of the ASF $A(x)$ that the autocorrelation $\Gamma_{A}(\theta)$ is equal to one at $\theta=0$; and as it is also hermitian, that $\Gamma_{A}^{(1)}(0)$ is purely imaginary. Hence, in (81) the numerator always tend to zero for $\theta\rightarrow 0$, but the denominator only does so if the efficiency $\rho$ is exactly equal to one. Thus, for any $\rho<1$, the generalized 0-BinSPADE does not collect information for narrow separation angles : $\mathcal{I}_{\text{0-BinSPADE},\rho<1}\rightarrow 0$ near [math].

On the other hand, with the expansion (26) one can develop the limit of (82) as :

[TABLE]

Thus, the generalized 1-BinSPADE is resilient to leakage in the deep sub-Rayleigh limit.