Sub-Rayleigh resolution of incoherent sources by array homodyning

Conventional incoherent imaging based on measuring the spatial intensity distribution in the image plane faces the resolution hurdle described by the Rayleigh diffraction criterion. Here, we demonstrate theoretically that quadrature statistics measured by means of array homodyne detection enables estimation of the distance between incoherent point sources well below the Rayleigh limit for sufficiently high signal-to-noise ratio. This capability is attributed to the availability of spatial coherence information between individual detector pixels acquired using the coherent detection technique.

While the optical band enables insightful observations of physical, chemical, and biological systems, resolving their spatial characteristics is often hindered by diffractive limitations of imaging instruments described by the fundamental Rayleigh criterion [1]. The Rayleigh limit follows from the direct detection of the spatial intensity distribution of the optical field in the image plane. The purpose of this Letter is to identify theoretically the capability of array homodyne detection [2][3][4] to resolve sub-Rayleigh features in optical imaging. The crucial benefit of array homodyning is the availability of the spatial coherence information for the optical field detected in the image plane. Such coherence is introduced by the transfer function of the imaging system between the source and the image planes even if contributions from individual points constituting the source are mutually incoherent. This observation underlies currently explored approaches to overcome the Rayleigh limit by detecting the optical field in a carefully selected basis of spatial modes in the image plane [5][6][7][8][9][10][11]. The scenario considered here will assume that individual pixels of the homodyne detector have dimension much smaller than the spatial variation of the transfer function, but contribute noise that is independent of their size. This noise model includes the important case of homodyne detection operated at the ultimate shot-noise limit. The sub-Rayleigh sensitivity of array homodyning will be demonstrated for the canonical example of a binary source, where light is emitted by two equally bright and mutually incoherent points.
For simplicity, we shall consider a one-dimensional model of the imaging system characterized by a normalized transfer function u(x), where x parametrises the transverse spatial coordinate in the image plane. Let For concreteness, a field produced by a binary source with two components u1(x) and u2(x) separated by 2d = x1 − x2 and centered at xc = (x1 + x2)/2 is schematically shown. The signal field is superposed with a strong, flat-wavefront local oscillator LO on a 50:50 beam splitter BS with outputs monitored by array detectors. Subtracting photocurrents in between pairs of matching individual pixels yields after rescaling the quadrature vector q = (. . . , qi−1, qi, qi+1, . . .).
the source be in general composed of a finite number of points emitting quasi-monochromatic light described by thermal statistics. The electromagnetic field in the image plane is then represented by the complex signal E (x) = l α l u l (x), where the lth source contributes the displaced transfer function u l (x) = u(x − x l ) multiplied by an amplitude α l characterized by a complex normal distribution α l ∼ CN (0, Pw l ) with a zero mean. For convenience, the variance of α l is expressed as a product of the total optical power P = ∞ −∞ E[|E (x)| 2 ] dx reaching the image plane and the relative weights w l of contributions from individual sources that add up to one, l w l = 1. As shown in Fig. 1, the field in the image plane is detected by means of balanced array homodyne detection with uniform both the efficiency and the noise figure across all the pixels. The ith detector pixel pair is cen-arXiv:2005.08693v1 [quant-ph] 18 May 2020 tered at x i and their size ∆x is assumed to be much smaller than the spatial variation of the transfer function. The differential photocurrent measured between the ith pair of pixels is a sum of contributions from the optical signal I s i and the detector noise I n i , the latter assumed to be Gaussian with a zero mean and uncorrelated between pixels. If the local oscillator phase is uniform across the array and set to zero, the signal contribution reads The detection noise variance Var[I n i ] = N is taken as independent of the pixel size. In such a case it is convenient to rescale the quadrature measured at the ith pixel pair as q i = (I s i + I n i )/ √ N . The quadrature vector q = (. . . , q i−1 , q i , q i+1 , . . .) T is then characterized by a multivariate normal distribution with a zero mean, E(q) = 0, and the covariance matrix where S = P/N is the signal-to-noise ratio (SNR) and I stands for the identity matrix. For an infinitesimal pixel size, the elements of the matrix Γ can be written as is the coherence function normalized by the total signal power in the image plane, The second expression given above follows directly from the absence of coherence between individual point sources. In the case of shot-noise-limited homodyning, the SNR reads twice the total average photon number reaching the image plane from the source. The precision 1/Var[θ] of estimating a parameter θ from a sample of N quadrature vectors q using an unbiased estimatorθ is upper bounded by the product N F θ , where F θ is the Fisher information (FI). For a multivariate normal distribution with a zero mean F θ takes the form [12,13]: In the following, we will consider a binary source comprising two equally bright points producing contributions u 1 (x) = u(x − x 1 ) and u 2 (x) = u(x − x 2 ) in the image plane located respectively at x 1 = x c +d and x 2 = x c −d, as shown in Fig. 1. Here d specifies the half-separation between the points and x c is the centroid of this binary source. Two models of the normalized real transfer function will be used in numerical examples: [hard aperture] (4) In both cases, the first derivative u (x) of the transfer function is used to characterize its spatial spread as which in turn defines for regular transfer functions the Rayleigh limit below which the resolution of conventional direct imaging is lost [14]. Fig. 2(a) depicts the FI F d for estimating the halfseparation d from the array homodyne measurement of field quadratures for the SNR S = 25, 100, and 400, calculated using Eq. (3) in the limit ∆x → 0. Details of the calculations are presented in Appendix. For d σ the Fisher information approaches the value S/[(1+2S −1 )σ 2 ] independently of the model of the transfer function selected in Eq. (4). In the sub-Rayleigh region, when d σ, a non-trivial feature appears in the form of a peak whose maximum shifts towards lower d with increasing SNR. Interestingly, for high SNR the peak shape does not depend noticeably on the model of the transfer function.
The dependence of the Fisher information F d on the half-separation d observed above can be understood with the help of the KarhunenLoéve decomposition [15] of the coherence function Γ(x, x ) defined in Eq. (2), which for the binary source takes the form In the case of two equally bright points and a real-valued transfer function the two eigenvalues read γ ± = (1±χ)/2, where χ = ∞ −∞ u 1 (x)u 2 (x) dx is the overlap between the two displaced transfer functions corresponding to individual points in the binary source. The respective normalized eigenmodes are e ± (x) = [u 1 (x) ± u 2 (x)]/ 2(1 ± χ).
For an infinitesimal pixel size ∆x, the KarhunenLoéve expansion of the coherence function directly facilitates the principal component analysis of the quadrature vector q. In the case of a binary source, the covariance matrix C has two distinguished eigenvectors e ± that can be written in a normalized form using the values of the eigenmodes e ± (x) on the pixel grid, These vectors define two principal components q ± = e T ± q of the multivariate quadrature distribution that are characterized by a zero mean and respective variances: The variances of all other components of q that are orthogonal to e ± are equal to one. As shown in Appendix, the FI F d can be written in the current case as a sum of two contributions, , that stem respectively from the dependence on d of the eigenvectors e ± and the eigenvalues V ± of the covariance matrix C. The change of the eigenvectors e ± with d results in the Rayleigh term F with the proportionality factor of the order of one and for d σ saturates at F For high SNR it is located in the region d σ and its shape turns out to be determined primarily by the dependence of V − on d. An approximate description of F (SR) d can be obtained in a closed analytical form for regular transfer functions by applying the Taylor series expansion up to the quadratic term in d to u(x − x c ± d). This yields χ ≈ 1 − 2d 2 /σ 2 with σ defined in Eq. (5). Consequently, according to Eq. (7) one has V − ≈ Sd 2 /σ 2 + 1. Using this expression in Eq. (3) specialized to the univariate case of the variance V − yields where we have introduced f (t) = 2t 2 /(1 + t 2 ) 2 . Fig. 2(b) compares the second approximate expression in Eq. (8) with the rescaled FI σ 2 F d /S for the soft aperture model. It is seen that the approximation given in Eq. (8) reproduces rather accurately the shape of the sub-Rayleigh feature for high SNR. An elementary analysis of the function f (t) yields the maximum of the sub-Rayleigh peak at d = σ/ √ S and the endpoints of the half-maximum interval located at ( √ 2 ± 1)σ/ √ S. Physically, the principal component q − = e T − q whose variance carries most information about the source separation in the sub-Rayleigh regime corresponds to the quadrature of the spatial mode given in the image plane by e − (x). For small separations, d σ, this mode function can be approximated using a straightforward Taylor series expansion by that does not depend explicitly on d. Furthermore, the variance V − that serves as the basis for estimating d effectively measures the optical power carried in the spa- . This relates the estimation recipe for the source separation emerging from the principal component analysis to the currently explored approaches to superresolution imaging based on spatial mode demultiplexing (SPADE) [5] and similar techniques, where the separation of a binary source is inferred from the fraction of the optical power directed to carefully defined spatial modes in the image plane [16]. For small separations, d σ, the relevant information is contained predominantly in the optical power measured for the mode v(x − x c ) [14]. While SPADE and similar techniques require a careful alignment of the detection apparatus hardware with respect to the field in the image plane in order to avoid a systematic error [17][18][19], the advantage of array homodyning is the ability to reconstruct the quadrature statistics for any relevant spatial mode through digital postprocessing of the quadrature vectors obtained from the pixelated measurement [20]. In particular, prior knowledge of the source centroid is not required to align the array homodyne detector. This is in contrast to SPADE-type techniques, where the need to determine the centroid results in an overhead in terms of the required signal [21]. However, it needs to be verified whether the array homodyne data are sufficient on their own to estimate the separation in the sub-Rayleigh region.
As the eigenmodes e ± (x) playing the central role when estimating d in the sub-Rayleigh region depend explicitly on the source centroid, let us examine first the precision of centroid estimation from the quadrature vectors. Fig. 2(d) depicts the FI F xc for estimating the centroid x c of a binary source as a function of the half-separation d calculated using Eq. (3). As detailed in Appendix, the dip seen for d σ stems from the fact that F xc is given by a sum of two contributions F xc = F can be explained intuitively as follows. For the infinitesimal pixel size assumed here, the spatial intensity distribution obtained from quadrature variances measured at individual pixels is overwhelmed by the detection noise and cannot be used standalone for reliable estimation of the centroid. Instead, the information about the centroid is primarily obtained from coherences between pixels that induce non-trivial covariances between individual quadratures according to Eq. (1). These coherences turn out to be most informative either when d ≈ 0, i.e. one is effectively dealing with a single point source, or when the two points constituting the binary source are well separated.
The above observation raises the issue whether array homodyning can be indeed used in the sub-Rayleigh region to determine with adequate precision the source separation without prior knowledge of the source centroid. This question is answered positively by the following algorithm. Use the measured quadrature vectors q to determine the variance V xr = Var[q xr ] of quadratures √ ∆x defined for a one-parameter family of spatial modes obtained by displacing v(x) specified in Eq. (9) by an arbitrary distance x r , as depicted in Fig. 3(a). For a general discrete source, the variance V xr is given in the limit ∆x → 0 by an integral expression As shown in Fig. 3(b) for the soft aperture model, in the case of a binary source the graph of V xr as a function of x r exhibits a two-lobe structure on top of the detection noise pedestal. The local minimum between the lobes can serve as an estimate for the centroidx c . This facilitates the estimation of the source separation from the gap between Vx c and the detection noise level by inverting Eq. (11) with Vx c used on the left hand side and x c inserted in lieu of x r on the right hand side. For d σ, the estimation formula takes an approximate form d ≈ σ (Vx c − 1)/S independently of the selected model of the transfer function.
The algorithm outlined above has been applied to Monte Carlo data generated by simulating for given source parameters 1000 realizations of the array homodyning experiment with the soft aperture model. In each realization, a sample of N = 500 quadrature vectors has been drawn for 1000 pixels of width 0.008σ each. This sample was used to compute the variance V xr as a function of x r and subsequently to determine from its local minimum the estimates for the centroidx c and the separationd. The precision of these estimates shown in Fig. 2(b, e) is given by the squared inverse of the standard deviations for the histograms ofd andx c determined from individual realizations, after being rescaled by N −1 . It is seen that for high SNR, the attainable precision follows the sub-Rayleigh part of the FI. However, it should be noted that the estimatord for the half-separation exhibits a minor negative bias that can be observed in Fig. 2(c). This is easily explained by the observation that estimating the half-separation d from V xc would produce no bias, and replacing V xc by the local minimum value Vx c can only decrease the estimated value of d. In contrast, the estimatorx c has no noticeable bias, as illustrated with Fig. 2(f). Similar results have been obtained for Monte Carlo simulations using the hard aperture model.
Concluding, the analysis of estimating separation for a binary source demonstrates the potential of array homodyne detection to resolve spatial features of composite sources well below the Rayleigh diffraction limit, provided that sufficiently high SNR can be attained. It should be noted that equivalent information about the source is obtained from the statistics of quadratures measured for orthonormal sets of spatial modes, e.g. Gauss-Hermite modes, that can be separated with the help of multiplane light conversion [22]. Spatially multimode coherent detection is currently being developed as a technique to boost the capacity of optical communication links via spatial division multiplexing [23]. In a preliminary study following work presented here we have found that a joint measurement of both field quadratures enables one to handle complex transfer functions, albeit at a factor of two penalty for the SNR. Finally, an interesting extension of the presented work would be the reconstruction of properties of more intricate composite sources, also two-dimensional, using the spatial coherence information supplied by array homodyning. ment Fund and the US Department of Navy award no. N62909-19-1-2127 issued by the Office of Naval Research.
Disclosure. The authors declare no conflicts of interest.

APPENDIX
Consider a scenario where a sample of N vectors composed of real random variables q = (. . . , q i−1 , q i , q i+1 , . . .) T is used to determine values of parameters θ j with the help of estimatorsθ j that are unbiased, i.e. E[θ j ] = θ j . The estimation precision can be characterized by the covariance matrix with elements C jj = Cov[θ j ,θ j ]. This covariance matrix satisfies the Cramér-Rao bound where F is the Fisher information matrix. Note that the variance of individual estimators is lower bounded by The second inequality is not necessarily tight when the Fisher information matrix is not diagonal.
If the variables q follow a normal multivariate distribution with a covariance matrix C and a zero mean, E[q] = 0, individual elements of the Fisher information matrix F are given by [12,13]: We shall consider the covariance matrix C in the form: where the principal components corresponding to the orthonormal eigenvectors e µ are distinguished by the fact that their variances V µ = Var[e T µ q] differ from one, V µ = 1. The matrix P ⊥ = I − µ e µ e T µ is the projection onto the subspace orthogonal to the subspace spanned by the eigenvectors e µ . The inverse of the covariance matrix is given explicitly by C −1 = µ V −1 µ e µ e T µ + P ⊥ and its derivative with respect to the parameter θ j reads For the sake of brevity, we have introduced shorthand notation ∂ j = ∂/∂θ j . Using the above expressions, the Fisher information matrix defined in Eq. (13) can be conveniently written as a sum of three components F = F (1) +F (2) +F (3) with elements given by The sums over µ are restricted to µ = µ owing to the fact that e T µ (∂ j e µ ) = 0, which follows directly from the normalization of the eigenvectors e µ . Furthermore, the orthogonality of the eigenvectors implies that for any µ, µ one has e T µ (∂ j e µ ) + e T µ (∂ j e µ ) = ∂ j (e T µ e µ ) = 0. The three components of the Fisher information matrix F (1) , F (2) , and the scalar products appearing in Eqs. (17) and (18) can be effectively expressed as integrals, e.g. e T µ ∂ j e µ = ∞ −∞ e µ (x)∂ j [e µ (x)]dx, provided that the functions e µ (x) vary slowly over the scale defined by ∆x. The above formalism can be applied in a straightforward manner to the estimation of the half-separation d and the centroid x c of a binary source composed of two equally bright points and observed using an imaging system described by a real transfer function u(x). For notational simplicity, it is convenient to use d and x c as the indices of the two-dimensional Fisher information matrix. Moreover, it can be verified by a direct calculation that in the case considered here the off-diagonal element of the Fisher information matrix vanishes, F xcd = 0. Hence the estimation of the half-separation and the centroid can be treated as statistically independent. For a binary source, the summation index in Eqs. (16)- (18) takes two values µ = ±. The respective eigenmodes are with the corresponding eigenvalues V ± = S(1 ± χ)/2 + 1. Here χ is the overlap between the two displaced transfer functions, which for small separations can be approximated up to the quadratic order in d using the Taylor series expansion as: In the second line integration by parts has been carried out assuming that the first derivative u (x) vanishes in the limit x → ±∞ and the definition of σ = ∞ −∞ [u (x)] 2 dx −1/2 from Eq. (5) of the main text has been used. For clarity, the calculation in Eq. (21) has been presented for x c = 0. For simplicity, the diagonal elements of the Fisher information matrix have been denoted in the main text as F j ≡ F jj , where j = d, x c . When the covariance matrix given in Eq. (14) has only two distinguished eigenvectors, the third contribution F (3) jj to the diagonal elements of the Fisher information matrix given by Eq. (18) reduces to Specifically, for j = d a straightforward calculation yields This immediately implies that F For small d the leading contribution comes from the first term with (V + − 1) 2 /V + ≈ S/(1 + S −1 ) and the integral ∞ −∞ [∂ d e + (x)] 2 dx producing an expression quadratic in d, The factor within the large round brackets is equal to 2/σ 4 for the soft aperture model and 4/(5σ 4 ) for the hard aperture model. For large separations the overlap χ vanishes and the eigenvalues V ± ≈ S/2 + 1 become equal to each other. Hence one can equivalently take as the eigenmodes in Eq. (24) the displaced transfer functions u(x − x c ± d) instead of e ± (x). The result is: The physical interpretation is that in the regime of large separations, the half-separation is obtained through estimation of the locations x 1 and x 2 of the two peaks produces by individual points in the binary source and calculating d = (x 1 − x 2 )/2. The same argument holds in the limit of large separations for the contribution F (2) xcxc that plays the role of the Rayleigh part F (R) xc in the case of estimating the centroid x c = (x 1 + x 2 )/2. Hence a calculation analogous to the one leading to Eq. (26) yields also Furthermore, a direct expansion into a power series shows that for high signal-to-noise ratio, S 1, and small d the Rayleigh part has the leading-order term of the form F (R) xc ≡ F (2) xcxc ≈ S 2 d 4 /σ 6 . The remaining part, given by Eq. (22) specialized to j = x c , plays a non-trivial form in the sub-Rayleigh region and has been denoted as F (SR) xc ≡ F (3) xcxc . Its approximate form for high signal-to-noise ratio S 1 can be obtained by neglecting terms of the order of d 2 /σ 2 in comparison to one, but retaining terms of the order of Sd 2 /σ 2 . Under these assumptions one has ∞ −∞ e − (x)∂ xc [e + (x)]dx 2 ≈ 1/σ 2 , V + ≈ S + 1, and V − ≈ Sd 2 /σ 2 + 1. Consequently, [small d] The last expression is used in Eq. (10) of the main text. For illustration, Fig. 4 depicts the Fisher information for estimating the half-separation d and the source centroid x c as a sum of the respective sub-Rayleigh and Rayleigh terms. Results for both the soft and the hard aperture models are shown. The sub-Rayleigh terms are compared with analytical approximations given in Eqs. (8) and (10) of the main text.