Sub-Rayleigh characterization of a binary source by spatially demultiplexed coherent detection

We investigate theoretically coherent detection implemented simultaneously on a set of mutually orthogonal spatial modes in the image plane as a method to characterize properties of a composite thermal source below the Rayleigh limit. A general relation between the intensity distribution in the source plane and the covariance matrix for the complex field amplitudes measured in the image plane is derived. An algorithm to estimate parameters of a two-dimensional symmetric binary source is devised and verified using Monte Carlo simulations to provide super-resolving capability for high ratio of signal to detection noise (SNR). Specifically, the separation between two point sources can be meaningfully determined down to $\textrm{SNR}^{-1/2}$ in the units determined by the spatial spread of the transfer function of the imaging system. The presented algorithm is shown to make a nearly optimal use of the measured data in the sub-Rayleigh region.


Introduction
When composite incoherent light sources are imaged using a conventional optical system, the system transfer function introduces partial spatial coherence of the optical field in the image plane [1]. This coherence escapes detection with direct imaging, i.e. a spatially resolved measurement of the intensity distribution, which results in the well-known Rayleigh curse, where spatial resolution is lost below the scale defined by the characteristic spread of the transfer function [2]. As discussed in a number of recent works [3][4][5][6][7][8][9][10][11][12][13][14], the Rayleigh limit can be overcome by measuring the optical intensity carried by one or more spatially extended modes in the image plane. The specific form of these modes is determined by the shape of the transfer function. Such a measurement effectively utilizes the spatial coherence information, as filtering out a spatial mode can be viewed as a coherent combination, with certain weights, of the field amplitudes across the image plane. However, it is recognized that superresolution techniques that are based on spatially demultiplexed intensity measurements require prior knowledge of the source properties, e.g. its centroid, or, in the absence of thereof, suffer from trade-offs in the precision attainable when estimating multiple parameters characterizing the source [15][16][17]. Attention also needs to be paid to imperfections, such as intermodal crosstalk [18,19], misalignment [19,20] and dark counts [19,[21][22][23].
The spatial coherence of light can be also measured by means of homodyne or heterodyne coherent detection. Then, rather than with intensity measurements, one can consider spatial mode demultiplexing followed with coherent detection to access the coherence information, using, in the simplest scenario, a single mode [7,8], or most generally, a set of mutually orthogonal modes. For the latter, such a measurement strategy is becoming feasible owing to the developments in spatial light conversion technology [24,25] combined with coherent detection techniques for high-capacity optical communications [26].
The purpose of this paper is to investigate the capability of spatially demultiplexed coherent detection to characterize and determine the features of a composite source in the sub-Rayleigh regime. In particular, we establish a general result for the spatial coherence information obtainable with this technique, as summarized in the covariance matrix for the complex amplitudes of the demultiplexed modes. Using this information, we then present an algorithm designed to estimate the source centroid, orientation, and the separation between the source points, for the generic model of a two-dimensional binary source comprising two equally bright points whose separation is well below the Rayleigh criterion. Importantly, the set of demultiplexed modes can be quite arbitrary, and does not need to match the transfer function of the imaging system, as the latter enters only the estimation algorithm. Nevertheless, when some coarse knowledge of the source centroid and the transfer function is known, the number of modes that need to be detected can be reduced substantially, as most information about the source will be contained in these few modes. This makes the scheme presented here more practical compared to the previously analyzed one-dimensional model of spatially resolved homodyne detection [27] that would produce vast amounts of raw data in the case of two-dimensional array detectors. Furthermore, we verify that the performance of the estimation algorithm presented here is close to optimal, namely that it nearly saturates the Cramér-Rao bound for the parameters of interest in the sub-Rayleigh region.
This paper is organized as follows. In Sec. 2 we introduce the model of the imaging system with spatially demultiplexed coherent detection and derive a general relation between the intensity distribution in the source plane and the covariance matrix for the complex amplitudes of the demultiplexed modes detected in the image plane. This result is specialized to the case of a binary symmetric source in Sec. 3, where also the super-resolving capability of coherent detection is identified using the mathematical concept of the Fisher information. In Sec. 4 the estimation algorithm is presented and applied to Monte Carlo simulated data. Finally, Sec. 5 concludes the paper.

Model
Consider a scenario where individual points in the source plane emit mutually incoherent thermal radiation contributing to the optical field ℰ src (r), where the subscript src indicates the quantity at the source plane. Here r = ( , ) are the spatial coordinates parameterizing the plane perpendicular to the propagation axis. Under these assumptions, the coherence function of the source field can be written as where is the total optical power and (r) is the normalized intensity distribution in the source plane, ∫ 2 r (r) = 1. The imaging system will be modeled by a transfer function (r) satisfying the normalization condition ∫ 2 r | (r)| 2 = 1. Taking for simplicity unit magnification of the optical system, the signal field ℰ(r) in the image plane is expressed by As depicted in Fig. 1, the field in the image plane is demultiplexed into a basis of orthonormal modes (r), labeled with a double index , that are subsequently subject to a coherent measurement of both conjugate quadratures, using e.g. heterodyne detection.
The result of the measurement on the th mode is a complex amplitude whose real and imaginary parts correspond to the two conjugate quadratures of the detected field. The detection noise will be modeled by an additive component noise ∼ CN (0, ) taken to follow a symmetric complex Gaussian distribution with a zero mean and a variance equal for all , Fig. 1. An imaging system with spatially demultiplexed coherent detection. For concreteness, the source is taken as composed of two equally bright points separated by 2 , oriented at an angle , and centered at r = ( , ). The optical field in the image plane is separated into a set of orthogonal spatially extended modes (r) that are subject to coherent detection which yields complex amplitudes . The source properties are determined from the covariance matrix with elements assumed to be uncorrelated with the detected optical signal as well as between the demultiplexed modes. For notational simplicity the amplitude will be specified in units determined by the variance of the detection noise: Combining the assumption of mutually incoherent thermal radiation and others made about the source characteristics, the imaging system, and the detection scheme, it is straightforward to obtain [28] that the variables are described by a circularly-symmetric complex multivariate Gaussian distribution with a zero mean, E[ ] = 0, and its covariance matrix C has elements equal to (4) Here S = / is the signal-to-noise ratio and¯(r ) are projections of the transfer function displaced by r onto the th demultiplexed mode, One can reproduce the case of pixelated array detection by taking (r) equal to 1/ √︁ Δ Δ over a rectangle of dimensions Δ × Δ centered at ( Δ , Δ ) and zero outside, with , running through the entire integer range, provided that the dimensions Δ , Δ are much smaller than the spatial variation of the transfer function (r). Such a detection scheme can be viewed as a direct, albeit noisy, measurement of the spatial coherence function in the image plane, with the noise contribution given by the second term of Eq. (4). For future reference, we also specify here a concrete example of Hermite-Gaussian functions indexed with non-negative integers , ≥ 0 Here determines the spatial extent of the modes and HG ( ) stand for the Hermite-Gaussian functions given by where ( ) denotes the th Hermite polynomial. A demultiplexing technique for mode functions specified in Eq. (6) has been recently demonstrated experimentally [29]. If the set of the mode functions (r) is complete, knowledge of the covariance matrix elements , makes it possible to reconstruct the variance of the coherently detected complex amplitude for any normalized mode (r) in the image plane. The explicit expression for the variance Var[ ] reads where¯= Thus, we emphasize that in principle coherent detection implemented for any choice of the set of demultiplexing modes, as long as it forms a complete basis, yields the same information about the field quadratures. Let us also note that for shot-noise limited coherent detection the signal-to-noise ratio S is given by the mean number of photons carried by the signal [21,30].

Binary source
Optical coherence in the image plane introduced by the transfer function (r) of the imaging system can provide means to identify sub-Rayleigh features of the source [4,5]. We will study this capability using the canonical example of a symmetric binary source with the normalized intensity distribution (r) of the form where r 1 , r 2 denote locations of light source components. It will be convenient to write r 1 = r + e and r 2 = r − e , where r = ( , ) is the source centroid, specifies the half-separation between the two source points, and the unit vector e = (cos , sin ) determines the source angular orientation. Henceforth we will use the soft-aperture model for the imaging system with a Gaussian transfer function given by The above formula assumes that the spatial spread of the transfer function defines the unit of length, meaning that centroid coordinates as well as separation appearing in equations are dimensionless quantities. Other quantities, such as parameters of the source or the spatial extent of the demultiplexed modes used in Eq. (6) will be specified relative to this unit. The second expression in Eq. (11) is written using notation introduced in Eq. (7).
The objective now is to use data collected by means of spatially demultiplexed coherent detection to estimate the centroid r , half-separation , and orientation of a symmetric binary source in the sub-Rayleigh regime when 1. Motivated by the one-dimension study carried out in [27], we will consider a normalized reference mode r ; (r) of the form The above expression can be viewed as a product of a one-dimensional transfer function HG 0 in one direction and its normalized derivative HG 1 in the orthogonal direction, rotated by and displaced by r = ( , ). The variance (r ; ) of the coherently detected complex amplitude for the mode r ; (r) can be calculated from the covariance matrix elements , using Eqs. (8) and (9) with r ; (r) used in lieu of (r) in Eq. (9).
A careful inspection of the dependence of (r ; ) on its arguments r and allows one to identify the following relation to the source parameters. As shown in Fig. 2(a), for a given the graph of (r ; ) as a function of r has the form of two lobes surrounding symmetrically the source centroid. While the angular orientation of the lobes is determined by , the value of the variance at the midpoint between the lobes, i.e., r , exhibits dependence on both and as depicted in Fig. 2(b). The actual expression for the variance at the source centroid reads: It is seen that the location of the maximum of (r = r ; ) as a function of corresponds to the orientation of the source, whereas the source separation can be read out from the height of the maximum above the detection noise level, equal to one in the chosen units. In the sub-Rayleigh region, when 1, the expression given in Eq. (13) can be simplified to the leading order in , which yields (r = r ; ) ≈ S 2 cos 2 ( − ) + 1.
In order to gain a simple insight into how precisely the source separation can be inferred from the variance map (r ; ) let us recall the notion of the Fisher information. Generally, for any unbiased local estimator˜of a real continuous parameter , the Fisher information F provides an upper bound on the attainable precision (Δ 2˜)−1 , defined as the inverse of the estimator variance Var[˜] normalized by the sample size used for estimation: When a parameter is estimated from a univariate circularly symmetric complex normal distribution with a variance , the Fisher information reads [31,32] For the sake of simplicity, assume for a moment that the source centroid is known perfectly and that the angle is set to . Using the approximate formula for (r = r ; = ) given in Eq. (14) it is elementary to derive the simplified expression for the Fisher information: Importantly, F attains the value of the order of S for ∼ S −1/2 . This contrasts with direct imaging, where the corresponding precision level is reached only for ∼ 1 in the units determined by the spread of the transfer function [21]. The above observation indicates the super-resolving capability of coherent detection for high signal-to-noise ratio, when S 1. Note, however, that F still vanishes for = 0, which is due to the presence of shot noise characteristic for quadrature detection [21,27].

Estimation algorithm
In order to estimate simultaneously the full set of parameters of a two-dimensional binary symmetric source the following numerical algorithm is proposed. First, within the interval 0 ≤ ≤ , select equidistant sampling points = / , where = 1, . . . , . For each , find locations of the two maxima of the variance (r ; = ) as a function of r and take the midpoint r = ( , ) between these maxima. Fitting a function given on the right hand side of Eq. (13) to the set of points , (r = r ; = ) yields estimates for the source half-separation˜and orientation˜. The estimate for the source centroid is given as an arithmetic average of all the midpoints,r = (˜,˜) = (1/ ) r . Importantly, note that in a practical scenario the variance (r ; ) is calculated from a finite set of detected modes (r). Because of this limitation one is able to access spatial coherence information only in a certain area of the image plane. Therefore one should possess at least a coarse knowledge of the source centroid and the transfer function in order to choose a set of demultiplexed modes (r) that covers the image. The better these parameters are known the lower number of modes needs to be detected.
The above estimation algorithm has been applied to a Monte Carlo simulated data for a binary source characterized by r = (−0.05, 0.1) and oriented at = /4. We assumed heterodyne detection implemented on a set of Hermite-Gauss modes defined in Eq. (6) with = 0.8 and the indices , restricted to the range , ≤ 4. The latter condition has been chosen to ensure that the detected modes carry nearly entire signal power. This requirement can be verified by calculating the sum , ≤4 ∫ 2 r (r )| (r )| 2 that exceeds 0.9994 for all values of the half-separation used in the Monte Carlo simulations. One realization of the experiment was taken to comprise = 500 elementary measurements of field quadratures for the chosen set of Hermite-Gauss modes. The statistical properties of the estimators were analyzed using 1000 of such realizations. The estimation quality can be further improved by taking into consideration more modes in the numerical postprocessing, which, in view of their extremely low signal amplitude, can be essentially approximated by just carrying only detection noise. Note that in the case of larger separations these higher order modes carry more information, however, in such instances one is outside the sub-Rayleigh regime and direct imaging techniques can be used instead. Specifically, the results presented below have been obtained for a heterodyne experiment simulated with Hermite-Gauss modes up to , ≤ 4 and then calculating the variance (r ; ) including modes 0 ≤ , ≤ 9, where for 4 < ≤ 9 or 4 < ≤ 9 complex amplitudes ∼ CN (0, 1) have been taken in the estimation algorithm. Furthermore, the interval 0 ≤ ≤ has been sampled with = 10 equidistant points.  Fig. 3(a) represents the simplified form of the Fisher information F derived in Eq. (17). It is seen that this elementary expression characterizes rather well the precision achieved by the devised algorithm for estimation of half-separation . The solid lines in Fig. 3(a-d) depict the Cramér-Rao bound on attainable precision that takes into account information available in the entire covariance matrix C. This bound is given by the Fisher information which for data distributed according to a circularly-symmetric complex normal distribution parameterized by a parameter entering through the covariance matrix C is equal to [32] It can be seen in Fig. 3(a,b) that coherent detection in the sub-Rayleigh regime outperforms precision offered by direct imaging [3], the latter indicated by dotted curves. Importantly, precision for orientation estimation offered by the latter detection scheme also exhibits Rayleigh limit. Unlike in the one dimensional case [17] this fact therefore prevents one from using various hybrid strategies based on first estimating centroid and orientation using direct imaging and then performing estimation of separation with the appropriately adjusted SPADE measurement [17]. Note that the above expression can be viewed as a multivariate generalization of Eq. (16). It is seen that the proposed estimation strategy is nearly optimal in the sub-Rayleigh region 1. The loss of precision for larger value of compared to the Cramér-Rao bound can be attributed to the need of analyzing statistics of more than just one spatial mode in the image plane [27]. It can be noticed in Fig. 3(e-h) that the estimates for the half-separation and the centroid coordinates exhibit a minor bias, which however remains much below the absolute values of the respective quantities and can be further reduced by including a larger number of modes.

Conclusions
We have investigated theoretically the capability of spatially demultiplexed coherent detection in the image plane to characterize a two-dimensional binary source composed of two mutually incoherent points. It has been found that the source separation can be estimated for values down to ∼ S −1/2 in units determined by spread of the transfer function of the imaging system. Here S is the signal-to-noise ratio of coherent detection, which for shot-noise limited scenario is given by the mean number of photons reaching the image plane. This provides a super-resolving capability when S 1. We presented and tested with Monte Carlo simulations an algorithm to estimate the source parameters: its half-separation, orientation and centroid location with precision approaching the Cramér-Rao bound in the sub-Rayleigh region. Central to the algorithm is the extraction of information about the source parameters from the variances of the quadrature measurements of the complex amplitude carried by a reference mode, for different orientation and displacement relative to the source parameters. Crucially, the choice of the reference mode can be decided later in the algorithm, and we need not to perform the experiments with the demultiplexed modes matching the imaging system transfer function-which in principle need not be known precisely at the time of the measurement-as variance in any mode can be obtained from covariance of another set of modes. Still, when coarse knowledge about the centroid and transfer function spread is given a priori, such that the set of demultiplexed modes covers the image, the estimation algorithm performs well with just a limited number of modes, with the reference mode to use identified accordingly. This makes spatial demultiplexing quadrature measurement more efficient than spatially resolved coherent detection which would require recording quadratures for a large number of pixels. Importantly, preliminary calculations indicate that the proposed estimation algorithm is capable of estimating source parameters with superresolving precision not only for two equally bright sources but also in more complicated scenarios involving unequal brightness.
A potentially useful aspect of spatially demultiplexed coherent detection is that in principle it yields equivalent information about the spatial coherence even if it is implemented in a plane other than the focal plane of the imaging system. Hence an out-of-focus coherent measurement should also have a super-resolving capability for high signal-to-noise ratio. In this scenario, however, one needs to consider a complex transfer function describing the imaging system. Another interesting direction of further research is characterization of sources extended along the axis of the imaging system [33].