# Gaussian quantum metrology for mode-encoded parameters

Giacomo Sorelli  Laboratoire Kastler Brossel, Sorbonne Université, ENS-
Université PSL, CNRS, Collège de France, 4 Place Jussieu, F-75252 Paris,
France  Fraunhofer IOSB, Ettlingen, Fraunhofer Institute of Optronics, System
Technologies and Image Exploitation, Gutleuthausstr. 1, 76275 Ettlingen,
Germany  Manuel Gessner  Departament de Física Teòrica, IFIC, Universitat de
València, CSIC, C/ Dr. Moliner 50, 46100 Burjassot (València), Spain  Nicolas
Treps  Laboratoire Kastler Brossel, Sorbonne Université, ENS-Université PSL,
CNRS, Collège de France, 4 Place Jussieu, F-75252 Paris, France  Mattia
Walschaers  Laboratoire Kastler Brossel, Sorbonne Université, ENS-Université
PSL, CNRS, Collège de France, 4 Place Jussieu, F-75252 Paris, France

###### Abstract

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.

##  1 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
$\Delta\vartheta$ of the estimator of a parameter $\vartheta$ is inversely
proportional to the quantum Fisher information of the quantum state
$\hat{\rho}_{\vartheta}$ 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, 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 populations (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 quantum-enhanced 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.

##  2 Preliminaries

###  2.1 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 $\left\\{u_{k}({\bf r},t)\right\\}_{k=1}^{N}$ of solutions of Maxwell’s
equations, orthonormal with respect to the inner product

$\left(u_{k}|u_{l}\right)=\int d^{3}{\bf r}u_{k}^{*}({\bf r},t)u_{l}({\bf
r},t)=\delta_{kl},$ (1)

and associating with each mode $u_{k}({\bf r},t)$ a pair of quadrature
operators $\hat{q}_{k}=(\hat{a}_{k}+\hat{a}_{k}^{\dagger})$ and
$\hat{p}_{k}=i(\hat{a}_{k}^{\dagger}-\hat{a}_{k})$, where
$\hat{a}_{k}^{\dagger}$ and $\hat{a}_{k}$ are standard creation and
annihilation operators. If we group all quadratures in the $2N-$dimensional
vector $\hat{\bf
x}=(\hat{q}_{1},\hat{p}_{1},\dots,\hat{q}_{n},\hat{p}_{n})^{\top}$, from the
canonical commutation relations
$[\hat{a}_{k},\hat{a}^{\dagger}_{l}]=\delta_{kl}$ for annihilation operators,
we obtain

$[\hat{x}_{j},\hat{x}_{k}]=2i\Omega_{jk},$ (2)

with the symplectic form $\Omega=\bigoplus_{k=1}^{N}\omega_{k}$, and
$\omega_{k}=i\sigma_{y}$, where we have introduced the notation
$\sigma_{i=x,y,z}$ for the standard $2\times 2$ Pauli matrices. Our preferred
phase space representation of an $N-$mode CV quantum state with density matrix
$\hat{\rho}$, and characteristic function $\chi({\bf y})=\tr[\exp\left(-i{\bf
y}^{\top}\Omega\hat{\bf x}\right)\hat{\rho}]$, is the Wigner function

$W({\bf x})=\int\frac{d{\bf y}^{2N}}{(2\pi)^{2N}}e^{-i{\bf y}^{\top}\Omega{\bf
x}}\chi({\bf y}).$ (3)

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)

$W({\bf x})=\frac{\exp\left[-({\bf x}-\bar{\bf x})^{\top}\sigma^{-1}({\bf
x}-\bar{\bf x})/2\right]}{(2\pi)^{N}\sqrt{\det\sigma}},$ (4)

which are completely determined by the displacement vector $\bar{\bf
x}=\langle\hat{\bf x}\rangle$, and the covariance matrix

${\bf\sigma}_{jk}=\frac{1}{2}\left\langle\left\\{(\hat{x}_{j}-\bar{x}_{j}),(\hat{x}_{k}-\bar{x}_{k})\right\\}\right\rangle,$
(5)

where $\\{\cdot,\cdot\\}$ denotes the anticommutator, and
$\langle\cdot\rangle$ the expectation value
$\langle\cdot\rangle=\tr\left[\cdot\hat{\rho}\right]$. For every physical
state, the covariance matrix satisfies the uncertainty inequality (Simon  _et
al._ , 1994)

${\bf\sigma}+i\Omega\geq 0.$ (6)

A useful property of the covariance matrix is that, according to Williamson
theorem (Williamson, 1936), it can be decomposed as

${\bf\sigma}=S{\bf\nu}S^{\top},\quad{\rm
with}\quad{\bf\nu}=\bigoplus_{k=1}^{N}\nu_{k}\mathds{1}_{2},$ (7)

where $S$ is a symplectic matrix, i.e. $S\Omega S^{\top}=\Omega$, we
introduced the notation $\mathds{1}_{n}$ for the $n-$dimensional identity
matrix, and ${\bf\nu}$ 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. $\nu_{k}\geq 1$.

###  2.2 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

$\displaystyle u_{k}({\bf r},t)$ $\displaystyle=\sum_{l=1}^{N}U_{kl}v_{l}({\bf
r},t),$ (8a) $\displaystyle v_{k}({\bf r},t)$
$\displaystyle=\sum_{l=1}^{N}U_{lk}^{*}u_{l}({\bf r},t),$ (8b)

with $U_{kl}=\left(v_{l}|u_{k}\right)$. 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

$\hat{\bf x}^{\prime}=O\hat{\bf x},$ (9)

with $O$ an orthogonal symplectic matrix, i.e
$OO^{\top}=O^{\top}O=\mathds{1}_{2N}$ and $O^{\top}\Omega O=\Omega$, with
elements $O_{2k-1,2l-1}=\real\left[U_{kl}\right]$,
$O_{2k,2l-1}=-\imaginary\left[U_{kl}\right]$,
$O_{2k-1,2l}=\imaginary\left[U_{kl}\right]$ and
$O_{2k,2l}=\real\left[U_{kl}\right]$.

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)

$\displaystyle\bar{\bf x}^{\prime}$ $\displaystyle=\mathcal{T}\bar{\bf
x}+\bar{\bf z},$ (10a) $\displaystyle{\bf\sigma}^{\prime}$
$\displaystyle=\mathcal{T}{\bf\sigma}\mathcal{T}^{\top}+\mathcal{N},$ (10b)

where $\bar{\bf z}$ is a real $2N-$dimensional vector, while $\mathcal{T}$ and
$\mathcal{N}$ are $2N\times 2N$ real matrices with
$\mathcal{N}=\mathcal{N}^{\top}$ and satisfying the positivity condition
$\mathcal{N}+i\mathcal{T}\Omega\mathcal{T}^{\top}\geq i\Omega$. From Eq. (9)
is easy to see that a mode basis change is a Gaussian channel (see Eqs. (10))
with $\mathcal{N}=0$ and $\mathcal{T}=O$.

###  2.3 Quantum estimation theory

Let us now assume that we want to estimate a parameter $\vartheta$ encoded in
a quantum state $\hat{\rho}_{\vartheta}$ from $M$ independent measurements of
a given positive operator-valued measure (POVM) defined by the operators
$\hat{K}_{\mu}$. Using classical post-processing techniques, from the
measurements’ results, we can extract an unbiased estimator
$\tilde{\vartheta}$ of the parameter as well as its standard deviation
$\Delta\tilde{\vartheta}$. 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)

$\Delta\tilde{\vartheta}\geq\frac{1}{\sqrt{M\mathcal{F}_{\vartheta,\hat{K}_{\mu}}}}$
(11)

with the Fisher information $\mathcal{F}_{\vartheta,\hat{K}_{\mu}}$ defined by

$\mathcal{F}_{\vartheta,\hat{K}_{\mu}}=\sum_{\mu}p(\mu|\vartheta)\left[\partial_{\vartheta}\log\left(p(\mu|\vartheta)\right)\right]^{2}$
(12)

where $p(\mu|\vartheta)=\Tr\left[\hat{K}_{\mu}\hat{\rho}_{\vartheta}\right]$
is the conditional probability of obtaining the result $\mu$ for a given value
of $\vartheta$, and we introduced the compact notation
$\partial_{\vartheta}\cdot=\partial\cdot/\partial\vartheta$ for the
derivative. The Fisher information optimized over all possible POVMs

$F_{\vartheta}=\max_{\hat{K}_{\mu}}\mathcal{F}_{\vartheta,\hat{K}_{\mu}},$
(13)

is the quantum Fisher information (QFI), and establishes the ultimate
metrological sensitivity (Braunstein and Caves, 1994). In general, the QFI can
be computed as

$F_{\vartheta}=\Tr\left[\hat{\mathcal{L}}^{2}_{\vartheta}\hat{\rho}_{\vartheta}\right],$
(14)

where $\hat{\mathcal{L}}_{\vartheta}$ 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)

$2\partial_{\vartheta}\hat{\rho}_{\vartheta}=\hat{\mathcal{L}}_{\vartheta}\hat{\rho}_{\vartheta}+\hat{\rho}_{\vartheta}\hat{\mathcal{L}}_{\vartheta}.$
(15)

When $\rho_{\vartheta}$ is an $N-$mode Gaussian state defined by the
displacement vector $\bar{\bf x}$ and the covariance matrix $\sigma$, the SLD
is quadratic in the quadratures (Monras, 2013; Šafránek  _et al._ , 2015)

$\hat{\mathcal{L}}_{\vartheta}=L^{(0)}_{\vartheta}+{\bf
L}^{(1)\top}_{\vartheta}\hat{\bf x}+\frac{1}{2}\hat{\bf
x}^{\top}L^{(2)}_{\vartheta}\hat{\bf x}$ (16)

with

$\displaystyle L^{(0)}_{\vartheta}$ $\displaystyle=-\frac{1}{2}\Tr\left[\sigma
L_{\vartheta}^{(2)}\right]-{\bf L}^{(1)\top}_{\vartheta}\bar{\bf
x}-\frac{1}{2}\bar{\bf x}^{\top}L^{(2)}_{\vartheta}\bar{\bf x}$ (17a)
$\displaystyle{\bf L}^{(1)}_{\vartheta}$
$\displaystyle=\sigma^{-1}\left(\partial_{\vartheta}\bar{\bf
x}\right)-L^{(2)}_{\vartheta}\bar{\bf x}$ (17b) $\displaystyle
L^{(2)}_{\vartheta}$
$\displaystyle=\frac{1}{2}\sum_{l=0}^{3}\sum_{jk=1}^{N}\frac{a_{jk}^{(l)}}{\nu_{j}\nu_{k}-(-1)^{l}}\left(S^{\top}\right)^{-1}A_{jk}^{(l)}S^{-1}$
(17c)

where $S$ and $\nu$ are the symplectic matrix and the symplectic eigenvalues
obtained from the Williamson decomposition of $\sigma$ as introduced in Eq.
(7), while
$a_{jk}^{(l)}=\Tr\left[A_{jk}^{(l)}S^{-1}(\partial_{\vartheta}\sigma)(S^{\top})^{-1}\right]$,
with $A_{jk}^{(l)}$ a set $2n\times 2n$ matrices that are zero everywhere
except in the $jk$ block where they are given by $\omega/\sqrt{2}$,
$\sigma_{z}/\sqrt{2}$, $\mathds{1}_{2}/\sqrt{2}$ and $\sigma_{x}/\sqrt{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

$F_{\vartheta}=F_{\sigma}+F_{\bar{\bf x}},$ (18)

with

$\displaystyle F_{\sigma}$
$\displaystyle=\Tr\left[L^{(2)}_{\vartheta}(\partial_{\vartheta}\sigma)\right]$
(19a)
$\displaystyle=\frac{1}{2}\sum_{l=0}^{3}\sum_{jk=1}^{N}\frac{\left(a_{jk}^{(l)}\right)^{2}}{\nu_{j}\nu_{k}-(-1)^{l}},$
$\displaystyle F_{\bar{\bf x}}$ $\displaystyle=(\partial_{\vartheta}\bar{\bf
x})^{\top}\sigma^{-1}(\partial_{\vartheta}\bar{\bf x}).$ (19b)

where $F_{\sigma}$ and $F_{\bar{\bf x}}$ are the contribution to the QFI
coming from variations of the covariance matrix $\sigma$ and the displacement
vector $\bar{\bf x}$, respectively.

##  3 Mode-encoded parameter estimation

Let us consider the estimation of a parameter $\vartheta$ encoded into a
Gaussian state $\hat{\rho}_{\vartheta}$ expressed into an $n-$dimensional mode
basis $u_{k}[\vartheta]({\bf 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 $\mathcal{H}_{n}={\rm span}\left(\\{u_{k}[\vartheta]({\bf
r},t)\\}\right)$. Since every basis of the mode Hilbert space
$\mathcal{H}_{n}$ would provide a description of the quantum state
$\hat{\rho}_{\vartheta}$ in terms of the smallest number $n$ of modes, the
choice of the mode basis $u_{k}[\vartheta]({\bf r},t)$ is not unique. Despite
this freedom of choice, in the most general parameter estimation scenarios,
every mode basis $u_{k}[\vartheta]({\bf r},t)$ will be parameter-dependent.
The latter fact implies that the Gaussian state $\hat{\rho}_{\vartheta}$
depends on $\vartheta$ not only explicitly through the displacement vector
$\bar{\bf x}_{\vartheta}$ and the covariance $\sigma_{\vartheta}$, but also
implicitly through the mode functions $u_{k}[\vartheta]({\bf r},t)$. Our goal
is this section is to calculate the QFI (18) taking into account both these
dependences.

###  3.1 Separation of mode and state parameter dependence

Our first step is to make the parameter dependence coming from the modes
$u_{k}[\vartheta]({\bf r},t)$ explicit in the covariance matrix and
displacement vector of the quantum state $\hat{\rho}_{\vartheta}$. To this
goal, we express them into a parameter-independent basis $v_{k}({\bf r},t)$:
Using Eq. (9), we get

$\displaystyle\sigma_{I}$ $\displaystyle=O\sigma_{\vartheta}O^{T},$ (20a)
$\displaystyle\bar{\bf x}_{I}$ $\displaystyle=O\bar{\bf x}_{\vartheta},$ (20b)

where we introduced the subscripts $I$ and $\vartheta$ to identify quantities
in the parameter-independent and parameter-dependent bases, respectively.
Naturally, the choice of the parameter-independent basis $v_{k}({\bf 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 $\hat{\rho}_{\vartheta}$, the parameter independent basis $v_{k}({\bf
r},t)$ must have dimension $N\geq 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
$\sigma_{\vartheta}$ in block diagonal form as

$\sigma_{\vartheta}=\begin{pmatrix}V_{n}&0\\\
0&\mathds{1}_{2(N-n)}\end{pmatrix},$ (21)

and the displacement vector as $\bar{\bf x}_{\vartheta}=(\bar{\bf
x}^{\top}_{n},0,\cdots,0)^{\top}$. To isolate the action of $O$ on the $n$
initially populated modes, it is convenient to rewrite it as a $1\times 2$
block matrix

$O=\begin{pmatrix}O_{n}&O_{N-n}\end{pmatrix},$ (22)

with $O_{n}$ and $O_{N-n}$ a $2N\times 2n$ and a $2N\times 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

$\sigma_{I}=O_{n}(V_{n}-\mathds{1}_{2n})O^{\top}_{n}+\mathds{1}_{2N}.$ (23)

Analogously, using Eq. (22) into Eq. (20b), we can rewrite the displacement
vector as

$\bar{\bf x}_{I}=O_{n}\bar{\bf x}_{n}.$ (24)

Equations (23) and (24) provide a description of the Gaussian state
$\hat{\rho}_{\vartheta}$ where the parameter dependence is fully expressed in
the covariance matrix $\sigma_{I}$ and the displacement vector ${\bf 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).

###  3.2 Covariance matrix contribution to the quantum Fisher information

We start with the calculation of $F_{\sigma}$ (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
$\sigma_{I}$ in the parameter independent basis with respect to the parameter

$\displaystyle\partial_{\vartheta}\sigma_{I}$
$\displaystyle=\left(\partial_{\vartheta}O_{n}\right)(V_{n}-\mathds{1}_{2n})O^{\top}_{n}$
(25)
$\displaystyle\;+O_{n}(V_{n}-\mathds{1}_{2n})\left(\partial_{\vartheta}O^{\top}_{n}\right)+O_{n}\left(\partial_{\vartheta}V_{n}\right)O^{\top}_{n}.$

To compute the quadratic term of the SLD $L^{(2)}_{\vartheta}$ (see Eq.
(17c)), we need the Williamson decomposition
$\sigma_{I}=S_{I}\nu_{I}S^{T}_{I}$ of the covariance matrix $\sigma_{I}$.
Using Eqs. (21) and (22), we can connect it to the Williamson decomposition
$V_{n}=S_{n}\nu S^{T}_{n}$ of the covariance of the $n$ initially populated
modes in the parameter dependent basis $u_{n}[\vartheta]({\bf r},t)$, and
obtain

$\displaystyle S_{I}$
$\displaystyle=\begin{pmatrix}O_{n}S_{n}&O_{N-n}\end{pmatrix},$ (26a)
$\displaystyle\nu_{I}$ $\displaystyle=\nu\oplus\mathds{1}_{2(N-n)}.$ (26b)

Accordingly, using the properties of the matrix $O$ (see App. A for details),
we can write

$S^{-1}_{I}\left(\partial_{\vartheta}\sigma_{I}\right)\left(S^{\top}_{I}\right)^{-1}=\begin{pmatrix}B_{n}&B_{\partial}^{\top}&0\\\
B_{\partial}&0&0\\\ 0&0&0\end{pmatrix},$ (27)

with

$\displaystyle B_{n}=\>$ $\displaystyle
S_{n}^{-1}D^{\top}_{n}\left(V_{n}-\mathds{1}_{2n}\right)\left(S_{n}^{\top}\right)^{-1}$
$\displaystyle+S_{n}^{-1}\left(V_{n}-\mathds{1}_{2n}\right)D_{n}\left(S_{n}^{\top}\right)^{-1}$
(28a)
$\displaystyle+S_{n}^{-1}\left(\partial_{\vartheta}V_{n}\right)\left(S_{n}^{\top}\right)^{-1},$
$\displaystyle B_{\partial}=\>$ $\displaystyle
D_{\partial}^{\top}\left(V_{n}-\mathds{1}_{2n}\right)\left(S_{n}^{\top}\right)^{-1}.$
(28b)

Here, $D_{n}$ is a $2n\times 2n$ matrix and $D_{\partial}$ is a $2n\times 2m$
matrix, constructed using, respectively, the coefficients $c_{kl}[\vartheta]$
and $c^{\prime}_{kl}[\vartheta]$ of the expansion

$\displaystyle\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)$
$\displaystyle=\sum_{l=1}^{n}c_{kl}[\vartheta]u_{l}[\vartheta]({\bf r},t)$
(29)
$\displaystyle\quad+\sum_{l=1}^{m}c^{\prime}_{kl}[\vartheta]u^{\prime}_{l}[\vartheta]({\bf
r},t),$

where the modes $u^{\prime}_{l}[\vartheta]({\bf r},t)$ form an $m(\leq
n)-$dimensional basis of the mode Hilbert space $\mathcal{H}_{\partial}={\rm
span}\left(\\{\partial_{\vartheta}u_{k}[\vartheta]({\bf
r},t)\\}\right)\setminus\mathcal{H}_{n}$ (See App. A). Accordingly, the
diagonal block $B_{n}$ contains a mode contribution (first two terms in Eq.
(28a)) due to the portion of the derivatives
$\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)$ within the space of the
initially populated modes $\mathcal{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_{\partial}$ and $B_{\partial}^{\top}$
only contain the mode contribution due to the leakage of the derivatives
$\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)$ from $\mathcal{H}_{n}$ to
$\mathcal{H}_{\partial}$.

Using Eq. (27), we can calculate the coefficients $a_{jk}^{(l)}$ in Eq. (17c),
which result in

$a_{jk}^{(l)}=\begin{cases}\Tr\left[A_{jk}^{(l)}B_{n}\right]&1\leq j,k\leq
n\\\ \Tr\left[\tilde{A}_{jk}^{(l)\top}B_{\partial}\right]&1\leq j\leq
n,n<k\leq(n+m)\\\
\Tr\left[\tilde{A}_{jk}^{(l)}B_{\partial}^{\top}\right]&n<k\leq(n+m),1\leq
j\leq n\\\ 0&j,k>(n+m)\end{cases}$ (30)

where $\tilde{A}_{jk}^{(l)}$ are $m\times n$ blocks of the matrices
$A_{ij}^{(l)}$. Finally, using Eq. (19a) and Eq. (26b), we can write the
covariance matrix contribution to the QFI as

$\displaystyle F_{\sigma}$
$\displaystyle=\frac{1}{2}\sum_{l=0}^{3}\sum_{jk=1}^{n}\frac{\left(a_{jk}^{(l)}\right)^{2}}{\nu_{j}\nu_{k}-(-1)^{l}}$
(31)
$\displaystyle\quad+\frac{1}{2}\sum_{l=0}^{3}\sum_{j=1}^{n}\sum_{k=1}^{m}\frac{\left(a_{j,k+n}^{(l)}\right)^{2}+\left(a_{k+n,j}^{(l)}\right)^{2}}{\nu_{j}-(-1)^{l}}.$

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}[\vartheta]({\bf r},t)$. On the other hand, the second term in Eq. (31)
contains a sum over the $n$ initially populated modes $u_{n}[\vartheta]({\bf
r},t)$ and another over their $m$ orthonormalized derivatives
$u^{\prime}_{n}[\vartheta]({\bf 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 $\hat{\rho}_{\vartheta}$ in the $n$ initially populated modes
$u_{k}[\vartheta]({\bf 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}({\bf r},t)$.

###  3.3 Displacement vector contribution to the quantum Fisher information

We now move on to compute $F_{\bar{\bf x}}$, as given by Eq. (19b), which
takes into account the contribution to the QFI coming from variations of the
displacement vector $\bar{x}_{I}$. To compute this term, we need the
derivative of Eq. (24)

$\partial_{\vartheta}\bar{\bf
x}_{I}=(\partial_{\vartheta}O_{n})\bar{x}_{n}+O_{n}(\partial_{\vartheta}\bar{\bf
x}_{n}),$ (32)

and the inverse of the covariance matrix $\sigma_{I}$ that, thanks to Eq.
(22), we can write as

$\sigma_{I}^{-1}=O_{n}V_{n}^{-1}O^{\top}_{n}+O_{N-n}O^{\top}_{N-n}.$ (33)

Finally, combining Eqs. (32) and (33), and using the properties of the
matrices $O_{n}$ and $O_{N-n}$ (see App. A), we obtain

$\displaystyle F_{\bar{\bf x}}$ $\displaystyle=(\partial_{\vartheta}\bar{\bf
x}_{n})^{\top}V^{-1}_{n}(\partial_{\vartheta}\bar{\bf x}_{n})$ (34)
$\displaystyle\;+(\partial_{\vartheta}\bar{\bf
x}_{n})^{\top}V^{-1}_{n}D_{n}^{\top}\bar{\bf x}_{n}+\bar{\bf
x}^{\top}_{n}D_{n}V^{-1}_{n}(\partial_{\vartheta}\bar{\bf x}_{n})$
$\displaystyle\;+\bar{\bf
x}^{\top}_{n}\left(D_{n}V^{-1}_{n}D^{\top}_{n}+D_{\partial}D^{\top}_{\partial}\right)\bar{\bf
x}_{n}.$

Similarly to what we observed for $F_{\sigma}$, $F_{\bar{\bf x}}$ only depends
on the displacement vector $\bar{\bf x}_{n}$ in the $n$ initially populated
modes $u_{k}[\vartheta]({\bf r},t)$ and their shapes. Moreover, we note that
the first term in in Eq. (34) only depends on variations of the displacement
vector $\bar{\bf x}_{n}$, while the last term only depends on changes of the
shapes of the $n$ initially populated modes $u_{k}[\vartheta]({\bf r},t)$. On
the other hand, in the two middle terms appear both
$(\partial_{\vartheta}\bar{\bf x}_{n}[\vartheta])$ 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

####  4.1.1 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]({\bf r})=u({\bf r}-{\bf r}_{0})$ with ${\bf 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]({\bf r})$ to have a well-defined
parity, s.t. it is orthogonal to its derivative:
$(\partial_{d}u_{0}|u_{0})=0$. Under these assumptions, we have $D_{n}=0$ and
$D_{\partial}=\eta\mathds{1}_{2}$, with $\eta=||\partial_{d}u({\bf r}-{\bf
r}_{0})||$. In this context, $\eta$ quantifies the spatial extent of the beam
we want to localize, e.g. for a Gaussian mode $u({\bf r})=\exp(-|{\bf
r}|^{2}/2w^{2})/\sqrt{\pi w^{2}}$, we have $\eta^{2}=1/2w^{2}$.

This estimation problem is fully defined by the mode $u_{0}[d]({\bf r})$. As a
consequence, the mean field contribution to the QFI (34) simplifies to

$F_{\bar{\bf x}}=(\partial_{d}{\bf x}_{0})^{\top}V_{0}^{-1}(\partial_{d}{\bf
x}_{0})+\eta^{2}||{\bf x}_{0}||^{2},$ (35)

and we can write the covariance contribution to the QFI (31) as

$F_{v}=\frac{1}{2}\sum_{l=0}^{3}\left(\frac{a_{l}^{2}}{\nu_{0}^{2}-(-1)^{l}}+\frac{b_{l}^{2}}{\nu_{0}-(-1)^{l}}\right),$
(36)

where we defined the coefficients

$\displaystyle a_{l}^{2}$
$\displaystyle=\tr\left[A_{l}S_{0}^{-1}(\partial_{d}V_{0})(S_{0}^{\top})^{-1}\right]^{2},$
(37a) $\displaystyle b_{l}^{2}$
$\displaystyle=\eta^{2}\tr\left[A_{l}\left(V_{0}-\mathds{1}_{2}\right)\left(S_{0}^{\top}\right)^{-1}\right]^{2}$
(37b)
$\displaystyle\quad+\eta^{2}\tr\left[A_{l}S_{0}^{-1}\left(V_{0}-\mathds{1}_{2}\right)\right]^{2},$

with

$\displaystyle A_{0}$
$\displaystyle=i\sigma_{y}/\sqrt{2};\;A_{1}=\sigma_{z}/\sqrt{2};$ (38)
$\displaystyle A_{2}$
$\displaystyle=\mathds{1}_{2}/\sqrt{2};\;A_{3}=\sigma_{x}/\sqrt{2},$

where we recall $\sigma_{x,y,z}$ are Pauli matrices.

We can now evaluate Eqs. (35) and (36) for different states of the mode
$u_{0}[d]({\bf r})$. Let us start by considering a coherent state
$\ket{\alpha}$, defined by the complex amplitude $\alpha$ that can be
parameter dependent. Accordingly, we have $\bar{\bf
x}_{0}=2(\real[\alpha],\imaginary[\alpha])$ and $V_{0}=\mathds{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

$F_{d,{\rm coh}}=|\partial_{d}\alpha|^{2}+4\eta^{2}N_{0},$ (39)

where we introduced the mean photon number $N_{0}=|\alpha|^{2}$. The second
term in Eq. (39) presents a shot-noise 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 $\alpha$
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
$\bar{\bf x}_{0}=0$ and $V_{0}=(2N_{0}+1)\mathds{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) are $a_{2}^{2}=8(\partial_{d}N_{0})^{2}$ and
$b_{2}^{2}=16N_{0}^{2}\eta^{2}$, resulting in

$F_{d,{\rm
th}}=\frac{(\partial_{d}N_{0})^{2}}{N_{0}(N_{0}+1)}+4\eta^{2}N_{0}.$ (40)

The $4\eta^{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}=|\alpha|^{2}$ to rewrite this term in function of
the mean photon number $N_{0}$ also in the coherent case. Accordingly, we get
$|\partial_{d}\alpha|^{2}=(\partial_{d}N_{0})^{2}/N_{0}$, 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.

####  4.1.2 Populating 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
$\partial_{d}u_{0}[d]({\bf r})$. In the following, we will see how our
formalism recovers this result, to extend it to different states of the mode
$u_{0}[d]({\bf r})$ and to take into account losses in the squeezed derivative
mode.

When populating the derivative mode, the mode Hilbert space $\mathcal{H}_{n}$,
as introduced in Sec. 3, is spanned by $u_{0}[d]({\bf r})=u({\bf r}-{\bf
r}_{0})$ and its normalized derivative $u_{1}[d]({\bf
r})=\partial_{d}u_{0}({\bf r}-{\bf r}_{0})/\eta$. On the other hand, the mode
Hilbert space $\mathcal{H}_{\partial}$ (see Sec. 3) only contains the second
derivative mode $u_{2}[d]({\bf r})=\left(\partial_{d}u_{1}[d]({\bf r})-\xi
u_{0}[d]({\bf r})\right)/\zeta$, with $\xi=(\partial_{d}u_{1}|u_{0})$ and
$\zeta=||\partial_{d}u_{1}[d]({\bf r})-\xi u_{0}[d]({\bf r})||$. Accordingly,
we have

$\displaystyle D_{n}$ $\displaystyle=\begin{pmatrix}0&\eta\mathds{1}_{2}\\\
\xi\mathds{1}_{2}&0\end{pmatrix},\quad D_{\partial}=\begin{pmatrix}0\\\
\zeta\mathds{1}_{2}\end{pmatrix}.$ (41)

Furthermore, we assume that the derivative mode $u_{1}[d]({\bf r})$ has no
mean field so that the mean field vector can be written as $\bar{\bf
x}^{\top}=(q_{0},p_{0},0,0)$. Therefore, the mean field term of the QFI (34)
results in

$F_{\bar{\bf x}}={\bf y}^{T}V^{-1}{\bf y}=(\partial_{d}\bar{\bf
x}_{0}^{\top},\eta\bar{\bf
x}_{0}^{\top})V^{-1}\begin{pmatrix}\partial_{d}\bar{\bf x}_{0}\\\ \eta\bar{\bf
x}_{0}\\\ \end{pmatrix}.$ (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 $V_{v}^{-1}$

$F_{\bar{\bf x}}=||{\bf y}||^{2}(V^{-1}_{v})_{0,0},$ (43)

where $(V_{v})_{0,0}$ is the variance of the $q-$quadrature of mode

$\displaystyle v_{0}[d]({\bf r},t)$
$\displaystyle=\frac{(\partial_{d}q_{0}+i\partial_{d}p_{0})}{||{\bf
y}||}u_{0}[d]({\bf r},t)$ (44)
$\displaystyle\quad+\frac{\eta(q_{0}+ip_{0})}{||{\bf y}||}u_{1}[d]({\bf r},t)$

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]({\bf 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.
$\partial_{d}q_{0}=\partial_{d}p_{0}=0$, the mode $v_{0}[d]({\bf r})$ equals
(up to a global phase) the derivative mode $u_{1}[d]({\bf r})$, which is
orthogonal to the mode $u_{0}[d]({\bf 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]({\bf r})$ which is partially
overlapping with $u_{0}[d]({\bf r})$.

Let us now discuss how the covariance matrix term of the QFI is modified by
population in the derivative mode $u_{1}[d]({\bf r})$. For simplicity, we will
focus on the case where the population of mode $u_{1}[d]({\bf r})$ is fully
uncorrelated with that of mode $u_{0}[d]({\bf r})$, therefore, the covariance
matrix takes the block diagonal form

$V=\begin{pmatrix}V_{0}&0\\\ 0&V_{1}\end{pmatrix}.$ (45)

Under these assumptions, the covariance matrix contribution to the QFI (31)
takes the form

$\displaystyle F_{v}$
$\displaystyle=\frac{1}{2}\sum_{l=0}^{3}\left(\frac{a_{l}^{2}}{\nu_{0}^{2}-(-1)^{l}}\right.$
(46)
$\displaystyle\quad+\left.\frac{\tilde{b}_{l}^{2}}{\nu_{0}\nu_{1}-(-1)^{l}}+\frac{d_{l}^{2}}{\nu_{1}-(-1)^{l}}\right),$

where we have defined the coefficients

$\displaystyle a_{l}^{2}$
$\displaystyle=\tr\left[A_{l}S_{0}^{-1}(\partial_{\vartheta}V_{0})(S_{0}^{\top})^{-1}\right]^{2},$
(47a) $\displaystyle\tilde{b}_{l}^{2}$
$\displaystyle=\left(\eta\tr\left[A_{l}S_{1}^{-1}\left(V_{0}-\mathds{1}_{2}\right)\left(S_{0}^{\top}\right)^{-1}\right]\right.$
(47b)
$\displaystyle\quad+\left.\xi\tr\left[A_{l}S_{1}^{-1}\left(V_{1}-\mathds{1}_{2}\right)\left(S_{0}^{\top}\right)^{-1}\right]\right)^{2}$
$\displaystyle\quad+\left(\eta\tr\left[A_{l}S_{0}^{-1}\left(V_{0}-\mathds{1}_{2}\right)(S_{1}^{-1})^{\top}\right]^{2}\right.$
$\displaystyle\quad+\left.\xi\tr\left[A_{l}S_{0}^{-1}\left(V_{1}-\mathds{1}_{2}\right)(S_{1}^{-1})^{\top}\right]\right)^{2},$
$\displaystyle d_{l}^{2}$
$\displaystyle=\zeta^{2}\tr\left[A_{l}S_{1}^{-1}\left(V_{1}-\mathds{1}_{2}\right)\right]^{2}$
(47c)
$\displaystyle\quad+\zeta^{2}\tr\left[A_{l}\left(V_{1}-\mathds{1}_{2}\right)(S_{1}^{-1})^{\top}\right]^{2}.$

To further illustrate how to use Eq. (46) in practice, let us now consider the
localisation of a thermal beam $u_{0}[d]({\bf r})$ aided by a squeezed vacuum
state in the derivative mode $u_{1}[d]({\bf r})$. Accordingly, we have

$\displaystyle V_{0}=(2N_{0}+1)\mathds{1}_{2},\quad
V_{1}=\begin{pmatrix}e^{-2r}&0\\\ 0&e^{2r}\end{pmatrix},$ (48)

which corresponds to

$\displaystyle S_{0}=\mathds{1}_{2},\quad{\rm with}\quad
S_{1}=\begin{pmatrix}e^{-r}&0\\\ 0&e^{r}\end{pmatrix},$ (49)

with $\nu_{0}=2N_{0}+1$ and $\nu_{1}=1$. Under these assumptions, the only
nonzero coefficients are $a^{2}_{2}=2(2N_{0}^{\prime}+1)^{2}$,
$\tilde{b}^{2}_{1}=16(N_{0}\eta-\xi)^{2}\sinh^{2}r$,
$\tilde{b}^{2}_{2}=16N_{0}^{2}\eta^{2}\cosh^{2}r$, and
$d^{2}_{1}=2\zeta^{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)

$\displaystyle F_{d,{\rm th-sq}}$
$\displaystyle=\frac{(\partial_{d}N_{0})^{2}}{N_{0}(N_{0}+1)}+\frac{4(N_{0}\eta-\xi)^{2}N_{1}}{N_{0}+1}$
(50) $\displaystyle\quad+4N_{0}\eta^{2}(N_{1}+1)+4\zeta^{2}N_{1},$

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}\gg N_{1}\gg 1$ limit,
where we have $F_{d,{\rm th-sq}}\sim 4N_{0}(2N_{1}+1)\eta^{2}\sim
2N_{0}\eta^{2}e^{2r}\sim e^{2r}F_{d,{\rm th}}/2$.

It is interesting to compare this result, with the quantum enhancement
achievable with a coherent state in mode $u_{0}[d]({\bf r})$. For simplicity,
let us consider the case where the mean field $\bar{\bf x}_{0}$ 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

$F_{d,{\rm coh-sq}}=4N_{0}\eta^{2}e^{2r}+4(\xi^{2}+\zeta^{2})N_{1}.$ (51)

The second term is negligible for $N_{0}\gg N_{1}$, and we obtain $F_{d,{\rm
coh-sq}}\sim e^{2r}F_{d,{\rm 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]({\bf r})$ when the mode
$u_{0}[d]({\bf 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. (43) and (50)
the average QFI of a coherent state in mode $u_{0}[d]({\bf r})$ combined with
a squeezed vacuum state in the normalized derivative mode $u_{1}[d]({\bf r})$
with random, uniformly distributed squeezing directions

$\displaystyle F_{d,{\rm avg}}$
$\displaystyle=\frac{4N_{0}\eta^{2}}{2\pi}\int(e^{2r}\cos^{2}\phi+e^{-2r}\sin^{2}\phi)d\phi$
$\displaystyle\quad+4(\xi^{2}+\zeta^{2})N_{1}$ (52)
$\displaystyle=4N_{0}\eta^{2}\cosh 2r+4(\xi^{2}+\zeta^{2})N_{1}.$

While in general, the convexity of the QFI ensures $F_{d,{\rm avg}}\geq
F_{d,{\rm th-sq}}$, for $N_{0}\gg N_{1}\gg 1$ we have $F_{d,{\rm avg}}\sim
2N_{0}\eta^{2}e^{2r}\sim F_{d,{\rm 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]({\bf 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]({\bf r})$ and the derivative mode $u_{1}[d]({\bf r})$
populated with an arbitrary zero-mean Gaussian state, i.e. a squeezed thermal
state. Accordingly, we have

$V_{1}=(2N_{T}+1)\begin{pmatrix}e^{-2r}&0\\\ 0&e^{2r}\end{pmatrix},$ (53)

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 $\nu_{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)

$\displaystyle F_{d,{\rm
th-g}}=\frac{(N_{0}^{\prime})^{2}}{N_{0}(N_{0}+1)}+\frac{4N_{S}\left(\xi(N_{T}+1)-\eta
N_{0}\right)^{2}}{2N_{0}N_{T}+N_{0}+N_{T}+1}$
$\displaystyle+\frac{4(N_{S}+1)\left(\eta^{2}N_{0}^{2}+N_{0}N_{T}\left(2\eta\xi+\zeta^{2}(2N_{T}+1)\right)\right)}{2N_{0}N_{T}+N_{0}+N_{T}}$
$\displaystyle+\frac{4N_{T}^{2}(N_{S}+1)\left(\zeta^{2}+\xi^{2}\right)}{2N_{0}N_{T}+N_{0}+N_{T}}+4N_{S}\zeta^{2}(N_{T}+1),$
(54)

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

$\displaystyle F_{d,{\rm th-th}}$
$\displaystyle=\frac{(N_{0}^{\prime})^{2}}{N_{0}(N_{0}+1)}+4\frac{(\eta
N_{0}+\xi N_{1})^{2}}{2N_{0}N_{1}+N_{0}+N_{1}}$
$\displaystyle\quad+4\zeta^{2}N_{1},$ (55)

and it is not hard to show that $F_{d,{\rm th-th}}$ (55) is always smaller
than the $F_{d,{\rm 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,{\rm th-th}}$ (55) is even
smaller than that for an unpopulated derivative mode $F_{d,{\rm th}}$ (40). To
better illustrate this interplay between squeezing and thermal noise, we
introduce the following parametrisation

$\displaystyle N_{S}=\chi N_{1}\quad{\rm and}\quad
N_{T}=\frac{(1-\chi)N_{1}}{1+2\chi N_{1}},$ (56)

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 $\chi=1$ the derivative mode is purely squeezed, while for
$\chi=0$ it is purely thermal, so that we can refer to $\chi$ as the squeezing
fraction. If we further assume that the beam we are trying to localize is
Gaussian, i.e. $u({\bf r})=\exp(-|{\bf r}|^{2}/2w^{2})/\sqrt{\pi w^{2}}$, we
show that for $\chi\geq 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 $\chi<1/2$ and small $N_{1}$ we obtain a worse
sensitivity compared to that when the derivative mode is in vacuum.

Figure 1: QFI (normalized by its maximum value when $u_{1}[d]({\bf 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 $\chi$.

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]({\bf r})$, and an arbitrary zero-mean Gaussian state of the
derivative mode $u_{1}[d]({\bf r})$, the QFI maintains the form (54) even
after losses. In fact, it is sufficient to perform the following substitutions

$\displaystyle N_{0}$ $\displaystyle=N_{0}^{\rm in}\kappa_{0}$ (57)
$\displaystyle N_{T}$
$\displaystyle=\frac{1}{2}\left(\left[\left(2\kappa_{1}(2N_{S}^{\rm
in}N_{T}^{\rm in}+N_{T}^{\rm in}+N_{S}^{\rm in})+1\right)^{2}\right.\right.$
(58) $\displaystyle\quad\left.\left.-4N_{S}^{\rm in}(2N_{T}^{\rm
in}+1)^{2}(N_{S}^{\rm in}+1)\right]^{1/2}-1\right)$ $\displaystyle\sinh(2r)$
$\displaystyle=\kappa_{1}\frac{2N_{T}^{\rm in}+1}{2N_{T}+1}\sinh(2r^{\rm
in}),$ (59)

where $\kappa_{0}$ and $\kappa_{1}$ are the attenuation coefficients of the
two modes $u_{0}[d]({\bf r})$, and $u_{1}[d]({\bf r})$, respectively; while
$N_{0}^{\rm in}$, $N_{T}^{\rm in}$ and $N_{S}^{\rm in}=\sinh^{2}r^{\rm in}$
are the populations of the mode $u_{0}[d]({\bf r})$, and the thermal and
squeezing components of the population of the mode $u_{1}[d]({\bf r})$,
respectively. Finally, in some applications, the attenuation coefficients
$\kappa_{0}$ and $\kappa_{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

$\displaystyle\qquad\frac{1}{2}\sum_{l=0}^{3}\frac{\tr\left[A_{l}S_{1}^{-1}\partial_{d}V_{1}\left(S_{1}^{-1}\right)^{\top}\right]^{2}}{\nu_{1}^{2}-(-1)^{l}}$
(60) $\displaystyle\quad=\frac{\left((2N_{T}+1)(\partial_{d}r)\cosh
r-2(\partial_{d}N_{T})\sinh r\right)^{2}}{4N_{T}(N_{T}+1)+2}$
$\displaystyle\qquad+\frac{\left((2N_{T}+1)(\partial_{d}r)\sinh
r-2(\partial_{d}N_{T})\cosh r\right)^{2}}{4N_{T}(N_{T}+1)}.$

###  4.2 Temporal separation between pulses

As a second example, we consider the estimation of the time delay $\tau$
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 $\tau$ 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-\tau/2)$ and $u(t+\tau/2)$ (see Fig. 2)

$\delta=\int u(t-\tau/2)u(t+\tau/2)dt.$ (61)

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

$\displaystyle u_{0}[\tau](t)$
$\displaystyle=\frac{u(t-\tau/2)+u(t+\tau/2)}{\sqrt{2(1+\delta)}},$ (62a)
$\displaystyle v_{0}[\tau](t)$
$\displaystyle=\frac{u(t-\tau/2)-u(t+\tau/2)}{\sqrt{2(1-\delta)}}.$ (62b)

We are interested in computing the QFI for the estimation of $\tau$, 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):

$\displaystyle u_{1}[\tau](t)$
$\displaystyle=\partial_{\tau}u_{0}[\tau](t)/\eta_{u}$ (63a) $\displaystyle
v_{1}[\tau](t)$ $\displaystyle=\partial_{\tau}v_{0}[\tau](t)/\eta_{v}$ (63b)
$\displaystyle u_{2}[\tau](t)$
$\displaystyle=(\partial_{\tau}u_{1}[\tau](t)-\xi_{u}u_{0}[\tau](t))/\zeta_{u}$
(63c) $\displaystyle v_{2}[\tau](t)$
$\displaystyle=(\partial_{\tau}v_{1}[\tau](t)-\xi_{v}v_{0}[\tau](t))/\zeta_{v},$
(63d)

where $\eta_{u}=||\partial_{\tau}u_{0}[\tau](t)||$,
$\eta_{v}=||\partial_{\tau}v_{0}[\tau](t)||$,
$\xi_{u}=(\partial_{\tau}u_{1}|u_{0})$,
$\xi_{v}=(\partial_{\tau}v_{1}|v_{0})$,
$\zeta_{u}=||\partial_{\tau}u_{1}[\tau](t)-\xi_{u}u_{0}[\tau](t)||$ and
$\zeta_{u}=||\partial_{\tau}v_{1}[\tau](t)-\xi_{v}v_{0}[\tau](t)||$. The
shapes of the modes $u_{i}[\tau](t)$ and $v_{i}[\tau](t)$ for the specific
case of Gaussian pulses $u(t)=e^{-t^{2}/2w^{2}}/(\pi w^{2})^{1/4}$ are
presented in Fig. 2. Using the modes (62) and (63), we can express the
matrices $D_{n}$ and $D_{\partial}$ (see Sec. 3.2) as

$\displaystyle D_{n}$ $\displaystyle=\begin{pmatrix}0&D_{\eta}\\\
D_{\xi}&0\end{pmatrix},\quad
D_{\partial}=\left(0,0,D_{\zeta}\right)^{\top},\quad{\rm with}$ $\displaystyle
D_{\eta}$ $\displaystyle=\begin{pmatrix}\eta_{u}\mathds{1}_{2}&0\\\
0&\eta_{v}\mathds{1}_{2}\end{pmatrix},\;D_{\xi}=\begin{pmatrix}\xi_{u}\mathds{1}_{2}&0\\\
0&\xi_{v}\mathds{1}_{2}\end{pmatrix},$ $\displaystyle D_{\zeta}$
$\displaystyle=\begin{pmatrix}\zeta_{u}\mathds{1}_{2}&\zeta_{v}\mathds{1}_{2}\end{pmatrix}.$
(64)

Figure 2: (top) Two Gaussian pulses temporally separated by $\tau=w$, with
their overlap $\delta$ 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 $\tau$. All modes have definite parity, in
particular the modes $u_{i}[\tau](t)$ (left - red) are even functions, while
the modes $v_{i}[\tau](t)$ (right - blue) are odd functions.

We have now specified all the mode-related quantities needed to compute the
QFI for the estimation of the temporal separation $\tau$ between pulses. To
proceed, we now make some further assumptions on the quantum state of the
pulses. In particular, we consider the modes $u_{0}[\tau](t)$ and
$v_{0}[\tau](t)$ to be in a general Gaussian state, and we allow for auxiliary
population of the orthogonalized first derivative modes $u_{1}[\tau](t)$ and
$v_{1}[\tau](t)$ with no mean field111Note that we allow for correlations
(classical or quantum) between the modes $u_{0}[\tau](t)$ and
$v_{0}[\tau](t)$. To focus on the role of squeezing in the derivative modes,
in the examples contained this paper, we do not consider such correlations.
However, we studied their role in diffraction-limited imaging in Sorelli  _et
al._ (2022).. Accordingly, we can write the mean field vector as $\bar{\bf
x}^{\top}=(\bar{\bf x}^{\top}_{0},0)$, with $\bar{\bf
x}_{0}=(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

$\displaystyle F_{\bar{\bf x}}=(\partial_{\tau}\bar{\bf x}_{0},D_{n}\bar{\bf
x}_{0})^{\top}V^{-1}(\partial_{\tau}\bar{\bf x}_{0},D_{n}\bar{\bf x}_{0})$
(65)

Similarly to what we discussed in Sec. 4.1, for every state with a nonzero
mean field, i.e. with $||\bar{\bf x}_{0}||\neq 0$, there always exists a mode
basis where $(\partial_{\tau}\bar{\bf x}_{0},D_{n}\bar{\bf x}_{0})$ has only
one nonzero component. In other words, there always exists an orthogonal
transformation $O$ s.t. $O(\partial_{\tau}\bar{\bf x}_{0},D_{n}\bar{\bf
x}_{0})=(\sqrt{||\partial_{\tau}\bar{\bf x}_{0}||^{2}+||D_{\eta}\bar{\bf
x}_{0}||^{2}},0,0,0)$. Accordingly the QFI only depends on the inverse
covariance matrix element $\left(O^{T}V^{-1}O\right)_{0,0}$. 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}[\tau](t)$ and $v_{1}[\tau](t)$ is uncorrelated
with that of the symmetric and antisymmetric superpositions $u_{0}[\tau](t)$
and $v_{0}[\tau](t)$ of the pulses we want to separate, so that we can write
the covariance matrix in block diagonal form

$V=\begin{pmatrix}V_{0}&0\\\ 0&V_{1}\end{pmatrix}.$ (66)

Under these assumptions, the covariance matrix contribution to the QFI takes
the form

$\displaystyle F_{V}=$
$\displaystyle\frac{1}{2}\sum_{l=0}^{3}\sum_{jk=0,1}\left(\frac{\left(a_{l}^{jk}\right)^{2}}{\nu_{0}^{j}\nu_{0}^{k}-(-1)^{l}}+\frac{\left(b_{l}^{jk}\right)^{2}}{\nu_{0}^{j}\nu_{1}^{k}-(-1)^{l}}\right.$
(67)
$\displaystyle\left.\quad+\frac{\left(c_{l}^{jk}\right)^{2}}{\nu_{1}^{j}\nu_{1}^{k}-(-1)^{l}}+\frac{\left(d_{l}^{jk}\right)^{2}}{\nu_{1}^{j}-(-1)^{l}}\right),$

where we introduced the coefficients

$\displaystyle\left(a_{l}^{jk}\right)^{2}$
$\displaystyle=\left(\tr\left[A^{(jk)}_{l}\partial_{\tau}V_{0}\right]\right)^{2},$
(68a) $\displaystyle\left(b_{l}^{jk}\right)^{2}$
$\displaystyle=\left(\tr\left\\{A_{l}^{(jk)}S_{0}^{-1}\left[D_{\xi}^{T}(V_{1}-\mathds{1}_{4})\right.\right.\right.$
(68b)
$\displaystyle\quad+\left.\left.\left.(V_{0}-\mathds{1}_{4})D_{\eta}\right]\left(S_{1}^{-1}\right)^{\top}\right\\}\right)^{2}$
$\displaystyle\quad+\left(\tr\left\\{A_{l}^{(jk)}S_{1}^{-1}\left[D_{\xi}^{T}(V_{1}-\mathds{1}_{4})\right.\right.\right.$
$\displaystyle\quad+\left.\left.\left.(V_{0}-\mathds{1}_{4})D_{\eta}\right]\left(S_{0}^{-1}\right)^{\top}\right\\}\right)^{2},$
$\displaystyle\left(c_{l}^{jk}\right)^{2}$
$\displaystyle=\left(\tr\left[A^{(jk)}_{l}\partial_{\tau}V_{1}\right]\right)^{2},$
(68c) $\displaystyle\left(d_{l}^{jk}\right)^{2}$
$\displaystyle=\left(\tr\left\\{A_{l}^{(jk)}D_{\zeta}(V_{1}-\mathds{1}_{4})\left(S_{1}^{-1}\right)^{\top}\right\\}\right)^{2},$
(68d)

with the matrices $A_{l}^{(jk)}$ as defined in Sec. 2.3 and $\nu_{0}^{j}$
($\nu_{1}^{j}$) 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_{l}^{jk}$, describes the
contribution of the population of the symmetric $u_{0}[\tau](t)$ and
antisymmetric $v_{0}[\tau](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. $\partial_{\tau}V_{0}\neq 0$. Similarly, the
third group of addends, depending on the coefficients $c_{l}^{jk}$, takes into
account the population of the derivative modes $u_{1}[\tau](t)$ and
$v_{1}[\tau](t)$, and is nonzero if and only if $\partial_{\tau}V_{1}\neq 0$.
The second group of addends, containing the coefficients $b_{l}^{jk}$, takes
into account how variations of the temporal separation $\tau$ leads to
coupling between the modes $u_{0}[\tau](t)$, $v_{0}[\tau](t)$ and their
derivative $u_{1}[\tau](t)$, $v_{1}[\tau](t)$. Finally, the addends containing
the coefficients $d_{l}^{jk}$ take into account how due to variations of
$\tau$ the derivative modes $u_{1}[\tau](t)$, $v_{1}[\tau](t)$ couple to the
second derivative modes $u_{2}[\tau](t)$ and $v_{2}[\tau](t)$.

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

$\displaystyle I(t)$
$\displaystyle=\langle\hat{E}^{\dagger}(t)\hat{E}(t)\rangle$
$\displaystyle=N_{0}(|u(t-\tau/2)|^{2}+|u(t+\tau/2)|^{2}),$ (69)

where we introduced the mean number of photons per pulse $N_{0}$, and the
electric field operator
$\hat{E}(t)=\sum_{j}\left(\hat{a}_{j}u_{j}[\tau](t)+\hat{b}_{j}v_{j}[\tau](t)\right)$,
with $\hat{a}_{j}$ and $\hat{b}_{j}$ the annihilation operators associated
with the even and odd modes $u_{j}[\tau](t)$ and $v_{j}[\tau](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}[\tau](t)$ and $v_{0}[\tau](t)$, with mean photon numbers
$N_{u}=N_{0}(1+\delta)$ and $N_{v}=N_{0}(1-\delta)$, respectively. Such a
state has no mean field $\bar{\bf x}_{0}=0$, so that its QFI is fully
determined by Eq. (67), and has a covariance matrix

$V_{0}=\begin{pmatrix}\left(2N_{0}(1+\delta)+1\right)\mathds{1}_{2}&0\\\
0&\left(2N_{0}(1-\delta)+1\right)\mathds{1}_{2}\end{pmatrix}.$ (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}[\tau](t)$
and $v_{1}[\tau](t)$ to be populated by two independent, equally-squeezed
vacuum states, described by the covariance matrix

$V_{1}=\begin{pmatrix}e^{-2r}&0&0&0\\\ 0&e^{2r}&0&0\\\ 0&0&e^{-2r}&0\\\
0&0&0&e^{2r}\\\ \end{pmatrix}.$ (71)

For such a quantum state, the QFI (67) takes the form (see App. C for the
explicit calculation of the coefficients (68))

$\displaystyle F_{\tau,{\rm th-sq}}$
$\displaystyle=\frac{2N_{0}\left[1+N_{0}(1+\delta^{2})\right]\left(\partial_{\tau}\delta\right)^{2}}{(1-\delta^{2})\left[(1+N_{0})^{2}-(N_{0}\delta)^{2}\right]}$
(72) $\displaystyle\quad+2(\zeta_{u}^{2}+\zeta_{v}^{2})\sinh^{2}r$
$\displaystyle\quad+4N_{0}\left[\eta^{2}_{u}(1+\delta)+\eta_{v}^{2}(1-\delta)\right]\cosh^{2}r$
$\displaystyle\quad+4\frac{\left(N_{0}(1+\delta)\eta_{u}-\xi_{u}\right)^{2}}{1+N_{0}(1+\delta)}\sinh^{2}r$
$\displaystyle\quad+4\frac{\left(N_{0}(1-\delta)\eta_{v}-\xi_{v}\right)^{2}}{1+N_{0}(1-\delta)}\sinh^{2}r.$

The behavior of the QFI (72), for Gaussian pulses, is plotted as red lines in
Fig. 3.

Figure 3: QFI (Normalized to its maximum value for thermal states in modes
$u_{0}[\tau](t)$ and $v_{0}[\tau](t)$) for the estimation of the temporal
separation $\tau$ between two thermal (red) or coherent pulses, either in
phase (blue) or out of phase (green), as a function of the temporal separation
$\tau$ 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}}/(\pi w^{2})^{1/4}$ for all panels.

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 $\phi$ between the coherent pulses

$\displaystyle I(t)$
$\displaystyle=\langle\hat{E}^{\dagger}(t)\hat{E}(t)\rangle$
$\displaystyle=N_{0}(|u(t-\tau/2)|^{2}+|u(t+\tau/2)|^{2})$ (73)
$\displaystyle\quad+2N_{0}u(t-\tau/2)u(t+\tau/2)\cos\phi.$

Such an intensity distribution can be obtained by populating the modes
$u_{0}[\tau](t)$ and $v_{0}[\tau](t)$ with coherent states, whose covariance
matrix is the identity $V_{0}=\mathds{1}_{4}$, and whose mean field is given
by $\bar{\bf x}_{0}=(\bar{\bf x}_{u},\bar{\bf x}_{v})^{\top}$, with

$\displaystyle\bar{\bf x}_{u}$
$\displaystyle=\sqrt{2N_{0}(1+\delta)}(1+\cos\phi,\sin\phi)$ (74a)
$\displaystyle\bar{\bf x}_{v}$
$\displaystyle=\sqrt{2N_{0}(1-\delta)}(1-\cos\phi,-\sin\phi).$ (74b)

In particular, from Eqs. (74), we can see that for in-phase ($\phi=0$)
coherent pulses, the mean field is fully determined by the $q$ quadrature of
mode $u_{0}[\tau](t)$. Similarly, when the two coherent pulses are out of
phase ($\phi=\pi$) the mean field is fully determined by the $q$ quadrature of
mode $v_{0}[\tau](t)$. As we did for thermal sources, we are going to consider
the derivative modes $u_{1}[\tau](t)$ and $v_{1}[\tau](t)$ by two independent
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

$F_{\tau,{\rm coh-sq}}=F_{\bar{\bf x},{\rm coh-sq}}+\left.F_{\tau,{\rm th-
sq}}\right|_{N_{0}=0},$ (75)

where the displacement term $F_{\bar{\bf x},{\rm coh-sq}}$ can be computed
from Eq. (65), and reads

$\displaystyle F_{\bar{\bf x},{\rm coh-sq}}$
$\displaystyle=(\partial_{\tau}\bar{\bf
x}_{0})^{\top}V_{0}^{-1}(\partial_{\tau}\bar{\bf x}_{0})+\bar{\bf
x}_{0}^{\top}D_{\eta}V_{1}^{-1}D_{\eta}\bar{\bf x}_{0}$
$\displaystyle=\frac{2N_{0}(1-\delta\cos\phi)(\partial_{\tau}\delta)^{2}}{1-\delta^{2}}$
(76) $\displaystyle\quad+2N_{0}e^{2r}\eta_{u}^{2}(1+\delta)(1+\cos\phi)^{2}$
$\displaystyle\quad+2N_{0}e^{2r}\eta_{v}^{2}(1-\delta)(1-\cos\phi)^{2}$
$\displaystyle\quad+2N_{0}e^{-2r}\sin^{2}\phi[\eta_{u}(1+\delta)+\eta_{v}(1-\delta)].$

We can see that, when the two coherent pulses are either in phase ($\phi=0$)
or out of phase ($\phi=\pi$), the last line in Eq. (76) vanishes and only the
squeezing-enhanced 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}[\tau](t)$ and $v_{1}[\tau](t)$, and for
$\phi=0,\pi$ the mean field (74) has a vanishing $p$ quadrature. The QFI for
the separation $\tau$ 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 $\tau\to 0$. In this regime (see App. B), we have
$(\partial_{\tau}\delta)^{2}/(1-\delta^{2})\sim(\Delta k)^{2}$, with

$(\Delta k)^{2}=\int[\partial_{t}u(t)]^{2}dt,$ (77)

and $\eta_{u}\sim\eta_{v}\sim\xi_{u}\sim\xi_{v}\sim\zeta_{u}\sim\zeta_{v}\sim
0$, which implies

$\displaystyle F_{\tau,{\rm th-sq}}$ $\displaystyle\xrightarrow{\tau\to
0}2N_{0}(\Delta k)^{2},$ (78a) $\displaystyle F_{\tau,{\rm coh-sq}}$
$\displaystyle\xrightarrow{\tau\to 0}2N_{0}(\Delta k)^{2}(1+\cos\phi).$ (78b)

Accordingly, independently of the squeezing value $r$, the QFI for in phase
coherent pulses vanishes for $\tau\to 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 $(\Delta k)^{2}=1/2w^{2}$).

To better understand this behavior, let us recall that the quantum state of
the finite overlap $\delta$ between the two pulses induces a $\tau-$dependent
population of the symmetric and antisymmetric modes $u_{0}[\tau](t)$ and
$v_{0}[\tau](t)$. It is this dependence on temporal separation, which enters
the QFI through $\partial_{\tau}\bar{\bf x}_{0}$ (in the coherent case) and
$\partial_{\tau}V_{0}$ (in the incoherent case), that dominates the QFI
behavior for $\tau\to 0$. This implies that the population of the derivative
modes $u_{1}[\tau](t)$ and $v_{1}[\tau](t)$, and in the particular the
squeezing thereof, has no impact on the $\tau\to 0$ behavior of the QFI.

On the contrary, for separations much larger than the pulses’ width, i.e.
$\tau\Delta k\gg 1$, the overlap $\delta$ tends to zero, and the populations
of the modes $u_{0}[\tau](t)$ and $v_{0}[\tau](t)$ become parameter
independent. The QFI is then dominated by the noise in derivative modes
$u_{1}[\tau](t)$ and $v_{1}[\tau](t)$. In particular, we have

$\displaystyle F_{\tau,{\rm th-sq}}\xrightarrow{\tau\to\infty}2(\Delta
k)^{2}\left(2\sinh^{2}r+N_{0}\cosh 2r\right),$ (79a) $\displaystyle
F_{\tau,{\rm coh-sq}}|_{\phi=0,\pi}\xrightarrow{\tau\to\infty}2(\Delta
k)^{2}\left(2\sinh^{2}r+N_{0}e^{2r}\right).$ (79b)

Accordingly, for large temporal separations $\tau$ 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 $\tau$ in the coherent
case is at most a factor two larger than that in the thermal one.

##  5 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 $\sim 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 $\tau$.
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).

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 carried out during the tenure of an ERCIM ‘Alain
Bensoussan’ Fellowship Programme. This work was funded by
MCIN/AEI/10.13039/501100011033 and 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 - QUANTUM 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+).

##  Appendix 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}[\vartheta]({\bf r},t)\\}$ and their derivatives
$\\{\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)\\}$. Let us start by
recalling the following Hilbert space definitions:

$\displaystyle\mathcal{H}_{n}$ $\displaystyle={\rm
span}\left(\\{u_{k}[\vartheta]({\bf r},t)\\}\right),$ (80)
$\displaystyle\mathcal{H}_{\partial}$ $\displaystyle={\rm
span}\left(\\{\partial_{\vartheta}u_{k}[\vartheta]({\bf
r},t)\\}\right)\setminus\mathcal{H}_{n}.$ (81)

We now assume that $m$ is the number of derivatives
$\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)$ that are linearly
independent from the $n$ initially populated modes, i.e. ${\rm
dim}(\mathcal{H}_{\partial})=m$. Accordingly, up to a reordering of the basis
$\\{\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)\\}$, we can always
construct a basis of $\mathcal{H}_{\partial}$ using the orthonormalized
version $u_{k}^{\prime}[\vartheta]({\bf r},t)$ of the derivatives of the first
$m$ initially populated modes $\partial_{\vartheta}u_{k}[\vartheta]({\bf
r},t)$.

We can now choose the modes $u_{k}^{\prime}[\vartheta]({\bf 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

$O_{N-n}=\begin{pmatrix}O_{\partial}&O_{N-n-m}\end{pmatrix},$ (82)

where $O_{\partial}$ and $O_{N-n-m}$ are matrices of dimensions $2N\times 2m$
and $2N\times 2(N-n-m)$, respectively. From the orthogonality of $O$, we can
obtain the following relations:

$\displaystyle O_{k}^{\top}O_{l}$ $\displaystyle=\delta_{kl}\mathds{1}_{{\rm
dim}(O_{k})},$ (83a) $\displaystyle\sum_{k}O_{k}O_{k}^{\top}$
$\displaystyle=\mathds{1}_{2N},$ (83b)

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\times 2$ blocks containing the mode
overlaps. Accordingly, it is sufficient to specify the derivative of the $kl$
block of $O_{n}$, which reads

$\partial_{\vartheta}\left(O_{n}\right)_{kl}=\begin{pmatrix}\real\left[(v_{l}|\partial_{\vartheta}u_{k}[\vartheta])\right]&-\imaginary\left[(v_{l}|\partial_{\vartheta}u_{k}[\vartheta])\right]\\\
\imaginary\left[(v_{l}|\partial_{\vartheta}u_{k}[\vartheta])\right]&\real\left[(v_{l}|\partial_{\vartheta}u_{k}[\vartheta])\right]\end{pmatrix}$
(84)

Combining Eqs. (84) and (29), we obtain

$\partial_{\vartheta}O_{n}=O_{n}D^{\top}_{n}+O_{\partial}D^{\top}_{\partial},$
(85)

where $D_{n}$ and $D_{\partial}$ are a $2n\times 2n$ and a $2n\times 2m$
matrices, respectively. Their $kl$ blocks are given by

$\displaystyle\left(D_{n}\right)_{kl}=\begin{pmatrix}\real\left(c_{kl}[\vartheta]\right)&-\imaginary\left(c_{kl}[\vartheta]\right)\\\
\imaginary\left(c_{kl}[\vartheta]\right)&\real\left(c_{kl}[\vartheta]\right)\end{pmatrix},$
(86)
$\displaystyle\left(D_{\partial}\right)_{kl}=\begin{pmatrix}\real\left(c^{\prime}_{kl}[\vartheta]\right)&-\imaginary\left(c^{\prime}_{kl}[\vartheta]\right)\\\
\imaginary\left(c^{\prime}_{kl}[\vartheta]\right)&\real\left(c^{\prime}_{kl}[\vartheta]\right)\end{pmatrix}.$
(87)

Using Eqs. (85) and (83), we then obtain

$\displaystyle O^{\top}_{n}(\partial_{\vartheta}O_{n})$
$\displaystyle=D_{n}^{\top},$ (88a) $\displaystyle
O^{\top}_{(N-n)}(\partial_{\vartheta}O_{n})$
$\displaystyle=(D_{\partial},0)^{\top}.$ (88b)

Let us conclude this appendix with few words on the coefficients
$c_{kl}[\vartheta]$ and $c^{\prime}_{kl}[\vartheta]$. The first are simply
given by the overlaps of the initially populated modes with their derivatives
$c_{kl}[\vartheta]=(u_{l}[\vartheta]|\partial_{\vartheta}u_{k}[\vartheta])$.
On the other hand, there exist several orthonormalization methods that can be
used to construct the modes $u_{k}^{\prime}[\vartheta]({\bf r},t)$, leading to
different expressions for the coefficients $c^{\prime}_{kl}[\vartheta]$. For
example, using the Gram-Schmidt procedure, the modes
$u_{k}^{\prime}[\vartheta]({\bf r},t)$ can be constructed iteratively as
$u_{k}^{\prime}[\vartheta]({\bf r},t)=\tilde{u}_{k}^{\prime}[\vartheta]({\bf
r},t)/\sqrt{(\tilde{u}_{k}^{\prime}[\vartheta]|\tilde{u}_{k}^{\prime}[\vartheta])}$
with

$\displaystyle\tilde{u}_{k}^{\prime}[\vartheta]({\bf r},t)$
$\displaystyle=\partial_{\vartheta}u_{k}[\vartheta]({\bf r},t)$ (89)
$\displaystyle\quad-\sum_{j=1}^{n}(u_{j}[\vartheta]|\partial_{\vartheta}u_{k}[\vartheta])u_{j}[\vartheta]({\bf
r},t)$
$\displaystyle\quad-\sum_{j=1}^{k-1}(u^{\prime}_{j}[\vartheta]|\partial_{\vartheta}u_{k}[\vartheta])u_{j}^{\prime}[\vartheta]({\bf
r},t).$

Accordingly, the coefficients $c^{\prime}_{kl}[\vartheta]$ are given by

$c^{\prime}_{kl}[\vartheta]=\begin{cases}\sqrt{(\tilde{u}_{k}^{\prime}[\vartheta]|\tilde{u}_{k}^{\prime}[\vartheta])}&{\rm
for}\;k=l\\\
(u^{\prime}_{j}[\vartheta]|\partial_{\vartheta}u_{k}[\vartheta])&{\rm
for}\;k<l\\\ 0&k>l\end{cases},$ (90)

resulting in a lower-triangular block matrix $D_{\partial}$.

##  Appendix B Derivative modes for temporal separation estimation

###  B.1 General expressions

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 $\tau$:

$\displaystyle\partial_{d}u_{0}[\tau](t)$
$\displaystyle=\frac{-\partial_{t}u(t-\tau/2)+\partial_{x}u(x+\tau/2)}{2\sqrt{2(1+\delta)}}$
$\displaystyle\quad-\delta^{\prime}\frac{u_{0}[\tau](t)}{2(1+\delta)},$ (91a)
$\displaystyle\partial_{d}v_{0}[\tau](t)$
$\displaystyle=\frac{-\partial_{t}u(t-\tau/2)-\partial_{t}u(t+\tau/2)}{2\sqrt{2(1-\delta)}}$
$\displaystyle\quad+\delta^{\prime}\frac{v_{0}[\tau](t)}{2(1-\delta)},$ (91b)

where we introduced $\delta^{\prime}=\partial_{\tau}\delta$. We assumed that
our pulses are symmetric, i.e. $u(t)=u(-t)$. We thus have
$(\partial_{t}u|u)=0$, which, combined with

$\displaystyle\delta^{\prime}$
$\displaystyle=-\int\partial_{t}u(t-\tau/2)u(t+\tau/2)dt$ (92)
$\displaystyle\;+\int u(t-\tau/2)\partial_{t}u(t+\tau/2)dt,$

implies that $(\partial_{\tau}v_{0}|u_{0})=(\partial_{\tau}u_{0}|v_{0})=0$.
Consequently, the orthogonalised first derivative modes are simply given by

$\displaystyle u_{1}[\tau](t)$
$\displaystyle=\partial_{\tau}u_{0}[\tau](t)/\eta_{u},$ (93a) $\displaystyle
v_{1}[\tau](t)$ $\displaystyle=\partial_{\tau}v_{0}[\tau](t)/\eta_{v},$ (93b)

with

$\displaystyle\eta_{u}^{2}$
$\displaystyle=||\partial_{\tau}u_{0}[\tau](t)||^{2}=\frac{(\Delta
k)^{2}-\beta}{4(1+\delta)}-\frac{(\delta^{\prime})^{2}}{4(1+\delta)^{2}},$
(94a) $\displaystyle\eta_{v}^{2}$
$\displaystyle=||\partial_{\tau}v_{0}[\tau](t)||^{2}=\frac{(\Delta
k)^{2}+\beta}{4(1-\delta)}-\frac{(\delta^{\prime})^{2}}{4(1-\delta)^{2}},$
(94b)

where we introduced

$\displaystyle(\Delta k)^{2}$
$\displaystyle=\int[\partial_{t}u(t)]^{2}dt,\quad{\rm and}$ (95)
$\displaystyle\beta$
$\displaystyle=\int\partial_{t}u(t-\tau/2)\partial_{t}u(t+\tau/2)dt.$ (96)

We now move to the construction of the orthonormalized second derivatives. The
fact that the modes $u_{1}[\tau](t)$ and $v_{1}[\tau](t)$ have, by
construction, opposite parity implies that
$(\partial_{\tau}u_{1}|u_{1})=(\partial_{\tau}v_{1}|v_{1})=0$, and
$(\partial_{\tau}u_{1}|v_{1})=-(\partial_{\tau}v_{1}|u_{1})$. Let us then
evaluate $(\partial_{d}u_{1}|v_{1})$ explicitly

$\displaystyle(\partial_{\tau}u_{1}|v_{1})$
$\displaystyle=\frac{(\partial_{\tau}u_{1}|\partial_{\tau}v_{0})}{\eta_{v}}$
(97)
$\displaystyle=\frac{(\partial^{2}_{\tau}u_{0}|\partial_{\tau}v_{0})}{\eta_{u}\eta_{v}}-\frac{\eta^{\prime}_{u}}{\eta_{v}\eta_{u}^{2}}(\partial_{\tau}u_{0}|\partial_{\tau}v_{0})$
$\displaystyle=\frac{(\partial^{2}_{\tau}u_{0}|\partial_{\tau}v_{0})}{\eta_{u}\eta_{v}},$

where in the last step we used that $u_{0}[\tau](t)$ and $v_{0}[\tau](t)$ are
even and odd functions of $t$, respectively. The second derivative of
$u_{0}(x)$ and $v_{0}(x)$ can be rewritten as

$\displaystyle\partial_{t}^{2}u_{0}[\tau](t)$
$\displaystyle=\frac{f_{u}[\tau](t)}{4\sqrt{2(1+\delta)}}-\frac{\delta^{\prime}}{1+\delta}\partial_{d}u_{0}[\tau](t)$
$\displaystyle\quad+C_{u}u_{0}[\tau](t),$ (98a)
$\displaystyle\partial_{t}^{2}v_{0}[\tau](t)$
$\displaystyle=\frac{f_{v}[\tau](t)}{4\sqrt{2(1-\delta)}}+\frac{\delta^{\prime}}{1-\delta}\partial_{d}v_{0}[\tau](t)$
$\displaystyle\quad+C_{v}v_{0}[\tau](t),$ (98b)

with

$\displaystyle f_{u}[\tau](t)$
$\displaystyle=\partial^{2}_{t}u(t-\tau/2)+\partial^{2}_{t}u(t+\tau/2),$ (99)
$\displaystyle f_{v}[\tau](t)$
$\displaystyle=\partial^{2}_{t}u(t-\tau/2)-\partial^{2}_{t}u(t+\tau/2),$ (100)

and

$\displaystyle C_{u}$
$\displaystyle=\frac{(\delta^{\prime})^{2}-2(1+\delta)\delta^{\prime\prime}}{4(1+\delta)^{2}},$
(101) $\displaystyle C_{v}$
$\displaystyle=\frac{(\delta^{\prime})^{2}+2(1-\delta)\delta^{\prime\prime}}{4(1-\delta)^{2}},$
(102)

where we have introduced the second derivative of the overlap parameter
$\delta^{\prime\prime}=\partial_{\tau}^{2}\delta$. Using
$(\partial_{\tau}u_{0}|\partial_{\tau}v_{0})=(u_{0}|\partial_{\tau}v_{0})=0$,
from Eq. (98a) we have

$\displaystyle(\partial_{d}u_{1}|v_{1})\propto$ (103)
$\displaystyle\int\partial_{t}^{2}u(t+\tau)\partial_{t}u(t)dt+\int\partial_{t}^{2}u(-t-\tau)\partial_{t}u(-t)dt$
$\displaystyle=\int\partial_{t}^{2}u(t+\tau)\partial_{t}u(t)dt-\int\partial_{t}^{2}u(t+\tau)\partial_{t}u(t)dt$
$\displaystyle=0,$

where we used that the pulse shape $u(t)$ is an even function of $t$.
Therefore, the orthogonalized second derivative modes are given by

$\displaystyle u_{2}[\tau](t)$
$\displaystyle=(\partial_{d}u_{1}[\tau](t)-\xi_{u}u_{0}([\tau](t))/\zeta_{u},$
(104) $\displaystyle v_{2}[\tau](t)$
$\displaystyle=(\partial_{d}v_{1}[\tau](t)-\xi_{v}v_{0}[\tau](t))/\zeta_{v},$
(105)

with $\xi_{u}=(\partial_{\tau}u_{1}|u_{0})$,
$\xi_{v}=(\partial_{\tau}v_{1}|v_{0})$,
$\zeta_{u}=||\partial_{\tau}u_{1}[\tau(t)-\xi_{u}u_{0}[\tau](t)||$ and
$\zeta_{v}=||\partial_{\vartheta}v_{1}[\vartheta](t)-\xi_{v}v_{0}[\vartheta](t)||$.

Let us explicitly calculate $\xi_{u}$ and $\xi_{v}$. This can be achieved by
using Eqs. (98a), (98b) and by noting that

$\displaystyle\int$ $\displaystyle\partial_{t}^{2}\left[u(t-\tau/2)\pm
u(t+\tau/2)\right]$ (106) $\displaystyle\times\left[u(t-\tau/2)\pm
u(t-\tau/2)\right]dt$
$\displaystyle=-\int\left[\partial_{t}u(t-\tau/2)\pm\partial_{t}u(t+\tau/2)\right]^{2}dt$
$\displaystyle=-2(\Delta k)^{2}\mp 2\beta,$

where we used partial integration and made the reasonable assumption that the
pulse shape $u(t)$ goes to zero at infinity. We then get

$\displaystyle\xi_{u}$
$\displaystyle=\frac{1}{\eta_{u}}\left(\frac{(\delta^{\prime})^{2}-2(1+\delta)\delta^{\prime\prime}}{4(1+\delta)^{2}}-\frac{(\Delta
k)^{2}+\beta}{4(1+\delta)}\right)$ (107a)
$\displaystyle=\frac{(\delta^{\prime})^{2}-2(1+\delta)\delta^{\prime\prime}-(1+\delta)\left((\Delta
k)^{2}+\beta\right)}{(1+\delta)\sqrt{(1+\delta)\left((\Delta
k)^{2}-\beta\right)-(\delta^{\prime})^{2}}},$ $\displaystyle\xi_{v}$
$\displaystyle=\frac{1}{\eta_{v}}\left(\frac{(\delta^{\prime})^{2}+2(1-\delta)\delta^{\prime\prime}}{4(1-\delta)^{2}}-\frac{(\Delta
k)^{2}-\beta}{4(1-\delta)}\right)$ (107b)
$\displaystyle=\frac{(\delta^{\prime})^{2}+2(1-\delta)\delta^{\prime\prime}-(1+\delta)\left((\Delta
k)^{2}-\beta\right)}{(1-\delta)\sqrt{(1-\delta)\left((\Delta
k)^{2}+\beta\right)-(\delta^{\prime})^{2}}}.$

Let us now compute explicitly the normalization constants of the modes
$u_{2}[\tau](t)$ and $v_{2}[\tau](t)$

$\displaystyle\zeta_{u}^{2}$
$\displaystyle=||\partial_{\tau}u_{1}[\tau](t)-\xi_{u}u_{0}[\tau](t)||^{2}$
(108a)
$\displaystyle=||\partial_{\tau}u_{1}[\tau](t)||^{2}+\xi_{u}^{2}-2\xi_{u}(\partial_{\tau}u_{1}|u_{0})$
$\displaystyle=||\partial_{\tau}u_{1}[\tau](t)||^{2}-\xi_{u}^{2},$
$\displaystyle\zeta_{v}^{2}$
$\displaystyle=||\partial_{\tau}v_{1}[\tau](t)-\xi_{v}v_{0}[\tau](t)||^{2}$
(108b)
$\displaystyle=||\partial_{\tau}v_{1}[\tau](t)||^{2}+\xi_{v}^{2}-2\xi_{v}(\partial_{\tau}v_{1}|v_{0})$
$\displaystyle=||\partial_{\tau}v_{1}[\tau](t)||^{2}-\xi_{v}^{2}.$

We can then expand

$\displaystyle||\partial_{\tau}u_{1}[\tau](t)||^{2}$
$\displaystyle=\left|\left|\frac{\partial^{2}_{\tau}u_{0}[\tau](t)}{\eta_{u}}-\frac{(\partial_{\tau}\eta_{u})^{2}}{\eta_{u}^{2}}\partial_{\tau}u_{0}[\tau](t)\right|\right|^{2}$
(109a)
$\displaystyle=\frac{||\partial^{2}_{\tau}u_{0}||^{2}+(\partial_{\tau}\eta_{u})^{2}}{\eta_{u}^{2}}$
$\displaystyle\;-2\frac{\partial_{\tau}\eta_{u}}{\eta_{u}^{3}}(\partial_{\tau}^{2}u_{0}|\partial_{\tau}u_{0}),$
$\displaystyle||\partial_{\tau}v_{1}[\tau](t)||^{2}$
$\displaystyle=\left|\left|\frac{\partial^{2}_{\tau}v_{0}[\tau](t)}{\eta_{v}}-\frac{(\partial_{\tau}\eta_{v})^{2}}{\eta_{v}^{2}}\partial_{\tau}v_{0}[\tau](t)\right|\right|^{2}$
(109b)
$\displaystyle=\frac{||\partial^{2}_{\tau}v_{0}[\tau](t)||^{2}+(\partial_{\tau}\eta_{v})^{2}}{\eta_{v}^{2}}$
$\displaystyle\;-2\frac{\partial_{\tau}\eta_{v}}{\eta_{v}^{3}}(\partial_{\tau}^{2}v_{0}|\partial_{\tau}v_{0}).$

From Eqs. (98a) and (98b), we then have

$\displaystyle||\partial^{2}_{\tau}u_{0}[\tau](t)||^{2}$
$\displaystyle=\frac{||f_{u}[\tau](t)||^{2}}{32(1+\delta)}+\frac{(\delta^{\prime})^{2}\eta_{u}^{2}}{(1+\delta^{2})}+C_{u}^{2}$
(110a)
$\displaystyle\;-\frac{\delta^{\prime}(f_{u}|\partial_{\tau}u_{0})}{4\sqrt{2}(1+\delta)^{3/2}}-C_{u}\frac{2(\Delta
k)^{2}+2\beta}{\sqrt{2(1+\delta)}}$
$\displaystyle||\partial^{2}_{\tau}v_{0}[\tau](t)||^{2}$
$\displaystyle=\frac{||f_{v}[\tau](t)||^{2}}{32(1-\delta)}+\frac{(\delta^{\prime})^{2}\eta_{v}^{2}}{(1-\delta^{2})}+C_{v}^{2}$
(110b)
$\displaystyle\;+\frac{\delta^{\prime}(f_{v}|\partial_{\tau}v_{0})}{4\sqrt{2}(1-\delta)^{3/2}}-C_{v}\frac{2(\Delta
k)^{2}-2\beta}{\sqrt{2(1-\delta)}},$

with

$\displaystyle||f_{u}[\tau](t)||^{2}$ $\displaystyle=2(\sigma+\epsilon)$
(111a) $\displaystyle||f_{v}[\tau](t)||^{2}$
$\displaystyle=2(\sigma-\epsilon),$ (111b)

and

$\displaystyle(f_{u}|\partial_{d}u_{0})$ $\displaystyle=\frac{(\Delta
k)^{2}+\beta}{1+\delta}\delta^{\prime}+\frac{\beta^{\prime}}{\sqrt{2(1+\delta)}},$
(112a) $\displaystyle(f_{v}|\partial_{d}v_{0})$ $\displaystyle=\frac{(\Delta
k)^{2}-\beta}{1-\delta}\delta^{\prime}+\frac{\beta^{\prime}}{\sqrt{2(1-\delta)}},$
(112b)

where we defined

$\displaystyle\sigma$
$\displaystyle=\int\left|\partial_{t}^{2}u(t)\right|^{2}dt,$ (113a)
$\displaystyle\epsilon$
$\displaystyle=\int\partial_{t}^{2}u(t-\tau/2)\partial_{t}^{2}u(t-\tau/2)dx.$
(113b)

Since the expressions of the mode quantities (especially $\zeta_{u}$ and
$\zeta_{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}}/(\pi w^{2})^{1/4}$ in the following:

$\displaystyle\eta_{u}$
$\displaystyle=\frac{\sqrt{\tau^{2}+4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)}\text{sech}\left(\frac{\tau^{2}}{8w^{2}}\right)}{8w^{2}},$
(114) $\displaystyle\eta_{v}$
$\displaystyle=\frac{\sqrt{4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)-\tau^{2}}\text{csch}\left(\frac{\tau^{2}}{8w^{2}}\right)}{8w^{2}},$
(115) $\displaystyle\xi_{u}$
$\displaystyle=-\frac{\sqrt{\tau^{2}+4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)}\text{sech}\left(\frac{\tau^{2}}{8w^{2}}\right)}{8w^{2}},$
(116) $\displaystyle\xi_{v}$
$\displaystyle=-\frac{\sqrt{4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)-\tau^{2}}\text{csch}\left(\frac{\tau^{2}}{8w^{2}}\right)}{8w^{2}},$
(117) $\displaystyle\zeta^{2}_{u}$
$\displaystyle=\frac{2w^{2}\sinh\left(\frac{\tau^{2}}{2w^{2}}\right)-\tau^{2}\left(1+\cosh\left(\frac{\tau^{2}}{4w^{2}}\right)\right)}{\left(\tau^{2}+4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)\right)^{2}}$
(118)
$\displaystyle\quad+\frac{\left(\tau^{4}+16w^{4}\right)\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)}{4w^{2}\left(\tau^{2}+4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)\right)^{2}}$
$\displaystyle\zeta_{v}^{2}$
$\displaystyle=\frac{2w^{2}\sinh\left(\frac{\tau^{2}}{2w^{2}}\right)-\tau^{2}\left(1-\cosh\left(\frac{\tau^{2}}{4w^{2}}\right)\right)}{\left(\tau^{2}-4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)\right)^{2}}$
(119)
$\displaystyle\quad-\frac{\left(\tau^{4}+16w^{4}\right)\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)}{4w^{2}\left(\tau^{2}-4w^{2}\sinh\left(\frac{\tau^{2}}{4w^{2}}\right)\right)^{2}}.$

Note that the fact that $\eta_{u,v}^{2}=\xi_{u,v}^{2}$ is a peculiarity of
Gaussian pulses and it is not true in general.

###  B.2 Small $\tau$ behaviour

Arguably, the most interesting regime for temporal separation estimation is
that of small $\tau$. Therefore, in the following, we discuss the behaviour of
the quantities computed above for $\tau\to 0$. Let us start by considering the
following series expansions

$\displaystyle\delta$ $\displaystyle=1-(\Delta
k)^{2}\frac{\tau^{2}}{2}+\sigma\frac{\tau^{4}}{24}+\mathcal{O}(\tau^{6}),$
(120a) $\displaystyle\beta$ $\displaystyle=(\Delta
k)^{2}-\sigma\frac{\tau^{2}}{2}+\mathcal{O}(\tau^{4}),$ (120b)
$\displaystyle\epsilon$ $\displaystyle=\sigma+\mathcal{O}(\tau^{2}).$ (120c)

Using Eqs. (120), it is possible to show that for $\tau\sim 0$ we have
$\eta_{u}\sim\eta_{v}\sim\xi_{u}\sim\xi_{v}\sim\zeta_{u}\sim\zeta_{v}\sim\tau$.
For example, for Gaussian pulses we have

$\displaystyle\eta_{u}$
$\displaystyle=\frac{\tau}{4\sqrt{2}w^{2}}+\mathcal{O}\left(\tau^{2}\right),$
(121a) $\displaystyle\eta_{v}$
$\displaystyle=\frac{\tau}{4\sqrt{6}w^{2}}+\mathcal{O}\left(\tau^{2}\right),$
(121b) $\displaystyle\xi_{u}$
$\displaystyle=-\frac{\tau}{4\sqrt{2}w^{2}}+\mathcal{O}\left(\tau^{2}\right),$
(121c) $\displaystyle\xi_{v}$
$\displaystyle=-\frac{\tau}{4\sqrt{6}w^{2}}+\mathcal{O}\left(\tau^{2}\right),$
(121d) $\displaystyle\zeta_{u}$
$\displaystyle=\frac{\tau}{4\sqrt{3}w^{2}}+\mathcal{O}\left(\tau^{2}\right),$
(121e) $\displaystyle\zeta_{v}$
$\displaystyle=\frac{\tau}{4\sqrt{5}w^{2}}+\mathcal{O}\left(\tau^{2}\right).$
(121f)

This behaviour implies that the contribution to the QFI coming from the
coefficients $b_{l}^{jk}$ and $d_{l}^{jk}$ vanishes for $\tau\to 0$ (see
Eqs.(68) and App. C).

##  Appendix 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 $\tau$ between two
incoherent thermal pulses aided by two squeezed vacuum states in the
derivative modes $u_{1}[\tau](t)$ and $v_{1}[\tau](t)$ defined by the
covariance matrices $V_{0}$ (70) and $V_{1}~{}\eqref{V_1_tau}$.

First, we note that $V_{0}$ (70) is already in Williamson form. Therefore, the
symplectic matrix $S_{0}$ entering in Eqs. (68) is the identity
$S_{0}=\mathds{1}_{4}$, while for the symplectic eigenvalues we have
$\nu_{0}^{0}=2N_{0}(1+\delta)+1$ and $\nu_{0}^{1}=2N_{0}(1-\delta)+1$. On the
other hand, the Williamson decomposition $V_{1}$ (71) is achieved by the
squeezing matrix

$S_{1}=\begin{pmatrix}e^{-r}&0&0&0\\\ 0&e^{r}&0&0\\\ 0&0&e^{-r}&0\\\
0&0&0&e^{r}\\\ \end{pmatrix},$ (122)

with symplectic eigenvalues $\nu_{1}^{0,1}=1$.

To compute the coefficients $a_{l}^{jk}$, we need the derivative of the matrix
$V_{0}$ (70). The latter depends on the temporal separation $\tau$ only
through the overlap parameter $\delta$. Accordingly, we have

$\partial_{\tau}V_{0}=\begin{pmatrix}2N_{0}(\partial_{\tau}\delta)\mathds{1}_{2}&0\\\
0&-2N_{0}(\partial_{\tau}\delta)\mathds{1}_{2}\end{pmatrix}.$ (123)

Since $\partial_{\tau}V_{0}$ is diagonal, the only nonzero $a_{l}^{jk}$
coefficients are

$\left(a^{00}_{2}\right)^{2}=\left(a^{11}_{2}\right)^{2}=8N_{0}(\partial_{\tau}\delta).$
(124)

Using $S_{0}=\mathds{1}_{4}$, and the fact that all matrices $D_{\xi}$,
$D_{\eta}$, and $V_{0}$ are diagonal, we have
$\left(b_{l}^{jk}\right)^{2}=2\left(\tr\left\\{A_{l}^{(jk)}M_{b}\right\\}\right)^{2},$
where we introduced the matrix

$M_{b}=\begin{pmatrix}X_{u}+Y_{u}&0&0&0\\\ 0&X_{u}-Y_{u}&0&0\\\
0&0&X_{v}+Y_{v}&0\\\ 0&0&0&X_{v}-Y_{v}\end{pmatrix},$ (125)

with

$\displaystyle X_{u}$ $\displaystyle=N_{0}(1+\delta)\eta_{u}\cosh r,$ (126a)
$\displaystyle X_{v}$ $\displaystyle=N_{0}(1-\delta)\eta_{v}\cosh r,$ (126b)
$\displaystyle Y_{u}$
$\displaystyle=\left[N_{0}(1+\delta)\eta_{u}-\xi_{u}\right]\sinh r,$ (126c)
$\displaystyle Y_{v}$
$\displaystyle=\left[N_{0}(1-\delta)\eta_{v}-\xi_{v}\right]\sinh r.$ (126d)

Consequently, the only nonzero $b_{l}^{jk}$ coefficients are

$\displaystyle\left(b_{1}^{00}\right)^{2}$
$\displaystyle=8\left[N_{0}(1+\delta)\eta_{u}-\xi_{u}\right]^{2}\sinh^{2}r,$
(127a) $\displaystyle\left(b_{2}^{00}\right)^{2}$
$\displaystyle=8N_{0}^{2}(1+\delta)^{2}\eta_{u}^{2}\cosh^{2}r,$ (127b)
$\displaystyle\left(b_{1}^{11}\right)^{2}$
$\displaystyle=8\left[N_{0}(1-\delta)\eta_{v}-\xi_{v}\right]^{2}\sinh^{2}r,$
(127c) $\displaystyle\left(b_{2}^{11}\right)^{2}$
$\displaystyle=8N_{0}^{2}(1-\delta)^{2}\eta_{v}^{2}\cosh^{2}r.$ (127d)

We assumed the covariance matrix $V_{1}$ of the derivative modes
$u_{1}[\tau](t)$ and $v_{1}[\tau](t)$ to be parameter independent, i.e.
$\partial_{\tau}V_{1}=0$, which leads to $c_{l}^{(jk)}=0$ for all $l,j$ and
$k$.

Finally, we can write the $d_{l}^{jk}$ coefficients as

$\left(d_{l}^{jk}\right)^{2}=2\left(\tr\left\\{A_{l}^{(jk)}M_{d}\right\\}\right)^{2}$
(128)

with

$M_{d}=D_{\zeta}(V_{1}-\mathds{1}_{4})S_{1}^{-1}=2\sinh
r\begin{pmatrix}\zeta_{u}\sigma_{z}&0\\\ 0&\zeta_{v}\sigma_{z}\end{pmatrix}.$
(129)

As consequence, the only nonzero $d_{l}^{jk}$ coefficients are

$\displaystyle\left(d_{1}^{00}\right)^{2}$
$\displaystyle=8\zeta_{u}^{2}\sinh^{2}r,$ (130a)
$\displaystyle\left(d_{1}^{11}\right)^{2}$
$\displaystyle=8\zeta_{v}^{2}\sinh^{2}r.$ (130b)

Substituting the coefficients in Eqs. (124), (127) and (130) into Eq. (67), we
then obtain the QFI (72).

## References

  * Helstrom (1976) C. W. Helstrom,  _Quantum detection and estimation theory_ , Vol. 3 (Academic press, New York, 1976).
  * Holevo (2011) A. S. Holevo,  _Probabilistic and statistical aspects of quantum theory_ , Vol. 1 (Springer, Berlin, 2011).
  * Paris (2009) M. G. A. Paris, International Journal of Quantum Information 07, 125 (2009).
  * Pezzè and Smerzi (2014) L. Pezzè and A. Smerzi, in  _Proceedings of the International School of Physics ”Enrico Fermi”_, Course 188, Varenna, edited by G. M. Tino and M. A. Kasevich (IOS Press, Amsterdam, 2014) pp. 691 – 741.
  * Giovannetti  _et al._ (2011) V. Giovannetti, S. Lloyd, and L. Maccone, Nature Photonics 5, 222 (2011).
  * Tóth and Apellaniz (2014) G. Tóth and I. Apellaniz, Journal of Physics A: Mathematical and Theoretical 47, 424006 (2014).
  * Taylor and Bowen (2016) M. A. Taylor and W. P. Bowen, Physics Reports 615, 1 (2016).
  * Tsang (2019) M. Tsang, Contemporary Physics 60, 279 (2019).
  * Giovannetti  _et al._ (2001) V. Giovannetti, S. Lloyd, and L. Maccone, Nature 412, 417 (2001).
  * Zhuang  _et al._ (2017) Q. Zhuang, Z. Zhang, and J. H. Shapiro, Phys. Rev. A 96, 040304 (2017).
  * Huang  _et al._ (2021) Z. Huang, C. Lupo, and P. Kok, PRX Quantum 2, 030303 (2021).
  * Acernese et al. (2019) F. Acernese et al. (Virgo Collaboration), Phys. Rev. Lett. 123, 231108 (2019).
  * Tse et al. (2019) M. Tse et al. (LIGO Collaboration), Phys. Rev. Lett. 123, 231107 (2019).
  * Fabre and Treps (2020) C. Fabre and N. Treps, Rev. Mod. Phys. 92, 035005 (2020).
  * Treps  _et al._ (2002) N. Treps, U. Andersen, B. Buchler, P. K. Lam, A. Maître, H.-A. Bachor, and C. Fabre, Phys. Rev. Lett. 88, 203601 (2002).
  * Treps  _et al._ (2003) N. Treps, N. Grosse, W. P. Bowen, C. Fabre, H.-A. Bachor, and P. K. Lam, Science 301, 940 (2003).
  * Cai  _et al._ (2021) Y. Cai, J. Roslund, V. Thiel, C. Fabre, and N. Treps, npj Quantum Information 7, 82 (2021).
  * Gessner  _et al._ (2022) M. Gessner, N. Treps, and C. Fabre, arXiv:2201.04050 (2022).
  * Tsang  _et al._ (2016) M. Tsang, R. Nair, and X.-M. Lu, Phys. Rev. X 6, 031033 (2016).
  * Paúr  _et al._ (2016) M. Paúr, B. Stoklasa, Z. Hradil, L. L. Sánchez-Soto, and J. Rehacek, Optica 3, 1144 (2016).
  * Boucher  _et al._ (2020) P. Boucher, C. Fabre, G. Labroille, and N. Treps, Optica 7, 1621 (2020).
  * Ansari  _et al._ (2021) V. Ansari, B. Brecht, J. Gil-Lopez, J. M. Donohue, J. Řeháček, Z. c. v. Hradil, L. L. Sánchez-Soto, and C. Silberhorn, PRX Quantum 2, 010301 (2021).
  * De  _et al._ (2021) S. De, J. Gil-Lopez, B. Brecht, C. Silberhorn, L. L. Sánchez-Soto, Z. c. v. Hradil, and J. Řeháček, Phys. Rev. Res. 3, 033082 (2021).
  * Mazelanik  _et al._ (2022) M. Mazelanik, A. Leszczyński, and M. Parniak, Nature Communications 13, 691 (2022).
  * Lupo and Pirandola (2016) C. Lupo and S. Pirandola, Phys. Rev. Lett. 117, 190802 (2016).
  * Holevo (1975) A. Holevo, IEEE Transactions on Information Theory 21, 533 (1975).
  * Weedbrook  _et al._ (2012) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Reviews of Modern Physics 84, 621 (2012).
  * Adesso  _et al._ (2014) G. Adesso, S. Ragy, and A. R. Lee, Open Systems & Information Dynamics 21, 1440001 (2014).
  * Caves (1981) C. M. Caves, Phys. Rev. D 23, 1693 (1981).
  * Pezzé and Smerzi (2008) L. Pezzé and A. Smerzi, Phys. Rev. Lett. 100, 073601 (2008).
  * Pinel  _et al._ (2012) O. Pinel, J. Fade, D. Braun, P. Jian, N. Treps, and C. Fabre, Phys. Rev. A 85, 010101 (2012).
  * Monras (2013) A. Monras, arXiv:1303.3682 (2013).
  * Šafránek  _et al._ (2015) D. Šafránek, A. R. Lee, and I. Fuentes, New Journal of Physics 17, 073016 (2015).
  * Jiang (2014) Z. Jiang, Phys. Rev. A 89, 032128 (2014).
  * Nair and Tsang (2016) R. Nair and M. Tsang, Phys. Rev. Lett. 117, 190801 (2016).
  * Sorelli  _et al._ (2022) G. Sorelli, M. Gessner, M. Walschaers, and N. Treps, Phys. Rev. Res. 4, L032022 (2022).
  * Braunstein and van Loock (2005) S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).
  * Serafini (2017) A. Serafini,  _Quantum continuous variables: a primer of theoretical methods_ (Tayor & Francis, Milton Park, 2017).
  * Simon  _et al._ (1994) R. Simon, N. Mukunda, and B. Dutta, Phys. Rev. A 49, 1567 (1994).
  * Williamson (1936) J. Williamson, American Journal of Mathematics 58, 141 (1936).
  * Walschaers (2021) M. Walschaers, PRX Quantum 2, 030204 (2021).
  * Holevo and Werner (2001) A. S. Holevo and R. F. Werner, Phys. Rev. A 63, 032312 (2001).
  * Braunstein and Caves (1994) S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 3439 (1994).
  * Sorelli  _et al._ (2021a) G. Sorelli, M. Gessner, M. Walschaers, and N. Treps, Phys. Rev. Lett. 127, 123604 (2021a).
  * Sorelli  _et al._ (2021b) G. Sorelli, M. Gessner, M. Walschaers, and N. Treps, Phys. Rev. A 104, 033515 (2021b).
  * Van Trees (2002) H. L. Van Trees,  _Detection, estimation, and modulation theory, part I: detection, estimation, and linear modulation theory_ (John Wiley & Sons, New York, 2002).
  * Van Trees (2001) H. L. Van Trees,  _Detection, estimation, and modulation theory, part III: Radar-Sonar Signal Processing and Gaussian Signals in Noise_ (John Wiley & Sons, New York, 2001).
  * Zhuang and Shapiro (2022) Q. Zhuang and J. H. Shapiro, Phys. Rev. Lett. 128, 010501 (2022).
  * Nichols  _et al._ (2018) R. Nichols, P. Liuzzo-Scorpo, P. A. Knott, and G. Adesso, Phys. Rev. A 98, 012114 (2018).