{'fraction_non_alphanumeric': 0.04851786397507606, 'fraction_numerical': 0.01718075318381402, 'mean_word_length': 4.308317214700193, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Moment-based superresolution: Formalism and applications

Sensitivity limits are usually determined using the Cram\'er-Rao bound. Recently this approach has been used to obtain the ultimate resolution limit for the estimation of the separation between two incoherent point sources. However, methods that saturate these resolution limits, usually require the full measurement statistics, which can be challenging to access. In this work, we introduce a simple superresolution protocol to estimate the separation between two thermal sources which relies only on the average value of a single accessible observable. We show how optimal observables for this technique may be constructed for arbitrary thermal sources, and we study their sensitivities when one has access to spatially resolved intensity measurements (direct imaging) and photon counting after spatial mode demultiplexing. For demultiplexing, our method is optimal, i.e. it saturates the quantum Cram\'er-Rao bound. We also investigate the impact of noise on the optimal observables, their measurement sensitivity and on the scaling with the number of detected photons of the smallest resolvable separation. For low signals in the image plane, we demonstrate that our method saturates the Cram\'er-Rao bound even in the presence of noise.


I. INTRODUCTION
Resolving small angular separations through an optical imaging system is an important problem both in microscopy and in astronomy. The most traditional imaging technique is a spatially resolved intensity measurement, also known as direct imaging. The resolution of this approach, as pointed out already by Abbe [1] and Rayleigh [2] at the end of the nineteenth century, is limited by diffraction and noise. However, for a large enough signalto-noise ratio, the diffraction limit is not fundamental [3], and can be overcome by superresolution techniques based on fluorescence microscopy [4][5][6], homodyne measurements [7][8][9] or intensity measurements in an appropriate spatial mode basis [10,11].
Recently, superresolution imaging was analysed from the point of view of quantum metrology [12][13][14][15][16][17], and the Fisher information was used to determine how well the separation between two incoherent sources can be resolved. In this framework, the diffraction limit manifests itself through the vanishing of the direct-imaging Fisher information for small source separations [18]. On the other hand, the quantum Fisher information, i.e. the Fisher information optimized over all possible measurements in the image plane, stays constant for small distances [18], leaving room for superresolution. Furthermore, it was shown [18] that the ultimate resolution limit given by the quantum Fisher information can be approached by photon counting after spatial-mode demultiplexing. Several experiments used spatial light modulators to implemented a simplified version of this demultiplexing technique, which distinguishes only between the fundamental and the first excited modes [19][20][21][22] On the other hand distance estimation from the demultiplexing of multiple spatial modes, realized using multi-plane light conversion [23], was only recently reported [24].
In general, to reach the Cramér-Rao bound, e.g. via maximum likelihood estimation, requires to measure the full photon counting statistics, which can be practically challenging. Here, we demonstrate that this is not necessary in superresolution imaging, where the ultimate resolution can be obtained using a moment-based estimation technique, that requires to measure only the average value of a single measurement observable. For this estimation technique, we identify the optimal observables when different measurements, such as spatially resolved intensity measurements (direct imaging) or photon counting after spatial mode demultiplexing, are available.
In particular, we focus on the estimation of the transverse separation between two thermal sources of arbitrary, and different brightnesses, and we determine the sensitivity that can be reached with our optimized observables, and consequently the minimal resolvable distance between the sources. For demultiplexing measurements, we construct the optimal observables also in presence of experimental imperfections such as misalignment, measurement crosstalk, and detector noise. Therefore, our results are directly relevant for practical applications. Even in the presence of noise, for low brightnesses of the sources, we prove that our approach is sufficient to saturate the Cramér-Rao bound. Finally, for arbitrary brightnesses of the sources, we demonstrate that our optimized demultiplexing measurement allows to reach the quantum Crámer-Rao bound [25,26], in the noiseless scenario, if sufficiently many modes are measured.
The paper is structured as follows: First, in Sec. II, we present our model for thermal sources in the image plane of an imaging system. After recalling the method of moments for parameter estimation in Sec. III, in Sec. IV, we use it to construct the optimal observable to estimate the source separation via spatial mode demultiplexing, while in Sec. V, we employ it to bound the sensitivity of ideal direct imaging. Detailed studies of the performances of our moment-based approach for ideal and noisy demultiplexing are presented in Sec. VI, and Sec. VII respec-tively. In Sec. VIII, we discuss the smallest source separation which is resolvable by the different measurement techniques. Section IX concludes our work.

II. THERMAL STATES IN THE IMAGE PLANE
We want to estimate the transverse distance between two point thermal sources located at positions ±r 0 , with r 0 = (d cos θ/2, d sin θ/2). The sources emit a total mean photon number equal to 2N in the spatial modes associated with the field operatorsŝ 1 andŝ 2 . We denote witĥ ρ a (N ) a thermal state with mean photon number N in the mode associated with the field operatorâ. Accordingly, the sources are described by the quantum statê ρ 0 =ρ s1 [(1 − γ)N ] ⊗ρ s2 [(1 + γ)N ] where −1 < γ < 1 takes into account the possibly different (but finite) temperatures of the two sources. Such a state is described by the density matrix ρ 0 = d 2 α 1 d 2 α 2 P s1,s2 (α 1 , α 2 ) |α 1 , α 2 α 1 , α 2 | , (1) where α 1/2 are coherent states of the field operatorŝ s 1/2 , and P s1,s2 (α 1 , α 2 ) = P s1 (α 1 )P s2 (α 2 ) is the Glauber-Sudarshan P −function, with The evolution of the field operators through a diffraction-limited imaging system, with transmissivity κ, is given by [26,27] whereĉ 1/2 are the field operators associated with the images u 0 (r ± r 0 ) of the two sources, with u 0 (r) the point spread function (PSF) of the imaging system, which we assume to be real up to a global phase. The field operatorsv 1/2 are associated with auxiliary modes, that are in the vacuum state. This beam-splitter model for the propagation through an imaging system is illustrated in Fig. 1 (a). The modes u 0 (r ± r 0 ) are non-orthogonal, and therefore the operatorsĉ † 1 andĉ 2 do not commute (see Fig. 1  (a)). As a consequence, these modes cannot be used to properly represent the quantum state of the sources in the image plane. To obviate this problem, we introduce the orthonormal image modes where δ represents the overlap between the source images δ = d 2 ru * 0 (r + r 0 )u 0 (r − r 0 ).
FIG. 1. Equivalent beam-splitter models for the propagation of thermal states through a diffraction-limited imaging system. (a) The source's modesŝ 1/2 are populated with two thermal states with photon numbers N (1 ± γ), and are mixed with the non-orthogonal vacuum modesv 1/2 on a beam splitter with transmissivity κ, resulting in the non-orthogonal image modesĉ 1/2 . (b) The symmetric and antisymmetric modeŝ s± have equal mean photon number N , but are classically correlated with phase-insensitive correlations ŝ † ±ŝ ∓ = γN . The modesŝ± are mixed with the orthogonal vacuum modesv± on two beam splitters with transmissivities κ±, resulting in the orthogonal image modesb±.
The relation between the field operatorsb ± associated to the modes u ± (r) and the field operatorsŝ ± = (ŝ 1 ± s 2 )/ √ 2 in the object plane can be obtained from the sum and difference of Eqs.

III. THE METHOD OF MOMENTS
We estimate the distance d between the two sources with the method of moments [17,29]. Following this approach, given an observableX, an estimatord for the parameter d is extracted from the sample meanx µ = µ i=1 x i /µ of µ independent measurements ofX. The distance estimator is obtained by comparing the sample meanx µ with a calibration curve given by the expectation value X = tr[Xρ(d, θ)] of the measurement oper-atorX as a function of the source separation d, which could be known either from theory, or from a previous calibration experiment.
When sufficiently many measurements are performed (µ 1), it follows from the central limit theorem that x µ is normally distributed with mean value X and variance (∆X) 2 /µ, with (∆X) 2 = X 2 − X 2 . Accordingly, we choose as a separation estimator, the parameter valued at which X equates the measurements meanx µ . The associated estimation error is given by (∆d) 2 = χ 2 [d, θ,X]/µ, with Moment-based estimation of the separation between the sources: the measured mean values of the available observables are linearly combined with optimal coefficients and compared with a calibration curve.
In practical situations, one does not have access to all possible measurement observables. Therefore, we now assume that we can measure only a finite number K of observablesX k , as well as their linear combinationŝ Xm =m ·X, withX = (X 1 , . . . ,X K ) T andm the measurement coefficients vector. Under these assumptions, it is possible to perform an analytical optimization over all possible linear combinations. Such an optimization yields the measurement sensitivity [29] (15) where Γ[d, θ,X] is the covariance matrix, whose elements are given by ∂d is the derivative vector, where we labelled as X k = X k the expectation value ofX k . The sensitivity given by Eq. (15) is achieved by the measurement coefficientsm = m [29], with where η is a normalization constant. The procedure discussed in this section to obtain the separation estimatioñ d from the optimal coefficients m in Eq. (16) is illustrated in Fig. 2.

IV. CONSTRUCTION OF THE OPTIMAL OBSERVABLE FOR DEMULTIPLEXING
Fisher-information based studies [18,25] suggested that photon counting after spatial mode demultiplexing allows to approach the ultimate limit for the separation estimation. Therefore, in this Section we assume that we have access to K spatial modes {v k (r)} with associated field operatorsâ k , and that we can measure the photon number in each modeN k =â † kâ k . We use the formalism presented in Sec. III to derive their optimal linear combinationN withm determined by equation (16). Once this optimal observable is known, one only need to access its average value. The latter can be derived by measuring the mean photon number in each of the modes {v k (r)} and then combining them according to the coefficients m..

A. Covariance matrix and derivative vector
We now use the P −function (11) to calculate the photon-number covariance matrix and derivative vector which are needed to compute the measurement sensitivity (15) and coefficients (16). To this goal, it is convenient to define an auxiliary mode basis {w i (r)} obtained by extending w 0 (r) = u + (r) and w 1 (r) = u − (r) to a complete orthonormal basis. Accordingly, for the field operatorsb i associated with this basis, we haveb 0 =b + andb 1 =b − . The field operatorsâ k in the measurement basis {v k (r)} can be expanded in terms of the field operatorsb k asâ k = l g klbl , with g kl = d 2 rv * k (r)w l (r). A one-dimensional comparison between the sources' images and the first two modes of the basis {w k (r)} and {v k (r)} is presented in Fig. 3.
Using the definition of the mean photon number in the measurement modes, we have where we used that allb i≥2 are in the vacuum. Analogously, the covariance matrix is given by The normally ordered expectation values of the field operatorsb i in Eqs. (18) and (19) can be calculated from the P −function (11) according to Given that P (β + , β − ) (see Eq. (11)) is a Gaussian involving only the modes u ± (r), the only non-zero expectation values are It is useful to write the mean photon number and the covariance matrices in terms of the sources images. For this purpose, we note that the coefficients g k± can be expressed in terms of the overlap functions with the source images f ±,k = drv * k (r)u(r ± r 0 ) as g k± = (f +,k ± f −,k )/ 2(1 ± δ), and using N ± = N κ(1 ± δ). Accordingly, inserting Eqs. (21) into Eqs. (18) and (19), using N ± = N κ(1 ± δ), we get the mean photon numbers (22) and the covariance matrix The derivative vector is obtained differentiating Eq. (22) with respect to the parameter d, which gives In the case of two equally bright sources (γ = 0), we recover the expressions given in [28]. Analogously, setting γ = 0 in Eq. (23), we have Γ kl [d, θ,N] = Γ 0 kl which coincides with the covariance matrix reported in [28].

B. Noise sources
For the rest of this work, we focus on the case of a Gaussian PSF u 0 (r) = 2/(πw 2 ) exp −|r| 2 /w 2 . For this PSF, the quantum Cramér-Rao bound can be approached by demultiplexing Hermite-Gauss (HG) modes with width matching that of the PSF [18,25], i.e. v k (r) = u k (r) ≡ u nm (r) with k = (n, m), such that u 00 (r) = u 0 (r). The HG modes are defined, for r = (x, y), as where H n (x) are the Hermite polynomials, and N nm = ((π/2)w 2 2 n+m n!m!) −1/2 is the normalization constant. For ideal measurements in the HG mode basis, we have the overlap functions f ±,k = β nm (±r 0 ), with that fully determine the covariance matrix (23) and the derivative vector (25), and consequently the measurement sensitivity (15) and coefficients (16). A detailed discussion of the performances of our estimation strategy for ideal measurements is reported in Sec. VI. We now study the impact of different noise sources on the covariance matrix (23) and the derivative vector (25). In particular, we consider misalignment between the demultiplexing basis and the source centroid, as well as crosstalks between the detection modes. These two imperfections affect Eqs. (23) and (25) only through modifications of the overlap funtions f ±,k . Finally, we consider dark counts at the detection stage, that modify the diagonal of the covariance matrix (23). A schematic illustration of how the different noise sources enter a demultiplexing measurement is presented in Fig. 4.

Misalignment
The assumption that the demultiplexing basis is perfectly centered with respect to the centroid of the two sources is often not true in practice. In fact, the source centroid is in general a priori unknown and needs to be pre-estimated, possibly via direct imaging, to optimally align the demultiplexer. In the case of faint sources, it was observed [18] that an imperfect positioning of the demultiplexer makes the Fisher information go to zero for small sources separations. Methods to mitigate this effect have been proposed by adaptively switching between demultiplexing and direct imaging [32] or by optimizing the detection basis [33]. An alternative approach, that we will not consider here, consists in estimating the source's centroid simultaneously with the separation [34].
Within our model, a two-dimensional shift r s = (d s cos θ s , d s sin θ s ) T of the centroid of the sources with respect to the demultiplexer axis enters the measurement sensitivity (15), and the coefficients (16) through the overlap functions f ±,k = β nm (±r 0 − r s ). As a trick, to compute the overlap functions β nm (a), we can use an analogy with the quantum mechanical harmonic oscillator. In particular, we can interpret the integral (27) as the overlap between the (n, m)−excited state of a two-dimensional harmonic oscillator and a coherent state (displaced to the phase-space coordinate a) of the same harmonic oscillator. We therefore obtain

Crosstalk
A recent experiment [24] identified imperfections in the mode decompositions as an important limitation for resolving the distance between incoherent point sources. For sources with low brightness, it was reported that in the presence of crosstalk in the demultiplexing basis the Fisher information drops to zero for small separations [35].
To include the impact of crosstalk between the detection modes in our model, we follow [35], and we describe it as a unitary matrix c kl that maps the ideal HG modes u l (r) into the actual measurement basis v k = l c kl u l (r). Accordingly, the overlap functions f ±,k are given by linear combinations of the functions (28), To assess the impact of crosstalk on the sensitivity of our method, as well as on the shape of our optimal observable, we numerically generate K×K unitary matrices by sampling the uniformly distributed random real coefficients λ i (0 ≤ λ i ≤ 1 and i λ 2 i = 1), that multiply the generalized Gell-Mann matrices G i [36]. The positive parameter allows us to tune the crosstalk strength. In experimentally relevant scenarios, crosstalk is generally weak [24,35]. Namely, the diagonal elements of the crosstalk matrix C ii ( ) = c ii are larger when compared to the off-diagonal ones C i =j ( ) = c i =j . Accordingly, we select to ensure In the following, when discussing results in the presence of crosstalk, we will consider multiple random realizations of the matrix c ij , and quantify the average crosstalk probability via the ensemble average |c ij | 2 .

Dark counts
Electronic noise at the detection stage introduces additional photon counts, the dark counts, which contain no information on the parameter value d. Consequently, the signal-to-noise ratio in each detection mode is reduced by dark counts. The impact of this noise source on distance estimation via spatial-mode demultiplexing has been investigated in the recent literature [37][38][39]. Despite the different approaches, all these works obtained the same qualitative result: dark counts cause the Fisher information to drop to zero for small source separations.
To include dark counts in our model, we add to the quantum mechanical photon number operator in the measurement modesN k a classical random variable ξ k . In particular, ξ k are thermally distributed with mean value To quantify the strength of dark counts with respect to the number of photons received in the image plane, we define σ k = N dc k /2N κ. The measured mean photon number in each detection mode N k and the covariance matrix Γ [d, θ,N] are now obtained not only by computing quantum mechanical expectation values, but also classical averages over the probability distribution of ξ k . Following this procedure, we obtain and, given that N dc k is independent on d, the derivative vector D[d, θ,N] is untouched by dark counts. On the other hand, the covariance matrix becomes with the extra diagonal term describing the noise induced by dark counts.

V. DIRECT IMAGING
In Secs. VI and VII, we will discuss the performances of ideal and noisy demultiplexing, respectively, and we will compare them with those of direct imaging. In order to perform such comparison, we now consider an ideal direct imaging system, and we evaluate its optimized moment-based sensitivity (15) that bounds the Fisher information according to Eq. (14).
Direct imaging estimates the separation between the two sources from the intensity distribution in the image plane. The mean intensity is given by In the first step, we used the expansion of the electric field in the basis {w k (r)} of the symmetric w 0 (r) = u + (r) and antisymmetric modes w 1 (r) = u − (r) (see Sec. IV), and the fact that, in the quantum state (8), all higher order modes w k≥2 (r) are in the vacuum state. In the second step, we used b † ±b± = N ± and b † ±b∓ = −γ N + N − (see Eqs. (21)), as well as the relation (4) between u 0 (r ± r 0 ) and u ± (r).
Let us now assume that the intensity measurements are performed with an ideal pixelized detector, with unity quantum efficiency and noiseless number-resolved photon counting at each pixel. The pixels are defined by the area with 2x p (2y p ) the horizontal (vertical) size of the pixels and (x i , y i ) the coordinate of the center of the pixel. We further assume that the total area A of the detector is fixed, and large enough to collect the full intensity in the image plane. Accordingly, the pixel area is fully determined by the number of segments on each axis N p , i.e.
x p × y p = A/N 2 p . The mean photon count per pixel is given by the mean intensity integrated over the pixel area with The intensity coherence function is defined as where in the last step we used the commutation re- . Using the expectation values (21), we can finally write the intensity coherence function as with The covariance matrix for intensity measurements with the ideal pixelized detector described above is obtained by integrating the coherence function (41) over the pixel area according to The derivative vector is obtained from the mean intensity per pixel (37), where the dependence on the parameter d is only contained in the functions Φ ij and Ψ ij (38), For a Gaussian PSF, the integrals (38) and (45) can be evaluated analytically, and result in In Eqs. (47), we have defined and The covariance matrix (43) and the derivatives vector (46), together with Eqs. (47) and (49), generalize the results in the Supplementary Material of [25] to arbitrary two-dimensional arrangements and unequal brightnesses of the two sources.
Let us now note that the covariance matrix (43), can be rewritten in the form where I is a diagonal matrix whose elements are given by the mean intensities per pixels I k=ij as reported in Eq. (37), and U = (ν (1) , ν (2) , ν (3) ) T is N 2 p × 3 matrix whose columns are given by the vectors where we used that according to Eqs. (47) the functions Φ ij , Ψ ij and Ξ ij are all real. Then, using the Woodbury matrix identity [40], we can write the inverse of the covariance matrix as where the correction Λ = U (1 3 + U T U ) −1 U T , with 1 3 the 3 × 3 matrix identity, can be evaluated numerically.
In our numerical calculations, independently of the mean photon number of the two thermal sources, we did not observe appreciable changes in the behaviour of the measurement sensitivity M [d, θ, N p ] for N p 50 × 50. Accordingly, N p = 50 × 50 was used to obtain all direct imaging curves reported in this work.
The representation (52) allows to write the optimal direct imaging sensitivity as where the first term is of order N κ, while the second one is of order (N κ) 2 and negative. As a consequence, the latter tends to reduce the relative sensitivity when increasing the source brightness. This behaviour can be observed in Fig. 5, where the direct imaging sensitivity M [d, θ, N p ] is plotted for different mean photon numbers of the sources.
The limit of continuous direct imaging is obtained for N p → ∞. In this limit, the first term in Eq. 53 becomes which coincides with the direct imaging Fisher information calculated assuming Poissonian sources in the N κ 1 regime [18]. Accordingly, when the number of photons in the image plane is low (N κ 1), when the second term in Eq. (53) is negligible, the optimized moment-based sensitivity saturates the Cramér-Rao bound for continuous direct imaging.
In Fig. 5, we see that the reduction of the measurement sensitivity due the Λ term in Eq. this regime, we can approximate the direct imaging sensitivity by considering an expansion of the integral (54) for d/2w 1. Such an expansion can be performed analytically and results in (55) The approximation (55) is compared (for γ = 0) with the exact numerical results for different mean photon numbers in Fig. 5.
It is interesting to note that for γ = 0, Eq. (55) does not vanish for d = 0 (see Fig. 6). Accordingly, an asymmetry in the intensity distribution between the two sources eases the separation estimation. In particular, we see that for extreme intensity imbalances γ → 1, the direct imaging sensitivity for small separations tends to its maximal value M DI [d, θ] → 2N κ/w 2 . We can understand this result considering that in the γ → 1 limit the estimation of the sources separation becomes the localization of a single source.

VI. IDEAL DEMULTIPLEXING
In this Section, we present analytical expressions for the optimal moment-based sensitivity, according to Eq. (15), attainable with ideal demultiplexing measurements and for the optimal observable achieving it, as given by Eq. (16).

A. Analytical inversion of the covariance matrix
We now focus on the case of ideal intensity measurements of K = (Q + 1) 2 HG modes u nm (r) with 0 ≤ n, m ≤ Q. In this case, the overlap functions β nm (±r 0 ) are obtained setting x s = y s = 0 in Eq. (28), which implies the symmetry β nm (−r 0 ) = (−1) n+m β nm (r 0 ).
Equation (56) causes the mean photon number in the measurement modes N k (see Eq. (22)) to be independent of the photon-number imbalance γ: Substituting Eq. (56) into Eq. (23) also simplifies the covariance matrix, which becomes with diag(ν) a diagonal K × K matrix whose elements are given by the vector ν k=(n,m) = N nm (see Eq. (57)), while the other two terms are outer products of vectors with elements ξ k=(n,m) = 1 + γ 2 N nm and ζ k=(n,m) = (−1) n+m 1 − γ 2 N nm . This particular form of the covariance matrix allows for its analytical inversion by two successive applications of the Sherman-Morrison formula [40]. Accordingly, we obtain Equation (59) is valid only under the assumption that the covariance matrix is invertible. According to the Sherman-Morrison formula, this is the case if and only if diag(ν) is invertible, i.e. when β nm (r 0 ) = 0 for every n, m. Looking at Eq. (28), for x s = y s = 0, we see that this condition is not satisfied for d = 0 or θ = kπ/2, with k = 0, 1, 2, . . . . Let us first focus on d = 0. In this case, β nm (r 0 ) = 0 for all n, m. However, this case corresponds to having only one source and the estimation problem under study is not well posed. On the other hand, for θ = kπ/2 with k even (odd), we have β nm (r 0 ) = 0 ∀n(m) = 0. We can therefore make the covariance matrix (58) invertible by removing all modes with n(m) = 0 from our measurement basis, and carry out the estimation within this smaller set of modes. This fact has a clear physical interpretation: The HG modes (see Eq. (26)) have the form u nm (x, y) = h n (x)h m (y), accordingly when θ is an odd (even) multiple of π/2 the two sources are aligned with the x−(y−)axis and the only relevant modes are u n0 (x, y) = h n (x) (u 0m (x, y) = h m (y)). Equation (59) yields the correct results also for θ = kπ/2 by extracting only the relevant matrix elements.

B. Measurement sensitivity
The measurement sensitivity M [d, θ,N] (15), which is the maximal sensitivity achievable by any linear combination of the mean photon numbers in the measurement modes N k [29], can be obtained by combining the inverse covariance matrix (59) with the derivative vector with x = d/2w, and is given by and We immediately see that, when the number of received photons is low (N κ 1), the sensitivity M [d, θ,N] is dominated by the F term. This term coincides with the Fisher information calculated assuming Poissonian sources in the N κ 1 regime (compare Eq. (63a) with Eq. (21) in the Supplementary Material of [35]). Moreover, F does not depend on γ, i.e. the number of received photons is low, the sensitivity in the estimation of d does not depend on the relative brightness of the two sources.
For higher received photon numbers, the terms of order (N κ) 2 become relevant. Given that the functions δ 1 , δ 2 and δ 3 are positive, the quadratic terms in the photon number always reduce the measurement sensitivity. Moreover, they depend on γ through A ± . Accordingly, the loss of sensitivity appearing for higher received photon numbers depends on how the photons are distributed between the two sources.
Let us now consider the limiting case where the full HG basis is measured, i.e. Q → ∞. In this limit, F → 1, δ 2 and δ 3 vanish, while δ 1 remains finite and provides a correction to the measurement sensitivity, which takes the form The equal brightness case is obtained by setting γ = 0 in Eq. (65), which results in which coincides with the quantum Fisher information for two equally bright thermal sources calculated in [25,26]. Accordingly, our estimation strategy saturates the quantum Cramér-Rao bound in the asymptotic limit. The second term in Eqs. (65) and (66) reduces the estimation sensitivity for intermediate separations, i.e. a dip appears in Figs. 7 and 8. This dip gets deeper when increasing the mean photon number N κ (see also [25,26]). Physically, this is due to the fact, that for higher brightness of the sources there is a non negligible probability of detecting multiple photons in the same modes which reduces the information that can be extracted from each photon.
The behaviour of the measurement sensitivity (62) for γ = 0 and different values of Q is compared with the quantum Fisher information (66) in Fig. 7. While the asymptotic expression (66) is independent on the alignment angle θ, this is not true at finite Q. In particular, M [d, θ,N] is minimal when the sources are aligned with the axes (θ = 0, π/2), and maximal when the sources are aligned along a bisector (θ = π/4, 3π/4). This behaviour can be observed in Fig. 7, where in panel (a) (θ = 0), for large separations d, the measurement sensitivity is clearly distinguishable from the quantum Fisher information for all considered values of Q. On the other hand, in panel (b) (θ = π/4), the measurement sensitivity cannot be discerned from the quantum Fisher information for Q ≥ 9. A similar behaviour was reported for the Fisher information, for the regime of low photon numbers in the image plane, in the Supplementary Material of [35]. The optimality of θ = π/4 is due to the fact that we assumed that we have a limited and equal number of modes available in each direction. Under this assumption, aligning the sources in one direction, we can use only Q of the Q 2 available modes. This limits the resolution especially for large separations. On the other hand, aligning the sources along the bisector, is equivalent to rotating the basis in order to have all Q 2 available modes along the sources axis, which provides the optimal resolution.
The asymptotic behaviour (Q → ∞) of the measurement sensitivity (65) for different brightnesses of the two sources (γ = 0) is illustrated in Fig. 8. We can see that the dip at intermediate separations gets shallower when increasing γ, and it disappears in the limit γ → 1. Accordingly, as observed for direct imaging in Sec. V, an intensity unbalance between the sources makes it easier to estimate their separation.

C. Measurements coefficients
Combining Eqs. (59) and (61) according to Eq. (16), we obtain the expression for the measurement coefficients (dashed lines in Fig. 11 in Sec. VII) From Eq. (67), we note that the coefficients of the optimal observable only depend on the sum of the indices i + j. Accordingly, independently on the sources orientations or their relative intensities, the method of moments prescribes to measure with the same weight the photon number in all HG modes of the same order. In other words, in the absence of noise, it is optimal to measure a circularly symmetric intensity distribution, which achieves the same sensitivity for each alignment angle θ. The above observation does not apply when the sources are aligned along the x−(y−)axis. In these cases, as discussed in Sec. VI A, one needs only to measure modes with i(j) = 0. The correct coefficients can be still obtained from Eq. (67) by setting i = 0(j = 0). Further comments on the behaviour of the measurement coefficients (67) is postponed to Sec. VII, where we will compare them with their counterpart in presence of noise.

A. Measurement sensitivity
Let us start our discussion on the sensitivity of our estimation strategy in the presence of noise by commenting on its optimality. We have already observed in Sec. VI that for ideal demultiplexing our method approaches the quantum Cramér-Rao bound for arbitrary source separations and brightnesses, when a sufficiently large number of modes is measured. Moreover, when the number of received photons is low (N κ 1), the covariance matrix (23) is dominated by its diagonal terms, and we have (69) As discussed in Sec. IV B, crosstalk and misalignment only affect the overlap functions f ±,k , while dark counts modify the diagonal of the covariance matrix. As a consequence, when all noise sources are considered, Eq. (69) remains valid if we substitute the mean photon number N k with the mean photon number plus noise: N k = N k + N dc k . Equation (69) is equal to the Fisher information for demultiplexing calculated assuming Poissonian sources in the low brightness regime (N κ 1) [18]. Accordingly, for N κ 1 our estimation strategy saturates the Cramér-Rao bound even in the presence of noise.
For thermal sources of arbitrary brightnesses and demultiplexing (noisy or not) in a finite number K of spatial modes, the Fisher information is unknown. In this regime, the sensitivity of our method provides a lower bound for the Fisher information. The latter can be saturated in practice by applying the simple estimation strategy presented in Sec. III to the optimal linear combination of mean photon numbers in the measurement modes discussed in Sec. VII B. The measurement sensitivity M [d, θ,N] is presented in Fig. 9 for different numbers of measured modes and different noise levels. We can see that all noise sources reduce the sensitivity with respect to its ideal value (62) (red and blue dashed curves in Fig. 9). In particular, we have that for equally bright sources, M [d, θ,N] vanishes for d → 0, and accordingly it gets harder to resolve small source separations.
Despite this loss of sensitivity for small separations, even when all noise sources are considered at the same time, as in Fig. 9 (d), demultiplexing outperforms ideal direct imaging (green lines). The regime where demultiplexing provides an advantage over direct imaging significantly depends on the brightnesses of the sources, and the noise levels. In particular, for large mean photon numbers and low noise, there are multiple crossings between the demultiplexing and the ideal direct imaging curves. To illustrate this behaviour, in Fig. 10, we plot, together with the ideal direct imaging curve (green), two curves corresponding to fixed misalignment and dark count levels but different crosstalk strengths for N κ = 10. We see that at low crosstalk levels, (red curve in Fig. 10) the demultiplexing sensitivity crosses the ideal direct imaging curve three times. Therefore, in the low noise regime, there are regions of larger separations where demultiplexing outperforms direct imaging. On the other hand, larger noise levels (blue curve in Fig. 10) reduce the demultiplexing sensitivity and cancel the larger-distance region where demultiplexing outperforms direct imaging. In Fig. 3 of [28], we investigated how the minimal separation at which the direct imaging sensitivity crosses the demultiplexing sensitivity varies with the noise level for different source brightnesses.
To better understand how the different noise sources affect the sensitivity, we now perform a perturbative expansion of M [d, θ,N] for small separations (x = d/2w 1). For simplicity, we focus on the case of equally bright thermal sources (γ = 0). Moreover, in practically relevant scenarios, we can assume that misalignment is of the order of the separation, or smaller. Therefore, we assume x = d/2w 1 and x s = d s /2w 1 to be of order a = ||(x, x s )|| 1 . From Eq. (28), we then obtain Combining Eqs. (29) and (25) + sin(2θ) c nm,11 c * nm,00 + c nm,10 c * nm,01 + c nm,01 c * nm,10 + c nm,00 c * nm,11 where we note that up to first order in a, misalignment does not affect the derivative vector.
The inverse of the 3 × 3 matrix in Eq. (72) can be obtained analytically and used to determine the measurement sensitivity, which results in 2N κ(|c 00,00 | 4 + σ 2 ) + |c 00,00 | 2 + σ + sin 4 θ |c 01,01 | 4 + cos 4 θ |c 10,10 | where, under the above assumption that detection noise is the same in all modes, we have σ = N dc /2N κ. In the above expansion, we considered σ of order zero in . Accordingly, in Eq. (73), the behaviour of the sensitivity for small x is dominated by dark counts (only diagonal terms of the cross-talk matrix appears).
Another relevant scenario might be the one where crosstalk and dark counts give contributions of the same order to the photon counts N nm in a given HG mode, i.e. when 2N κ|c ij | 2 ∼ N dc , or in other words when σ = N dc /2N κ is of order 2 . In this case, the mea- In this case, the probabilities of scattering from the first order modes u 01 and u 10 to the fundamental mode u 00 , enters in the behaviour of the measurement sensitivity for small x. Setting σ = 0 (corresponding to no dark counts, N dc = 0) in Eq. (74), we obtain the leading order expansion of M [d, θ,N] when only crosstalk is present (see also the Supplementary Material of [35]). Equations (73) and (74) can be simplified by considering a uniform crosstalk model [35], namely by setting all diagonal entries of the crosstalk matrix to t and all the off-diagonal ones to r such that |t| 2 +(D −1)|r| 2 = 1. For weak crosstalk, we have in addition |t| 2 ≈ 1 and |r| 2 1. Accordingly, Eqs. (73) and (74) becomes respectively We labelled as M dc (M ct ) the dark-counts (cross talk) dominated sensitivity. Equation (75) remains valid also in the case when no crosstalk or misalignment are present, and we expand to leading order in σ. Interestingly enough, with the assumption that the misalignment is of the same order as the separation, in both Eqs. (73) and (74) the effect of misalignment is hidden by the other noise sources. However, in Fig. 9 (a), we know that misalignment alone also cause the measurement sensitivity to go to zero for x → 0. To illustrate how this happens, we repeat the procedure described above in the absence of crosstalk and dark counts. The measurement sensitivity, at leading order in a, results in sin 4 θ x 2 sin 2 θ + 4x 2 s sin 2 θ s + cos 4 θ x 2 cos 2 θ + 4x 2 s cos 2 θ s  The approximations of the measurement sensitivity in the different noise regimes described by Eqs. (73) -(77) are compared with the exact numerical results in Fig. 9.
Interestingly, for all noise sources (Eqs. (73) -(77)) the measurement sensitivity goes to zero as x 2 for small separations. We observed, in Eq. (55), that also the direct imaging sensitivity (for γ = 0) presents the same scaling: However, for typical noise levels (as reported in Fig. 9), demultiplexing has generally a more favourable coefficient.

B. Measurement coefficients
We now discuss the behaviour of the optimal observable that allows to reach the sensitivity bounds illustrated in the previous section. To this goal, in Fig. 11, we study the dependence of the coefficients m ij (see Eq. (16)) of the optimal linear combination of intensity measurementsX m = ij m ij N ij in the HG modes u ij (r) as a function of the separation d between the two sources for 0 ≤ i, j ≤ Q with Q = 1 ( Fig. 11 (a-d)) and Q = 2 ( Fig. 11 (f)).
Arguably, the most interesting feature of Fig. 11 is that when d is small compared to the diameter 2w of the PSF u 0 (r) (d/2w 0.2), the coefficients m ij depend very weakly on d. Accordingly, in the relevant regime of small separations, a fixed observable can be used to estimate a vast range of parameter values. An exception to this behaviour is observed in Fig. 11 (a) and in its zoom, Fig. 11 (d), where we show the coefficients in presence of misalignment only. In this case, we observe that for separations of the order of the misalignment, the measurement coefficients present some modulations. In particular, for d ∼ d s , the shifted centroid significantly impacts the image decomposition in the u 01 (r) and u 10 (r) modes. Accordingly, the coefficients for these modes are slightly depleted, while the coefficient of u 11 (r) mode is increased. However, when misalignment is combined with other imperfections ( Fig. 11 (f)) these modulations are washed out by the reduction in signal-to-noise ratio induced by crosstalk and dark counts in higher order modes. In Sec. VI, we observed that, in the absence of noise, the measurement coefficients m ij only depend on the sum i + j (see Eq. (67)). In Fig. 11, we show how different noise sources can break this degeneracy. This effect is especially clear for the coefficients m 01 (blue) and m 10 (green). In particular, we see that when the center of the measurement basis does not lie on the sources axis (θ s = θ, dot-dashed lines in Fig. 11 (a) and (d)) m 01 and m 10 deviate in opposite directions from the θ d = θ = π/4 curve favouring the mode where the signal increased because of the centroid shift. Also crosstalk, which affects randomly the different modes, can remove the coefficients' degeneracy, as can be observed in the small difference between the green and blue lines in Fig. 11 (b), or in the more evident separation of the yellow and pink curves (m 02 and m 20 ) in Fig. 11 (f). On the other hand, in Fig. 11 (c), where the dark count level was assumed to be the same in all modes, the coefficients m 01 and m 10 are perfectly degenerate. Of course, this would have not been the case, if the electronic noise level had not been the same in the different modes. Finally, in Fig. 11 (f), we notice that in the presence of weak noise, even though the degeneracy is removed, the coefficients m ij still move in groups (green-blue for i + j = 1, orange-yellow-pink for i + j = 2, brown-grey for i + j = 3 and purple for i + j = 4).
From Fig. 11 (f), we can also see how the different coefficients change with the separation d between the two sources. First of all, for small separations, the mode u 00 contains no information on d, accordingly m 00 = 0 for a vast range of separations. Moreover, in the absence of noise (dashed lines Fig. 11 (f)), all coefficients m ij>0 (67) are different from zero for all values of d, and higher-order modes have larger weights. In fact, even though, for small separations, the image of the two sources produces very low signals in the higher-order modes, these signals are noiseless. Accordingly, the optimal observable amplifies these small signals to extract the most information on the parameter d out of them. On the other hand, different noise sources introduce photon counts that contain no information on the separation. Accordingly, in the presence of noise (solid lines Fig. 11 (f)), our optimal observable prescribes to measure only those modes where the signal-to-noise ratio is high enough. In particular, for small values of d, all coefficients are zero (or slightly negative) except m 01 and m 10 . Therefore, demultiplexing into the HG modes u 01 (r) and u 10 (r) (blue and green) is sufficient to achieve the optimal resolution. For larger separations, the other available modes start to become relevant, at first u 20 (r), u 11 (r), and u 02 (r) (orange, yellow and pink), then u 12 (r), and u 21 (r) (brown and grey), and finally u 22 (r) (purple). Increasing the separation, the higher-order modes become dominant, and all coefficients tend to their values in the absence of noise (the latter effect is visible for m 01 and m 10 in Fig. 11 (f)).

VIII. MINIMAL RESOLVABLE DISTANCE
We conclude this work by discussing how the minimal distance, that can be resolved with our moment-based estimation strategy, scales with the number of detected photons. For the sake of simplicity, we focus here on thermal sources with equal brightness (γ = 0).
We consider the distance between the two sources to be resolvable as long as the estimation error ∆d is smaller than the value d of the parameter itself, namely when d min /∆d ≥ 1 [35]. Following our moment-based estimation approach, we have (∆d) 2 = 1/µM [d, θ,N], accordingly the minimal resolvable distance is defined by the condition Given two thermal sources each emitting on average N photons, an imaging system with a transmissivity κ, after µ measurements of our optimal observable, the mean number of detected photons is given by N det = µ2N κ [41]. Accordingly, for a fixed imaging system, N det is affected either by µ, which could be changed by increasing or decreasing the detection time, or by the brightness of the sources N . In Fig. 12, we study how for noisy demultiplexing, the solution of Eq. (78) changes when we vary N det either by increasing µ for N κ = 1 (blue curves) or by changing N κ for µ = 1 (red curves). While in the first case, every noise source induces the same scaling d min ∼ N 1/4 det , in the latter case, the non trivial dependence of M [d min , θ,N] on N leads to more complicated behaviours. This difference is due to the thermal statistics of our sources. On the other hand, in the case of coherent sources there would not be any difference, since measuring µ copies of a coherent state |α give the same results as measuring |µα .
To get an analytical understanding of these behaviours, we use the small-separation expansions of the measurement sensitivity M [d, θ,N]. For ideal demultiplexing, our moment-based estimation strategy saturates the quantum Fisher information (62). In the small d regime, the latter is dominated by the constant term 2N κ/w 2 , which defines the solution of Eq. (78) Accordingly, ideal demultiplexing provides a "shot noise" scaling of the minimal resolvable distance which is represented as black lines in Fig. 12.
Let us now consider the impact of the different noise sources on the minimal resolvable distance achievable with demultiplexing measurements. For small values of N det , the minimal resolvable distance is determined by the expansion of M [d, θ,N] around a finite value of d. Independently on the noise level, the leading term of such an expansion is independent on d. As a consequence, we have the d min ∼ 1/ √ N det scaling, we observed in the case of ideal demultiplexing. Increasing N det , the solution of Eq. (78) is determined by the quadratic behaviour of M [d, θ,N] around d ∼ 0 (see Eqs. (75) -(77)), which induces the changes of scaling observed in Fig. 12 (similar considerations can be found in [35]).
To evaluate the impact of misalignment on the large N det scaling minimal resolvable distance, we use Eq. (77), which we recall is valid for misalignment of the order of the separation d s ≈ d 1, and no crosstalk and dark counts (see Table I). This expression leads to which is represented as a blue dashed line in Fig. 12 (a). A similar √ d s N −1/4 det scaling was reported in [33] where misalignment was studied in the low brightness regime. We note that when N is increased for µ = 1 the red curve in Fig. 12 (a) presents an almost flat region before falling on a N −1/4 det line with a less favourable coefficient than the one predicted by Eq. (80). This different scaling is due to the second order terms in N κ which are present in the offdiagonal elements of the covariance matrix (23) whose inverse determines the measurement sensitivity M [d, θ,N]. The analytical expression (80) cannot capture this behaviour. In fact, it was derived from Eq. (77) which was obtained under the assumption that misalignment is of the same order of the separation d ∼ d s ∼ a. The latter condition, which does not hold for the minimal resolvable distance when N det is large, causes the higher terms in N κ to disappear from M [d, θ,N] when expanding to leading order in a. On the other hand, for low photon numbers N κ 1, Eq. (77) is accurate. Accordingly Eq. (80) perfectly captures the scaling with the detected photon number when N det is increased by changing the number of measurements µ with N κ = 1 (blue curve).
Let us now discuss the effect of crosstalk on the minimal resolvable distance. To do so, we consider the uniform crosstalk model, and we insert Eq. (76) with σ = 0 (no dark counts) into Eq. (78), which results in Equation (81) is presented as blue dashed lines in Fig. (12) (b), and coincides with the results obtained in [35] for Poissonian sources. We can see that the blue (fixed N κ = 1 and varying µ) curve approach Eq. (81) for large N d , while the red (fixed µ = 1 and varying N κ) lies slightly above it. On the other hand, bright thermal states are more noisy, and this is reflected in the broad red error bands for large values of N κ.
Finally, let us focus on the role of dark counts. In this case, the large N det scaling of the minimal resolvable distance can be obtained by inserting Eq. (75) into Eq. (78). Here, we considered the dark counts level N dc to be constant as we increase N det , which results in cos(4θ) + 3 4N dc (N dc + 1) + 1 N dc (N dc + 1) + 2N κ(2N κ + 1) From this expression, we can observe that when the number of detected photons N det is increased by increasing µ, the minimal resolvable distance scales as d min ∼ N −1/4 det . On the other hand, when N det is increased with the source brightness, we obtain the same scaling that we had in the ideal case N −1/2 det . This has a very clear physical interpretation: if we detect N det photons by accumulating several measurements (µ 1) at low brightnesses (N κ 1), we also accumulate dark counts. On the other hand, in a single shot measurement (µ 1) at high brightness (N κ 1) with the same detector dark counts become negligible, we obtain the scaling of an ideal measurement. These different behaviours can be observed by comparing the scaling of the red (µ = 1, variable N ) and blue lines (2N κ = 1, variable µ) in Fig. 12 (c).
Finally, let us mention that even though Eq. (83) predicts the correct scaling with N and µ, it tends to underestimate d min . This is due to the quadratic approximation (75) which, as visible in Fig. 9 (c), grows faster than the exact measurement sensitivity, and consequently d µM dc [d, θ,N] crosses 1 (see Eq. (78)) earlier than the exact curve d µM [d, θ,N]. This underestimation is particularly relevant for low dark counts levels, N dc 1, and in particular it leads to an unphysical better-than-ideal scaling for N dc < ( 4 + 2/µ − 2)/4. In Fig. 12 (d), we present the behaviour of the minimal resolvable distance, when all noise sources are present at the same time. As discussed above, when the number of detected photons N det is increased by accumulating many measurements µ (blue curves), we are also increasing the dark counts. In this case, dark counts are the dominant noise source (compare with the blue curve in Fig. 12 (c)), and the error bands due to crosstalk become negligible. On the other hand, when N det is increased with the brightness of the sources N (red curve in Fig. 12 (d)), dark counts quickly become negligible, and the dominant noise source is crosstalk (compare with the red curve in Fig. 12 (c)).
In conclusion, let us compare the noisy demultiplexing results with those of direct imaging. The large N d scaling of the minimal resolvable distance for noiseless direct imaging can be obtained from Eq. (55), which, for γ = 0, results in M [d, θ, N p ] ≈ 8x 2 , and leads to (see also [35]) Equation (84) is plotted as green lines in Fig. 12 and for a vast range of noise conditions it stays above the noisy demultiplexing curves. An exception is represented by the blue curves in Fig. 12 (c) and (d). However, in Fig. 12 (c) and (d), we chose N dc = 1 for which the approximation (83) gives sensible results. For large source brightnesses, this value corresponds to a reasonable ratio between dark counts and photons counts in each demultiplexed mode.
On the other hand, it is overpessimistic for 2N κ = 1, as it can be quickly understood by considering, for example, that N 01 /2N κ ≈ 10 −5 at x = d/2w ≈ 10 −2 . Therefore, we can conclude that, in most practical situations, spatial mode demultiplexing allows to resolve significantly smaller separations than direct imaging.

IX. CONCLUSION
We discussed in detail a method to extract the separation of two thermal sources, with arbitrary and possibly different brightnesses, from a single optimized measurement observable. For this imaging application, we considered different measurements such as direct imaging with pixel detectors and realistic spatial mode demultiplexing. Our results show how, even in presence of different relevant noise sources, demultiplexing allows for better resolutions than direct imaging. In several realistic scenarios, we constructed the optimal observable for demultiplexing measurements, which depends very weakly on the separation d, providing stability of the estimation procedure over a vast domain of separations. Moreover, in the limiting case of noiseless demultiplexing, for arbitrary received photon numbers and separations, the optimized observable approaches the quantum Cramér-Rao bound if sufficiently many modes are measured. Finally, for low photon numbers in the image plane, the sensitivity M [d, θ,N] of our method saturates the Fisher infor-mation even for noisy demultiplexing.
From our results is evident that the coefficients of the optimal linear combination of photon number measurements are severely affected by noise. This dependence must be taken into account if one wants to achieve optimality. This could be done by calculating the coefficients from experimental measurements of the covariance matrix and the derivative vector, obtained with the help of two test sources. Alternatively, one could characterize the different noise sources, and then compute the coefficients theoretically. In this latter case, a precise measure of the crosstalk matrix of a demultiplexer, and the dark count level of detectors are not problematic. On the other hand, misalignment errors are mainly due to an imperfect knowledge of the source centroid. These errors may be limited by scanning the multiplexer position, or by using adaptive methods such as the one proposed in [32].
Finally, let us point out that in experimental realizations (see for example [24]) all available observables can be measured at the same time which allows to look at arbitrary linear combinations in post processing. Therefore, the estimation technique presented here can be understood as the optimal post-processing technique that makes the best use of the available data. As a consequence, our approach fits particularly well for the estimation of a dynamically changing parameter. In fact, provided that the photon detectors are faster than the typical time scale of the parameter changes, the optimal measurement coefficients at each time can be selected in post-processing.