Gaussian quantum metrology for mode-encoded parameters

Quantum optical metrology aims to identify ultimate sensitivity bounds for the estimation of parameters encoded 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 displacement of a beam and to the temporal separation between two pulses. For these examples, we show how the estimation sensitivity can be enhanced by adding squeezing into specific modes.


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 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 (Helstrom, 1976;Holevo, 2011;Paris, 2009;Pezzè and Smerzi, 2014;Giovannetti et al., 2011;Tóth and Apellaniz, 2014).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 (Taylor and Bowen, 2016;Tsang, 2019), to remote sensing with lidars and radars (Giovannetti et al., 2001;Zhuang et al., 2017;Huang et al., 2021), to gravitational wave detection (Acernese et al., 2019;Tse et al., 2019).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 spatio-temporal distribution is conveniently described in terms of modes.i.e. normalized solutions of Maxwell's equations in vacuum (Fabre and Treps, 2020).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 context of mode-encoded parameter estimation focused on specific problems.For example, the case where the total light's intensity is not affected by the parameter, but its distribution among different modes is, was considered for the estimation of a small lateral beam displacement (Treps et al., 2002(Treps et al., , 2003)), or in the estimation of spectral parameters of a frequency comb (Cai et al., 2021).A general theory for this fixed-intensity scenario was recently presented by Gessner et al. (2022).Two different (mathematically equivalent) problems, that lately attracted a lot of attention, are the estimation of the separation between two point sources analysed through a diffraction-limited imaging system (Tsang et al., 2016;Paúr et al., 2016;Boucher et al., 2020) or the temporal separation between two pulses (Ansari et al., 2021;De et al., 2021;Mazelanik et al., 2022).For these problems, the parameter of interest is encoded in the shape of two (spatial or temporal) modes in the detection plane with separation-dependent pop-ulations (Lupo and Pirandola, 2016).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 (Holevo, 1975;Weedbrook et al., 2012;Adesso et al., 2014).
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 non-linear 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 (Caves, 1981), and is now a key ingredient of several quantumenhanced metrological schemes (Treps et al., 2003;Pezzé and Smerzi, 2008;Acernese et al., 2019;Tse et al., 2019).While previous studies of the quantum Fisher information for Gaussian states exist (Pinel et al., 2012;Monras, 2013;Šafránek et al., 2015;Jiang, 2014), 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 spatio-temporal 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 examples of such applications, we reconsider the estimation of the transverse displacement of a beam and the temporal separations between two pulses: For the former case, we extend the results for coherent beams of (Pinel et al., 2012) to thermal beams, we show that squeezing in the right mode provides a quantum enhancement also in this case, and we discuss how to include the effect of thermal noise and losses.For the latter, we confirm known results for thermal (Nair and Tsang, 2016;Lupo and Pirandola, 2016) and coherent pulses (Sorelli et al., 2022), and we investigate the possibility of a quantum enhancement populating additional modes with squeezed light.
Our paper is organized as follows: First, in Sec. 2, we recall 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. 3. Section 4 contains the application of our results to the estimation of the transverse displacement of a beam and the temporal separation of two pulses.Section 5 concludes our work.

Gaussian states
An N −mode continuous variable (CV) (Braunstein and van Loock, 2005;Serafini, 2017) 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 quadrature operators qk = (â k + â † k ) and pk = i(â † k − âk ), where â † k and âk are standard creation and annihilation operators.If we group all quadratures in the 2N −dimensional vector x = (q 1 , p1 , . . ., qn , pn ) , 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 (Holevo, 1975;Weedbrook et al., 2012;Adesso et al., 2014) which are completely determined by the displacement vector x = x , and the covariance matrix where {•, •} denotes the anticommutator, and • the expectation value For every physical state, the covariance matrix satisfies the uncertainty inequality (Simon et al., 1994) σ + iΩ ≥ 0. (6) A useful property of the covariance matrix is that, according to Williamson theorem (Williamson, 1936), 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.

Mode and state transformations
Let us recall that the space of solutions of Maxwell's equation is a Hilbert space (Fabre and Treps, 2020;Walschaers, 2021).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 (Fabre and Treps, 2020), which implies that the quadrature vector transforms according to with O an orthogonal symplectic matrix, i.e 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 (Holevo and Werner, 2001) where z 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.

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 operatorvalued measure (POVM) defined by the operators Kµ .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 (Helstrom, 1976;Holevo, 2011;Paris, 2009;Pezzè and Smerzi, 2014;Giovannetti et al., 2011) 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 (Braunstein and Caves, 1994).In general, the QFI can be computed as where Lϑ is the symmetric logarithmic derivative (SLD), implicitly defined by the equation (Helstrom, 1976;Holevo, 2011;Paris, 2009;Pezzè and Smerzi, 2014;Giovannetti et al., 2011) When ρ ϑ is an N −mode Gaussian state defined by the displacement vector x and the covariance matrix σ, the SLD is quadratic in the quadratures (Monras, 2013;Šafránek et al., 2015) where S and ν are the symplectic matrix and the symplectic eigenvalues obtained from the Williamson decomposition of σ as introduced in Eq. (7), while a 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 (Serafini, 2017), we can write the QFI for a Gaussian state as with where F σ and F x are the contribution to the QFI coming from variations of the covariance matrix σ and the displacement vector x, respectively.

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 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 vector xϑ 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) taking into account both these dependences.

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 does not 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 as xϑ = (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. (83b) 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 as 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).

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 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 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 Using Eq. ( 27), we can calculate the coefficients a (l) jk in Eq. (17c), which result in ij .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 state ρϑ 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).

Displacement vector contribution to the quantum Fisher information
We now move on to compute F x, as given by Eq. (19b), which takes into account the contribution to the QFI coming from variations of the displacement vector xI .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), we obtain Similarly to what we observed for F σ , F x only depends on the displacement vector xn 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 vector xn , 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 appear both (∂ ϑ xn [ϑ]) and D n .Accordingly, they combine mode variations with changes in the displacement vector.
4 Application to spatial and temporal resolution 4.1 Spatial beam positioning

A single populated mode
As a first example, we consider the estimation of the transverse displacement d of a light beam whose spatial profile is defined by the mode u 0 [d](r) = u(r − r 0 ) with r 0 = (d, 0), where, without loss of generality, we assumed the beam to be displaced along the x axis.Furthermore, we consider the mode u[d](r) to have a well-defined parity, s.t. it is orthogonal to its derivative: In this context, η quantifies the spatial extent of the beam we want to localize, e.g. for a Gaussian mode u(r) This estimation problem is fully defined by the mode u 0 [d](r).As a consequence, the mean field contribution to the QFI (34) simplifies to and we can write the covariance contribution to the QFI (31) as where we defined the coefficients where we recall σ x,y,z are Pauli matrices.
We can now evaluate Eqs.(35) and (36) for different states of the mode u 0 [d](r).Let us start by considering a coherent state |α , defined by the complex amplitude α that can be parameter dependent.Accordingly, we have x0 = 2(Re[α], Im[α]) and V 0 = 1 2 .In this case, is not hard to verify that the covariance matrix contribution (36) vanishes, F v = 0, and the QFI is fully determined by the displacement term (35), which reduces to where we introduced the mean photon number N 0 = |α| 2 .The second term in Eq. (39) presents a shotnoise scaling and is inversely proportional to the beam size: small displacements of a larger beam are harder to estimate.On the other hand, the first term in Eq. (39) takes into account how α depends on the transverse displacement of the beam.Such a dependence could be induced by position-dependent losses.
Let us now consider the localization of a thermal beam, for which we have x0 = 0 and V 0 = (2N 0 + 1)1 2 .As opposed to the coherent case discussed above, in this case the displacement contribution (35) vanishes, and the QFI is fully determined by the covariance matrix term (36).Since V 0 is proportional to the identity, the only nonzero coefficients in Eqs.(37 The 4η 2 N 0 term is identical to the one in Eq. (39).Accordingly, when the mean photon number N 0 does not depend explicitly on the transverse displacement, we have the same QFI for thermal and coherent beams.On the other hand, the explicit dependence of the mean photon number N 0 on the parameter induces a quite different dependence.To make this difference more explicit, we use N 0 = |α| 2 to rewrite this term in function of the mean photon number N 0 also in the coherent case.Accordingly, we get which is a factor N 0 +1 larger than the corresponding term in the thermal case.As a consequence, if the d−dependence of mean photon number dominates the QFI, such as in the case of strong displacement-dependent losses, coherent states provide a significant advantage over thermal states.This is due to the fact that for coherent states, a variation of the mean photon number consists in a change of mean field, while for thermal states it is a change of the covariance matrix, and the former is more efficient than the latter in making two Gaussian distributions distinguishable.In the following, we will see how our formalism recovers this result, to extend it to different states of the mode u 0 [d](r) and to take into account losses in the squeezed derivative mode.When populating the derivative mode, the mode Hilbert space H n , as introduced in Sec. 3, is spanned by u 0 [d](r) = u(r − r 0 ) and its normalized derivative u 1 [d](r) = ∂ d u 0 (r − r 0 )/η.On the other hand, the mode Hilbert space H ∂ (see Sec. 3) only contains the second derivative mode Furthermore, we assume that the derivative mode u 1 [d](r) has no mean field so that the mean field vector can be written as x = (q 0 , p 0 , 0, 0).Therefore, the mean field term of the QFI (34) results in (42) As noted by Pinel et al. (2012), the QFI (42) can be rewritten as a function of a unique element of the inverse covariance matrix where (V v ) 0,0 is the variance of the q−quadrature of mode Accordingly, for states with a nonzero mean field and a parameter-independent covariance matrix (e.g.coherent states), it is necessary and sufficient to squeeze the q quadrature of mode v 0 [d](r) to quantum enhance our beam positioning capability.It is interesting to observe that, if the mean field vector does not depend explicitly on the beam displacement d, i.e.
which is orthogonal to the mode u 0 [d](r) that defines the beam shape.In this case, this effect has been exploited experimentally to enhance position estimation with a so called quantum laser pointer (Treps et al., 2003).In a more general scenario, e.g. in presence of position-dependent losses, Eq. (44) prescribes to squeeze a mode v 0 [d](r) which is partially overlapping with u 0 [d](r).
Let us now discuss how the covariance matrix term of the QFI is modified by population in the derivative mode u 1 [d](r).For simplicity, we will focus on the case where the population of mode u 1 [d](r) is fully uncorrelated with that of mode u 0 [d](r), therefore, the covariance matrix takes the block diagonal form Under these assumptions, the covariance matrix contribution to the QFI (31) takes the form where we have defined the coefficients To further illustrate how to use Eq.(46) in practice, let us now consider the localisation of a thermal beam u 0 [d](r) aided by a squeezed vacuum state in the derivative mode u 1 [d](r).Accordingly, we have which corresponds to with ν 0 = 2N 0 + 1 and ν 1 = 1.Under these assumptions, the only nonzero coefficients are 16N 2 0 η 2 cosh 2 r, and d 2 1 = 2ζ 2 sinh 2 r.Substituting into Eq.(46), we obtain the following expression the QFI (for a zero mean state, the contribution in Eq. (43) vanishes) where we have introduced the number of photons N 1 = sinh 2 r in the squeezed derivative mode.Given that a thermal state has no preferred direction in phase space, we find that the result in Eq. ( 50) remains valid if we modify the squeezing direction.Furthermore, we can see that the QFI (50) is always larger than the one in Eq. (40) for a thermal state alone.This becomes particularly evident if we assume that N 0 does not explicitly depend on the parameter, and we consider the N 0 N 1 1 limit, where we have It is interesting to compare this result, with the quantum enhancement achievable with a coherent state in mode u 0 [d](r).For simplicity, let us consider the case where the mean field x0 does not depend explicitly on the transverse beam displacement d.In such a case, combining Eq. (43) with Eq. (50) (setting the number of thermal photons to zero), we obtain The second term is negligible for N 0 N 1 , and we obtain F d,coh−sq ∼ e 2r F d,coh .Accordingly, for the positioning a bright thermal beam aided with a squeezed state in the derivative mode, we have a quantum enhancement which is just a factor two smaller than that we obtain adding squeezing in the derivative mode u 1 [d](r) when the mode u 0 [d](r) is in a coherent state.We can understand this result by considering a thermal state as an ensemble average over coherent states with Gaussian distributed amplitudes and uniformly distributed phases.Accordingly, when adding squeezing in the derivative mode, the relative orientation between the coherent states in the ensemble and the squeezing will result sometimes in an enhancement and sometimes in a reduction of the sensitivity (see Eq. ( 43)).To make this statement more quantitative, we compute from Eqs.
While in general, the convexity of the QFI ensures th−sq , which supports our interpretation that the the quantum advantage enabled by squeezing for thermal states can be seen as an average over the sensitivity enhancements/diminutions obtained for coherent states.
We demonstrated above how squeezing in the normalized derivative mode u 1 [d](r) can lead to a sensitivity enhancement in the estimation of the displacement of a Gaussian beam.However, in practical situations it is hard to get a squeezed state which is not corrupted by noise.To illustrate what happens in these more practical scenarios, let us consider a thermal state in mode u 0 [d](r) and the derivative mode u 1 [d](r) populated with an arbitrary zero-mean Gaussian state, i.e. a squeezed thermal state.Accordingly, we have where r quantifies the squeezing strength, while N T quantifies the amout of thermal noise.Accordingly, the matrix S 1 in Eq. (49) remains the same, while the symplectic eigenvalue become ν 1 = 2N T + 1.The total photon number in the derivative mode for such a state is given by N 1 = N T + N S + 2N T N S , with the squeezing contribution given by N S = sinh 2 r.Following the same steps as above, we now obtain the following expression for the covariance matrix contribution to the QFI (46) which reduces to Eq. (50) when N T = 0. On the other hand, when N S = 0 and the population of the derivative mode becomes purely thermal, we obtain and it is not hard to show that F d,th−th (55) is always smaller than the F d,th−sq (50): unsurprisingly, populating the derivative mode with squeezing is always better than populating that with thermal noise.
In fact, for small values of N 1 , the QFI F d,th−th (55) is even smaller than that for an unpopulated derivative mode F d,th (40).To better illustrate this interplay between squeezing and thermal noise, we introduce the following parametrisation which allows to vary the amount of squeezing and thermal noise while keeping constant the total number of photons N 1 in the derivative mode.In particular, for χ = 1 the derivative mode is purely squeezed, while for χ = 0 it is purely thermal, so that we can refer to χ as the squeezing fraction.If we further assume that the beam we are trying to localize is Gaussian, i.e. u(r) = exp −|r| 2 /2w 2 / √ πw 2 , we show that for χ ≥ 1/2 and N 1 > 0, the QFI (54) is always larger than that for unpopulated derivative mode, i.e. for N 1 = 0. On the other hand, as presented in Fig. 1, for χ < 1/2 and small N 1 we obtain a worse sensitivity compared to that when the derivative mode is in vacuum.
Optical metrology protocols are generally very sensitive to photon losses, it is therefore important to illustrate how such losses can be taken account.Accordingly, it is useful to note that for a thermal state of mode u 0 [d](r), and an arbitrary zero-mean Gaussian state of the derivative mode u 1 [d](r), the QFI maintains the form (54) even after losses.In fact, it is sufficient to perform the following substitutions where κ 0 and κ 1 are the attenuation coefficients of the two modes u 0 [d](r), and u 1 [d](r), respectively; while N in 0 , N in T and N in S = sinh 2 r in are the populations of the mode u 0 [d](r), and the thermal and squeezing components of the population of the mode u 1 [d](r), respectively.Finally, in some applications, the attenuation coefficients κ 0 and κ 1 can be parameter dependent.In those cases, not only N 0 depends on the transverse displacement d (as taken into account by the first term in Eq. ( 54)), but also N T and r.This leads to an additional term in the QFI which takes the form

Temporal separation between pulses
As a second example, we consider the estimation of the time delay τ between two light pulses with the same temporal profile defined by the mode u(t), which for simplicity, we will assume to be real and even, i.e. u(t) = u(−t).From a parameter estimation point of view, this problem is most interesting when the separation τ between the pulses is smaller than (or comparable to) the pulse width.In this context, there is a finite overlap between the modes u(t − τ /2) and u(t + τ /2) (see Fig. 2) Accordingly, as discussed by Lupo and Pirandola (2016); Sorelli et al. (2021a,b) for the spatial domain, it is convenient to describe the problem in terms of the two orthonormal modes We are interested in computing the QFI for the estimation of τ , when the two modes (62), and eventually their derivatives, are populated.Accordingly, we complement the modes (62) with their orthonormalized first and second derivatives (see App. B for detailed calculations): where The shapes of the modes u i [τ ](t) and v i [τ ](t) for the specific case of Gaussian pulses u(t) = e −t 2 /2w 2 /(πw 2 ) 1/4 are presented in Fig. 2. Using the modes (62) and (63), we can express the matrices D n and D ∂ (see Sec. 3.2) as We have now specified all the mode-related quantities needed to compute the QFI for the estimation of the temporal separation τ between pulses.To proceed, we now make some further assumptions on the -1.0 -0.5 0 0.5 quantum state of the pulses.In particular, we consider the modes u 0 [τ ](t) and v 0 [τ ](t) to be in a general Gaussian state, and we allow for auxiliary population of the orthogonalized first derivative modes u 1 [τ ](t) and v 1 [τ ](t) with no mean field1 .Accordingly, we can write the mean field vector as x = (x 0 , 0), with x0 = (q u 0 , p u 0 , q v 0 , p v 0 ).The mean field term of the QFI (34) can then be expressed as Similarly to what we discussed in Sec.4.1, for every state with a nonzero mean field, i. .Therefore, the use of quantum resources, such as squeezing, to increase such matrix element can lead to an enhanced sensitivity (Pinel et al., 2012).
Let us now have a look at the covariance matrix contribution to the QFI (31).To this goal, we will make the simplifying assumption that the population of the derivative modes u 1 [τ ](t) and v 1 [τ ](t) is uncorrelated with that of the symmetric and antisymmetric superpositions u 0 [τ ](t) and v 0 [τ ](t) of the pulses we want to separate, so that we can write the covariance matrix in block diagonal form Under these assumptions, the covariance matrix contribution to the QFI takes the form where we introduced the coefficients with the matrices A (jk) l as defined in Sec.2.3 and ν j 0 (ν j 1 ) the symplectic eigenvalues of the covariance matrix V 0 (V 1 ).Accordingly, we have four groups of addends in the QFI (67): The first one, depending on the coefficients a jk l , describes the contribution of the population of the symmetric u 0 [τ ](t) and antisymmetric v 0 [τ ](t) superpositions of the two pulses.These terms are nonzero if and only if the covariance matrix V 0 explicitly depends on the temporal separation, i.e. ∂ τ V 0 0. Similarly, the third group of addends, depending on the coefficients c jk l , takes into account the population of the derivative modes u 1 [τ ](t) and v 1 [τ ](t), and is nonzero if and only if ∂ τ V 1 0. The second group of addends, containing the coefficients b jk l , takes into account how variations of the temporal separation τ leads to coupling between the modes u 0 Finally, the addends containing the coefficients d jk l take into account how due to variations of τ the derivative modes Let us now evaluate the QFI (67) for a specific quantum state of the two pulses.In particular, we are interested in two equally-bright fully-incoherent pulses whose intensity distribution is given by where we introduced the mean number of photons per pulse N 0 , and the electric field operator Ê(t) = j âj u j [τ ](t) + bj v j [τ ](t) , with âj and bj the annihilation operators associated with the even and odd modes u j [τ ](t) and v j [τ ](t), respectively (see Fig. 2).It is not hard to see that the intensity distribution I(t) (69) is achieved by a thermal state of the modes u 0 [τ ](t) and v 0 [τ ](t), with mean photon numbers Such a state has no mean field x0 = 0, so that its QFI is fully determined by Eq. ( 67), and has a covariance matrix (70) In Sec.4.1, we have seen that adding squeezing to the derivative mode improve the sensitivity, even for the spatial localization of an incoherent thermal beam.To verify, whether this is the case also for the temporal separation between two thermal pulses, we assume the derivative modes u 1 [τ ](t) and v 1 [τ ](t) to be populated by two independent, equally-squeezed vacuum states, described by the covariance matrix For such a quantum state, the QFI (67) takes the form (see App. C for the explicit calculation of the coefficients (68)) The behavior of the QFI (72), for Gaussian pulses, is plotted as red lines in Fig. 3.
For comparison, we will now also evaluate the QFI for the temporal separation of two equally bright fully coherent pulses.As opposed to Eq. ( 69), in this case the intensity distribution also contains an interference term depending on the relative phase φ between the coherent pulses Such an intensity distribution can be obtained by populating the modes u 0 [τ ](t) and v 0 [τ ](t) with coherent states, whose covariance matrix is the identity V 0 = 1 4 , and whose mean field is given by x0 = (x u , xv ) , with In particular, from Eqs. (74), we can see that for inphase (φ = 0) coherent pulses, the mean field is fully determined by the q quadrature of mode u 0 [τ ](t).
Similarly, when the two coherent pulses are out of phase (φ = π) the mean field is fully determined by the q quadrature of mode v 0 [τ ](t).As we did for thermal sources, we are going to consider the derivative modes u 1 [τ ](t) and v 1 [τ ](t) by two independent for the estimation of the temporal separation τ between two thermal (red) or coherent pulses, either in phase (blue) or out of phase (green), as a function of the temporal separation τ in units of the pulse width w.For each panel, we considered a mean photon number of N 0 = 1 per pulse, and different levels of squeezing in the derivative modes, as quantified by the parameter r = 0 (top), r = 0.5 (middle), and r = 1 (bottom).The pulse shape is assumed Gaussian u(t) = e −t 2 /2w 2 /(πw 2 ) 1/4 for all panels.
squeezed vacuum states (see Eqs. (66) and ( 71)).Under these assumptions, the covariance matrix contribution to the QFI can be obtained simply by setting N 0 = 0 into Eq.(72).Accordingly, we have where the displacement term F x,coh−sq can be com-puted from Eq. (65), and reads We can see that, when the two coherent pulses are either in phase (φ = 0) or out of phase (φ = π), the last line in Eq. (76) vanishes and only the squeezingenhanced term proportional to e 2r survives.This is consistent with the fact that the covariance matrix V 1 (71) presents squeezing along the q quadrature of modes u 1 [τ ](t) and v 1 [τ ](t), and for φ = 0, π the mean field (74) has a vanishing p quadrature.The QFI for the separation τ between in phase and out of phase coherent pulses are presented as blue and green lines in Fig. 3.
Let us now compare the expressions for the QFI for thermal and coherent pulses aided by squeezing in the derivative modes reported in Eqs. ( 72) and (75), respectively.We start by comparing the behaviours for vanishingly small separations τ → 0. In this regime (see App. B), we have and Accordingly, independently of the squeezing value r, the QFI for in phase coherent pulses vanishes for τ → 0, while that for out of phase coherent pulses is twice the one for incoherent pulses (see Fig. 3 where for Gaussian pulses we have (∆k) 2 = 1/2w 2 ).
To better understand this behavior, let us recall that the quantum state of the finite overlap δ between the two pulses induces a τ −dependent population of the symmetric and antisymmetric modes u 0 [τ ](t) and v 0 [τ ](t).It is this dependence on temporal separation, which enters the QFI through ∂ τ x0 (in the coherent case) and ∂ τ V 0 (in the incoherent case), that dominates the QFI behavior for τ → 0. This implies that the population of the derivative modes u 1 [τ ](t) and v 1 [τ ](t), and in the particular the squeezing thereof, has no impact on the τ → 0 behavior of the QFI.
On the contrary, for separations much larger than the pulses' width, i.e. τ ∆k 1, the overlap δ tends to zero, and the populations of the modes u 0 [τ ](t) and v 0 [τ ](t) become parameter independent.The QFI is then dominated by the noise in derivative modes u 1 [τ ](t) and v 1 [τ ](t).In particular, we have Accordingly, for large temporal separations τ we have a squeezing enhancement.Such an enhancement is always larger for coherent pulses than for thermal pulses.However, similarly to what we observed for the spatial localization of a beam, the QFI enhancement for large τ in the coherent case is at most a factor two larger than that in the thermal one.

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 two paradigmatic problems: the estimation of the transverse displacement of a beam, and of the temporal separation between two pulses.
In the study of the transverse displacement we showed that if the mean photon number of the beam is independent of its transverse position, the displacement of a coherent and thermal beam can be estimated with the same sensitivity.On the other hand, if the mean number of photons N 0 in the beam depends on its transverse displacement, e.g. because of position dependent losses, this dependence adds an additional term to the QFI which is ∼ N 0 times larger for coherent beams than for thermal ones.Furthermore, we showed that the sensitivity in the estimation of a transverse displacement can be enhanced by adding squeezing to a mode shaped like the derivative of the beam.Such a squeezing-enabled quantum enhancement is at most a factor two larger for coherent beams than for thermal ones.
We then moved to the time domain and considered the estimation of the temporal separation between two coherent or thermal pulses.Such pulses are described by two temporal modes (the symmetric and anti-symmetric superpositions of the pulses) whose shape and populations depend on the separation parameter τ .We showed that the interplay between these two dependences plays a fundamental role in the choice of which modes one needs to squeeze to achieve a quantum enhancement.For large temporal separations, when the pulses have a negligible overlap, they are are most sensitive to the changes in the mode shapes.Accordingly, in this regime a quantum enhancement is possible by adding squeezing to the derivatives of the symmetric and anti-symmetric superpositions of the pulses.As for the transverse displacement estimation, the quantum enhancement achieved for coherent pulses is at most a factor two larger than the one obtained for thermal ones.On the other hand, for small temporal separations, when the pulses have a significant overlap, the QFI is dominated by how photons redistributes among the symmetric and anti-symmetric superpositions of the pulses.As a consequence, populating the derivative modes has no effect on the sensitivity in this regime.
Our approach could be readily applied to other mode-encoded parameter estimation scenarios in various field of science and technology ranging from astronomy to microscopy (Gessner et al., 2022).Moreover, parameters encoded into time-frequency modes appears in the characterization of frequency combs (Cai et al., 2021), or in radars that estimate the distance of a reflecting target from the temporal profile of chirped pulses (Van Trees, 2002, 2001) (recent studies have addressed this problem in the quantum regime (Zhuang and Shapiro, 2022;Gessner et al., 2022)).Finally, the applicability of our approach could be further broadened by considering the simultaneous estimation of multiple parameters (Nichols et al., 2018).source of inspiration for this work.We also thank Ilya Karuseichyk for useful discussion.This work was partially funded by French ANR under COS-MIC 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 carried out during the tenure of an ERCIM 'Alain Bensoussan' Fellowship Programme.

A Properties of the basis-change matrices
Here we derive some useful properties of the matrix O, and its blocks O n and O N from the properties of the n initially populated modes {u k [ϑ](r, t)} and their derivatives {∂ ϑ u k [ϑ](r, t)}.Let us start by recalling the following Hilbert space definitions: We now assume that m is the number of derivatives We can now choose the modes u k [ϑ](r, t) as the first m among the N −n auxiliary vacuum modes that we use to describe the quantum state of the system in the parameter dependent basis.In light of this, it is convenient to further decompose the matrix O N −n as where O ∂ and O N −n−m are matrices of dimensions 2N ×2m and 2N ×2(N −n−m), respectively.From the orthogonality of O, we can obtain the following relations: where the sum in Eq. (83b) runs over the total number of column blocks we decomposed the matrix O into.
Let us now compute the derivative of the matrix O n .In Sec.2.2, we have seen that that O is composed by 2 × 2 blocks containing the mode overlaps.Accordingly, it is sufficient to specify the derivative of the kl block of O n , which reads ) Combining Eqs.(84) and (29), we obtain where D n and D ∂ are a 2n × 2n and a 2n × 2m matrices, respectively.Their kl blocks are given by Using Eqs.(85) and (83), we then obtain Let us conclude this appendix with few words on the coefficients c kl [ϑ] and c kl [ϑ].The first are simply given by the overlaps of the initially populated modes with their derivatives c kl ).On the other hand, there exist several orthonormalization methods that can be used to construct the modes u k [ϑ](r, t), leading to different expressions for the coefficients c kl [ϑ].For example, using the Gram-Schmidt procedure, the modes u k [ϑ](r, t) can be constructed iteratively as Accordingly, the coefficients c kl [ϑ] are given by where we introduced δ = ∂ τ δ.We assumed that our pulses are symmetric, i.e. u(t) = u(−t).We thus have (∂ t u|u) = 0, which, combined with Consequently, the orthogonalised first derivative modes are simply given by where we introduced We now move to the construction of the orthonormalized second derivatives.The fact that the modes u 1 [τ ](t) and v 1 [τ ](t) have, by construction, opposite parity implies that where in the last step we used that u 0 [τ ](t) and v 0 [τ ](t) are even and odd functions of t, respectively.The second derivative of u 0 (x) and v 0 (x) can be rewritten as with where we have introduced the second derivative of the overlap parameter δ where we used that the pulse shape u(t) is an even function of t.Therefore, the orthogonalized second derivative modes are given by Let us explicitly calculate ξ u and ξ v .This can be achieved by using Eqs.(98a), (98b) and by noting that where we used partial integration and made the reasonable assumption that the pulse shape u(t) goes to zero at infinity.We then get Let us now compute explicitly the normalization constants of the modes u 2 [τ ](t) and v 2 [τ ](t) We can then expand From Eqs. (98a) and (98b), we then have and where we defined Since the expressions of the mode quantities (especially ζ u and ζ v ) computed above for a generic pulse shape u(t) are fairly complicated, we present their explicit expressions for a Gaussian pulse u(t) = e −t 2 /2w 2 /(πw 2 ) 1/4 in the following: Note that the fact that η 2 u,v = ξ 2 u,v is a peculiarity of Gaussian pulses and it is not true in general.

B.2 Small τ behaviour
Arguably, the most interesting regime for temporal separation estimation is that of small τ .Therefore, in the following, we discuss the behaviour of the quantities computed above for τ → 0. Let us start by considering the following series expansions This behaviour implies that the contribution to the QFI coming from the coefficients b jk l and d jk l vanishes for τ → 0 (see Eqs.(68) and App.C).

C Calculation of the QFI (72)
In this Appendix, we explicitly compute the coefficients (68) that lead to the QFI (72) for the estimation of the temporal separation τ between two incoherent thermal pulses aided by two squeezed vacuum states in the derivative modes u 1 [τ ](t) and v 1 [τ ](t) defined by the covariance matrices V 0 (70) and V 1 (71).
First, we note that V 0 (70) is already in Williamson form.Therefore, the symplectic matrix S 0 entering in Eqs. with symplectic eigenvalues ν 0,1 1 = 1.To compute the coefficients a jk l , we need the derivative of the matrix V 0 (70).The latter depends on the temporal separation τ only through the overlap parameter δ.Accordingly, we have Since ∂ τ V 0 is diagonal, the only nonzero a jk l coefficients are , where we introduced the matrix ) with We assumed the covariance matrix V 1 of the derivative modes u 1 [τ ](t) and v 1 [τ ](t) to be parameter independent, i.e. ∂ τ V 1 = 0, which leads to c (jk) l = 0 for all l, j and k. with (129) As consequence, the only nonzero d jk l coefficients are Substituting the coefficients in Eqs.(124), ( 127) and (130) into Eq.(67), we then obtain the QFI (72).
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 4.1.2Populating the derivative mode It was demonstrated by Pinel et al. (2012); Gessner et al. (2022), that the QFI (39) can be enhanced by adding squeezing to the derivative mode ∂ d u 0 [d](r).

Figure 1 :
Figure 1: QFI (normalized by its maximum value when u 1 [d](r) is in vacuum) for the estimation of the transverse shift d of a thermal Gaussian beam with mean photon number N 0 = 10 (top) and N 0 = 1 (bottom), assisted by a thermal squeezed state with mean photon number N 1 , as function of the squeezing fraction χ.

Figure 2 :
Figure 2: (top) Two Gaussian pulses temporally separated by τ = w, with their overlap δ represented as a shaded area.(bottom) Set of orthonormal modes constructed from the two Gaussian pulses and their first and second derivatives with respect to τ .All modes have definite parity, in particular the modes u i [τ ](t) (leftred) are even functions, while the modes v i [τ ](t) (right -blue) are odd functions.

Figure 3 :
Figure3: QFI (Normalized to its maximum value for thermal states in modes u 0 [τ ](t) and v 0 [τ ](t)) for the estimation of the temporal separation τ between two thermal (red) or coherent pulses, either in phase (blue) or out of phase (green), as a function of the temporal separation τ in units of the pulse width w.For each panel, we considered a mean photon number of N 0 = 1 per pulse, and different levels of squeezing in the derivative modes, as quantified by the parameter r = 0 (top), r = 0.5 (middle), and r = 1 (bottom).The pulse shape is assumed Gaussian u(t) = e −t 2 /2w 2 /(πw 2 ) 1/4 for all panels.
This work was funded by MCIN/AEI/10.13039/501100011033and the European Union "NextGenerationEU" PRTR fund [RYC2021-031094-I].This work has been founded by the Ministry of Economic Affairs and Digital Transformation of the Spanish Government through the QUANTUM ENIA project call -QUAN-TUM SPAIN project, by the European Union through the Recovery, Transformation and Resilience Plan -NextGenerationEU within the framework of the Digital Spain 2026 Agenda, and by the CSIC Interdisciplinary Thematic Platform (PTI+) on Quantum Technologies (PTI-QTEP+).

Finally
r, t) that are linearly independent from the n initially populated modes, i.e. dim(H ∂ ) = m.Accordingly, up to a reordering of the basis {∂ ϑ u k [ϑ](r, t)}, we can always construct a basis of H ∂ using the orthonormalized version u k [ϑ](r, t) of the derivatives of the first m initially populated modes ∂ ϑ u k [ϑ](r, t).
In this appendix, we construct the orthogonalized first and second derivatives of the modes u 0 [d](t) and u 0 [d](t).Let us start by computing the derivatives of Eqs.(62) with respect to the parameter τ :