{'fraction_non_alphanumeric': 0.05600876375163155, 'fraction_numerical': 0.016432034309155324, 'mean_word_length': 3.993598696461825, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, '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}}
A modern description of Rayleigh's criterion

Rayleigh's criterion states that it becomes essentially difficult to resolve two incoherent optical point sources separated by a distance below the width of point spread functions (PSF), namely in the subdiffraction limit. Recently, researchers have achieved superresolution for two incoherent point sources with equal strengths using a new type of measurement technique, surpassing Rayleigh's criterion. However, situations where more than two point sources needed to be resolved have not been fully investigated. Here we prove that for any incoherent sources with arbitrary strengths, a one- or two-dimensional image can be precisely resolved up to its second moment in the subdiffraction limit, i.e. the Fisher information (FI) is non-zero. But the FI with respect to higher order moments always tends to zero polynomially as the size of the image decreases, for any type of non-adaptive measurement. We call this phenomenon a modern description of Rayleigh's criterion. For PSFs under certain constraints, the optimal measurement basis estimating all moments in the subdiffraction limit for 1D weak-source imaging is constructed. Such basis also generates the optimal-scaling FI with respect to the size of the image for 2D or strong-source imaging, which achieves an overall quadratic improvement compared to direct imaging.


I. INTRODUCTION
Rayleigh's criterion, as a long-standing textbook theorem, puts a fundamental limit on the power of optical resolution [1,2]. It states that when two points are separated from each other by a distance smaller than the width of pointspread function (PSF) of the optical system, namely in the subdiffraction limit, it becomes essentially difficult to distinguish them. Recently however, researchers made a breakthrough towards surpassing Rayleigh's criterion using a new type of measurement technique, by looking at the imaging problem from the perspective of quantum metrology [3][4][5][6][7][8][9][10].
In metrology, Fisher information (FI) characterizes the ultimate precision of parameter estimation through Cramér-Rao bound [11][12][13]. When estimating the separation between two equal strength incoherent sources, it can be shown that FI tends to zero as they become closer when we use direct imaging approach (i.e. counting photons at different positions on the imaging plane). However, the quantum Fisher information (QFI, equal to the maximum FI over all possible quantum measurements) remains a constant, implying the possibility of superresolution [3]. In fact, many types of measurement have been proposed to achieve this kind of superresolution [3][4][5][6][14][15][16][17][18] and some of these approaches have already been demonstrated experimentally [19][20][21][22]. For example, when the PSF is Gaussian, it is possible to achieve the highest estimation precision by projecting the optical field onto Hermite-Gaussian modes [3,16,17].
While this new approach appears to be a promising candidate to substantially improve imaging resolution, many questions are yet to be answered: (1) What is the ultimate precision one can achieve, in a general imaging scenario, given experimentalists access to all types of measurement? (2) Which type of measurement achieves such precision? In this paper, we tackle these two questions by conducting a comprehensive Fisher information analysis in the general scenario where the incoherent source distribution on the object plane is arbitrary.
A direct way to parametrize an image is to use positions and intensities of each point as parameters to be estimated. However, it may not be the perfect choice because the position of one specific point does not tell much about the structure of the whole image. Instead, we can use moments to characterize an image which has wide applications in image analysis [23]. Since the difficulty involved in calculating QFI increases significantly as the number of sources increase, we only consider the limiting values of QFIs as the size of the image tends to zero (much smaller than the width of PSF) which we call "the subdiffraction limit".
In this paper, we choose normalized moments (normalized so that it has dimension of length) as parameters to be estimated, where detailed calculations for Gaussian PSFs and the spatial-mode demultiplexing (SPADE) measurement scheme are contained in Refs. [7,17]. We obtain the fundamental precision limit of estimating moments in the subdiffraction limit which formulated a modern description of Rayleigh's criterion, as opposed to the traditional Rayleigh's criterion restricted by direct imaging. We find that the FI with respect to (wrt) second moment remains a positive value in the subdiffraction limit, in accordance with previous work on estimating the separation between two coherent source. However, the FI wrt higher order moments always vanishes in the subdiffraction limit for nonadaptive measurements, answering question (1). This result shows the capability of going beyond direct imaging will not provide unlimited power and only push image resolution one step forward -from the first moment (the centroid of the image) to the second moment. To be specific, if we use s to represent the size of an image, the FI wrt to the K-th order moment vanishes as O(s K−2 ) (O(s K−1 )) when K is even (odd), compared with O(s 2K−2 ) using direct imaging.
Based on the FI analysis, we also obtain optimal quantum measurements (in the subdiffraction limit) corresponding to the optimal FI. It is shown in this paper that when PSF is under certain constraints, the optimal measurement basis is strongly related to its derivatives. Roughly speaking, the probability from projecting the optical field onto the K-th order derivative of the PSF provides information of the 2K-th order moment of the image. And choosing derivatives as the measurement basis successfully classifies information of different moments into different measurement outcomes, which will provide optimal FIs wrt these moments in the subdiffraction limit. In this paper, we partially answer question (2) by first providing optimal quantum measurement scheme for second moment. For higher order moments, we prove the optimality of this scheme for 1D weak-source imaging. For 2D imaging or for strong-source imaging, such scheme only provides the optimal scaling of FI wrt s, but the coefficient may be further improved.

II. SUMMARY OF RESULTS
Here we briefly summarize our results on Fisher information analysis for incoherent optical imaging.
• In Sec. III, we provide the formalism of the far-field imaging of incoherent optical sources, where we use P representation of optical states to express the Fisher information matrix (FIM).
• In Sec. IV, we consider imaging for weak incoherent sources in one-dimensional imaging. We show that the Fisher information (FI) with respect to normalized moments decreases polynomially as the size of the image decrease, by order-of-magnitude analysis. To be specific, the FI wrt second moments remains a constant as the size of the image tends to zero, and the FI wrt to higher order moments drops to zero.
• In Sec. V, we generalize the statement in Sec. IV to sources with arbitrary strength, again by order-of-magnitude analysis.
• In Sec. VI, we detail the FI analysis wrt to second moments by providing the exact value of FI and corresponding optimal measurements, as FI wrt second moment is not influenced by Rayleigh's criterion.
• In Sec. VII, we generalize all discussions about one-dimensional imaging to two-dimensional imaging, including calculating FI wrt to second moments in 2D.
• In Sec. VIII, we detail the FI analysis wrt to all moments and show how the optimal scaling of FI can be achieved wrt all moments, which is improved quadratically when compared to direct imaging. Sec. VIII also serves as a justification of the order-of-magnitude analysis in Sec. IV and Sec. V.
We also summarize the contents of each appendix here: • Appendix (A) discusses the condition under which the series expansion of probabilities and FIs. For a wellbehaved point spread function, the series expansion of probability converges uniformly and therefore the FIs can also be expanded wrt different orders of the size of the image. We also point out that our analysis can only be applied to non-adaptive measurements in order for the series expansion to be valid.
• Appendix (B) provides the first three terms in the series expansion of measurement probability for arbitrary incoherent sources, which is not explicitly given in Sec. V.
• Appendix (C) provides an alternative way to parametrize second moments in 2D imaging, as opposed to the one in Sec. VII.
• Appendix (D), Appendix (E) and Appendix (F) complement discussions in Sec. VIII in terms of optimizing FI wrt odd moments for weak incoherent sources in 1D imaging, 2D imaging, generalization to arbitrary strengths.
• Appendix (G) discusses the pre-estimation of the centroid. We provide a measurement scheme which is optimal for weak sources and at least 96.4% efficient for strong sources.
The main results in this paper are also summaried in Table I and Table II for further reference.

III. FORMALISM
Consider a one-dimensional object composed of J points on the object plane. The original field on the object plane can be expressed using P representation [24], where α = (α 1 , . . . , α J ) T is the column vector of complex field amplitudes for J optical spatial modes and where |0 is the vacuum state, a † j and a j are the canonical creation and annihilation operators at position x j . Suppose the fields are uncorrelated at different points on the object plane, then P Γ0 (α) is the independent Gaussian distribution of the J modes: where (Γ 0 ) j ≥ 0 is the average photon number emitted at the jth point and Γ 0 = ((Γ 0 ) 1 , . . . , (Γ 0 ) J ) T .
The imaging system maps the source operators a j , a † j into the image operators ψ j , ψ † j with an attenuation factor η: Here η is the transmission probability.
x is the canonical creation operator at position x and v † j is the creation operator of the auxiliary environmental modes [6]. Moreover, we assume the PSF satisfies the following assumption which will later be used in determining the optimal measurement basis. This assumption is easily satisfied, for example, when PSFs are real (real PSFs can be implemented by a two-lens system [25]), e.g. ψ PSF (x) ∝ e −x 2 /4σ 2 ; or when they are even, e.g. ψ PSF (x) ∝ e ikx 2 /2z sinc(x/σ). The field on the image plane is expressed as where and Γ := ηΓ 0 is the average photon number received from each mode. We also define the average photon number on the image plane := J j=1 Γ j (which is usually a small number) and the relative source strength γ j := Γ j / for later use. We can see that after integrating all phases in α, only those photon number diagonal terms will survive and we may write where π m is the probability of having m photons in the state and ρ m is an m-photon multimode Fock state.
Our goal is to extract information of the image from ρ. We use a set of positive operators {E(n)} satisfying n E(n) = I to represent the positive-operator valued measure (POVM) performed on ρ [12,26]. The resultant probability distributions are where E[·] represents expectation values under Gaussian distribution P Γ (α).
The Cramér-Rao bound [11] Σ F −1 (10) provides the ultimate precision limit in terms of parameter estimation, where " " means the LHS minus the RHS is positive semi-definite, Σ k is the error covariance matrix wrt parameters {M k } k≥1 and is the corresponding Fisher information matrix (FIM). M k are some functions of {x j , Γ j }, later chosen to be the normalized moments.

IV. THE ULTIMATE RESOLUTION LIMIT FOR WEAK INCOHERENT SOURCES
The probability of measurement outcome n is where ψ = (ψ 1 , . . . , ψ J ) T is the column vector of annihilation operators ψ j . In the limit where the average photon number on the image plane is small (the value of is considered known because it is easy to measure), we can expand it as a series in : where p(n) := J j=1 γ j 0|ψ j E(n)ψ † j |0 . Since the first term contains no information of the object, the FIM will be dominated by the second term, which corresponds to the situation where only one photon is detected. To study the behavior of FIM in the subdiffraction limit, we expand ψ j around its centroidX. One should be careful with the convergence radius of the series expansion though, which has a lower bound independent of the measurement E(n) (see Appendix (A)). The second term in Eq. (13) becomes example, if the object contains only two points, there are only three degrees of freedom -the positions of two points and the ratio of their strengths, then we choose the first three moments as the parameters to be measured. We use s = max i,j |x j − x i | to characterize the size of the image and conduct FI analysis in the subdiffraction limit when s → 0. Here we assume the centroid of the imageX = J j=1 γ j x j is known accurately either based on existing telescopic data or pre-estimation. In this case, we have M 1 = 0. In Appendix (G), we provide a measurement scheme for pre-estimation ofX. In 1D imaging, the scheme is optimal for weak sources and at least 96.4% efficient for strong sources. The methodology behind this scheme is not clear until Sec. VIII. Therefore we are not going to explain it in detail here.
Since any converging power series is dominated by its first non-zero term as s → 0, we have Note that when the terms of lower order than k in P (n; {x j , Γ j } does not vanish, should be bounded by a power of s with higher order than O(s −1 ). From Eq. (16), the FI for k ≥ 2 would be which indicates the following theorem: Theorem 1 (Modern Rayleigh's criterion for one-dimensional imaging): For imaging of incoherent point sources in the subdiffraction limit, the estimation variance of moment M k>2 increases inverse-polynomially as s decreases; meanwhile, the estimation variance of the second moment M 2 is bounded by a constant independent of s.
Note that we only need to bound the diagonal element of the FIM because the variance in estimation M k satisfies where the equality holds true when F is diagonal. A simple schematic illustration of above theorem is shown in Fig. 1. Further justifications are contained in Sec. V, Sec. VI and Sec. VIII. We discuss the validity of this order-of-magnitude analysis in Appendix (A). We emphasize here that the measurement is assumed to be non-adaptive in this paper and our analysis does not include the case where measurement can be adaptively modified (Appendix (A)) assuming prior knowledge on the moments to be estimated. And the adaptive measurement is excluded because it requires demanding experimental techniques. A more general analysis through direct calculation of quantum Fisher information, which can be applied to all type of measurement, can be found in Ref. [27]. Consider two point sources with equal source strengths. The distance between them equal to 2M2 can be measured precisely, therefore it shall be easy to distinguish (a1) and (a2). (b) Images (b1) and (b2) have the same M2 but different M4. Consider four point sources with equal source strengths. It is difficult to estimate the third and higher moments to distinguish the two images from each other.

V. THE ULTIMATE RESOLUTION LIMIT FOR INCOHERENT SOURCES WITH ARBITRARY STRENGTHS
In this section, we generalize the above discussion in weak source limit to sources with arbitrary strengths. In Eq. (12), we replace ψ † α with its expansion According to Wick's theorem (Isserlis' theorem) [28], any moment of Gaussian distributions can be calculated using the values of second order moments Here ] vanishes when integrating wrt phases of α.
We observe that P (n; {x j , Γ j }) can be decomposed into a power series of O(s), like in Eq. (14), where which is the probability of outcome n when all J points are located at the centroidX with thermal average 'excitation' number . Hence, we have shown that order-of-magnitude analysis is still valid. Specially, for 1, the expansion of Q k (n) depends solely on p k (n) and M k : and Eq. (21) simplifies to Eq. (13) for weak incoherent sources. We also notice that Q 2 (n)/Q 0 (n) = O( s 2 ) (see Appendix (B)), which means the subdiffraction limit (requiring Q 2 (n) Q 0 (n)) needs smaller s as increases.

VI. FI WRT SECOND MOMENT AND CORRESPONDING OPTIMAL MEASUREMENT
In Sec. IV, we have shown that there is a possibility to obtain a non-zero FI wrt M 2 . We are now going to find the exact value of the optimal FI wrt second moment and corresponding measurement basis. First, let's consider the weak-source scenario, As s → 0, P (n; {x j , Γ j }) and will be dominated by its first non-zero term, therefore according to Eq. (14), where we define a set of 0-null measurement outcomes We also note that p 0 (n) = 0 implies p 1 (n) = 0. Since E(N w 0 ) is Hermitian and non-negative, its eigenstates corresponding to non-vanishing eigenvalues must be orthogonal to ψ †X |0 and Re[ 0|ψ where the first equality is achieved when ψ (1) † X |0 is an eigenstate of E(N w 0 ) with an eigenvalue equal to one. For example, is optimal, in accordance with the optimality of the SPADE measurement scheme for Gaussian PSFs [3]. Furthermore, if ψ PSF (x) is an even function, its derivative will be odd and we can also choose E(N w 0 ) to be I−P

2
where P is the parity operator satisfying P · f (x) = f (−x), which is the so-called SLIVER measurement scheme [14]. This type of measurement does not depend on the specific expressions of the point-spread functions.
We emphasize that above discussions are only applicable in the subdiffraction limit and the optimal measurement should be modified for finite s. When we consider the special case where there are only two equal strength point sources, however, Eq. (27) remains optimal even when s is large [3].
When we use direct imaging approach, i.e. {E(n)} = {a † x a x dx}, the 0-null measurement outcomes have zero measure and lim s→0 F 22 = 0, because the probability density of the photon position x is which explains the traditional Rayleigh's criterion.
For an arbitrary source strength where the 0-null measurement outcome We also note that Q 0 (n) = 0 implies Q 1 (n) = 0 (see Appendix (B)). A detailed calculation of Eq. (12) shows that when n ∈ N 0 , and hence It has the exact same expression as Eq. (26), meaning FI wrt the second moment grows linearly as the source strength grows, following the standard quantum limit [29]. Our results agree with previous work on estimating the separation between two incoherent sources for arbitrary source strengths [5,6].
The measurement is optimal when (ψ †X ) k ψ (1) † X |0 are all eigenstates of E(N 0 ) with eigenvalues equal to one. For example, is optimal, in accordance with the optimality of fin-SPADE and pix-SLIVER [5].

VII. GENERALIZATION TO TWO-DIMENSIONAL IMAGING
Results in previous sections can be directly generalized to two-dimensional imaging. Suppose there are J point sources at positions (x j , y j ). The normalized moments are redefined as following: which fully characterizes the object configuration. Also, the size of the image s : We can expand the creation and annihilation operators around the centroid (∂ k X denotes the k-th order derivative wrtX) and calculate the probability distribution P (n; Similar to Eq. (16), we have the following order-of-magnitude analysis: and similar to Eq. (17), the diagonal elements of the FI matrix is thus extending the modern description of Rayleigh's criterion to 2D imaging: Theorem 2 (Modern Rayleigh's criterion for two-dimensional imaging): For imaging of incoherent point sources in the subdiffraction limit, the estimation variance of any moment M k with k + > 2 increases inversepolynomially as s decreases; however, the estimation variance of the second moment M 20 , M 11 and M 02 are bounded by a constant independent of s.
A simple schematic illustration above theorem is shown in Fig. 2. We are now going to find the exact values of FI wrt M 20 , M 11 and M 02 and corresponding optimal measurements. For simplicity we consider the weak source scenario. For arbitrary source strength, the FIs are still the same and the optimal measurements E(n) should be replaced with This assumption is satisfied, for example, when the PSF is real. The second order term of P (n; {x j , y j , Γ j }) is We only consider 0-null measurement outcome n ∈ N w 0 = {n| 0|E(n)|0 = 0|ψXȲ E(n)ψ †XȲ |0 = 0, ∀k} because for n / ∈ N w 0 , the zeroth order term of P (n; {x j , y j , Γ j }) would be positive and does not contribute to the FI as s → 0. Furthermore, we assume E(n) = ΠE(n)Π where Π is the projection onto the space span{ψ (10) † XȲ |0 , ψ (01) † XȲ |0 } because any component of E(n) perpendicular to it does not contribute to Q 2 (n) in the first order expansion of and consequently only affects the value of the FI in higher order terms of .
Then we can write every operator as a two-dimensional matrix in basis where ∆k 2 x := 0|ψ where Note that ρ 2 depends not only on the PSF via (∆k x , ∆k y , r) but also on the second moments. The FIM can be then be calculated using Q 2 (n) for any specific POVM {E(n)}. One way to parametrize the second moments is to define M 20 = X 2 , M 02 = Y 2 and M 11 = βXY , where X, Y is the standard deviation along x-and y-axis and β is the correlation between the distributions along x-and y-axis.
If we approximate the image by a Gaussian distribution P (x, y) = the contour lines of P (x, y) will be ellipses described by Different distributions can be distinguished from each other if we can precisely estimate the values of (X, Y, β). Another way to parametrize the second moments is to use The major and minor length of the ellipses correspond to the square root Λ 1,2 of the eigenvalues of C and the orientation θ is associated with the direction of its eigenvectors. Estimation wrt (Λ 1 , Λ 2 , θ) is discussed in Appendix (C). First let's consider the singular case where β = 1, |M 2 11 | = M 2 20 M 2 02 and ρ 2 is pure. It happens when all points sources are aligned on the same line, e.g. when there are only two point sources [4]. The optimal measurement can be determined by calculating quantum Fisher information matrix (QFIM) wrt X and Y : where the Hermitian operator L µ is the symmetric logarithmic derivative of ρ 2 wrt µ defined via ∂ µ ρ 2 = 1 2 (L µ ρ 2 +ρ 2 L µ ) [13]. The QFIM derived from Eq. (42) is The optimal measurement can be chosen to be any rank-one projection onto an orthonormal basis of the real space For 2D PSF satisifying the following more strict assumption (generalized from Eq. (5)): } where the parity operators P 1(2) satisfies P 1 (2) f (x, y) = f (−x, y) (f (x, −y)) (the same as 2D-SLIVER [4]). This type of measurement does not depend on the specific expressions of PSFs. In fact, any measurement E(n) = µν=+,− m µ,ν |e µ e ν | can be transformed into a PSF-independent version by replacing |e + e + | with When M 20 , M 11 and M 02 are indepedent parameters, β < 1, |M 2 11 | < M 2 20 M 2 02 and ρ 2 is a mixed state. The QFIM wrt (X, Y, β) is However, the QFIM is not simultaneously achievable for (X, Y, β), meaning the quantum Cramér-Rao bound Σ J −1 it not attainable. The optimal measurement for (X, Y ) is {ΠE(n 1 )Π = |e + e + | , ΠE(n 2 )Π = |e − e − |} and the optimal measurement for β is {ΠE(n 1 )Π = |e 1 e 1 | , ΠE(n 2 )Π = |e 2 e 2 |}, where and . We note that when β = 0, the optimal measurement basis for (X, Y ) and β are mutually unbiased. In fact, any three parameters characterizing ρ 2 can never be measured simultaneously using projection-valued measurement (PVM) because ρ 2 is only rank two. In practice, we can switch between different types of measurements during the measurement process. The resultant FIM will be the average of FIMs wrt each measurement.

VIII. ESTIMATION OF ALL MOMENTS IN THE SUBDIFFRACTION LIMIT
Even though the information of normalized moments with an order higher than two is jeopardized in the subdiffraction limit, it is worth figuring out the maximum FI achievable and the optimal measurement corresponding to it as one may still need to measure the high-order normalized moments even when the FI is low and the estimation cost is expensive. In this section, we will assume all moments are inpedendent variables and we only consider weak source scenario here. Generalization to sources with arbitrary strengths is contained in the Appendix (F). Ref. [17] contains a detailed discussion on the special case where the source is weak and the PSF is Gaussian, but the optimality was not proved there.
Eq. (5)  For simplicity, let's first look at the one-dimensional case with weak sources ( 1), (The analysis for arbitrary source strengths is detailed in Appendix (F).) According to Eq. (16) and Eq. (17), the lowest power of s F k can attain is max{k, } − 2 if and only if there is an E(n) such that P (n; {x j , Γ j }) is zero until the min{k, }-th order of s. However, this condition is not necessarily satisfiable for each moments.
In order for p 0 (n) = 0| ψX E(n)ψ †X |0 = 0, E(n) has to be orthogonal to ψ †X |0 (ψ †X |0 is not in the support of E(n)). Similarly, according to Eq. (30), in order for Q 2 (n) (up to the first order of ) to be zero, E(n) has to be orthogonal to ψ (1) † X |0 . We define -null measurement outcomes and we have N w ⊆ N w −1 for all , that is, -null measurement is ( − 1)-null. Then for all ≥ 0, Q 2 = O( 2 ) requires n ∈ N w . Suppose n ∈ N w −1 , then and We assume derivatives of the PSF {∂ k X ψ PSF (x −X), k ≥ 0} form a linear independent subset in L 2 (C). An orthonormal set {b (k) (x), k ≥ 0} can be generated via Gram-Schmidt process such that b ( ) (x) is orthogonal to every ∂ k X ψ PSF (x −X) with k ≤ − 1 and For example, when the PSF is Gaussian, {b (k) (x), k ≥ 0} are the Hermite-Gaussian modes; when the PSF is a sinc function, {b (k) (x), k ≥ 0} are the spherical Bessel functions of the first kind. We also notice that, according to Eq. (5), , k is odd} and they are orthogonal subspaces.
Then F 2 2 is maximized when b ( ) † X |0 is an eigenstate of E(n) with an eigenvalue equal to one. The resultant FI is For example when = 1, b (1) (x) = 1 ∆k ∂X ψ PSF (x −X) and Eq. (55) gives Eq. (26). We can show that it is possible for the FI to attain the lowest power of s (the highest precision) for even moments. To be specific, if k = 2 is even, by projecting quantum states on the image plane onto basis {b , F kk is maximized and proportional to the (2 − 2)-th power of s, as indicated in Eq. (17). Moreover, according to the Cramér-Rao bound (Eq. (10)), The estimation precision of M 2 is lower bounded by the value of (F 2 2 ) −1 . Meanwhile, the choice of measurement basis {b ( ) † X |0 , ≥ 0} not only minimizes the value of (F 2 2 ) −1 but also makes F diagonal, which means that the second equality in Eq. (56) holds true. Therefore, we conclude that {b ( ) † X |0 , ≥ 0} is an optimal basis for estimation of even moments for weak incoherent sources. Note that {b ( ) † X |0 , ≥ 0} may not be a complete basis, but any POVM is optimal as long as it contains projections onto them and other terms E(n) contained in {E(n)} is irrelevant because they do not affect the FIM in the lowest order approximation. We do not write out the irrelevant part of POVM in our discussion.
For odd moments, however, the above arguments do not apply. If we require n ∈ N w to satisfy then E(n) is not supported by ψ (k) † X |0 for all k ≤ . Consequently, we have which implies negligible contribution from weak sources. Therefore, in order to take odd moments into account, we need to relax Eq. (57) by choosing n ∈ N w −1 \N w +1 to keep the O( ) term in Q 2 +1 (n). The coefficient of (M 2 +1 ) 2 +1 can be non-zero when E(n) is supported by both ψ ( +1) † X |0 and ψ ( ) † X |0 . Meanwhile, would be non-zero at O( ) too. In the subdiffraction limit (s → 0), the denominator in Eq. (17) is dominated by Q 2 when n ∈ N w −1 \N w +1 . As shown in Appendix (D), we can maximize F 2 +1 2 +1 and in the meantime make the estimation of odd moments independent from the estimation of even moments (by letting F 2 +1 2 = F 2 2 +1 = O(s 2 )).
Then analogous to Eq. (56), F 2 +1 2 +1 fully characterizes the estimation precision of M 2 +1 . It is maximized when |0 }. Up to the lowest order of s and , In the meantime, we can also calculate F 2 2 which is exactly its optimal value as in Eq. (55). Therefore, |0 } achieves the optimal precision for both M 2 and M 2 +1 simultaneously.
To conclude, we can use the following two subsets of measurement basis: , is even} and |0 , is odd} (divided into two subsets so that they don't overlap) to estimate {M k |k = 4k or 4k + 1, k ≥ 1} and {M k |k = 4k + 2 or 4k + 3, k ≥ 0}, respectively. Each moment can be measured with the optimal precision and independently from other moments (the FIM is diagonal). However, each one of B w 1,2 can only extract half of the whole moment information: B w 1 estimates moments with orders equal to multiples of 4 plus 0 or 1; B w 2 estimates moments with orders equal to multiples of 4 plus 2 or 3. If one only needs to estimate even moments, , ≥ 0} is optimal. Now let's consider the two-dimensional case. Similar to the one-dimensional case, we define Q 2K is derived from Taylor expansion of Eq. (12). We notice that Q 2K can be written as E[ Ψ K | E(n) |Ψ K ] for some unnormalized state |Ψ K . Hence Q 2K is always non-negative and is equal to zero (up to the first order of ) if and only if n ∈ N w K . Based on the method of induction, we conclude that Q 2K = O( 2 ) if and only if n ∈ N w K . Therefore, by choosing proper measurement basis for n ∈ N w K−1 \N w K , one can estimate M L 2K−L with an FI up to O(s 2K−2 ) for all 0 ≤ L ≤ 2K. In general, the optimal measurement basis depends on the value of each moments.
For M L 2K+1−L , consider the O(s 2K+1 ) term in P (n; {x j , Γ j }): Clearly, if n ∈ N w K , Q 2K+1 (n; {M k , k + ≤ 2K + 1}) = 0. Therefore we should focus on measurement E(n) such that n ∈ N w K−1 \N w K+1 . Similar to 1D imaging, the optimal scaling we can obtained for Suppose ψ PSF (x, y) is separable and ψ PSF (x, y) = ψ 1,PSF (x)ψ 2,PSF (y). One can generate two orthonormal sets {b Similar to 1D imaging, one can project ρ onto these basis to extract information of moments (see Table III) and achieve the optimal scaling of s (but not necessarily the optimal coefficients). As before, one type of measurement can only estimate part of all the moments (1/4 to be specific) and by combining different types of measurements one Types of measurement Measurement basis L Moments estimated can grasp information of all moments. In practice, combining {B w i } 6 i=1 will be enough to extract all the information of moments from ρ. For further justifications and calculations of FIs see Appendix (E).
In the case of sources with arbitrary strenghs, we show in Appendix (F) that the same scaling wrt s is still achievable by replacing every E(n) with k=0 give the same FIs as in the weak source scenario. However, the coefficient may be further improved using other basis, due to the fact that information of high order moments can be obtained by detecting several low order derivative operators simultaneously, which is neglectable when the source is weak. In contrast to estimation of the second moment, when estimating higher order moments, the optimal precision increases superlinearly (instead of linearly) as the source strength grows in the subdiffraction limit.

IX. CONCLUSION
We have performed a comprehensive Fisher information analysis on general imaging scenarios in the subdiffraction limit, where the improvement of image resolution is considered difficult due to the positive width of point spread functions. We conclude that, for any incoherence sources, a 1D or 2D image can be precisely estimated up to its second moment and the higher order moments are difficult to estimate in the sense that the error increase inversepolynomially as the size of image decreases. The imaging situation considered in the paper is quite general where both the number of point sources and source strengths can be arbitrary. The problem of pre-estimation of centroid is also worked out.
For real point spread functions, we put forward a measurement scheme which provides the optimal Fisher information in the subdiffraction limit. The measurement basis is constructed based on the derivates of the point spread function, which are closely related to moments of an image. The optimal measurement scheme for second moment is discussed in detail. For higher order moments, compared with direct imaging approach, our measurement scheme guarantees at least a quadratic improvement of Fisher information in terms of the scaling wrt the size of the image. The coefficient of Fisher information is also optimal for weak sources, but can be further improved for strong sources. It is not clear, though, which measurement basis is optimal in terms of the exact value of Fisher information for strong sources.
The generality of our results has a cost though -the Fisher information is only calculated in the limiting case where the size of the image tends to zero. Direct calculations for a positive size can be difficult and it remains unsolved how to identify the optimal measurement scheme when the size is not too small (in the subdiffraction limit) and also not too large (the point spread function can be viewed as a delta function). Our results, however, is an important theoretical result towards the ultimate resolution limit for incoherent optical imaging.
Note added.-Recently, Ref. [27] appeared, which directly calculates the quantum Fisher information wrt moments for subdiffraction incoherent optical imaging. This approach can be applied to all types of measurements, without the non-adaptivity restriction in our analysis. Our results on arbitrary source strength, generalization to two-dimensional imaging and optimal scaling achieving measurement, however, are not covered in Ref. [27]. In this section, we justify the series expansion of the probability P (n; {x j , Γ j }) around its centroid. For simplicity, we only consider weak sources in 1D imaging. For single-photon measurement, We want to know when the following series will converge uniformly to P (n; {x j , Γ j }): where Let the radius of convergence R = (lim sup k→∞ |P k (n)| 1/k ) −1 , then Eq. (A2) converges uniformly as long as s < R [30].
Next we show that the radius of convergence R ≥ R 0 where R 0 independent of E(n).
where ψ ( ) PSF represents the -th order derivative of f and ψ ( ) Then Therefore when s < R 0 ≤ R, Eq. (A2) uniformly converges. For example, for a Gaussian PSF we see that R 0 ≥ σ from Eq. (A4). Therefore in the subdiffraction limit (s σ), the series expansion is always valid. However, things may break down when s > R 0 which may happen if ψ PSF (x) has complex sub-wavelength structure.
When s < R 0 , the diagonal element of the Fisher information matrix is where we assume where N + = {n : P k (n) ≥ 0}, N − = {n : P k (n) < 0} and E(N ± ) = n∈N± E(n).
The order-of-magnitude analysis above is valid only when is reasonably small when s is small. We argue that this is usually true for non-adaptive measurements: • Consider first the case when P k (n)(M k ) k k >k P k (n)(M k ) k , then clearly may be large. However, the contribution to F kk is negligible.
• When P k (n)(M k ) k ≈ k >k P k (n)(M k ) k and (when P k (n) = 0 for all k ≤ k) the first and second terms in cancel each other out, up to the lowest order of s. Above analysis could become invalid. However, it requires a special design of measurement based on prior knowledge of the moments. We exclude this type of adaptive measurement in our discussion.
Appendix B: First three terms in the series expansion of measurement probability for arbitrary incoherent sources We aim to expand P (n, {x j , Γ j }) in series of O(s k ) where s is the size of the image. To do this we replace ψ † α with ∞ k=0 First of all, we calculate the value of denominator which gives Therefore, The zeroth order term is where we use E[e −|A (0) | 2 |A (0) | 2k ] = k! k (1+ ) k+1 . The first order term is For n ∈ N 0 , Q 0 (n) = Q 1 (n) = 0 and only the first term in Q 2 (n) survives, which gives Eq. (30). The second term Re[ 0|(ψX ) k ψ (2) X E(n)(ψ †X ) k+1 |0 ] in Eq. (B5) vanishes for n ∈ N 0 because E(n) is Hermitian and non-negative and its eigenstates corresponding to non-vanishing eigenvalues must be orthogonal to (ψ †X ) k |0 for all k.
for estimation of θ. We note here that Eq. (C2) and Eq. (C3) are mutually unbiased.
The measurement basis |Φ ± also leads to F 2 +1 2 = F 2 2 +1 = O(s 2 ) which means F is effectively diagonal and the second equality in the above inequality holds, because up to the lowest order of s we have Appendix E: Measurement basis and corresponding FIs for weak incoherent sources in 2D imaging According to Eq. (62), by choosing measurement basis xy , one can achieve the optimal scaling of s (but not necessarily the optimal coefficients) for FIs wrt M 2L 2K−2L for all K and L ≤ K: By choosing measurement basis , one can achieve the optimal scaling of s for FIs wrt M 2L+1 2K−(2L+1) for all K is even (or odd) and L < K: 2K as a parameter can be estimated simultaneously with precision O(s −2K ) and independently of M 2L+1 2K−(2L+1) . Here we have used the property that b (k ) (x) is orthogonal to ∂ k X ∂ Ȳ ψ PSF (x −X, y −Ȳ ) as long as k and k (or and ) do not have the same parity (i.e. are not both even or odd). To conclude, B w 1,2 cover the estimation of moments whose orders on x-and y-axis are both even or both odd. The optimal FI scaling is O(s 2K−2 ) in this case, where 2K is the sum of orders on x-and y-axis.
For moments who have different parities on x-and y-axis, we can use basis Based on Eq. (63), we can calculate the following FIs (up to the lowest order of s): and M 2L 2K−2L (M 2K−2L 2L ) can be estimated simultaneously and independently with M 2L+1 2K−2L (M 2K−2L 2L+1 ): which are exactly their optimal values as in Eq. (E2).
Appendix F: Estimation of higher order moments with arbitrary source strengths Here we only consider 1D imaging, the discussion can be easily generalized to 2D imaging. As already shown in Sec. VI. Only 0-null measurement n ∈ N 0 = {n|Q 0 (n) = 0|(ψX ) k E(n)(ψ †X ) k |0 = 0, ∀k} contributes to the FI wrt M 2 . Using the method of induction, we define -null measurement where When n ∈ N −1 , M 2 first apprears in Q 2 (n, {M k , k ≤ 2 }). Up to the lowest order of s, When n ∈ N −1 , Q 2 (n; {M k , k ≤ 2 }) has the following form: where the remainder term Q R 2 (n; {M k , k ≤ 2 − 1}) contains only moments with orders lower than 2 . We note that Q 2 (n, {M k , k ≤ 2 }) contains only terms like (summing over k ≥ max{K + 0 , K − 0 }) ± m . From Wick's theorem and Eq. (20), it is clear that the only term dependent on M 2 corresponds to K + 0 = K − 0 = 1, m ± = 1 and m± = . When m ± = 1 and m± = , we have proving Eq. (F4). Therefore, by choosing the modified measurement Q R 2 (n; {M k , k ≤ 2 − 1}) = 0 and the same FI (Eq. (55)) wrt M 2 is recovered using the modified measurement B 0 . We can also modify other basis analogously by allowing multi-photon detection of ψ †X and it will provide the same HereM 1 andM 2 is redefined using X R as the centroid. According to Appendix (D), the optimal measurement in terms of estimatingM 1 =X − X R can be an arbitrary projection onto two orthonormal basis in the real span of is optimal. The corresponding FI is which is the same as Eq. (26). Therefore, if we want to estimate both the second moment M 2 and the centroidX, a straightforward method is to first use half of the whole resource to locateX such that δX s and then use the rest half to estimate M 2 as described in Sec. VI. The effective FIM would be half of the optimal ones Eq. (G3) and Eq. (26), When we want to estimate even higher order moments, the resource required to locateX is neglectable. Now we show Eq. (G4) is the optimal precision we can get (in the subdiffraction limit) and the QFIM Eq. (G5) is not achievable. For any POVM {E(n)}, the only case when P (n, {x j , Γ j }) does not lead to a zero-FI wrt M 2 is when there is an E(n) such that and we use the relationM 2 2 = M 2 2 +M 2 1 . Note that A 2 (n) ≈ 0| ψ (1) and A 1 (n) 2 ≤ 4A 2 (n)A 0 (n), we have When P (n, {x j , Γ j }) is dominated by Q 0 (n), we also have (G13) Therefore, any achievable FIM must satisfies and Clearly the last three equalities are simultaneously satisfied when FIM is Eq. (G4), implying the optimality of our measurement scheme. The situation becomes a bit more complicated for arbitrary source strengths. First, we expand P (n; {x j , Γ j }) around X R up to O(s) Since the quantum state is photon number diagonal, the optimal measurement estimatingM 1 must also be photon number diagonal [13], that is, {E(n)} should contains {E(n k ), k ≥ 1} where E(n k ) = Π k E(n k )Π k and Π k is projection onto k-photon subspace. In this case, we shall write where the equality holds when {E(n k )} is an arbitrary projection onto two orthonormal basis in the real span of {(ψ † X R ) k−1 ψ (1) † X R |0 , (ψ † X R ) k |0 } as long as Q 0 (n) Q 1 (n) Q 2 (n) is satisfied. For example, E(n k,± ) = 1 2 is optimal. Therefore, if we want to estimate both the second moment M 2 and the centroidX =M 1 + X R , a straightforward method is to first use half of the whole resource to locateX such that δX ≤ s and then use the rest half to estimate M 2 as described in Sec. VI. Note that to achieve the optimal precision wrt M 1 , one has to count the number of detected photons by projecting the quantum state ontõ The resource required to locateX when we want to estimate even higher order moments is still neglectable as in the weak source scenario. Now we consider the possiblity of further improving Eq. (G20), here we show that above scheme is at least 96.4% efficient. According to Appendix (B), we have, up to O(s 2 ), P (n, {x j , Γ j }) = A 0 (n) + A 1 (n)M 1 + A 2 (n)M 2 2 + A 3 (n)M 2 1 .
For different measurement outcome n, there are only two situations: • P (n, {x j , Γ j }) = O(s 2 ), then Other terms can be ignored in the subdiffraction limit.
• P (n, {x j , Γ j }) = O(1), then and A 2 (n) and A 3 (n) can be ignored in the subdiffraction limit. For simplicity we can assume Eq. (G25) and Eq. (G26) are also true.
One important property derived from this relation is that n A 2 (n) = n A 3 (n) = ∆k 2 .
We can easily generalize above measurement scheme to 2D imaging when the PSF is separable. ψ with optimal FIs equal to F 10 10 = 4 ∆k 2 x , F 01 01 = 4 ∆k 2 y .
We won't discuss simultaneous estimation of the centroidM 10 ,M 01 and the second moments M 20 , M 11 , M 02 here.