Imaging arbitrary incoherent source distributions with near quantum-limited resolution

We demonstrate an approach to obtaining near quantum-limited far-field imaging resolution of incoherent sources with arbitrary distributions. Our method assumes no prior knowledge of the source distribution, but rather uses an adaptive approach to imaging via spatial mode demultiplexing that iteratively updates both the form of the spatial imaging modes and the estimate of the source distribution. The optimal imaging modes are determined by minimizing the estimated Cramér-Rao bound over the manifold of all possible sets of orthogonal imaging modes. We have observed through Monte Carlo simulations that the manifold-optimized spatial mode demultiplexing measurement consistently outperforms standard imaging techniques in the accuracy of source reconstructions and comes within a factor of 2 of the absolute quantum limit as set by the quantum Cramér-Rao bound. The adaptive framework presented here allows for a consistent approach to achieving near quantum-limited imaging resolution of arbitrarily distributed sources through spatial mode imaging techniques.

www.nature.com/scientificreports/ Optimal imaging methodology. The setup of the adaptive spatial mode imager, shown in Fig. 1, consists of a two-dimensional source region F(R) emitting light which is imaged in the far-field through a hard aperture. At the image plane, the field is sorted into a set of orthonormal spatial modes � = {φ j (r)|j = 0, . . . , J − 1} , and measured by individual photon counting detectors. In this work, we assume a photon source emitting spatially incoherent light with detection events that can be modeled as a Poisson process 6 . The photon counts in each channel are then used to reconstruct a source estimate, as well as to update the optimal set of imaging spatial modes for the next measurement iteration. We note that restriction to incoherent sources is not fundamental to the approach, and extension to coherent sources is also possible.
To optimize the imaging for arbitrary incoherent source distributions, we parametrize the source by decomposing it into a set of orthogonal functions: where F(R) is the source brightness in the object plane, and f k (R) is the k th function in the source decomposition basis set. We will refer to these functions as source modes, to distinguish them from the spatial imaging modes. To make the problem computationally tractable we restrict the number of source modes to a finite value K, such that F(R) can be adequately approximated in its subspace. The imaging problem then becomes one of accurately estimating the multiparameter coefficient vector c = (c 1 , c 2 , . . . , c K ) in the presence of noise. We restrict our attention to the case where the imaging noise is dominated by the photon shot noise and other noise sources are negligible, as is often the case for photon detection in the visible spectrum.
To properly frame this measurement problem, we use the statistical tools of the Fisher information matrix and the corresponding Cramér-Rao bound, which sets the lowest achievable mean squared error (MSE) of the coefficient estimates of any unbiased estimator. The classical Fisher information per photon measured with the orthonormal imaging modes is constructed as 6 where P j is the probability, conditioned on a detection event, of detecting a photon in imaging mode φ j . For incoherent source distributions this probability is given by where |ψ PSF R � is the field point spread function (PSF) in the image plane originating from source point R . Replacing F(R) with its orthogonal decomposition allows us to directly evaluate the partial derivative terms. The covariance matrix of the coefficient vector for any unbiased estimator satisfies the inequality known as the Cramér-Rao bound. Specifically, the MSE of any unbiased estimator of the coefficient c k is bounded by the corresponding diagonal entry of the Fisher information matrix inverse: Adaptive spatial mode imaging system. Light from the far-field source distribution F(R) is captured by a hard aperture and imaged onto a spatial mode sorter of imaging modes . At each measurement iteration, the number of photons detected in each channel n j is used to generate a source estimate, as well as an updated set of optimal spatial imaging modes for the next measurement iteration. www.nature.com/scientificreports/ The optimal imaging problem can then be cast into one of minimizing the MSE of the source coefficients by minimizing the diagonal elements of the Fisher information matrix inverse. The objective function to be minimized is where W is a weighting matrix that can be applied depending on the nature of the imaging task. In the results shown here, we choose W to be equal to the identity matrix, corresponding to minimizing the mean of the diagonal elements of the Cramér-Rao bound, CRB = 1/K K k=1 [I] −1 kk . The objective function L(�) is minimized over the manifold of all possible orthonormal imaging modes using numerical manifold optimization methods 35 . Our specific implementation uses the Scaled Gradient Projection Method (SGPM) designed for the minimization of functions over the set of orthonormal matrices, known as the Stiefel manifold 36 .

Results
Optimized imaging of point sources and extended sources. Before demonstrating the full adaptive imaging method, we first characterize the maximum imaging improvement attainable assuming the true source distribution is known for the purposes of computing the Fisher information matrix in Eq. (2). We first apply the MO-SPADE imaging method to one-dimensional imaging of a line of incoherent point sources, where the locations of the sources are known a priori and the task is to estimate the magnitudes of the point sources. This scenario applies, for instance, in estimating the brightness of poorly resolved astronomical objects (e.g., exoplanets), whose locations are known from indirect methods. The source modes are then Dirac delta functions located at the point sources. For simplicity we assume a shift-invariant Gaussian field PSF of where we have normalized the image plane spatial coordinate x by the magnification factor of the imaging system. The width σ of the Gaussian PSF is approximately related to the radius r c of a corresponding Airy disk 37 by σ ≈ r c /3. Figure 2 displays the results of MO-SPADE imaging of five equally spaced point sources using mode sorting detection with five spatial imaging modes. We find that for point source spacings x in the sub-Rayleigh regime, MO-SPADE imaging attains a lower CRB compared to direct imaging, with an advantage that grows larger as x → 0 . In terms of resolution, MO-SPADE allows for smaller x point source spacings before reaching the equivalent CRB value of direct imaging, showing an increase in effective resolution by over a factor of 1.5 in the deep sub-Rayleigh regime. The spatial imaging modes which attain this performance are shown in Fig. 2b. Their profiles are influenced by the source distribution and PSF, and become increasingly complex for larger numbers of imaging modes and sub-Rayleigh source modes. We confirmed the performance enhancement via Monte Carlo simulations for �x = 0.3σ , collecting a mean value of N photons for each trial (Fig. 2d). The non-negative least squares (NNLS) estimates of the point source amplitudes follow the CRB limits as expected for large values of N , and confirms that MO-SPADE imaging achieves an order of magnitude reduction in the MSE of the source reconstruction. For low photon numbers, the MSE of the NNLS estimates falls below the CRB due to its nonzero bias, saturating at a maximum level where the estimation error is of approximately the same magnitude as the source amplitudes.
MO-SPADE imaging also shows enhanced performance when imaging arbitrary extended objects. The source region is approximately discretized into a set of non-overlapping rectangle functions with width a. We have chosen rectangle functions for simplicity, but other basis choices may be more appropriate depending on the imaging task. The optimization of the spatial imaging modes then proceeds in an identical way as with the point source demonstration. Figure 3 shows the results of MO-SPADE imaging for three different example source distributions of finite extent, including a uniform distribution as well as sources with smoothly varying amplitude variations. MO-SPADE imaging achieves more than a 10-fold reduction in the CRB of the extended www.nature.com/scientificreports/ source amplitude estimation in the deep sub-Rayleigh regime, analogous to the point source imaging task. In both of these imaging scenarios, the direct imaging CRB converges towards the optimal CRB for large a or x , in agreement with previous studies on two point sources.
Comparison to the quantum limit. We also compare our results to the QCRB, which sets a lower bound on the CRB that can be achieved for any possible measurement of a quantum system. We compute the QCRB on an extended source distribution by taking the inverse of the quantum Fisher information matrix is the symmetric logarithmic derivative (SLD) of density matrix ρ computed as 6 : where q and |e q � are the eigenvalues and eigenvectors of ρ . This version of the quantum Cramér-Rao bound is sometimes referred to as the SLD-CRB, and additional details of the computation are provided in the Supplementary Information. The multi-parameter QCRB is in general not guaranteed to be attainable by a physically realizable measurement 38 . The conditions under which this bound is attainable have been studied in several recent works [39][40][41][42][43][44] . If the SLD operators commute, that is [L i , L j ] = 0 for all i, j, then the bound is attainable with a projective measurement over the common eigenbasis of the SLD operators. If the SLD operators do not commute, then the CRB of all the parameters cannot be saturated simultaneously with a single projective measurement www.nature.com/scientificreports/ basis. However, if the system satisfies the condition tr(ρ[L i , L j ]) = 0 , known as weak commutativity, then the quantum bound is still attainable, although collective measurements over multiple copies of the system may be necessary 40 . In the case of far-field imaging of incoherent sources, the weak commutativity condition holds and the bound is always attainable in principle with the aid of collective measurements 24 . Figure 2c as well as Fig. 3(d-f) show the comparison between the mean CRB of the MO-SPADE imaging and the mean QCRB. We find that in the well-resolved regime the MO-SPADE CRB follows closely the quantum CRB, while in the sub-Rayleigh regime the ratio between the MO-SPADE and quantum CRB approaches, but never exceeds, a factor of 2, as shown in the upper panels of Fig. 3(d-f), which plots the ratio between the two CRBs. These results suggest that for far-field incoherent imaging, collective measurements can provide at most a factor of 2 improvement in the CRB compared to the optimal projective measurements as provided by MO-SPADE imaging, even in the deep sub-Rayleigh regime.
Adaptive imaging on unknown source distributions. So far in this work, the CRB and optimal imaging modes have been derived under the assumption of a known source distribution. In most practical cases, the true source distribution is unknown, and it is the purpose of the imaging measurement itself to estimate the source distribution. To solve this problem, an iterative, adaptive approach to imaging is proposed, where the measurement duration is divided into multiple measurement periods, and the source distribution and optimal MO-SPADE imaging modes are alternately estimated. At each adaptive iteration i the estimated Fisher information is calculated using the source estimate c est i−1 from the previous iteration along with both previous and current imaging modes, weighted according to the relative number of collected photons N i , The set of orthonormal imaging modes, i , are optimized by minimizing the estimated CRB: To initialize the adaptive iterations, we set 0 as the direct imaging basis to take the first measurement. At each iteration, all previous measurement data is used in calculating the estimated source coefficients. By following this iterative scheme, we obtain progressively better estimates of the source distribution, and simultaneously approach the optimal imaging modes for the true source distribution. We note that since previous measurement information collected with different imaging bases is included in the source reconstructions, the number of adaptive imaging modes in each iteration can be less than the total number of sources modes. This can potentially allow for large improvements in the imaging while only using a relatively small number of adaptive imaging modes. Additional details on the adaptive algorithm are provided in the Supplementary Information.  www.nature.com/scientificreports/ To validate the adaptive approach, we performed a Monte Carlo simulation of the reconstruction performance using the NNLS estimator for the one-dimensional source shown in Fig. 3b, discretized into 28 source mode rectangles with width a = 0.8σ . In this simulation, only 8 imaging modes are used for each adaptive step, and the measurement duration is split equally among the initial direct imaging step and two adaptive MO-SPADE measurements-a configuration that we have empirically found to work well. The results of this simulation are shown in Fig. 4. At low photon levels, the bias of the NNLS estimator allows the MSE to fall below the respective CRB bounds, and saturates at a maximum error value, as was the case for the point source estimation problem. As the number of photons increases, the MSE begins to follow the CRB, where this trend continues until the approximation of the source by discrete rectangular modes imposes a floor on the reconstruction error. In the CRB-following region, the adaptive MO-SPADE imaging attains a performance that is more than an order of magnitude better than direct imaging and is within a small factor of the optimal (prior-knowledge) MO-SPADE imaging.
An example of the spatial imaging modes which achieve this performance are shown in Fig. 4(b,c). While many of the modes are relatively unchanged between the two iterations, some of the modes exhibit altered oscillation patterns. The small differences between the iterations agrees with our empirical findings that additional iterations of the adaptive imaging procedure beyond the first two tend to produce only marginal improvements in image reconstruction quality.
Finally, we demonstrate the performance of the adaptive imaging on arbitrary two-dimensional source distributions. For these results, we assume a two-dimensional Gaussian PSF with width σ , corresponding to apodized imaging through a circular aperture. As in the previous one-dimensional example of adaptive imaging, the measurement time is divided evenly among the initial direct imaging step and the two adaptive modal measurements. Figure 5 shows the results of the adaptive imaging measurement for a source that exhibits spatial frequency chirp and varying contrast, as well as for a Siemens star target. Both sources are placed on top of large uniform backgrounds, simulating the case of low-contrast imaging. The case of low-contrast imaging is particularly advantageous for MO-SPADE imaging, as the bias of the non-negative estimators is lower in the relevant regime where the MSE of the image reconstruction is on the order of the image contrast. Comparing the NNLS reconstruction results, adaptive MO-SPADE imaging outperforms direct imaging and is able to resolve features below the classical Rayleigh resolution limit, with a MSE several times smaller than the direct imaging result. For these plots, the number of imaging modes used in the adaptive modal measurements is equal to the number of source modes plus one, which in both cases equates to 577 imaging modes. However, it is likely that far fewer modes could be used for each adaptive iteration while still obtaining close to optimal results. www.nature.com/scientificreports/

Discussion
Our results demonstrate a general framework for the imaging of arbitrary distributions with near quantumlimited resolution and accuracy. Further improvements to the adaptive imaging method described here could be made by incorporating the bias of estimators into the optimization framework 45 . This would allow for improved imaging performance with estimators that are biased at low photon numbers, such as NNLS and maximum likelihood estimators (e.g. Richardson-Lucy deconvolution). The adaptive imaging method presented here for incoherent thermal sources can be applied to many fields such as multi-emitter fluorescence microscopy and astronomical imaging. Spatial mode imaging is a rapidly developing field, and methods which can incorporate the large number of spatial modes required for adaptive modal imaging of arbitrary 2D sources are being actively investigated 20 . Based on these promising trends, it may be possible to experimentally demonstrate the modal imaging results presented here in the near future.