Giacomo Sorelli

Manuel Gessner

Mattia Walschaers

Nicolas Treps

Moment-based superresolution: Formalism and applications

Sensitivity limits are usually determined using the Cramér-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ér-Rao bound. We also investigate the impact of noise on the optimal observables, their measurement sensitivity, and 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ér-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 by Abbe [START_REF] Abbe | Beiträge zur theorie des mikroskops und der mikroskopischen wahrnehmung[END_REF] and Rayleigh [START_REF] Rayleigh | Investigations in optics, with special reference to the spectroscope[END_REF] at the end of the nineteenth century, is limited by diffraction and noise. However, for a large enough signal-to-noise ratio, the diffraction limit is not fundamental [START_REF] Goodman | Statistical Optics[END_REF] and can be overcome by superresolution techniques based on fluorescence microscopy [START_REF] Hell | Breaking the diffraction resolution limit by stimulated emission: Stimulated-emission-depletion fluorescence microscopy[END_REF][START_REF] Klar | Fluorescence microscopy with diffraction resolution barrier broken by stimulated emission[END_REF][START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF], homodyne measurements [START_REF] Hsu | Optimal optical measurement of small displacements[END_REF][START_REF] Delaubert | Quantum measurements of spatial conjugate variables: Displacement and tilt of a Gaussian beam[END_REF][START_REF] Pinel | Ultimate sensitivity of precision measurements with intense Gaussian quantum light: A multimodal approach[END_REF], or intensity measurements in an appropriate spatial-mode basis [START_REF] Helstrom | Resolution of point sources of light as analyzed by quantum detection theory[END_REF][START_REF] Tsang | Resolving starlight: A quantum perspective[END_REF].

Superresolution imaging has been analyzed from the point of view of quantum metrology [START_REF] Helstrom | Quantum Detection and Estimation Theory[END_REF][START_REF] Braunstein | Statistical Distance and the Geometry of Quantum States[END_REF][START_REF] Holevo | Probabilistic and Statistical Aspects of Quantum Theory[END_REF][START_REF] Paris | Quantum estimation for quantum technology[END_REF][START_REF] Giovannetti | Advances in quantum metrology[END_REF][START_REF] Pezzè | Quantum Theory of Phase Estimation[END_REF] 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 [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF]. 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 [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF], leaving room for superresolution. Furthermore, it was shown [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF] 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 implement a simplified version of this demultiplexing technique, which distinguishes only between the fundamental and the first excited modes [START_REF] Paúr | Achieving the ultimate optical resolution[END_REF][START_REF] Tang | Fault-tolerant and finiteerror localization for point emitters within the diffraction limit[END_REF][START_REF] Yang | Far-field linear optical superresolution via heterodyne detection in a higher-order local oscillator mode[END_REF][START_REF] Tham | Beating Rayleigh's Curse by Imaging Using Phase Information[END_REF]. On the other hand, distance estimation from the demultiplexing of multiple spatial modes, realized using mul-tiplane light conversion [START_REF] Morizur | Programmable unitary spatial mode manipulation[END_REF], was only recently reported [START_REF] Boucher | Spatial optical mode demultiplexing as a practical tool for optimal transverse distance estimation[END_REF].

In general, reaching the Cramér-Rao bound, e.g., via maximum-likelihood estimation, requires the measurement of 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, which requires one 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 the 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 us to reach the quantum Cramér-Rao bound [START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF][START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF], 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 Secs. VI and VII, respectively. In Sec. VIII we discuss the smallest source separation which is resolvable by the different measurement techniques. Section IX summarizes our work.

II. THERMAL STATES IN THE IMAGE PLANE

We want to estimate the transverse distance between two thermal point 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 by ρ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 state ρ0

= ρs 1 [(1 -γ )N] ⊗ ρs 2 [(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 s 1 ,s 2 (α 1 , α 2 )|α 1 , α 2 α 1 , α 2 |, (1)
where |α 1 [START_REF] Rayleigh | Investigations in optics, with special reference to the spectroscope[END_REF] are coherent states of the field operators ŝ1 [START_REF] Rayleigh | Investigations in optics, with special reference to the spectroscope[END_REF] and P s 1 ,s 2 (α 1 , α 2 ) = P s 1 (α 1 )P s 2 (α 2 ) is the Glauber-Sudarshan P function, with

P s 1 (2) (α 1 (2) ) = 1 π N (1 ∓ γ ) e -|α 1 (2) | 2 /(1∓γ )N . ( 2
)
The evolution of the field operators through a diffractionlimited imaging system, with transmissivity κ, is given by [START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF][START_REF] Shapiro | The quantum theory of optical communications[END_REF]

ĉ1 = √ κ ŝ1 + √ 1 -κ v1 , (3a) ĉ2 = √ κ ŝ2 + √ 1 -κ v2 , ( 3b
)
where ĉ1 [START_REF] Rayleigh | Investigations in optics, with special reference to the spectroscope[END_REF] 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 operators v1 [START_REF] Rayleigh | Investigations in optics, with special reference to the spectroscope[END_REF] are associated with auxiliary modes, which 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 nonorthogonal 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

u ± (r) = u 0 (r + r 0 ) ± u 0 (r -r 0 ) √ 2(1 ± δ) , ( 4
)
where δ represents the overlap between the source images The relation between the field operators b± associated with the modes u ± (r) and the field operators ŝ± = (ŝ 1 ± ŝ2 )/ √ 2 in the object plane can be obtained from the sum and difference of Eqs. (3) [START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF],

δ = d 2 r u * 0 (r + r 0 )u 0 (r -r 0 ). (5)
b± = √ κ ± ŝ± + 1 -κ ± v± , (6)
with v± associated with auxiliary modes, which are in the vacuum state, and κ ± = κ (1 ± δ).

We now use Eq. ( 6) to propagate the quantum state ρ0 of the sources to the image plane. First, we note that the transformation to the modes ŝ± can be interpreted as a 50:50 beam splitter

ŝ+ ŝ- = 1 √ 2 1 1 1 -1 ŝ1 ŝ2 . ( 7
)
According to Eq. ( 7), the coherent states |α 1 |α 2 correspond to the coherent states

| α 1 +α 2 √ 2 | α 1 -α 2 √ 2
≡ |α + |α -of the field operators ŝ± [see Fig. 1(b) for a schematic illustration]. The quantum state of the sources can therefore be rewritten in terms of the latter coherent states as

ρ0 = d 2 α + d 2 α -P s + ,s -(α + , α -)|α + , α -α + , α -|, (8)
where

P s + ,s -(α + , α -) = P s 1 ,s 2 ( α + +α - √ 2 , α + -α - √ 2 )
. Accordingly, the modes u ± (r) have both mean photon numbers ŝ † ± ŝ± = N and the photon-number imbalance γ appears in the classical phase-insensitive correlations ŝ † ± ŝ∓ = Nγ [see orange lines in Fig. 1(b)]. Going through the imaging system according to Eq. ( 6), we have

|α + , α -→ | √ κ + α + , √ κ -α -≡ |β + , β -,
with |β ± coherent states of the field operators b± . We can therefore write the quantum state in the image plane as

ρ(d, θ ) = d 2 β + d 2 β -P b + ,b -(β + , β -)|β + , β -β + , β -|, (9)
with

P b + ,b -(β + , β -) = 1 κ + κ - P s + ,s - β + √ κ + , β - √ κ - . ( 10
)
Combining Eqs. ( 10) and ( 2), we can write the P function as

P b + ,b -(β + , β -) = 1 π 2 det V e -β † V β , ( 11
)
where we have defined β = (β + , β -) T and

V = N + γ √ N + N - γ √ N + N - N - , (12)
with N ± = Nκ ± . For equally bright sources (γ = 0), the offdiagonal elements of the matrix V vanish and Eq. ( 11) reduces to the product of two Gaussian functions corresponding to ρ(θ, d ) = ρb + (N + ) ⊗ ρb -(N -), as reported in Ref. [START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF]. Finally, the sources in the image plane are described by Eqs. ( 9) and ( 11), with the information on the parameter d contained in the shape of the modes u ± (r) and the mean photon numbers N ± .

III. METHOD OF MOMENTS

We estimate the distance d between the two sources with the method of moments [START_REF] Pezzè | Quantum Theory of Phase Estimation[END_REF][START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF]. Following this approach, given an observable X , an estimator d for the parameter d is extracted from the sample mean xμ = μ i=1 x i /μ of μ independent measurements of X . The distance estimator is obtained by comparing the sample mean xμ with a calibration curve given by the expectation value X = tr[ X ρ(d, θ )] of the measurement operator X 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 value d at which X equates the measurement mean xμ . The associated estimation error is given by ( d

) 2 = χ 2 [d, θ, X ]/μ, with χ 2 [d, θ, X ] = ( X ) 2 ∂ X ∂d 2 . ( 13
)
The quantity χ 2 [d, θ, X ] determines the sensitivity of the method-of-moments estimation strategy for the quantum state ρ(d, θ ) and the observable X . It obeys the chain of inequalities

χ -2 [d, θ, X ] F[d, θ, X ] F Q [d, θ], ( 14
)
where the second is saturable by an optimal operator X [START_REF] Braunstein | Statistical Distance and the Geometry of Quantum States[END_REF][START_REF] Pezzè | Quantum Theory of Phase Estimation[END_REF][START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF][START_REF] Kholevo | A generalization of the Rao-Cramer inequality[END_REF][START_REF] Fröwis | Tighter quantum uncertainty relations following from a general probabilistic bound[END_REF]. Here F[d, θ, X ] is the Fisher information that bounds the achievable sensitivity when estimating d from measurements of X according to the Cramér-Rao lower bound ( d ) 2 (μF[d, θ, X ]) -1 [START_REF] Helstrom | Quantum Detection and Estimation Theory[END_REF][START_REF] Braunstein | Statistical Distance and the Geometry of Quantum States[END_REF]. Finally,

F Q [d, θ] = max X (F[d, θ, X ]
) is the quantum Fisher information which defines the ultimate metrological sensitivity [START_REF] Braunstein | Statistical Distance and the Geometry of Quantum States[END_REF].

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 observables

X 1 X 2 X 3 . . . X K m 1 m 2 m 3 . . . m K xµ d X ΔX d ± Δd FIG. 2.
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.

Xk , as well as their linear combinations X m = m • X, with X = ( X1 , . . . , XK ) T and m 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 [START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF] [START_REF] Paris | Quantum estimation for quantum technology[END_REF] where [d, θ, X] is the covariance matrix, whose elements are given by k,l [d, θ, X] = Xk Xl -X k X l , and D[d, θ, X] = ∂X k ∂d is the derivative vector, where we defined as X k = Xk the expectation value of Xk . The sensitivity given by Eq. ( 15) is achieved by the measurement coefficients m = m [START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF], with

M[d, θ, X] = max m (χ -2 [d, θ, X m]) = D[d, θ, X] T [d, θ, X] -1 D[d, θ, X],
m = η -1 [d, θ, X]D[d, θ, X], ( 16
)
where η is a normalization constant. The procedure discussed in this section to obtain the separation estimation 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 [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF][START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF] suggested that photon counting after spatial-mode demultiplexing allows one 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 mode Nk = â † k âk . We use the formalism presented in Sec. III to derive their optimal linear combination

Nm = m • N = K i=1 m i Ni , ( 17
)
with m determined by Eq. [START_REF] Giovannetti | Advances in quantum metrology[END_REF]. Once this optimal observable is known, we 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.

u 0 (r -r 0 )

u 0 (r + r 0 ) Sources' images w 0 (r) = u + (r) w 1 (r) = u -(r) Auxiliary modes v 0 (r) = u 0 (r) v 1 (r) = u 1 (r) Measurement modes FIG.
3. The images of the two sources u 0 (r ± r 0 ) [red (left)] are compared with the first two modes of the auxiliary mode basis {w i (r)} [blue (middle)] and the ideal HG measurement basis {v i (r) = u i (r)} [green (right)]. The modes {v i (r)} are modified by misalignment and crosstalk.

A. Covariance matrix and derivative vector

We now use the P function [START_REF] Tsang | Resolving starlight: A quantum perspective[END_REF] to calculate the photonnumber covariance matrix and derivative vector which are needed to compute the measurement sensitivity [START_REF] Paris | Quantum estimation for quantum technology[END_REF] and coefficients [START_REF] Giovannetti | Advances in quantum metrology[END_REF]. For 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 operators bi associated with this basis, we have b0 = b+ and b1 = b-. The field operators âk in the measurement basis {v k (r)} can be expanded in terms of the field operators bk as âk = l g kl bl , with g kl = d 2 r v * 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

N k = Nk = â † k âk = i j=± g * ki g k j b † i b j , ( 18
)
where we used that all bi 2 are in the vacuum. Analogously, the covariance matrix is given by kl

[d, θ, N] = Nk Nl - N k N l , with Nk Nl = â † k â † l âk âl + δ kl â † k âk = mnpq=± g * mk g * nl g pk g ql b † m b † n bp bq + δ kl N k . ( 19
)
The normally ordered expectation values of the field operators bi in Eqs. ( 18) and ( 19) can be calculated from the P function [START_REF] Tsang | Resolving starlight: A quantum perspective[END_REF]

according to b † m bn = β * m β n P(β + , β -)d 2 β + d 2 β -, (20a) b † m b † n bp bq = β * m β * n β p β q P(β + , β -)d 2 β + d 2 β -. (20b)
Given that P(β + , β -) [see Eq. [START_REF] Tsang | Resolving starlight: A quantum perspective[END_REF]] is a Gaussian involving only the modes u ± (r), the only nonzero expectation values are b

† ± b± = N ± , (21a) b † ± b∓ = -γ N + N -, (21b) b † ± b † ± b± b± = 2N 2 ± , (21c) b † ± b † ∓ b± b∓ = N + N -(1 + γ 2 ), (21d) b †
± b † ± b∓ b∓ = 2γ 2 N + N -, (21e) b † ± b † ∓ b∓ b∓ = -2γ N ± N + N -. (21f)
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 = dr v * k (r)u(r ± r 0 ) as

g k± = ( f +,k ± f -,k )/ √
2(1 ± δ) and using N ± = Nκ (1 ± δ). Accordingly, inserting Eqs. [START_REF] Yang | Far-field linear optical superresolution via heterodyne detection in a higher-order local oscillator mode[END_REF] into Eqs. ( 18) and ( 19), using N ± = Nκ (1 ± δ), we get the mean photon numbers [START_REF] Tham | Beating Rayleigh's Curse by Imaging Using Phase Information[END_REF] and the covariance matrix

N k = Nκ (| f +,k | 2 + | f -,k | 2 ) -γ Nκ (| f +,k | 2 -| f -,k | 2 )
kl [d, θ, N] = 0 kl + γ 1 kl + γ 2 2 kl , ( 23
)
with 0 kl = (Nκ ) 2 (| f -,k | 2 | f -,l | 2 + | f +,k | 2 | f +,l | 2 +2 Re[ f -,k f * -,l f +,l f * +,k ]) +δ kl Nκ (| f +,k | 2 + | f -,k | 2 ), (24a)
1 kl = -δ kl Nκ (| f +,k | 2 -| f -,k | 2 ), (24b)
2 kl = (Nκ ) 2 (| f -,k | 2 | f -,l | 2 + | f +,k | 2 | f +,l | 2 -2 Re[ f -,k f * -,l f +,k f * +,l ]). ( 24c
)
The derivative vector is obtained by differentiating Eq. ( 22) with respect to the parameter d, which gives

D k [d, θ, N] = 2Nκ Re f * +,k ∂ f +,k ∂d + f * -,k ∂ f -,k ∂d -γ Re f * +,k ∂ f +,k ∂d -f * -,k ∂ f -,k ∂d . ( 25
)
In the case of two equally bright sources (γ = 0), we recover the expressions given in Ref. [START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF]. Analogously, setting γ = 0 in Eq. ( 23), we have kl [d, θ, N] = 0 kl , which coincides with the covariance matrix reported in Ref. [START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF].

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 [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF][START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF], 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

u nm (x, y) = N nm H n √ 2x w H m √ 2y w e -(x 2 +y 2 )/w 2 , ( 26
)
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

β nm (a) ≡ d 2 r u * nm (r)u 00 (r -a), (27)
which fully determine the covariance matrix [START_REF] Morizur | Programmable unitary spatial mode manipulation[END_REF] and the derivative vector [START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF] and consequently the measurement sensitivity [START_REF] Paris | Quantum estimation for quantum technology[END_REF] and coefficients [START_REF] Giovannetti | Advances in quantum metrology[END_REF]. 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 [START_REF] Morizur | Programmable unitary spatial mode manipulation[END_REF] and the derivative vector [START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF]. In particular, we consider misalignment between the demultiplexing basis and the source centroid, as well as crosstalk between the detection modes. These two imperfections affect Eqs. ( 23) and ( 25) only through modifications of the overlap functions f ±,k . Finally, we consider dark counts at the detection stage, which modify the diagonal of the covariance matrix [START_REF] Morizur | Programmable unitary spatial mode manipulation[END_REF]. 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 preestimated, possibly via direct imaging, to optimally align the demultiplexer. In the case of faint sources, it was observed [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF] that an imperfect positioning of the demultiplexer makes the Fisher information go to zero for small source separations. Methods to mitigate this effect have been proposed by adaptively switching between demultiplexing and direct imaging [START_REF] Grace | Approaching quantum-limited imaging resolution without prior knowledge of the object location[END_REF] or by optimizing the detection basis [START_REF] De Almeida | Discrimination and estimation of incoherent sources under misalignment[END_REF]. An alternative approach, which we will not consider here, consists in estimating the source's centroid simultaneously with the separation [START_REF] Řehaček | Multiparameter quantum metrology of incoherent point sources: Towards realistic superresolution[END_REF].

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 [START_REF] Paris | Quantum estimation for quantum technology[END_REF] and the coefficients ( 16) through the overlap functions f ±,k = β nm (±r 0r 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 [START_REF] Shapiro | The quantum theory of optical communications[END_REF] 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

β nm (±r 0 + r s ) = 1 √ n!m! exp - d s cos θ s ± d 2 cos θ 2 + d s sin θ s ± d 2 sin θ 2 2w 2 × d s w cos θ s ± d 2w cos θ n d s w sin θ s ± d 2w sin θ m . ( 28
)

Crosstalk

A recent experiment [START_REF] Boucher | Spatial optical mode demultiplexing as a practical tool for optimal transverse distance estimation[END_REF] 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 [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF].

To include the impact of crosstalk between the detection modes in our model, we follow [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF] 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 [START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF],

f ±,mn = pq c nm,pq β pq (±r 0 -r s ). ( 29
)
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

C( ) = exp -i K 2 -1 i=1 λ i G i ( 30
)
by sampling the uniformly distributed random real coefficients λ i (0 λ i 1 and i λ 2 i = 1), which multiply the generalized Gell-Mann matrices G i [START_REF] Bertlmann | Bloch vectors for qudits[END_REF]. The positive parameter allows us to tune the crosstalk strength. In experimentally relevant scenarios, crosstalk is generally weak [START_REF] Boucher | Spatial optical mode demultiplexing as a practical tool for optimal transverse distance estimation[END_REF][START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF]. 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 that

|c i j | 2 = 1 K (K -1) K k =l=1 |c kl | 2 1. (31)
In the following, when discussing results in the presence of crosstalk, we will consider multiple random realizations of the matrix c i j and quantify the average crosstalk probability via the ensemble average |c i j | 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 (dc). The impact of this noise source on distance estimation via spatial-mode demultiplexing has been investigated in the recent literature [START_REF] Len | Resolution limits of spatial mode demultiplexing with noisy detection[END_REF][START_REF] Lupo | Subwavelength quantum imaging with noisy detectors[END_REF][START_REF] Oh | Quantum Limits of Superresolution in a Noisy Environment[END_REF]. 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 modes Nk a classical random variable ξ k . In particular, ξ k are thermally distributed with mean value

ξ k = N dc k . ( 32
)
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 by computing not only quantum-mechanical expectation values, but also classical averages over the probability distribution of ξ k . Following this procedure, we obtain

N k = N k + N dc k (33)
and, given that N dc k is independent of d, the derivative vector D[d, θ, N] is untouched by dark counts. On the other hand, the covariance matrix becomes

kl [d, θ, N] = kl [d, θ, N] + δ kl N dc k N dc k + 1 , ( 34
)
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 a 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

I (r) = Ê (+) † (r)E (+) (r) = i, j=± u * i (r)u j (r) b † i b j = Nκ[(1 + γ )|u 0 (r + r 0 )| 2 + (1 -γ )|u 0 (r -r 0 )| 2 ]. (35)
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

A i j = {(x, y) : |x -x i | x p , |y -y i | y p }, (36)
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

I i j = Nκ (1 + γ ) i j + Nκ (1 -γ ) i j , ( 37
)
with i j = A i j |u 0 (r + r 0 )| 2 d 2 r, ( 38a
) i j = A i j |u 0 (r -r 0 )| 2 d 2 r. (38b)
The intensity coherence function is defined as

(r, r ) = Ê (+) † (r)E (+) (r) Ê (+) † (r )E (+) (r ) -I (r)I (r ), (39)
with

Ê (+) † (r) Ê (+) (r) Ê (+) † (r ) Ê (+) (r ) = i jkl=± u * i (r)u j (r)u * k (r )u l (r ) b † i b j b † k bl = i jkl=± u * i (r)u j (r)u * k (r )u l (r ) b † i b † k b j bl + δ(r -r )I (r), (40)
where in the last step we used the commutation relation [ bi , b † j ] = δ i j and the completeness relation i w * i (r)w i (r ) = δ(rr ). Using the expectation values [START_REF] Yang | Far-field linear optical superresolution via heterodyne detection in a higher-order local oscillator mode[END_REF], we can finally write the intensity coherence function as

(r, r ) = 0 (r, r ) + γ 1 (r, r ) + γ 2 2 (r, r ), ( 41
) with 0 (r, r ) = N 2 κ 2 {|u 0 (r + r 0 )| 2 |u 0 (r + r 0 )| 2 + |u 0 (r -r 0 )| 2 |u 0 (r -r 0 )| 2 +2 Re[u * 0 (r + r 0 )u 0 (r -r 0 )u * 0 (r -r 0 )u 0 (r + r 0 )]} +Nκδ(r -r )[|u 0 (r + r 0 )| 2 + |u 0 (r -r 0 )| 2 ], (42a) 1 (r, r ) = Nκδ(r -r )[|u 0 (r -r 0 )| 2 -|u 0 (r -r 0 )| 2 ], (42b) 2 (r, r ) = N 2 κ 2 {|u 0 (r + r 0 )| 2 |u 0 (r + r 0 )| 2 + |u 0 (r -r 0 )| 2 |u 0 (r -r 0 )| 2 -2 Re[u * 0 (r + r 0 )u 0 (r -r 0 )u * 0 (r -r 0 )u 0 (r + r 0 )]}. ( 42c
)
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

i jkl [d, θ, N p ] = A i j d 2 r A kl d 2 r (r, r ) = 0 i jkl + γ 1 i jkl + γ 2 2 i jkl , ( 43
)
with 0 i jkl = N 2 κ 2 ( i j kl + i j kl + 2 Re{ * i j kl }) + δ ik δ jl Nκ ( i j + i j ), ( 44a
) 1 i jkl = δ ik δ jl Nκ ( i j -i j ), ( 44b
) 2 i jkl = N 2 κ 2 ( i j kl + i j kl -2 Re{ * i j kl }), (44c)
where we have introduced

i j = A i j u * 0 (r + r 0 )u 0 (r -r 0 )d 2 r. (45)
The derivative vector is obtained from the mean intensity per pixel [START_REF] Len | Resolution limits of spatial mode demultiplexing with noisy detection[END_REF], where the dependence on the parameter d is only contained in the functions i j and i j [START_REF] Lupo | Subwavelength quantum imaging with noisy detectors[END_REF],

D i j [d, θ, N p ] = Nκ (1 + γ ) ∂ i j ∂d + Nκ (1 -γ ) ∂ i j ∂d . ( 46
)
For a Gaussian PSF, the integrals ( 38) and ( 45) can be evaluated analytically and result in

i j = 1 4 {erf[C -(x -)] -erf[C -(x + i )]}{erf[S -(y - i )] -erf[S -(y + i )]}, ( 47a
) i j = 1 4 {erf[C + (x - i )] -erf[C + (x + i )]}{erf[S + (y - i )] -erf[S + (y + i )]}, ( 47b
) i j = 1 4 e -d 2 /2w 2 (erf[x - i ] -erf[x + i ])(erf[y - i ] -erf[y + i ]). (47c)
In Eqs. (47) we have defined

C ± (x) = d cos θ √ 2w ± x, ( 48a
)
S ± (y) = d sin θ √ 2w ± y, (48b)
and x ± i = √ 2(x i ± x p )/w and y ± i = √ 2(y i ± y p )/w. Analogously, for the derivatives

∂ i j
∂d and

∂ i j ∂d we obtain ∂ i j ∂d = cos θ 2 √ 2π w (e -C 2 -(x - i ) -e -C 2 -(x + i ) ){erf[S -(y - i )] -erf[S -(y + i )]} + sin θ 2 √ 2π w (e -S 2 -(y - i ) -e -S 2 -(y + i ) ){erf[C -(x - i )] -erf[C -(x + i )]}, ( 49a
)
∂ i j ∂d = cos θ 2 √ 2π w (e -C 2 + (x - i ) -e -C 2 + (x + i ) ){erf[S + (y - i )] -erf[S + (y + i )]} + sin θ 2 √ 2π w (e -S 2 + (y - i ) -e -S 2
+ (y + i ) ){erf[C + (x - i )] -erf[C + (x + i )]}. ( 49b
)
The covariance matrix (43) and the derivative vector (46), together with Eqs. ( 47) and (49), generalize the results in the Supplemental Material of [START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF] 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

[d, θ, N P ] = I + UU T , ( 50
)
where I is a diagonal matrix whose elements are given by the mean intensities per pixels I k=i j as reported in Eq. ( 37) and U = (ν (1) , ν (2) , ν (3) ) T is N 2 p × 3 matrix whose columns are given by the vectors ν (1) k=i j = Nκ 1 + γ 2 i j , (51a)

ν (2) k=i j = Nκ 1 + γ 2 i j , ( 51b
)
ν (3) k=i j = Nκ 2(1 -γ 2 ) i j , ( 51c
)
where we used that, according to Eqs. (47), the functions i j , i j , and i j are all real. Then, using the Woodbury matrix identity [START_REF] Hager | Updating the inverse of a matrix[END_REF], we can write the inverse of the covariance matrix as

-1 i jkl [d, θ, N P ] = δ ik δ jl I i j -i jkl , ( 52
)
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 behavior 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 us to write the optimal direct-imaging sensitivity as

M[d, θ, N p ] = N p i j=1 1 I i j ∂I i j ∂d 2 - N p i jkl=1 D i j [d, θ, N p ] i jkl D kl [d, θ, N p ], (53)
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 behavior 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

F DI [d, θ] = 1 I (r) ∂I (r) ∂d 2 d 2 r, (54)
which coincides with the direct-imaging Fisher information calculated assuming Poissonian sources in the Nκ 1 regime [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF]. 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. ( 53) is relevant only for intermediate source separations. In fact, the behavior of the measurement sensitivity M[d, θ, N p ] [Eq. ( 53)] for small distances between the sources is dominated by the diagonal part of the covariance matrix (43). Accordingly, in 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

M DI [d, θ] = 2Nκ w 2 [γ 2 + 4x 2 (2 -5γ 2 + 3γ 4 )] + O(x 4 ).
(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 for d = 0 (see 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 by considering that in the γ → 1 limit the estimation of the source 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 by setting x s = y s = 0 in Eq. [START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF], which implies the symmetry

β nm (-r 0 ) = (-1) n+m β nm (r 0 ). ( 56
)
Equation ( 56) causes the mean photon number in the measurement modes N k [see Eq. ( 22)] to be independent of the photon-number imbalance γ :

N nm = 2Nκβ 2 nm (r 0 ). ( 57
)
Substituting Eq. (56) into Eq. ( 23) also simplifies the covariance matrix, which becomes

= diag(ν) + ξξ T /2 + ζζ T /2, ( 58
)
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 [START_REF] Hager | Updating the inverse of a matrix[END_REF]. Accordingly, we obtain

-1 mnm n = δ mm δ nn 2Nκβ 2 mn (r 0 ) - (-1) m+n+m +n A + -B[(-1) m+n + (-1) m +n ] + A - A + A --B 2 , ( 59
)
with

A ± = 2 1 ± γ 2 + 2Nκ Q mn=0 β 2 mn (r 0 ), (60a)
B = 2Nκ Q mn=0 (-1) m+n β 2 mn (r 0 ). (60b)
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 and 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 and 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] [Eq. ( 15)], which is the maximal sensitivity achievable by any linear combination of the mean photon numbers in the measurement modes N k [START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF], can be obtained by combining the inverse covariance matrix (59) with the derivative vector

D mn [d, θ, N] = 2κN wx (m + n -x 2 )β 2 mn (r 0 ), (61)
with x = d/2w, and is given by

M[d, θ, N] = 2Nκ w 2 [F -(2Nκ )(δ 1 + δ 2 + δ 3 )], (62)
where

F = Q mn=0 (n + m -x 2 ) 2 x 2 β 2 mn (r 0 ), (63a)
δ 1 = A + A + A --B 2 S 2 1 , ( 63b
)
δ 2 = - 2B A + A --B 2 S 1 S 2 , ( 63c
)
δ 3 = A - A + A --B 2 S 2 2 , ( 63d
)
with

S 1 = Q p,q=0 (-1) p+q p + q -x 2 x β 2 pq (r 0 ), (64a)
S 2 = Q p,q=0 p + q -x 2 x β 2 pq (r 0 ). ( 64b
)
We immediately see that, when the number of received photons is low (Nκ 1), the θ, N] is dominated by the F term. This term coincides with the Fisher information calculated by assuming Poissonian sources in the Nκ 1 regime [compare Eq. (63a) with Eq. ( 21) in the Supplemental Material of [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF]]. 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, and δ 1 remains finite and provides a correction to the measurement sensitivity, which takes the form

M inf [d, θ, N] = lim Q→∞ M[d, θ, N] = 2Nκ w 2 1 - 4(1 -γ 2 )Nκe -4x 2 x 2 [(γ 2 + 1)Nκ + 1] (1 -γ 4 )N 2 κ 2 (1 -e -4x 2 ) + 2Nκ + 1 . ( 65
)
The equal-brightness case is obtained by setting γ = 0 in Eq. ( 65), which results in

M inf [d, θ] = 2Nκ w 2 - 8N 2 κ 2 x 2 (Nκ + 1)e -4x 2 w 2 [(Nκ + 1) 2 -N 2 κ 2 e -4x 2 ]
, (66) which coincides with the quantum Fisher information for two equally bright thermal sources calculated in Refs. [START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF][START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF]. 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 and8. This dip gets deeper when increasing the mean photon number Nκ (see also [START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF][START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF]). 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 behavior 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 of 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 behavior can be observed in Fig. 7. In Fig. 7 values of Q considered. On the other hand, in Fig. 7(b) (θ = π/4), the measurement sensitivity cannot be discerned from the quantum Fisher information for Q 9. Similar behavior was reported for the Fisher information, for the regime of low photon numbers in the image plane, in the Supplemental Material of [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF]. The optimality of θ = π/4 is due to the fact that we have 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 behavior (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 imbalance 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)

w η m i j = i + j -x 2 x - 2κN A + A --B 2 × {[(-1) i+ j A + -B]S + -[(-1) i+ j B -A -]S -} (67)
with the normalization constant

η[d, θ, Q] = Q i, j=0 m 2 i j -1/2 . ( 68
)
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 measurement, 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 comment on the behavior of the measurement coefficients (67) is postponed to Sec. VII, where we will compare them with their counterpart in the presence of noise.

VII. NOISY DEMULTIPLEXING

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

M[d, θ, N] ≈ k=1 1 N k ∂N k ∂d 2 . ( 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 replace the mean photon number N k with the mean photon number plus noise: 69) is equal to the Fisher information for demultiplexing calculated by assuming Poissonian sources in the low-brightness regime (Nκ 1) [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF]. 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.

N k = N k + N dc k . Equation (
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 behavior, 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 (upper) solid 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, higher-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 [START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF], 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

β 00 (r 0 -r s ) = 1 - x 2 2 -2x 2 s -2xx s cos(θ -θ s ) + O(a 3 ), (70a)
β 10 (r 0 -r s ) = 2x s cos θ s + x cos θ + O(a 3 ), (70b)
β 01 (r 0 -r s ) = 2x s sin θ s + x sin θ + O(a 3 ), (70c)
β 11 (r 0 -r s ) = (2x s sin θ s + x sin θ ) × (2x s cos θ s + x cos θ ) + O(a 3 ), ( 70d
)
β 20 (r 0 -r s ) = (2x s cos θ s + x cos θ ) 2 √ 2 + O(a 3 ), ( 70e
)
β 02 (r 0 -r s ) = (2x s sin θ s + x sin θ ) 2 √ 2 + O(a 3 ), ( 70f
)
β nm (r 0 -r s ) = O(a 3 ), n + m 3. (70g)
Combining Eqs. ( 29) and ( 25) with the expansions (70), we get the expression for the derivative vector

w 2Nκ D[d, θ, N] nm = x[-2|c nm,00 | 2 + cos 2 θ (2|c nm,10 | 2 + √ 2c nm,20 c * nm,00 + √ 2c nm,00 c * nm,20 ) + sin 2 θ (2|c nm,01 | 2 + √ 2c nm,02 c * nm,00 + √ 2c nm,00 c * nm,02 ) + 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 )] + O(a 2 ), (71)
where we note that up to first order in a, misalignment does not affect the derivative vector.

Let us now consider a generic weak crosstalk distribution, which corresponds to setting 1 in Eq. [START_REF] Kholevo | A generalization of the Rao-Cramer inequality[END_REF]. In fact, for 1, we can write the crosstalk matrix as C( ) ≈ 1i G, where the elements of the matrix G are of order unity. Accordingly, the off-diagonal elements of the crosstalk matrix, c nm,kl with n = k and m = l, are of order , with 0 < 1. Under this assumption, we can restrict ourselves to the smallest square covariance matrix containing all terms of order . with = N dc (N dc + 1) the dark-count term, which we assume to be weak and the same for all modes. This truncated covariance matrix contains terms up to order 3 . In particular, the leading order in is 0 for 00,00 ; for 01,00 , 10,00 , 00,01 , 00,10 , 10,10 , and 01,01 ; and 3 for 10,01 and 01,10 .

The inverse of the 3 × 3 matrix in Eq. ( 72) can be obtained analytically and used to determine the measurement sensitivity, which results in

M[d, θ, N] = 2Nκ w 2 |c 00,00 | 4 2Nκ (|c 00,00 | 4 + σ 2 ) + |c 00,00 | 2 + σ + sin 4 θ |c 01,01 | 4 + cos 4 θ |c 10,10 | 4 2Nκσ 2 + σ + O( ) x 2 + O(a 3 ), (73)
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 behavior of the sensitivity for small x is dominated by dark counts (only diagonal terms of the crosstalk matrix appear).

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 i j | 2 ∼ N dc , or, in other words, when σ = N dc /2Nκ is of order 2 . In this case, the measurement sensitivity would be approximated by

M[d, θ, N] = 2Nκ w 2 sin 4 θ |c 01,01 | 4 |c 01,00 | 2 + σ + cos 4 θ |c 10,10 | 4 |c 10,00 | 2 + σ + O( -1 ) x 2 + O(a 3 ). (74)
In this case, the probabilities of scattering from the first order modes u 01 and u 10 to the fundamental mode u 00 enter in the behavior 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 Supplemental Material of [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF]). Equations ( 73) and (74) can be simplified by considering a uniform crosstalk model [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF], 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) become, respectively,

M dc [d, θ, N] ≈ 2Nκ w 2 cos(4θ ) + 3 8Nκσ 2 + 4σ + 1 2Nκ (σ 2 + 1) + σ + 1 x 2 , ( 75
) M ct [d, θ, N] ≈ 2Nκ w 2 cos(4θ ) + 3 4(|r| 2 + σ ) x 2 . ( 76
)
We denote by M dc (M ct ) the dark-count (crosstalk)-dominated sensitivity. Equation (75) remains valid also in the case when no crosstalk or misalignment is 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 causes 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

M mis [d, θ, N] ≈ 2Nκ w 2 sin 4 θ x 2 sin 2 θ + 4x 2 s sin 2 θ s + cos 4 θ x 2 cos 2 θ + 4x 2 s cos 2 θ s x 2 . ( 77
)
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:

w 2 M DI [d, θ]/2Nκ ∼ 8x 2 .
However, for typical noise levels (as reported in Fig. 9), demultiplexing has generally a more favorable coefficient.

B. Measurement coefficients

We now discuss the behavior of the optimal observable that allows us to reach the sensitivity bounds illustrated in the preceding section. For this goal, in Fig. 11 we study the dependence of the coefficients m i j [see Eq. ( 16)] of the optimal linear combination of intensity measurements Xm = i j m i j N i j in the HG modes u i j (r) as a function of the separation d between the two sources for 0 i, j Q with Q = 1 [Figs. 11(a 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 i j 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 behavior is observed in Fig. 11(a) and in its close-up in Fig. 11(d), where we show the coefficients in the 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(e)] 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 i j 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 = θ , dash-dotted lines in Figs. 11(a cients' 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(e). On the other hand, in Fig. 11(c), where the dark-count level is 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(e) we notice that in the presence of weak noise, even though the degeneracy is removed, the coefficients m i j still move in groups (green and blue for i + j = 1; orange, yellow, and pink for i + j = 2; brown and gray for i + j = 3; and purple for i + j = 4).

From Fig. 11(e) 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(e)], all coefficients m i j>0 [Eq. ( 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(e)], our optimal observable prescribes the measurement of 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 gray), 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(e)].

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 [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF]. Following our moment-based estimation approach, we have ( d ) 2 = 1/μM[d, θ, N]; accordingly, the minimal resolvable distance is defined by the condition

d min μM[d min , θ, N] = 1. ( 78
)
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κ. 1 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κ coherent sources there would not be any difference, since measuring μ copies of a coherent state |α gives the same results as measuring |μα . To get an analytical understanding of these behaviors, 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),

d min = w √ N det . ( 79
)
Accordingly, ideal demultiplexing provides a "shot-noise" scaling of the minimal resolvable distance, which is represented by 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 of the noise level, the leading term of such an expansion is independent of 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 behavior 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 Ref. [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF]).

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

d min = √ 2d s w N 1/4 det (cos 4 θ sec 2 θ s + sin 4 θ csc 2 θ s ) 1/4 , ( 80
)
which is represented as a blue dashed line in Fig. 12(a). A similar √ d s N -1/4 det scaling was reported in Ref. [START_REF] De Almeida | Discrimination and estimation of incoherent sources under misalignment[END_REF], where misalignment was studied in the low-brightness regime. We note that when N is increased for μ = 1 the red (lighter) curve in Fig. 12(a) presents an almost flat region before falling on a N -1/4 det line with a less favorable 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 off-diagonal elements of the covariance matrix [START_REF] Morizur | Programmable unitary spatial mode manipulation[END_REF] whose inverse determines the measurement sensitivity M[d, θ, N]. The analytical expression (80) cannot capture this behavior. 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

d min = w N 1/4 det |r| 2 3 + cos(4θ ) 1/4 . ( 81
)
Equation ( 81) is presented as blue dashed lines in Fig. 12(b) and coincides with the results obtained in Ref. [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF] for Poissonian sources. We can see that the blue (darker) curve (fixed Nκ = 1 and varying μ) approaches Eq. ( 81) for large N d , while the red (lighter) one (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 consider the dark counts level N dc to be constant as we increase N det , which results in

d min = √ 2w (N 2 κ 2 μ) 1/4 cos(4θ ) + 3 4N dc (N dc + 1) + 1 N dc (N dc + 1) + 2Nκ (2Nκ + 1) -1/4 . ( 82
)
For large values of Nκ, Eq. ( 82) simplifies to

d min = √ 2w √ Nκμ 1/4
cos(4θ ) + 3 N dc (N dc + 1)

-1/4 . ( 83
)
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 becoming negligible, we obtain the scaling of an ideal measurement. These different behaviors 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 behavior 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 [START_REF] Gessner | Superresolution Limits from Measurement Crosstalk[END_REF])

d min = w (N d ) 1/4 1 2 1/4 .
(84) 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 Figs. 12(c) and 12(d). However, in Figs. 12(c) and 12(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 us to resolve significantly smaller separations than direct imaging.

IX. CONCLUSION

We have 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 the 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 information even for noisy demultiplexing.

From our results it 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 the latter case, precise measures 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 Ref. [START_REF] Grace | Approaching quantum-limited imaging resolution without prior knowledge of the object location[END_REF].

Finally, let us point out that in experimental realizations (see, for example, [START_REF] Boucher | Spatial optical mode demultiplexing as a practical tool for optimal transverse distance estimation[END_REF]) all available observables can be measured at the same time, which allows us to look at arbitrary linear combinations in postprocessing. Therefore, the estimation technique presented here can be understood as the optimal postprocessing 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 the photon detectors are faster than the typical timescale of the parameter changes, the optimal measurement coefficients at each time can be selected in postprocessing.

FIG. 1 .

 1 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 nonorthogonal vacuum modes v1/2 on a beam splitter with transmissivity κ, resulting in the nonorthogonal image modes ĉ1/2 . (b) The symmetric and antisymmetric modes ŝ± have equal mean photon number N, but are classically correlated with phase-insensitive correlations ŝ † ± ŝ∓ = γ N. The modes ŝ± are mixed with the orthogonal vacuum modes v± on two beam splitters with transmissivities κ ± , resulting in the orthogonal image modes b±

FIG. 4 .

 4 FIG.[START_REF] Hell | Breaking the diffraction resolution limit by stimulated emission: Stimulated-emission-depletion fluorescence microscopy[END_REF]. Graphical illustration of a noisy demultiplexing measurement. The image of the two sources enters a demultiplexer, which performs mode sorting affected by crosstalk. Photon-counting measurements affected by electronic noise are performed at each demultiplexer's output.

1 FIG. 5 .

 15 FIG.5. Direct-imaging sensitivity M[d, θ, N p ] for two equally bright thermal sources (γ = 0) with a different mean number of received photons Nκ. We considered an alignment angle θ = π/4 and a square detector of side 6w divided into N p = 50 × 50 pixels. The black dashed line represents the analytical approximation for small source separation given by Eq. (55).

1 FIG. 6 .

 16 FIG. 6. Direct-imaging sensitivity M[d, θ, N p ] for thermal sources with fixed total mean photon number Nκ = 1.5 and different brightness imbalances: γ = 0 [red (bottom) line], γ = 0.25 [blue (middle) line], and γ = 0.5 [green (top) line]. We considered an alignment angle θ = π/4 and a square detector of length 6w divided into N p = 50 × 50 pixels. The dashed lines represent the analytical approximation for small source separation given by Eq. (55).

1 FIG. 7 .

 17 FIG. 7. Measurement sensitivity (62) for equally bright sources (γ = 0) and for different values of Q (solid lines) and alignment angle (a) θ = 0 and (b) θ = π/4. The black dashed line represents the quantum Fisher information (66). We assumed Nκ = 1.5 for both panels.

1 FIG. 8 .

 18 FIG. 8. Asymptotic measurement sensitivity (65) for fixed Nκ = 1.5 and different brightness imbalances γ .

5 FIG. 9 .

 59 FIG. 9. Measurement sensitivity M[d, θ, N] for demultiplexing into HG modes u nm (r) with n, m Q [K = (Q + 1) 2 ] with Q = 1 [red (upper) solid line] and Q = 2 [blue (middle) solid line] with different noise sources: (a) misalignment (d s /2w = 0.01 and θ s = π/4), (b) crosstalk ( |c i j | 2 = 0.0017), (c) dark counts (σ k = 0.001 ∀ k), and (d) all three noise sources at the same time. Black dashed lines starting from the zero represent short-distance approximations; in particular, we used (a) Eq. (77), (b) Eq. (76) with |t| 2 = |c i j | 2 and σ = 0, (c) Eq. (75) with σ = 0.001, and (d) Eq. (76) with |t| 2 = |c i j | 2 and σ = 0.001. Dashed lines starting from one show the results for ideal measurements (62) for Q = 1 [red (upper) line] and Q = 2 [blue (lower) line]. The green (lower) solid line describes direct-imaging results (53). For all plots, we assumed Nκ = 1.5, θ = π/4, and γ = 0.

1 FIG. 10 .

 110 FIG. 10. Demultiplexing sensitivity M[d, θ, N] in the presence of dark counts σ k = 0.001 ∀ k, misalignment d s /2w = 0.01 and θ s = π/4, and two different crosstalk levels |c i j | 2 = 0.0004 [red (upper) solid line] and |c i j | 2 = 0.004 [blue (lower) solid line]. Solid lines and bands represent means and standard deviations computed from 500 realizations of the crosstalk matrices. The green dashed line represents the ideal direct-imaging sensitivity. All curves correspond to θ = π/4, Nκ = 10, and γ = 0.

Following

  )-11(d)] and Q = 2 [Fig. 11(e)].

FIG. 11 .

 11 FIG. 11. Dependence on the source separation d of the measurement coefficients m i j for intensity measurements in the HG modes basis u i j (r) with (a)-(d) 0 i, j 1 (K = 4) and (e) 0 i, l 2 (K = 9). Different noise sources are considered: (a) and (d) misalignment [d s /2w = 0.01 and θ s = π/4 (solid) and θ s = π/3 (dotted)], (b) crosstalk ( |c i j | 2 = 0.0017), (c) dark counts (σ = 0.001), and (e) all three combined (d s /2w = 0.01, θ s = π/4, |c i j | 2 = 0.0017, and σ = 0.001). Dashed lines represents the coefficients in the noiseless case. Solid lines and bands in (b) and (d) represents the mean and one standard deviation computed over 500 crosstalk matrices. All plots correspond to Nκ = 1.5 and θ = π/4.

= 1 [

 1 FIG. 12. Scaling of the minimal resolvable distance d min with the number of detected photons N det = μ(2Nκ ) changed by varying either the number of measurements μ [blue (darker gray) line] or the brightness of the sources N [red (lighter gray) line]. We considered demultiplexing measurements into K = (Q + 1) 2 = 9 HG modes u i j (r) with 0 i, j Q = 2 and different noise sources: (a) misalignment (d s /2w = 0.01 and θ s = π/4), (b) crosstalk ( |c i j | 2 = 0.0017), (c) dark counts (N dc = 1), and (d) all three combined. Black (lower) and green (upper) straight lines represent the N -1/2 det scaling of ideal demultiplexing and the N -1/4 det scaling of ideal direct imaging, respectively. Dashed lines represent analytical approximations for large N det given by (a) Eq. (80), (b) Eq. (81), and (c) and (d) Eq. (82).

TABLE I .

 I Tabular representation of the role of the different noise sources in the perturbative expansions reported in Sec. VII A: √ indicates the dominant noise source, ( √ ) indicates a noise source that is present but not dominant (in the case where the noise source affects the equation, its order is also included in the parentheses), and ✗ denotes a noise source that is not included. Equation numbers given in square brackets correspond to the uniform crosstalk model.

	Equation	Misalignment		Crosstalk	Dark counts
	Eq. (73) [Eq. (75)]	(	√	)	(	√	,	σ )			√
	Eq. (74) [Eq. (76)] Eq. (77)	(	√ √	)			√ ✗		(	√	, σ	✗	2 )

To be more precise, μ2Nκ is the number of photons received in the image plane, which is always larger than or equal to the number of detected photons. The main cause of undetected photons is the finite number K of modes in the demultiplexing basis. However, at the level of the minimal resolvable distance, the population of the high-order modes is usually negligible. It is therefore justified to refer to N det = μ(2Nκ ) as the number of detected photons.

ACKNOWLEDGMENTS G.S. acknowledges financial support from ONERA, the French aerospace laboratory. M.G. acknowledges funding from the LabEx ENS-ICFP Grants No. ANR-10-LABX-0010 and No. ANR-10-IDEX-0001-02 PSL*. This work was partially funded by French ANR under COSMIC Project No. ANR-19-ASTR-0020-01. This work received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 899587. This work was supported by the European Union's Horizon 2020 research and innovation program under the QuantERA program through the project ApresSF.