{'fraction_non_alphanumeric': 0.04453334610962054, 'fraction_numerical': 0.017200076657723266, 'mean_word_length': 4.435546875, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gaussian quantum metrology for mode-encoded parameters: general theory and imaging applications

Quantum optical metrology aims to identify ultimate sensitivity bounds for the estimation of parameters econded into quantum states of the electromagnetic field. In many practical applications, including imaging, microscopy, and remote sensing, the parameter of interest is not only encoded in the quantum state of the field, but also in its spatio-temporal distribution, i.e. in its mode structure. In this mode-encoded parameter estimation setting, we derive an analytical expression for the quantum Fisher information valid for arbitrary multimode Gaussian fields. To illustrate the power of our approach, we apply our results to the estimation of the transverse separation between two point sources observed with a diffraction-limited imaging system. In this context, we determine the ultimate resolution provided by mutual coherence and squeezing of the sources.


I. INTRODUCTION
A fundamental task in quantum metrology is to identify the ultimate sensitivity limit in the estimation of a parameter encoded into a quantum state. Even under ideal conditions, when all technical noise sources are removed, quantum noise poses unavoidable limitations to such estimation. In spite of that, quantum parameter estimation theory provides the tools to reduce quantum noise by optimizing the output measurements. This optimization leads to the quantum Cramér Rao lower bound, which states that the minimal uncertainty ∆ϑ of the estimator of a parameter ϑ is inversely proportional to the quantum Fisher information of the quantum stateρ ϑ where the parameter is encoded [1][2][3][4][5][6]. This bound can be further optimized by finding quantum states that, for a given parameter, maximize the value of the quantum Fisher information.
Electromagnetic fields play a privileged role as metrological probes in a variety of branches of science and technology, ranging from imaging and microscopy [7,8], to remote sensing with lidars and radars [9][10][11], to gravitational wave detection [12,13]. In several of these applications, the parameter of interest does not only modify the quantum state of the probe light, but also its spatio-temporal distribution. Such a spatiotemporal distribution is conveniently described in terms of modes. i.e. normalized solutions of Maxwell's equations in vacuum [14]. For example, spatial modes of light describe the different components of an image, while the properties of an optical pulse are encoded into frequency-time modes.
Previous works in this contests of mode-encoded parameter estimation focused on specific problems. For example, it was considered the case when the total light's intensity is not affected by the parameter, but its distribution among different modes is; this happens in the estimation of a small lateral beam displacement [15,16], or in the estimation of spectral parameters of a frequency comb [17]. A general theory for this fixed-intensity scenario was recently presented [18]. A different problem, that lately attracted a lot of attention, is the estimation of the separation between two point sources analysed through a diffraction limited imaging system [19][20][21][22].
For this problem, the parameter of interest is encoded in the shape of two transverse modes in the image plane [19]. Moreover, the imaging system induces separation-dependent photon losses that affect the quantum state of the light in the image plane. Pushed by the need to go beyond these case studies, in this work, we study mode-encoded parameter estimation with arbitrary multimode Gaussian states, i.e. photonics quantum states fully defined by the first two moments of their quadratures [23][24][25].
Gaussian states play a central role in quantum optics: they describe important classical states such as coherent states, representing lasers operating above threshold, and thermal states, describing fully incoherent light. Furthermore, non-classical Gaussian states can be produced deterministically in nonlinear optical processes. Among the latter states, there are squeezed states, whose reduced quantum noise has been proposed as a useful resource since the early days of quantum parameter estimation [26], and is now a key ingredient of several quantum-enhanced metrological schemes [12,13,16,27]. While previous studies of the quantum Fisher information for Gaussian states exist [28][29][30][31], they focused on the estimation of parameters defining the first two moments of the quadratures, e.g. mean field, phase, and squeezing.
The aim of this work is to overcome these limitations, to study the estimation of parameters encoded in the spatiotemporal profile of the electromagnetic field, and therefore to broaden the applicability of Gaussian quantum metrology to new fields of technology such as imaging, microscopy, and temporal (or spectral) beam profiling. As an example of such an application, we reconsider the estimation of the transverse distance between two point sources with our approach. We confirm previously known results [20,21], and extends them investigating the impact of different Gaussian resources such as mutual coherence and squeezing of the sources. For the former, we study the difference between coherence due to thermal correlations and that induced by displacement. For the latter, we differentiate between those states that are entangled in the image plane, and those that are not.
Our paper is organized as follows: First, in Sec. II, we recall arXiv:2202.10355v1 [quant-ph] 21 Feb 2022 some basic facts about Gaussian states and quantum parameter estimation. We then derive the analytical expression of the quantum Fisher information for mode-encoded parameter estimation with Gaussian states, in Sec. III. Section IV contains the application of our results to the estimation of the separation between two point sources; and Sec. V concludes our work.

A. Gaussian states
An N−mode continuous variable (CV) [32,33] quantum system can be described by choosing a mode basis, i.e. a set {u k (r, t)} N k=1 of solutions of Maxwell's equations, orthonormal with respect to the inner product and associating with each mode u k (r, t) a pair of quadra- k andâ k are standard creation and annihilation operators. If we group all quadratures in the 2N−dimensional vector x = (q 1 ,p 1 , . . . ,q n ,p n ) , from the canonical commutation relations [â k ,â † l ] = δ kl for annihilation operators, we obtain with the symplectic form Ω = N k=1 ω k , and ω k = iσ y , where we have introduced the notation σ i=x,y,z for the standard 2 × 2 Pauli matrices. Our preferred phase space representation of an N−mode CV quantum state with density matrixρ, and characteristic function χ(y) = tr exp −iy Ωx ρ , is the Wigner function In this work, we restrict ourselves to the study of Gaussian states. N−mode Gaussian states are CV states with Gaussian Wigner function [23][24][25] which are completely determined by the displacement vector x = x , and the covariance matrix where {·, ·} denotes the anticommutator, and · the expectation value · = tr ·ρ . For every physical state, the covariance matrix satisfies the uncertainty inequality [34] σ + iΩ ≥ 0.
A useful property of the covariance matrix is that, according to Williamson theorem [35], it can be decomposed as where S is a symplectic matrix, i.e. S ΩS = Ω, we introduced the notation 1 n for the n−dimensional identity matrix, and ν is a diagonal matrix whose element are known as symplectic eigenvalues. The uncertainty inequality (6) implies that the symplectic eigenvalues must be larger than unity, i.e. ν k ≥ 1.

B. Modes and states transformations
Let us recall that the space of solutions of Maxwell's equation is a Hilbert space [14,36]. Accordingly, different mode bases are connected via unitary transformations with U kl = (v l |u k ). Under mode basis changes, creation operators follow the same transformation rules as the modes [14], which implies that the quadrature vector transforms according tox A mode basis change is a particular case of a Gaussian channel: a completely positive, trace preserving map transforming Gaussian states into Gaussian states. Such channels are completely determined by their transformation rules for the displacement vector and the covariance matrix [37] x = Tx +z, wherez is a real 2N−dimensional vector, while T and N are 2N × 2N real matrices with N = N and satisfying the positivity condition N + iT ΩT ≥ iΩ. From Eq. (9) is easy to see that a mode basis change is a Gaussian channel (see Eqs. (10)) with N = 0 and T = O.

C. Quantum estimation theory
Let us now assume that we want to estimate a parameter ϑ encoded in a quantum stateρ ϑ from M independent measurements of a given positive operator-valued measure (POVM) defined by the operatorsK µ . Using classical post-processing techniques, from the measurements results, we can extract an unbiased estimatorθ of the parameter as well as its standard deviation ∆ϑ. The latter is bounded according to the Cramér-Rao inequality [1][2][3][4][5] ∆ϑ ≥ 1 with the Fisher information F ϑ,K µ defined by where p(µ|ϑ) = Tr K µρϑ is the conditional probability of obtaining the result µ for a given value of ϑ, and we introduced the compact notation ∂ ϑ · = ∂ · /∂ϑ for the derivative. The Fisher information optimized over all possible POVMs is the quantum Fisher information (QFI), and establishes the ultimate metrological sensitivity [38]. In general, the QFI can be computed as where the symmetric logarithmic derivative (SLD) is implicitly defined by the equation [1][2][3][4][5] 2∂ ϑρϑ =L ϑρϑ +ρ ϑLϑ .
When ρ ϑ is an N−mode Gaussian state defined by the displacement vectorx and the covariance matrix σ, the SLD is quadratic in the quadratures [29,30] with where S and ν are the symplectic matrix and the symplectic eigenvalues obtained from the Williamson decomposition of σ as introduced in Eq. (7), while a (l) jk a set 2n × 2n matrices that are zero everywhere except in the jk block where they are given by ω/ √ 2, σ z / √ 2, 1 2 / √ 2 and σ x / √ 2 for l = 0, 1, 2 and 3, respectively. Substituting Eqs. (16) and (17) into Eq. (14) and using the properties of the characteristic function [33], we can write the QFI for a Gaussian state as with where F σ and Fx are the contribution to the QFI coming from variations of the covariance matrix σ and the displacement vectorx, respectively.

III. MODE-ENCODED PARAMETER ESTIMATION
Let us consider the estimation of a parameter ϑ encoded into a Gaussian stateρ ϑ expressed into an n−dimensional mode basis u k [ϑ](r, t), with n the smallest number of modes necessary to describe the system. We will refer to the Hilbert space spanned by these modes as H n = span ({u k [ϑ](r, t)}). Since every basis of the mode Hilbert space H n would provide a description of the quantum stateρ ϑ in the terms of the smallest number n of modes, the choice of the mode basis u k [ϑ](r, t) is not unique. Despite this freedom of choice, in the most general parameter estimation scenarios, every mode basis u k [ϑ](r, t) will be parameter-dependent. The latter fact implies that the Gaussian stateρ ϑ depends on ϑ not only explicitly through the displacement vectorx ϑ and the covariance σ ϑ , but also implicitly through the mode functions u k [ϑ](r, t). Our goal is this section is to calculate the QFI (18) in the scenario described above.

A. Separation of mode and state parameter dependence
Our first step is to make the parameter dependence coming from the modes u k [ϑ](r, t) explicit in the covariance matrix and displacement vector of the quantum stateρ ϑ . To this goal, we express them into a parameter independent basis v k (r, t): using Eq. (9), we get where we introduced the subscripts I and ϑ to identify quantities in the parameter-independent and parameter-dependent bases, respectively. Naturally, the choice of the parameter independent basis v k (r, t) is not unique. However, since this basis doesn't contain any information on the parameter, its choice does not affect the final expression for the QFI, as will become clear at the end of our calculation. Given that n is the smallest number of modes necessary to represent the stateρ ϑ , the parameter independent basis v k (r, t) must have dimension N ≥ n. To take into account this change in dimension, we complement the state in the parameter dependent mode basis with N − n vacuum modes, so that we can write the covariance matrix σ ϑ in block diagonal form as and the displacement vector asx ϑ = (x n , 0, · · · , 0) . To isolate the action of O on the n initially populated modes, it is convenient to rewrite it as a 1 × 2 block matrix with O n and O N−n a 2N × 2n and a 2N × 2(N − n) matrices, respectively. Some useful properties of these matrices and their derivatives are reported in App. A. Substituting Eqs. (22) and (21) into Eq. (20a), and using the properties of the matrices O n and O N−n (See Eq. (A4b) in App. A), we can rewrite the covariance matrix in the mode independent basis as Analogously, using Eq. (22) into Eq. (20b), we can rewrite the displacement vector asx Equations (23) and (24) provide a description of the Gaussian stateρ ϑ where the parameter dependence is fully expressed in the covariance matrix σ I and the displacement vector x I . In particular, the transformation properties of the n initially populated modes appear now explicitly through the matrix O n . In the following, we are going to use these expressions to compute the two terms in Eq. (18).

B. Covariance matrix contribution to the quantum Fisher information
We start with the calculation of F σ (see Eq. (19a)) which describes the contribution to the sensitivity due to variations of the covariance matrix. Let us start by taking the derivative of the covariance matrix σ I in the parameter independent basis with respect to the parameter To compute the quadratic term of the SLD L (2) ϑ (see Eq. (17c)), we need the Williamson decomposition σ I = S I ν I S T I of the covariance matrix σ I . Using Eqs. (21) and (22), we can connect it to the Williamson decomposition V n = S n νS T n of the covariance of the n initially populated modes in the parameter dependent basis u n [ϑ](r, t), and obtain Accordingly, using the properties of the matrix O (see App. A for details), we can write with Here, D n is a 2n × 2n matrix and D ∂ is a 2n × 2m matrix, constructed using, respectively, the coefficients c kl [ϑ] and c kl [ϑ] of the expansion where the modes u l [ϑ](r, t) form an m(≤ n)−dimensional basis of the mode Hilbert space Accordingly, the diagonal block B n contains a mode contribution (first two terms in Eq. (28a)) due to the portion of the derivatives ∂ ϑ u k [ϑ](r, t) within the space of the initially populated modes H n , and a contribution given by the explicit the dependence of the covariance matrix V n on the parameter. On the other hand, the off diagonal blocks B ∂ and B ∂ only contain the mode contribution due to the leakage of the derivatives whereÃ (l) jk are m × n blocks of the matrices A (l) i j . Finally, using Eq. (19a) and Eq. (26b), we can write the covariance matrix contribution to the QFI as The sum in the first term in Eq. (31) only runs over the n initially populated modes. Accordingly, it describes the contribution to the QFI given by variations of the state within the n initially populated modes u n [ϑ](r, t). On the other hand, the second term in Eq. (31) contains a sum over the n initially populated modes u n [ϑ](r, t) and another over their m orthonormalized derivatives u n [ϑ](r, t). Therefore, it takes into account the contribution to the QFI due to the coupling between the initially populated modes and their derivatives induced by parameter variations. Finally, let us note that Eq. (31) is completely determined by the covariance matrix V n of the statê ρ ϑ in the n initially populated modes u k [ϑ](r, t), and by the shape of the modes themselves, but, as anticipated, it does not depend on the choice of the auxiliary parameter-independent basis v k (r, t).

C. Displacement vector contribution to the quantum Fisher information
We now move to compute Fx, as given by Eq. (19b), which takes into account the contribution to the QFI coming from variations of the displacement vectorx I . To compute this term, we need the derivative of Eq. (24) and the inverse of the covariance matrix σ I that, thanks to Eq. (22), we can write as Finally, combining Eqs. (32) and (33), and using the properties of the matrices O n and O N−n (see App. A), weobtain Similarly, to what we observed for F σ , Fx only depends on the displacement vectorx n in the n initially populated modes u k [ϑ](r, t) and their shapes. Moreover, we note that the first term in in Eq. (34) only depends on variations of the displacement vectorx n , while the last term only depends on changes of the shapes of the n initially populated modes u k [ϑ](r, t). On the other hand, in the two middle terms both (∂ ϑxn [ϑ]) and D n appear, accordingly they combine modes variations with changes in the displacement vector.

IV. A CASE STUDY: SOURCE SEPARATION ESTIMATION
The formalism developed in Sec. III is fully general and can be used to compute the QFI for the estimation of a parameter encoded in n modes of an arbitrary Gaussian state. To illustrate how to use this framework in practice, we now apply it to the estimation of the transverse separation between two point sources observed through a diffraction-limited imaging system.

A. Modes and states in the image plane
We start by describing how two point sources transversally separated by a distance d in the object plane propagate through a diffraction-limited imaging system characterized by a point spread function u 0 (r), which (up to a global phase factor) we assume to be real [21]. Moreover, without loss of generality, we can assume the sources to be aligned along the x−axis, so that their positions are given by ±r 0 = (±d/2, 0).
The two point sources are associated with quadrature operatorsx (s) . Propagation through the imaging system transforms the source operators into the image operatorsx 1 = (q 1 ,p 1 ) andx 2 = (q 2 ,p 2 ) associated with the modes u 0 (r ± r 0 ), according tô where κ is the attenuation factor of the imaging system, andx (v) 1 andx (v) 2 are quadrature operators of auxiliary nonorthogonal vacuum modes that describe the effect of diffraction. We remark that this model is well defined only for κ ≤ 0.5, which is in accordance with the fact that we are considering only light propagating forward from the point sources to the far field of the imaging system [21,39].
To obviate to the diffraction-induced loss of orthogonality discussed above, it is convenient to introduce the images of the symmetric and antisymmetric superposition of the source modes From Eq. (35), we can then obtain the displacement vectorx and the covariance matrix V of the quadraturesx = (x + ,x − ) from those of the source modesx 0 = (x (s) 1 ,x (s) 2 ) by using the Gaussian channel (See Eq. (10)) where describes the mode transformation to the symmetric and antisymmetric superposition of the source modes (which we labelled with the subscript BS since it is equivalent to a 50 − 50 beam splitter), while describe two independent lossy channels with attenuation factors κ ± = κ(1 ± δ) depending on the parameter d through the mode overlap δ. From Eqs. (38), we can see that in the object plane the parameter d is fully encoded in the source modes while the covariance matrix V 0 and the displacement vector x 0 only contain information on the sources' photon statistics. On the other hand, because of the separation-dependent attenuation factors κ ± , for the quantum state in the image plane the dependence on the parameter we want to estimate is both explicit inx and V and implicit in the mode functions in Eq. (36).
To compute the QFI using the general theory we presented in Sec. III, we first identify the modes u ± [d](r) as the initially populated modes forming a parameter dependent basis of the Hilbert space H n=2 , which in this particular case, we are going to relabel as H ± [40]. Then, in order to compute the matrices D ± and D ∂ , we need to construct the basis {u k [d](r)} of the Hilbert space H ∂ = span({∂ d u ± [d](r)}) \ H ± (see Eq. (29)). It is not difficult to verify (See App. B) that this basis is given by the derivative modes ∂ d u ± [d](r) themselves divided by the parameter-dependent normalization factors with δ = ∂ d δ, and Equations (37), (41), and (42) depend on the particular shape of the PSF u 0 (r) of the imaging system. In the following, we will derive expressions for the QFI valid for arbitrary PSFs. However, to produce figures that illustrate the behaviour of the QFI for different quantum states of the sources, we will specify to a Gaussian PSF u 0 (r) = 2/(πw 2 ) exp −|r| 2 /w 2 . The explicit expressions for δ, δ , (∆k) 2 and β for a Gaussian PSF are reported in App. C. For the coefficients in Eq. (29), the fact that the modes with k, l = ±. Accordingly, for the separation estimation problem, the matrix D ± = 0 , while D ∂ is a 4 × 4 diagonal matrix D ∂ = η + 1 2 ⊕ η − 1 2 . This greatly simplifies the calculation of the QFI, in fact in this case both B ± and B ∂ are 4 × 4 square matrices [41], and in particular, the matrix In a similar way, the displacement contribution to the QFI simplifies to We now have all the ingredients to compute the QFI for the estimation of the separation between two point sources in arbitrary Gaussian states. In the following, we will present the results for some specific examples.

B. Correlated thermal states
As a first example, we consider two sources in a correlated thermal state. These sources are characterized, in the object plane, by a two-mode zero mean Gaussian state with covariance matrix where N 0 is the mean photon number emitted by each source, 0 ≤ γ ≤ 1 quantifies the correlation strength, and R(φ) is a rotation matrix, with the angle φ determining the phase of the correlations. The quantities γ and φ also determines the degree of mutual coherence between the two sources (see App. D for a more in depth discussion), which is given by DMC = γe iφ for these particular states. Accordingly, uncorrelated thermal sources (γ = 0) are fully incoherent, while perfectly correlated thermal sources (γ = 1) are mutually coherent [42,43]. For every other value of 0 < γ < 1, we are in presence of partially mutual coherent sources.
Let us start by considering some simplified expressions that are available when the covariance matrix V becomes diagonal (C = 0 in Eq. (45)), i.e. when the two thermal sources are uncorrelated (γ = 0), or correlated either perfectly in phase (φ = 0) or in antiphase (φ = π). For uncorrelated thermal sources, we recover the known expression [20,21] which is plotted as dashed black lines in Fig. 1. For in phase and antiphase correlated thermal states instead we obtain The results in Eq. (49) (blue and purple lines in Fig. 1) were obtained in [21] using an expression for the QFI only valid for photon-number diagonal states. With our general approach for arbitrary Gaussian states, we generalized these results to arbitrary values of γ and φ. Even though the general expressions are too cumbersome to be reported here, they can be analytically computed, and their behaviour is presented in Fig. 1 for an imaging system with a Gaussian PSF. In Fig. 1, we can see that while for large separations all curves converge to F d = 2κN 0 /w 2 , the small separation behaviour is seriously affected by the degree of mutual coherence γe iφ of the states. In particular, which means that, for small separations, for arbitrary values of γ, when |φ| < π/2 correlated thermal states perform worse than uncorrelated (i.e. fully incoherent) thermal states. This is due to the fact that for |φ| > π/2, we have destructive interference in the image plane, which enhances the information on the separation contained in the total number of detected photons in the image plane. This is particularly true for the case of mutually coherent thermal sources (γ = 1), where for antiphase correlations (φ = π), we have completely destructive interference for d = 0, which leads to a factor 2 enhancement of the QFI over the fully incoherent case. On the other hand, for in phase mutually coherent thermal sources (φ = 0), we have fully constructive interference for d = 0 which causes the QFI to vanish when d → 0.

C. Displaced thermal states
To illustrate how to compute the displacement contribution to the QFI (43), and to further elucidate the role of partial coherence, we now study displaced thermal states. In particular, we consider the following quantum state of the sources in the object plane where ρ(N th ) is a thermal state with mean photon number N th , and D(α) and D αe iφ describe displacements of the same amplitude separated by an angle φ in phase space. Without loss of generality, we can assume the first source to be displaced along the q quadrature (i.e. we assume α ∈ R), and write the covariance matrix and the displacement vector of the state (51) as andx 0 = 2(α, 0, α cos φ, α sin φ). For this kind of states, each source mode has a mean photon number N 0 = N th + α 2 . For a fair comparison, we will assume N 0 to be fixed and we define N th = (1 − γ)N 0 and α 2 = γN 0 . With these definitions the degree of mutual coherence between the sources reads DMC = γe iφ , which is the same form we obtained for correlated thermal states, with the difference that now γ is not the correlation strength, but the ratio between the coherent contribution to the mean photon number α 2 (coming from displacement) and the total photon number N 0 ; and φ here is the relative phase of the two displacements. In particular, now for γ = 1 the sources are in a coherent state, while for γ = 0 we obtain uncorrelated thermal states as in Sec. IV B. In the following, we discuss how this different source of mutual coherence affects the QFI for the estimation of the separation between the sources. To obtain the QFI, we first compute the covariance matrix and the displacement vector in the image plane, which are We can immediately observe that the covariance matrix (53) can be obtained from the one in Eq. (45) by setting the correlation strength to zero and replacing N 0 with N th = (1 − γ)N 0 . Accordingly, the covariance matrix contribution to the QFI is obtained by doing the substitution N 0 → (1−γ)N 0 in Eq. (48): It remains to compute the contribution coming from the displacement. This is given by Eq. (43), which using together with Eqs. (41) and (42) results in Combining Eqs. (55) and (57) we obtain the complete expression for the QFI F disp Fig. 2). The QFI for the separation between two displaced thermal sources in Fig. 2 (a-c) behaves similarly to the one for correlated thermal sources in Fig. 1, where the same values of the degree of mutual coherence γe iφ , and the mean photon number N 0 have been considered. In particular, both figures show the same behaviour for very small and large separations: On the other hand, displaced thermal sources (solid lines with dots in Fig. 2) have a larger QFI than correlated thermal sources (lighter lines in Fig. 2) in the intermediate regime of separations of the order of the PSF width w. Accordingly, in this regime, especially for high values of the amplitude of the mutual coherence γ, whether coherence originates from displacement or classical correlations has an impact on the separation estimation. On the contrary, in the low photon number limit (see Fig. 2 (d-f)), the two sources of coherence lead to substantially indistinguishable behaviours of the QFI, as given by which for φ = π/2 coincides with the constant QFI for incoherent thermal sources (γ = 0) [19]. It is interesting to note that for γ = 1, Eq. (55) vanishes and Eq. (57) reduces exactly to Eq. (58) for arbitrary values of N 0 . In other words, when the two sources are both in coherent states the QFI per photon does not depend on the mean photon number N 0 of sources. Accordingly, in this case, we have a shot-noise scaling of the estimation sensitivity ∆d ∼ N −1/2 0 . Conversely, the decay of the QFI per photon on distances of the order of the PSF width for increasing N 0 , for displaced or correlated thermal sources, leads to a less favorable scaling for intermediate source separations.
Recently, several theoretical and experimental works debated the role of (partial) coherence in the low photon number regime considered in Eq. (58) [44][45][46][47][48][49][50][51][52][53]. The origin of the debate are two different models for the image plane density matrix: the first restricts to the single photon subspace (ρ =ρ 1 ), while the second considers a sum of a vacuum and single pho-ton term (ρ = p 0ρ0 + p 1ρ1 ). Removing the vacuum contributionρ 0 , one neglects the losses induced by the imaging system that, in presence of coherence, induce parameter-dependent fluctuations of the photon number in the image plane, and therefore can increase the QFI [45,49,53]. However, this information becomes unavailable when the total number of photon in the object plane is unknown, and in this case one recovers the results of the pure single-photon model. A rigorous description of this scenario can be obtained with a multiparameter estimation approach for the simultaneous estimation of the mean number of emitted photons N 0 and the separation d [54]. The Gaussian model presented here always includes a vacuum contribution, and in particular for correlated thermal states in the low photon limit, the image plane quantum state takes the formρ = p 0ρ0 + p 1ρ1 . On the other hand, displaced thermal states are not photon-number diagonal, and, even for low brightness of the sources, they cannot be described by the two models discussed above. However, as observed in Eq. (58), in the low photon limit, this difference between the two states does not affect the QFI.

D. Squeezed states
As a final example, to illustrate how our approach applies to a quantum scenario, we consider two squeezed sources. In particular, we study two vacuum sources equally squeezed along two directions separated by an angle θ in phase space. For simplicity, in the following, we will focus on two particular cases: two sources squeezed along the same quadrature (θ = 0) and two sources squeezed along two opposite quadratures (θ = π/2). Without loss of generality, we assume the first source to be squeezed along the q quadrature, so that, in the object plane, the two sources are characterized by a zero mean Gaussian Wigner function, with covariance matrix with S (ξ) = diag(e −ξ , e ξ ) a single-mode squeezing matrix, and R(θ) a rotation matrix. The squeezing parameter ξ is related to the mean photon number N 0 emitted by each source by the relation N 0 = sinh 2 ξ [24,33]. In the image plane, such squeezed states remain zero mean Gaussian states, with covariance matrix with The derivative of the covariance matrix (60) with respect to the source separation is given by with To complete our calculation of the QFI, we need the Williamson decomposition (7) of the covariance matrix (60). Let us start from two sources squeezed in the same direction. In this case, for θ = 0, the action of the mapping to the symmetric and antisymetric source modes (see Eq. (39)) leaves the state unchanged. This state is separable, and so is its image as can be seen by the diagonal structure of the covariance matrix V θ=0 [55,56]. The Williamson decomposition of such a diagonal matrix results in a squeezing transformation applied to a thermal state in each mode (see App. E). Using this decomposition, we find that the only non zero coefficients in the QFI expansion (31) are a 1 11 , a 1 22 , a 2 11 and a 2 22 that determine the QFI contribution due to changes in the covariance matrix V 0 , and a 1 13 = a 1 31 , a 1 24 = a 1 42 , a 2 13 = a 1 31 , and a 2 24 = a 1 42 that give the QFI contribution due to deformations of the modes. The resulting QFI (whose analytical expression is to long to be reported here) is plotted as dashed lines in Fig. 3. Due to the significant losses (recall that κ < 0.5), these uncorrelated squeezed states become substantially semiclassical in the image plane. Accordingly, their QFI per photon behaves similarly to that of uncorrelated (incoherent) thermal sources, with a dip for intermediate separations of the order of the PSF width that gets deeper for larger mean photon numbers [20,21].
Let us now move to two sources squeezed in opposite directions (θ = π/2). In this case, the quantum state in the symmetric and antisymetric source modes is a two-mode squeezed vacuum, which is entangled, and remains entangled in the image plane for every non zero separation d > 0 [56]. For these entangled states, Eqs. (60) and (62), together with the Williamson decomposition of a covariance matrix in standard form (see App. E) imply that the only non-zero coefficients in Eq. (31) are a (1) 12 , a (1) 21 , a (2) 11 and a (2) 22 that establish the contribution to the QFI given by changes of the covariance matrix V π/2 , and a (1) 13 = a (1) 31 , a (1) 24 = a (1) 42 , a (2) 14 = a (2) 41 and a (2) 23 = a (2) 32 giving the contribution to the QFI coming from the changes of the modes. Combining all the different contributions, we obtain the expressions for the QFI plotted as solid lines in Fig. 3. Focusing on Fig. 3 (a), we note that for low brightness of the sources (N 0 1), and high losses (κ 0.25) the QFI for these entangled states, despite a slightly deeper dip, is fairly similar to that for independent squeezed states (dashed lines). On the other hand, for reduced losses (orange curve) the QFI per photon features a modest increase for intermediate separations. Arguably more interesting is the fact that, (see Fig. 3 (b) and (c)), for these entangled states the QFI per photon is substantially unaffected by the mean photon number N 0 (i.e. the squeezing) of the sources. Accordingly, in this case, similarly to what we observed for coherent states in Sec. IV C, we have a shot-noise scaling of the estimation sensitivity, ∆d ∼ N −1/2 0 , while non-entangled squeezed states provide a less favourable scaling for intermediate separations. with those in Fig. 2, we observe that for small separations, squeezing does not enhance the sensitivity when compared to a vast class of displaced or correlated thermal states with phase |φ| > π/2. In particular, for small separations, coherent sources with opposite phases provide a QFI which is twice as large as that of entangled squeezed sources.

V. CONCLUSION
In this paper, we determined the ultimate sensitivity limit for the estimation of a parameter encoded into the quantum state as well as the mode structure of a multimode Gaussian state of the electromagnetic field. In particular, we presented an analytical expression for the QFI, bounding the estimation sensitivity through the Cramér Rao lower bound, which can be calculated from the first two moments of the states and the dependence on the parameter of the mode functions. Such an expression expands the field of use of Gaussian quantum metrology to the estimation of parameters encoded into the spatio-temporal distribution of an electromagnetic signal.
We illustrated how to apply our general formalism by studying a paradigmatic imaging problem: the estimation of the transverse separation between two sources observed through a diffraction-limited imaging system. For this imaging problem, we dealt with two-mode Gaussian states with the parameter of interest encoded in the shape of the two image-plane modes as well as in their covariance matrix and displacement vector. We investigated the effects of different Gaussian resources such as partial coherence and squeezing.
For partial coherence, we computed the QFI for arbitrary complex degrees of mutual coherence DMC = γe iφ for arbitrarily bright Gaussian sources. We observed that for phases |φ| ≥ π/2, destructive interference provides a resolution enhancement in the regime of small sources' separations. Moreover, we studied both mutual coherence originating from thermal correlations and from displacement, and showed that for bright sources and intermediate separations (of the order of the PSF width) these two sources of coherence leads to different sensitivities. This difference was missed by previous studies of partial coherence, which only considered the low photon limit [44][45][46][47][48][49][50][51][52][53]. This is particularly important for coherent sources, for which the QFI per photon becomes independent on the brightness of the sources leading to a shot-noise scaling of the estimation error, i.e. ∆d ∼ N −1/2 0 .
In the study of squeezed sources, we distinguished between separable and entangled states in the image modes basis. On the one hand, we showed that for both these states the sensitivity is highly affected by the losses induced by the imaging system which render squeezed states semiclassical. Accordingly, they behave similarly to uncorrelated thermal states. on the other hand, entangled states present a QFI per photon which is not significantly affected by the brightnesses of the sources, while for separable states the QFI per photon for intermediate separations significantly decays with the sources' mean photon number. Accordingly, for separation of the order of the point spread function width, entangled states provide a sensitivity scaling enhancement when compared to separable states. However, in opposition to other quantum metrology scenarios, e.g. phase estimation [3][4][5], in this context this scaling enhancement does not go beyond the shot noise limit ∆d ∼ N −1/2 0 , which we observe for coherent sources. Our approach could be readily applied to other spatial mode-encoded parameter estimation scenarios in various field of science and technology ranging from astronomy to microscopy [18]. Moreover, parameters encoded into timefrequency modes appears in the characterization of frequency combs [17], or in radars that estimate the distance of a reflecting target from the temporal profile of chirped pulses [57,58] (recent studies have addressed this problem in the quantum regime [18,59]). Finally, the applicability of our approach could be further broadened by considering the simultaneous estimation of multiple parameters [60].

ACKNOWLEDGMENTS
We are very grateful to Claude Fabre for our illuminating discussions on the role of modes in quantum optical metrology, which have been a main source of inspiration for this work. We also thank Ilya Karuseichyk for useful discussion. This work was partially funded by French ANR under COSMIC project (ANR-19-ASTR0020-01). This work received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 899587. This work was supported by the European Union's Horizon 2020 research and innovation programme under the QuantERA programme through the project ApresSF. This work received funding from Ministerio de Ciencia e Innovación (MCIN) / Agencia Estatal de Investigación (AEI) for Project No. PID2020-115761RJ-I00 and support of a fellowship from "la Caixa" Foundation (ID 100010434) and from the European Union's Horizon 2020 research and innovation program under Marie Skłodowska-Curie Grant Agreement No. 847648, fellowship code LCF/BQ/PI21/11830025. This work was partially funded by CEX2019-000910-S [MCIN/AEI], Fundació Cellex, Fundació Mir-Puig, and Generalitat de Catalunya through CERCA.

Correlated thermal states
The image plane covariance matrix for correlated thermal states V in Eq. (45) does not feature any squeezing. Accordingly, the Bloch-Messiah reduction [62,63] guarantees that its Williamson decomposition can be achieved using a symplectic orthogonal matrix S ± and its symplectic eigenvalues coincide with its standard eigenvalues ν ± = N + + N − ± Z + 1 with N ± and Z given in Sec. IV B. It should be noted that these eigenvalues change order depending on the degree of coherence of the state, as well as the source separation: Moreover, in the case of perfectly correlated thermal sources (γ = 1) one of the two eigenvalues ν ± becomes one, i.e. one of the two modes in the Williamson decomposition becomes pure. We recall that in this case, the SLD (see Eqs. (16) and (17)) is not well defined, and the QFI should be computed by taking the limit γ → 1 at the end of the calculations [33].
Let us now assume that ν ± appear as entries of the diagonal matrix ν in decreasing order, the symplectic orthogonal matrix S ± such that V ± = S ± νS ± can be written as where S ± = (N + −N − )/|N + −N − |. We can notice that, for γ = 0 or φ = 0, the matrix V ± becomes diagonal. Correspondingly, Z = |N + − N − | which implies that X = 1 and Y = 0, so that S ± becomes the identity.

Squeezed states
Here, we present the Williamson decomposition of the covariance matrix V θ (see Eq. (60)) for the two particular cases discussed in Sec. IV D: two sources squeezed along the same quadrature (θ = 0) and two sources squeezed along opposite quadratures (θ = π/2).

a. Sources squeezed in the same direction
For two sources squeezed in the same direction, the image plane covariance matrix (60) reduces to the diagonal matrix with Since V N 0 is not proportional to the identity, even though V θ=0 is diagonal, it is not in Williamson form. To bring it in Williamson form, we rewrite V θ=0 as the covariance matrix of a pair of squeezed thermal states where Equating Eqs. (E4) and (E6), we obtain the following expressions for the mean photon numbers of the thermal states and for the squeezing parameters defining the symplectic matrices S ξ ± .