# Thermometry of Gaussian quantum systems using Gaussian measurements

Marina F.B. Cenni  ICFO-Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
Ludovico Lami  Institut für Theoretische Physik und IQST, Universität Ulm,
Albert-Einstein-Allee 11, D-89069 Ulm, Germany  Antonio Acín  ICFO-Institut de
Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels (Barcelona), Spain  ICREA-Institució Catalana de Recerca i
Estudis Avançats, 08010, Barcelona, Spain  Mohammad Mehboudi  Département de
Physique Appliquée, Université de Genève, 1211 Genève, Switzerland

###### Abstract

We study the problem of estimating the temperature of Gaussian systems with
feasible measurements, namely Gaussian and photo-detection-like measurements.
For Gaussian measurements, we develop a general method to identify the optimal
measurement numerically, and derive the analytical solutions in some relevant
cases. For a class of single-mode states that includes thermal ones, the
optimal Gaussian measurement is either Heterodyne or Homodyne, depending on
the temperature regime. This is in contrast to the general setting, in which a
projective measurement in the eigenbasis of the Hamiltonian is optimal
regardless of temperature. In the general multi-mode case, and unlike the
general unrestricted scenario where joint measurements are not helpful for
thermometry (nor for any parameter estimation task), it is open whether joint
Gaussian measurements provide an advantage over local ones. We conjecture that
they are not useful for thermal systems, supported by partial analytical and
numerical evidence. We further show that Gaussian measurements become optimal
in the limit of large temperatures, while photo-detection-like measurements do
it for when the temperature tends to zero. Our results therefore pave the way
for effective thermometry of Gaussian quantum systems using experimentally
realizable measurements.

###### Contents

  1. I Introduction and motivation
  2. II Preliminaries: The Phase-space formalism
    1. II.1 Gaussian States and Measurements
      1. II.1.1 The overlap of two systems
  3. III Thermometry of Gaussian quantum systems: Optimal vs feasible measurements
    1. III.1 Ultimate bounds on thermometry of Gaussian systems
    2. III.2 Thermometry using Gaussian measurements
      1. III.2.1 Observation.—Pure Gaussian measurements are always better
      2. III.2.2 Classical Fisher information of Gaussian measurements
      3. III.2.3 Homodyne detection and equivalence with the error propagation
      4. III.2.4 Heterodyne measurement
      5. III.2.5 The optimal Gaussian measurement
      6. III.2.6 Optimal Gaussian measurements on multimode systems
      7. III.2.7 Multimode Gaussian states at thermal equilibrium: Local versus global Gaussian measurements
    3. III.3 Photo-Detection
      1. III.3.1 Multimode photo-detection
  4. IV Conclusions
  5. A More about Gaussian measurements
    1. A.1 General form of Gaussian measurements
    2. A.2 Local Gaussian measurements, reduced states and covariance matrices
  6. B The classical Fisher information of Gaussian probability distributions
  7. C Additivity of the Gaussian-measured Fisher information in a special case
  8. D Numerical optimisation, and the algorithm

##  I Introduction and motivation

The high precision required for the operation of quantum devices demands for
the estimation of their parameters with minimum possible error. In particular,
estimating temperature is a problem that has attained substantial popularity
in recent years Mehboudi  _et al._ (2019a); De Pasquale and Stace (2018),
since it is a parameter that highly impacts the technological usage of quantum
systems. Thermometry process, like any other parameter estimation problem, can
be seen as different stages of (i) preparation of a probe or thermometer, (ii)
interaction of the probe with the sample, and (iii) measurement of the probe
and data analysis Giovannetti  _et al._ (2011, 2006); Tóth and Apellaniz
(2014), see also Fig. 1. When there are no restrictions on any of these
stages, one finds the ultimate limits on thermometry precision Mehboudi  _et
al._ (2019a); De Pasquale and Stace (2018). These limits have deeply increased
our understanding of fundamental restrictions/possibilities in any thermometry
procedure Correa  _et al._ (2015). Nonetheless, one also needs to address
frameworks that are experimentally feasible, hence revisit the fundamental
limits subject to experimental restrictions on the different stages above. In
this work, we are in particular concerned with restrictions on the performed
measurements, stage (iii) above. To our knowledge, there have been so far few
works addressing measurement restrictions, such as limited resolution
Jorgensen  _et al._ (2020) or coarse grained measurements Hovhannisyan  _et
al._ (2021), in thermometry applications.

Our main goal is to study the use of experimentally feasible measurements in
thermometry processes using continuous variable (CV) Gaussian systems. CV
Gaussian models are very successful at describing various physical systems of
relevance, such as Bosonic gases, electromagnetic fields, and Josephson
junctions, to name a few. This makes them versatile candidates for several
quantum information processing as well as quantum simulation tasks Weedbrook
_et al._ (2012); Lewenstein  _et al._ (2007, 2012); Serafini (2017a); Diamanti
and Leverrier (2015); Navascués and Acín (2005). When there is no limitation
on the performed measurements, thermometry (and more generally metrology) of
Gaussian CV systems is characterized readily, in particular the optimal
measurement and the ultimate precision can be found analytically Monras
(2013); Jiang (2014); Nichols  _et al._ (2018). Nonetheless, a rigorous
analysis of quantum metrology with realistic measurements, namely Gaussian as
well as photo-detection measurements, is missing. Our work addresses this
question in the context of temperature estimation.

Several previous theoretical works have already focused on thermometry via
Gaussian probes, e.g., exploiting Bosonic impurities embedded in cold gases,
or light modes that interact with an external field Mehboudi  _et al._
(2019b); Planella  _et al._ (2020); Correa  _et al._ (2017); Hovhannisyan and
Correa (2018); Mukherjee  _et al._ (2019); Martín-Martínez  _et al._ (2013).
Unfortunately, in these systems the optimal measurement for thermometry can be
experimentally very demanding and effectively out of reach. For example,
consider a scenario where a light mode is used as thermometer. After
interacting with a sample with an unknown temperature $T$, the light mode can
relax to a thermal state with the same temperature. By performing measurements
on the light mode one aims to estimate the temperature. The maximally
informative measurement in this case is an energy measurement. This
corresponds to ideal photon-counting, which is extremely difficult. It is thus
advisable to investigate experimentally realisable measurements, such as
realistic photo-detection or Gaussian measurements and benchmark their
performance by comparing them with the known optimal measurement. In this
work, we follow this approach and first consider single-mode systems. We
analytically prove that the optimal Gaussian measurement corresponds to a
homodyne or heterodyne measurement, depending on the temperature. We
investigate the conditions under which the ultimate bound on thermometry
precision is attained. We find that at high temperatures Gaussian measurements
are almost optimal, while they perform poorly in the low-temperature regime,
where, in turn, photo-detection is close to optimal. We then move to the
multi-mode scenario. While in the case where measurements are not restricted
the additivity of the Fisher information implies that the optimal measurement
consist of the independent application of the optimal measurement for single
modes, whether the same holds true when restricting the measurements to be
Gaussian is open. We conjecture—supported by numerical evidence—that this is
indeed the case, and analytically prove it for special cases.

Figure 1: Thermometry protocols can be broken into three parts of (i)
preparation of a probe, (ii) interaction with the sample and parametrisation,
and (iii) measurement and data analysis. The focus of this work is on the last
part of the protocol (iii). In particular, we explore the precision bounds
restricted to experimentally feasible measurements, namely Gaussian or photo-
detection measurements.

##  II Preliminaries: The Phase-space formalism

We start by setting up the notation and introducing the necessary tools from
the phase-space formalism, such as Gaussian systems and Gaussian measurements.
For additional details, we refer the interested reader to the more complete
references Holevo (2011); Barnett and Radmore (2002); Serafini (2017b). The
basics of parameter estimation in Gaussian systems is discussed in the next
section, where we also present some of our results.

###  II.1 Gaussian States and Measurements

Every trace class111A linear bounded operator $T:\mathcal{H}\to\mathcal{H}$
acting on a separable Hilbert space $\mathcal{H}$ is said to be of trace class
if $\operatorname{Tr}|T|<\infty$, where
$|T|\coloneqq\mathchoice{{\hbox{$\displaystyle\sqrt{T^{\dagger}T\,}$}\lower
0.4pt\hbox{\vrule
height=8.77776pt,depth=-7.02223pt}}}{{\hbox{$\textstyle\sqrt{T^{\dagger}T\,}$}\lower
0.4pt\hbox{\vrule
height=8.77776pt,depth=-7.02223pt}}}{{\hbox{$\scriptstyle\sqrt{T^{\dagger}T\,}$}\lower
0.4pt\hbox{\vrule
height=6.17221pt,depth=-4.93779pt}}}{{\hbox{$\scriptscriptstyle\sqrt{T^{\dagger}T\,}$}\lower
0.4pt\hbox{\vrule height=4.80554pt,depth=-3.84445pt}}}$, and for a positive
semidefinite $A\geq 0$ the formal trace
$\operatorname{Tr}A\coloneqq\sum_{i}\langle i|A|i\rangle$ can be shown not to
depend on the chosen orthonormal basis $\\{|i\rangle\\}_{i}$ of $\mathcal{H}$.
operator defined in a Hilbert space has an equivalent representation in terms
of $\bm{S}$-ordered Phase-space Quasi-probability Distributions
[($\bm{S}$)-PQDs], which are distributions defined over a real symplectic
space that is referred to as phase-space. In order to characterize this
representation, we first define the frame operators. These operators are
denoted by $\Delta(\bm{z})$ and are often used to define the $\bm{S}$-ordered
quasi-probability distributions
$W({O},\bm{z})\coloneqq\operatorname{tr}[{O}\Delta(\bm{z})]$ corresponding to
the operator ${O}$, that can either be an observable or a density operator in
the Hilbert space222Throughout this paper we use a notation in which phase
space related quantities are in bold font, e.g., ${\bm{\sigma}}$ and
${\bm{d}}$ represent the covariance matrix and displacement vector. On the
contrary, we write the operators in the Hilbert space of the system in normal
font, e.g., ${\rho}$ represents the density matrix. For the positive operator-
valued measure (POVM) elements like ${M_{{\bm{a}}|{s}}}$ the operator itself
is in the Hilbert space and represented with normal font, but the outcome is
often a phase-space vector like ${\bm{a}}$ thus we use the bold font..

First we note that, for a given quantum state $\rho$ and any set of POVM
operators $\\{{M}_{{{\bm{a}}|s}}\\}$—with ${s}$ labeling the specific choice
of POVM and ${\bm{a}}$ the different outcomes for that given choice—the
outcome probabilities for the measurements are given by the Born rule,

$P(\bm{a}|\rho;s)=\operatorname{tr}[\rho{M}_{\bm{a}|{s}}].$ (1)

To show how any $m$-mode operator like $\rho$ can be represented by
($\bm{S}$)-PQDs in phase-space, we start by expanding $\rho$ in terms of the
canonical displacement operators,

$\displaystyle\rho$ $\displaystyle=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{y}\operatorname{tr}[\rho D(\bm{y})]D(-\bm{y})$ (2)

where $D(\bm{y})\coloneqq\exp(-i\bm{y}\bm{\Omega}R^{\intercal})$ is the
displacement (or Weyl) operator, $R=(q_{1},p_{1},\dots,q_{m},p_{m})$ is the
vector of canonical quadrature operators, $[q_{j},p_{k}]=i\delta_{j,k}$, and
$\bm{y},\bm{z}\in\mathbb{R}^{2m}$, and
$\bm{\Omega}=\bigoplus_{j=1}^{m}\left(\begin{smallmatrix}0&1\\\
-1&0\end{smallmatrix}\right)$ is the symplectic form enforcing the uncertainty
relations in phase-space.

The $\Delta(\bm{z})$ operators are defined as the Fourier transforms of the
displacement operator

$\Delta(\bm{z})\coloneqq\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,D(\bm{y})\,e^{-i\bm{z}\bm{\Omega}\bm{y}^{\intercal}},$
(3)

where the integral is intended to converge in the weak operator topology. The
Wigner function for the density operator $\rho$ is then given by

$\displaystyle W(\bm{z}|\rho)$
$\displaystyle=\operatorname{tr}\big{[}\rho\Delta(\bm{z})\big{]}=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,\operatorname{tr}[\rho
D(\bm{y})]\,e^{-i\bm{z}\bm{\Omega}\bm{y}^{\intercal}}.$ (4)

Since $\operatorname{tr}[D(\bm{y})]=(2\pi)^{m}\delta^{2m}(\bm{y})$, the
operator $\Delta(\bm{z})$ satisfies the following relations:

$\operatorname{tr}[\Delta(\bm{z})]=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,\operatorname{tr}[D(\bm{y})]\,e^{-i\bm{z}\bm{\Omega}\bm{y}^{\intercal}}=\frac{1}{(2\pi)^{m}}$
(5)

and

$\int
d^{2m}\bm{z}\,\Delta(\bm{z})=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,D(\bm{y})\,\int
d^{2m}\bm{z}e^{-i\bm{z}\bm{\Omega}\bm{y}^{\intercal}}=I,$ (6)

where we have used that

$\displaystyle\int
d^{2m}\bm{z}e^{-i\bm{z}\bm{\Omega}\bm{y}^{\intercal}}=(2\pi)^{2m}\delta^{2m}(\bm{y}).$
(7)

By taking the inverse Fourier transform of Eq. (4) and using Eq. (7), one
obtains

$\displaystyle\int
d^{2m}\bm{z}\,W(\bm{z}|\rho)e^{i\bm{z}\bm{\Omega}\bm{y}^{\intercal}}=\operatorname{tr}[\rho
D(\bm{y})].$ (8)

Substituting Eq. (8) into Eq. (2) gives

$\displaystyle\begin{split}\rho&=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{z}\,W(\bm{z}|\rho)\int
d^{2m}\bm{y}D(-\bm{y})\,e^{i\bm{y}\Omega\bm{y}^{\intercal}}\\\
&=(2\pi)^{2m}\int d^{2m}\bm{z}\,W(\bm{z}|\rho)\Delta(\bm{z}).\end{split}$ (9)

Finally, by using these results in Eq. (1), we obtain

$P(\bm{a}|\rho;s)=\operatorname{tr}[\rho{M}_{\bm{a}|{s}}]=(2\pi)^{M}\int
d^{\,2m}\\!\bm{z}\,W(\bm{z}|\rho)\,W(\bm{a}|{s},\bm{z}),$ (10)

where
$W(\bm{a}|{s},\bm{z})=\operatorname{tr}\big{[}{M}_{\bm{a}|{s}}\Delta(\bm{z})\big{]}$
is the Wigner function associated to the POVM measurement operators.333The
only measurement operators we will encounter will be of trace class. From this
relation, we see that the information contained in the probabilities given by
Eq. (1) is equivalently contained in the integral overlap of the corresponding
Wigner functions.

####  II.1.1 The overlap of two systems

For any two quantum states $\rho_{1}$ and $\rho_{2}$ we can use (2) to compute
their overlap

$\displaystyle\operatorname{tr}[\rho_{1}\rho_{2}]$
$\displaystyle=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{y}\operatorname{tr}[\rho_{1}D(\bm{y})]\operatorname{tr}[\rho_{2}D(-\bm{y})],$
(11)

Now, the term $\operatorname{tr}[D(\bm{y})\rho_{i}]$ is nothing but the
characteristic function of the $i$th system. That is, by definition,

$\displaystyle\chi(\bm{y}|\rho_{i})\coloneqq\operatorname{tr}[D(\bm{y})\rho_{i}].$
(12)

For Gaussian systems $\rho_{i}$ with covariance matrices ${\bm{\sigma}}_{i}$
and displacement operators ${\bm{d}}_{i}$, the characteristic function is
given by the Gaussian distribution

$\displaystyle\chi(\bm{y}|\rho_{i})=e^{-\frac{1}{2}\bm{y}\bm{\Omega}\bm{\sigma}_{i}\bm{\Omega}\bm{y}^{\intercal}+i\bm{d}_{i}\bm{\Omega}\bm{y}^{\intercal}},$
(13)

with the displacement vector and the covariance matrix (CM) defined as
$\bm{d}_{i}|_{k}=\left<R_{k}\right>_{\rho_{i}}$, and
$\bm{\sigma}_{i}|_{kl}=\left<\\{R_{k},R_{l}\\}_{+}\right>_{\rho_{i}}/2-\left<R_{k}\right>_{\rho_{i}}\left<R_{l}\right>_{\rho_{i}}$,
where $\left<A\right>_{\rho_{i}}\coloneqq\operatorname{tr}[A\rho_{i}]$ and
$\\{,\\}_{+}$ represents the anti-commutator. Having this form for the
Gaussian characteristic function, we can write

$\displaystyle\operatorname{tr}[\rho_{1}\rho_{2}]$
$\displaystyle=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{y}\operatorname{tr}[\rho_{1}D(\bm{y})]\operatorname{tr}[\rho_{2}D(-\bm{y})]$
$\displaystyle=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{y}e^{-\frac{1}{2}\bm{y}\bm{\Omega}(\bm{\sigma}_{1}+\bm{\sigma}_{2})\bm{\Omega}\bm{y}^{\intercal}+i(\bm{d}_{1}-\bm{d}_{2})\bm{\Omega}\bm{y}^{\intercal}}$
$\displaystyle=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{y}e^{-\frac{1}{2}\bm{y}{\tilde{\bm{\sigma}}}\bm{y}^{\intercal}+i{\tilde{\bm{d}}}\bm{y}^{\intercal}},$
(14)

where we define
${\tilde{\bm{\sigma}}}\coloneqq\bm{\Omega}(\bm{\sigma}_{1}+\bm{\sigma}_{2})\bm{\Omega}$,
and ${\tilde{\bm{d}}}\coloneqq(\bm{d}_{1}-\bm{d}_{2})\bm{\Omega}$. In order to
solve this last integral, we need to find the orthogonal basis of
${\tilde{\bm{\sigma}}}$, such that the integration variables decouple.
Performing this elementary transformation yields

$\displaystyle\operatorname{tr}[\rho_{1}\rho_{2}]$
$\displaystyle=\frac{1}{(2\pi)^{m}}\int
d^{2m}\bm{y}e^{-\frac{1}{2}\bm{y}{\tilde{\bm{\sigma}}}\bm{y}^{\intercal}+i{\tilde{\bm{d}}}\bm{y}^{\intercal}}=\frac{1}{(2\pi)^{m}}\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=15.41556pt,depth=-12.3325pt}}}{{\hbox{$\textstyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=11.9pt,depth=-9.52003pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=9.9pt,depth=-7.92003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=9.9pt,depth=-7.92003pt}}}e^{-\frac{1}{2}\tilde{\bm{d}}{\tilde{\bm{\sigma}}}^{-1}\tilde{\bm{d}}^{\intercal}}$
$\displaystyle=\frac{1}{\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}(\bm{\sigma}_{1}+\bm{\sigma}_{2})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}(\bm{\sigma}_{1}+\bm{\sigma}_{2})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}(\bm{\sigma}_{1}+\bm{\sigma}_{2})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}(\bm{\sigma}_{1}+\bm{\sigma}_{2})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}}e^{-\frac{1}{2}{(\bm{d}_{1}-\bm{d}_{2})}{(\bm{\sigma}_{1}+\bm{\sigma}_{2})}^{-1}{(\bm{d}_{1}-\bm{d}_{2})^{\intercal}}}.$
(15)

The Gaussian measurement operators of interest here can be written as
displaced Gaussian states $\rho_{s}^{M}$, with the displacement vector and the
covariance matrix given by ${\bm{d}}_{s}^{M}$ and ${\bm{\sigma}}_{s}^{M}$,
respectively.444When necessary we use the superscript $M$ to refer to
measurement operators, displacement vectors, and covariance matrices, to
prevent confusion with those of the system. We can write

$\displaystyle{M}_{\bm{a}|{s}}=D(\bm{a})\rho_{s}^{M}D^{\dagger}(\bm{a})/(2\pi)^{m}.$
(16)

Our aim is then to find the characteristic function of this Gaussian POVM
element. By using the cyclic property of the trace, and the properties of the
displacement operator
$D(\bm{y})D(\bm{a})=\exp[-i\bm{y}\bm{\Omega}\bm{a}^{\intercal}/2]D(\bm{y}+\bm{a})$,
and $D^{\dagger}(\bm{y})=D(-\bm{y})$, we have that

$\displaystyle\chi(\bm{y}|{M}_{\bm{a}|{s}})$
$\displaystyle=\operatorname{tr}[D(\bm{a})\rho_{s}^{M}D^{\dagger}(\bm{a})D(\bm{y})]/(2\pi)^{m}$
$\displaystyle=\frac{e^{i\bm{a}\bm{\Omega}\bm{y}^{\intercal}}}{(2\pi)^{m}}\chi(\bm{y}|\rho_{s}^{M}),$
(17)

which means that—up to a $1/(2\pi)^{m}$ constant—${M}_{\bm{a}|{s}}$ has a
Gaussian characteristic function with the same covariance matrix as
$\rho_{s}^{M}$, but with displacement vector increased by the amount $\bm{a}$
referring to the measurement outcome. By putting the pieces together, we have
that, for a Gaussian measurement performed on a Gaussian state $\rho$ with
displacement vector $\bm{d}$ and CM $\bm{\sigma}$, the outcome probability
distribution takes the (also Gaussian) form

$\displaystyle
P(\bm{a}|\rho;s)=\operatorname{tr}[\rho{M}_{\bm{a}|{s}}]=\frac{1}{(2\pi)^{m}\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}}e^{-\frac{1}{2}{(\bm{d}_{s}^{M}+\bm{a}-\bm{d})}{({\bm{\sigma}_{s}^{M}}+\bm{\sigma})}^{-1}{(\bm{d}_{s}^{M}+\bm{a}-\bm{d})^{\intercal}}},$
(18)

where ${\bm{d}}_{s}^{M}$ and ${\bm{\sigma}}_{s}^{M}$ are respectively the
displacement vector and the covariance matrix of ${\rho}_{s}^{M}$, i.e., the
Gaussian state of the measurement. The above is the most general form of the
probability distribution of the outcome of a Gaussian measurement on a
Gaussian state; however, in what follows we shall assume for simplicity that
both $\rho$ and $\rho_{s}^{M}$ have vanishing displacements, so that

$\displaystyle
P(\bm{a}|\rho;s)=\frac{1}{(2\pi)^{m}\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}}e^{-\frac{1}{2}\bm{a}{({\bm{\sigma}_{s}^{M}}+\bm{\sigma})}^{-1}\bm{a}^{\intercal}}.$
(19)

The normalisation of this probability can be checked easily

$\displaystyle\int d^{2m}\bm{a}P(\bm{a}|\rho;s)$
$\displaystyle=\frac{1}{(2\pi)^{m}\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}}\int
d^{2m}\bm{a}e^{-\frac{1}{2}\bm{a}{({\bm{\sigma}_{s}^{M}}+\bm{\sigma})}^{-1}\bm{a}^{\intercal}}.$
$\displaystyle=\frac{1}{(2\pi)^{m}\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}}\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})^{-1}}\,}$}\lower 0.4pt\hbox{\vrule
height=17.70667pt,depth=-14.16539pt}}}{{\hbox{$\textstyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})^{-1}}\,}$}\lower 0.4pt\hbox{\vrule
height=12.43332pt,depth=-9.9467pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})^{-1}}\,}$}\lower 0.4pt\hbox{\vrule
height=9.43332pt,depth=-7.54669pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{(2\pi)^{2m}}{{\rm
det}({\bm{\sigma}_{s}^{M}}+\bm{\sigma})^{-1}}\,}$}\lower 0.4pt\hbox{\vrule
height=9.43332pt,depth=-7.54669pt}}}=1.$ (20)

See the appendix for a generalisation of Gaussian measurements that
incorporates noise and auxiliary modes.

##  III Thermometry of Gaussian quantum systems: Optimal vs feasible
measurements

As already advanced, we are interested in the measurement part of thermometry
protocols. Let us first recall some basics of metrology in quantum systems.
Given a quantum system we aim to estimate a parameter like temperature $T$
imprinted on its density matrix denoted by $\rho(T)$. To this end, one
performs a POVM measurement with positive elements $\\{M_{{\bm{a}}|{s}}\\}$,
where $s$ labels the specific measurement and ${\bm{a}}$ the different
outcomes for that specific measurement. We have $\int
d{\bm{a}}M_{{\bm{a}}|{s}}=I$ for all $s$, with $I$ being the identity
operator. The outcomes of the measurement are then mapped into an estimate of
the true parameter by using an estimator function $T_{\rm est}({\bm{a}})$. One
can quantify the error on such an estimation process by means of the mean
square error $\delta^{2}T_{\rm est}(\rho;s)\coloneqq\int
d{\bm{a}}P(\bm{a}|\rho;s)(T_{\rm est}({\bm{a}})-T)^{2}$, with
$P(\bm{a}|\rho;s)=\operatorname{tr}[M_{{\bm{a}}|{s}}\rho(T)]$ being the
probability of observing the outcome ${\bm{a}}$ conditioned on measuring $s$.
Notice that the estimation error depends on the performed measurement. It is
well known that for any unbiased estimator, that is, for any estimator
satisfying $\int d{\bm{a}}P(\bm{a}|\rho;s)T_{\rm est}({\bm{a}})=T$, the mean
square error is lower bounded by the Cramér-Rao bound

$\displaystyle\delta^{2}T_{\rm est}(\rho;s)\geqslant\frac{1}{m{\cal F}^{\rm
C}(\rho;s)},$ (21)

with ${\cal F}^{\rm C}(\rho;s)$ being the (classical) Fisher information given
by

$\displaystyle{\cal F}^{\rm C}(\rho;s)\coloneqq\left<\left[\partial_{T}\log
P(\bm{a}|\rho;s)\right]^{2}\right>=\int
d{\bm{a}}P(\bm{a}|\rho;s)\left[\partial_{T}\log P(\bm{a}|\rho;s)\right]^{2}.$
(22)

The Cramér-Rao bound can be saturated in the asymptotic limit by choosing the
maximum likelihood estimator Newey and McFadden (1994). A good candidate to
quantify the precision of measurements is therefore given by the Fisher
information; the larger, the more precise the measurement.

It is a natural and fundamental question to ask: what is the ultimate bound on
precision? This is obtained by finding the measurement that minimises Eq.
(21). The corresponding precision is quantified by the so-called quantum
Fisher information (QFI), which can be formally defined as

$\displaystyle{\cal F}^{\rm Q}(\rho)\coloneqq\max_{s}{\cal F}^{\rm
C}(\rho;s)=\operatorname{tr}[\rho\Lambda^{2}].$ (23)

Here, $\Lambda$ is the symmetric logarithmic derivative (SLD), a Hermitian
operator defined by

$\displaystyle 2\partial_{T}\rho(T)=\Lambda\rho(T)+\rho(T)\Lambda.$ (24)

The SLD also characterizes the optimal measurement, which turns out to be
projective in the basis of $\Lambda$.

###  III.1 Ultimate bounds on thermometry of Gaussian systems

The problem of finding the ultimate bounds in parameter estimation of Gaussian
quantum systems was first addressed in Monras (2013), further explored in
Jiang (2014), and extended to multiple parameter estimation in Nichols  _et
al._ (2018). It was shown that the optimal measurement can be always written
as a second order operator in terms of the quadratures, that is

$\displaystyle\Lambda=\sum_{kl}{\bm{C}}^{(2)}_{kl}R_{k}R_{l}+\sum_{k}{\bm{C}}^{(1)}_{k}R_{k}+{\bm{C}}^{(0)},$
(25)

where ${\bm{C}}^{(2)}$ is a $2N\times 2N$ matrix of real numbers,
${\bm{C}}^{(1)}$ is a $2N$ dimensional real vector, and ${\bm{C}}^{(0)}$ is a
real number. These three quantities are given by the solution of the following
equations:

$\displaystyle\partial_{T}{\bm{\sigma}}$
$\displaystyle=2{\bm{\sigma}}{\bm{C}}^{(2)}{\bm{\sigma}}+\frac{1}{2}{\bm{\Omega}}{\bm{C}}^{(2)}{\bm{\Omega}},$
(26) $\displaystyle{\bm{C}}^{(1)}$
$\displaystyle={\bm{\sigma}}^{-1}\partial_{T}{\bm{d}}-2{\bm{C}}^{(2)}{\bm{d}},$
(27) $\displaystyle{\bm{C}}^{(0)}$
$\displaystyle=-{\bm{C}}^{(1)T}{\bm{d}}-\frac{1}{2}\operatorname{Tr}[{\bm{C}}^{(2)}{\bm{\sigma}}]-{\bm{d}}^{\intercal}{\bm{C}}^{(2)}{\bm{d}}.$
(28)

In turn, the quantum Fisher information reads

$\displaystyle{\cal F}^{\rm
Q}({\bm{\sigma}},{\bm{d}})=\partial_{T}{\bm{d}}^{T}{\bm{\sigma}}^{-1}\partial_{T}{\bm{d}}+2\operatorname{Tr}[{\bm{C}}^{(2)}{\bm{\sigma}}{\bm{C}}^{(2)}{\bm{\sigma}}+\frac{1}{4}{\bm{C}}^{(2)}{\bm{\Omega}}{\bm{C}}^{(2)}{\bm{\Omega}}].$
(29)

Despite its elegant and simple form, the optimal measurement (25) is,
unfortunately, often an experimentally challenging one to perform. Thus, we
will explore more feasible measurements in the next sections. It is worth
mentioning that the problem of phase estimation by means of Gaussian
measurements has been addressed in Oh  _et al._ (2019) for a single mode,
while using Gaussian measurements in the Bayesian approach to metrology was
considered in Morelli  _et al._ (2021); Kundu  _et al._ (2021).

###  III.2 Thermometry using Gaussian measurements

Gaussian measurements are an important family of experimentally realisable
measurements. They have been defined formally in section II.1; alternatively,
they can be characterized as all those measurements that are realizable by
appending ancillas in the vacuum state, applying arbitrary Gaussian unitaries,
and performing homodyne detections Giedke and Cirac (2002). We can
characterize them by means of their covariance matrix ${\bm{\sigma}}_{s}^{M}$.
Our first observation is that the optimal measurement should have a covariance
matrix representing a pure quantum state.

####  III.2.1 Observation.—Pure Gaussian measurements are always better

To prove this claim, we show that any non-pure Gaussian measurement can be
cast as classical post-processing of another pure Gaussian measurement,
therefore, not increasing its (classical Fisher) information.

The POVM elements of a Gaussian measurement can be written in the form

${M}_{\bm{a}|{s}}=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,D\big{(}\bm{y}\big{)}e^{{-\bm{y}{\bm{\sigma}_{s}^{M}}\bm{y}^{\intercal}/2-i\bm{a}\bm{\Omega}\bm{y}^{\intercal}}}.$
(30)

If the measurement is not pure, then there exists a pure covariance matrix
$\bm{\sigma}_{t}^{M}$ s.t. ${\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M}\geqslant
0$. This follows directly from Williamson’s theorem Williamson (1936) (see
also (Lami  _et al._ , 2017, Lemma 5) for a direct proof). The POVM
measurement elements corresponding to this pure CM are

${M}_{\bm{b}|t}=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,D\big{(}\bm{y}\big{)}e^{{-\bm{y}\bm{\sigma}_{t}^{M}\bm{y}^{\intercal}/2-i\bm{b}\bm{\Omega}\bm{y}^{\intercal}}}.$
(31)

If ${M}_{\bm{a}|{s}}$ can be reproduced by ${M}_{\bm{b}|t}$, then we should
find the probability measure $\pi({\bm{a}}|{\bm{b}})$ s.t.,

$\displaystyle{M}_{\bm{a}|{s}}=\int\frac{d^{\,2m}\\!\bm{b}\,}{(2\pi)^{m}}\pi({\bm{a}}|{\bm{b}}){M}_{\bm{b}|t}~{}~{}\forall{\bm{a}}.$
(32)

By combining (31) and (32) we obtain

$\displaystyle{M}_{\bm{a}|{s}}=\iint\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\frac{d^{\,2m}\\!\bm{b}\,}{(2\pi)^{m}}\pi({\bm{a}}|{\bm{b}})\,D\big{(}\bm{y}\big{)}e^{{-\bm{y}\bm{\sigma}_{t}^{M}\bm{y}^{\intercal}/2-i\bm{b}\bm{\Omega}\bm{y}^{\intercal}}}.$
(33)

If we compare with (30), we should have

$\displaystyle\int\frac{d^{\,2m}\\!\bm{b}\,}{(2\pi)^{m}}\pi({\bm{a}}|{\bm{b}})e^{-i\bm{b}\bm{\Omega}\bm{y}^{\intercal}}=e^{{-\bm{y}({\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M})\bm{y}^{\intercal}/2-i\bm{a}\bm{\Omega}\bm{y}^{\intercal}}}.$
(34)

Finally by taking the inverse Fourier transform we find

$\displaystyle\pi({\bm{a}}|{\bm{b}})$
$\displaystyle=\frac{e^{\frac{1}{2}({\bm{b}}-{\bm{a}}){\bm{\Omega}}({\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M})^{-1}{\bm{\Omega}}^{\intercal}({\bm{b}}-{\bm{a}})^{\intercal}}}{\mathchoice{{\hbox{$\displaystyle\sqrt{\det({\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M})\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{\det({\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M})\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{\det({\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M})\,}$}\lower
0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det({\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M})\,}$}\lower
0.4pt\hbox{\vrule height=3.75pt,depth=-3.00002pt}}}}.$ (35)

Note that the fact that the Fourier transform exists owes to the fact that
${\bm{\sigma}_{s}^{M}}-\bm{\sigma}_{t}^{M}\geqslant 0$.

In this manner, we have proven that the optimal Gaussian measurements should
be chosen from the pool of pure covariance matrices, which can always be
written as ${\bm{\sigma}}_{s}^{M}={\bm{S}}^{M}({\bm{S}}^{M})^{\intercal}$,
with ${\bm{S}}^{M}$ being a symplectic transformation.

####  III.2.2 Classical Fisher information of Gaussian measurements

As we saw in section II.1, when performing a Gaussian measurement on Gaussian
systems the resulting probability distribution, given by (18), is also
Gaussian. It is charachterized by its covariance matrix
${\bm{V}}\coloneqq{\bm{\sigma}}_{s}^{M}+{\bm{\sigma}}$, with
${\bm{\sigma}}_{s}^{M}$ and ${\bm{\sigma}}$ being the covariance matrices of
the measurement and the system, respectively. In what follows we focus on
systems with zero first-order moments. This is motivated by the fact that, at
thermal equilibrium, one can always change the basis with a temperature-
independent transformation to one with zero displacements. Moreover, in most
non-thermal scenarios the quadratic nature of the system-sample Hamiltonian is
such that the displacement vector of the system vanishes Mehboudi  _et al._
(2019b); Correa  _et al._ (2017); Planella  _et al._ (2020); Lampo  _et al._
(2017); Khan  _et al._ (2021). The Fisher information of classical Gaussian
distributions is well known, and we also present a simple proof for it in
appendix. The final expression is given by

$\displaystyle{\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{s}^{M})$
$\displaystyle=\frac{1}{2}\operatorname{Tr}\left[\left(({\bm{\sigma}}+{\bm{\sigma}}_{s}^{M})^{-1}\partial_{T}{\bm{\sigma}}\right)^{2}\right],$
(36)

where we explicitly include in the argument that the classical Fisher
information is obtained for measurement ${\bm{\sigma}}_{s}^{M}$ conditioned on
the covariance matrix of the system being ${\bm{\sigma}}$.

####  III.2.3 Homodyne detection and equivalence with the error propagation

Homodyne detection is a projective measurement in the basis of one of the
quadratures. For example, consider a single-mode Gaussian system, and assume
that we want to measure it in the position basis. The measurement covariance
matrix can be written as
${\bm{\sigma}}_{R_{1}}^{M}\coloneqq\lim_{r\to\infty}{\rm diag}[1/r,r]$. By
replacing in (36) one finds the classical Fisher Information for the position
measurement

$\displaystyle{\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{R_{1}}^{M})$
$\displaystyle=\lim_{r\to\infty}\frac{1}{2}\operatorname{Tr}\left[\left(({\bm{\sigma}}+{\rm
diag}[1/r,r])^{-1}\partial_{T}{\bm{\sigma}}\right)^{2}\right]=\frac{|\partial_{T}{\bm{\sigma}}_{11}|^{2}}{2{\bm{\sigma}}_{11}^{2}}.$
(37)

Often times, the inverse of the error propagation is an alternative way of
characterising the precision. The error propagation characterizes the
estimation error when using some observable $O$ and is defined as

$\displaystyle\delta^{2}T(\rho;O)_{\rm EP}\coloneqq\frac{{\rm
Var}(O)}{|\partial_{T}\left<O\right>|^{2}}.$ (38)

Regarding position measurement, if we replace the observable with the square
of the position operator, i.e., $O\to R_{1}^{2}$, then by using Wick’s theorem
and Eqs. (37) and (38) one finds that

$\displaystyle\delta^{2}T(\rho;R_{1}^{2})_{\rm EP}=\frac{1}{{\cal F}^{\rm
C}({\bm{\sigma}};{\bm{\sigma}}_{R_{1}}^{M})}.$ (39)

This equivalence between the error propagation of the squared quadratures and
the inverse of classical Fisher information holds for any parameter estimation
problem. Needless to say that, in case we were to measure the momentum, we
just need to replace $R_{1}\to R_{2}$ and
${\bm{\sigma}}_{11}\to{\bm{\sigma}}_{22}$. Note that this result holds
regardless of the basis, and in particular for rotated quadrature measurements
as well; one just needs to rewrite the covariance matrix of the system
${\bm{\sigma}}$ in the rotated basis. The equivalence between the error
propagation and the inverse of the classical Fisher information does not
necessarily hold in multimode scenario. Nonetheless, for the trivial case with
a tensor product state in a basis that does not depend on the parameter—where
the classical Fisher information becomes additive—the same result holds. That
is to say if the system covariance matrix satisfies
${\bm{\sigma}}=\oplus{\bm{\sigma}}^{k}$, and
$\partial_{T}{\bm{\sigma}}=\oplus_{k}\partial_{T}{\bm{\sigma}}^{k}$, then any
local measurement in the form of
${\bm{\sigma}}_{s}^{M}=\oplus_{k}{\bm{\sigma}}_{\rm k,HD}^{M}$—where
${\bm{\sigma}}_{\rm k,HD}^{M}$ is an arbitrary Homodyne detection on the $k$th
mode—has a classical Fisher information that is equivalent to the inverse of
the error propagation for $O=\otimes_{k}\left(R_{\rm k,HD}\right)^{2}$—where
$R_{\rm k,HD}$ is the local quadrature in the Hilbert space of the $k$th mode
corresponding to ${\bm{\sigma}}_{\rm k,HD}^{M}$.

####  III.2.4 Heterodyne measurement

A Heterodyne measurement can be represented—in the single-mode case—by the
covariance matrix ${\bm{\sigma}}_{\rm Het}^{M}={\rm diag}[1,1]$, which remains
the same under local rotation of the quadratures. One can find the
corresponding classical Fisher information by replacing this in Eq. (36). If
the system’s covariance matrix and its derivative can be diagonalised in the
same basis, i.e., ${\bm{\sigma}}={\rm
diag}[{\bm{\sigma}}_{11},{\bm{\sigma}}_{22}]$, and
$\partial_{T}{\bm{\sigma}}={\rm
diag}[\partial_{T}{\bm{\sigma}}_{11},\partial_{T}{\bm{\sigma}}_{22}]$ then the
classical Fisher information gets a simple form

$\displaystyle{\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{\rm
Het}^{M})=\frac{1}{2}\left(\left[\frac{\partial_{T}{\bm{\sigma}}_{11}}{1+{\bm{\sigma}}_{11}}\right]^{2}+\left[\frac{\partial_{T}{\bm{\sigma}}_{22}}{1+{\bm{\sigma}}_{22}}\right]^{2}\right).$
(40)

If the system is at thermal equilibrium, its covariance matrices is diagonal
in a basis that is independent of temperature, and therefore the above
equation is useful. Otherwise, if one deals with non-thermal states (such as,
e.g., steady states that are non-thermalMehboudi  _et al._ (2019b); Correa
_et al._ (2017); Planella  _et al._ (2020); Lampo  _et al._ (2017); Khan  _et
al._ (2021) or dynamical non-thermal states Correa  _et al._ (2015)), one
should avoid using this, and directly use Eq. (36).

####  III.2.5 The optimal Gaussian measurement

Given a covariance matrix and its derivative, it is a natural question to find
the Gaussian measurement with the maximum Fisher information

$\displaystyle{\bm{\sigma}}_{\rm max}^{M}$ $\displaystyle\coloneqq{\rm
arg}\max_{{\bm{\sigma}}_{s}^{M}}\frac{1}{2}\operatorname{Tr}\left[\left(({\bm{\sigma}}+{\bm{\sigma}}_{s}^{M})^{-1}\partial_{T}{\bm{\sigma}}\right)^{2}\right],$
(41) $\displaystyle{\bm{\sigma}}_{\rm min}^{M}$ $\displaystyle\coloneqq{\rm
arg}\min_{{\bm{\sigma}}_{s}^{M},{\rm
pure}}\frac{1}{2}\operatorname{Tr}\left[\left(({\bm{\sigma}}+{\bm{\sigma}}_{s}^{M})^{-1}\partial_{T}{\bm{\sigma}}\right)^{2}\right].$
(42)

As already advanced, ${\bm{\sigma}}_{\rm max}^{M}$ should represent a pure
state. We also defined ${\bm{\sigma}}_{\rm min}^{M}$ as the covariance matrix
of the pure measurement with minimum Fisher information. If the covariance
matrix of the system and its derivative are proportional to identity (e.g., if
we have a thermal state), then it is easy to find the optimal Gaussian
measurement in the single-mode case. Let ${\bm{\sigma}}=\nu{\bm{I}}_{2}$, and
$\partial_{T}{\bm{\sigma}}={\nu}^{\prime}{\bm{I}}_{2}$. Then, we also note
that the measurement covariance matrix can be cast as
${\bm{\sigma}}_{r}^{M}={\bm{O}}{\bm{r}}^{M}{\bm{O}}^{T}$, with ${\bm{O}}$
being an orthogonal transformation, and ${\bm{r}}^{M}={\rm diag}[{r},1/{r}]$
with the squeezing parameter $r\in[0,\infty)$. The classical Fisher
information does not depend on the orthogonal rotation ${\bm{O}}$, and reads

$\displaystyle{\cal F}^{\rm
C}({\bm{\sigma}};{\bm{\sigma}}_{r}^{M})=\frac{(\nu^{\prime})^{2}}{2}\left[\frac{1}{(\nu+r)^{2}}+\frac{1}{(\nu+1/r)^{2}}\right].$
(43)

Note that this does not change by the transformation $r\to 1/r$, and
therefore, we can limit the domain to $r\in[1,\infty)$. We can find the value
of $r$ corresponding to the minimum and maximum of the Fisher information by
simple algebra, which is summarised in Table 1. For example, the maximum is
given by the expression

${\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{\max}^{M})=\max_{r\geqslant
1}{\cal F}^{\rm
C}({\bm{\sigma}};{\bm{\sigma}}_{r}^{M})={\nu^{\prime}}^{2}\max\left\\{\frac{1}{2\nu^{2}},\,\frac{1}{(\nu+1)^{2}}\right\\}.$
(44)

It can be seen that the optimal measurement is either Homodyne or Heterodyne
detection. The freedom in the orthogonal rotation ${\bm{O}}$ means that one
can perform these measurements in an arbitrary basis.

If the state is at thermal equilibrium, then $\nu=\coth(\omega/2T)$, with
$\omega$ being the energy of the Bosonic mode and $T$ being the temperature.
The fact that the optimal Gaussian measurement depends on temperature is
interesting; on the contrary, the optimal non-Gaussian measurement is always
performed in the basis of the system Hamiltonian, regardless of the
temperature. Furthermore, our result shows not only that at high temperatures
the Heterodyne detection is the best Gaussian measurement, but also that its
performance approaches that of the optimal measurement. To see this, note that
the quantum Fisher information of a single Bosonic mode (that is a Harmonic
oscillator) at thermal equilibrium is given by Correa  _et al._ (2015)

$\displaystyle{\cal F}^{\rm Q}({\rm
HO})=\frac{\omega^{2}}{4T^{4}\sinh^{2}(\omega/2T)},$ (45)

and thus at $T\to\infty$ the leading order behaves as ${\cal F}^{\rm Q}({\rm
HO})=1/T^{2}$. Similarly, one can see that at the same asymptotic limit, the
leading order of the Fisher information of Heterodyne measurement behaves as
${\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{\rm max}^{M})={\cal F}^{\rm
C}({\bm{\sigma}};{\bm{\sigma}}_{\rm Het}^{M})=1/T^{2}$.

Consideration | $r_{\min}$ | ${\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{\rm min}^{M})$ | $r_{\max}$ | ${\cal F}^{\rm C}({\bm{\sigma}};{\bm{\sigma}}_{\rm max}^{M})$
---|---|---|---|---
$1\leqslant\nu\leqslant 2$ | $1$ | $\frac{\nu^{\prime 2}}{(1+\nu)^{2}}$ | $\infty$ | $\frac{\nu^{\prime 2}}{2\nu^{2}}$
$2\leqslant\nu\leqslant 1+\mathchoice{{\hbox{$\displaystyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}$ | $\frac{\nu(\nu^{2}-3)}{2}$ | $\frac{\nu^{\prime 2}(\nu^{2}-2)}{2(\nu^{2}-1)}$ | $\infty$ | $\frac{\nu^{\prime 2}}{2\nu^{2}}$
$1+\mathchoice{{\hbox{$\displaystyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{2\,}$}\lower 0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}\leqslant\nu$ | $\frac{\nu(\nu^{2}-3)}{2}$ | $\frac{\nu^{\prime 2}(\nu^{2}-2)}{2(\nu^{2}-1)}$ | 1 | $\frac{\nu^{\prime 2}}{(\nu+1)^{2}}$
Table 1: Maximal and minimal values of the CFI ${\cal F}^{\rm
C}({\bm{\sigma}};{\bm{\sigma}}_{s}^{M})$ for Gaussian measurements in one-mode
thermal states parametrized by $\nu$, and the measurements that achieve them
with the squeezing parameter $r$. Only pure measurements are considered. The
optimal measurement is either Homodyne detection (for the low-temperature
regime) or Heterodyne detection (for high temperatures). This is contrary to
the optimal [non-Gaussian] measurement, which is a projective measurement in
the energy basis regardless of the temperature. Figure 2: Temperature
estimation of single-mode thermal systems at thermal equilibrium. We compare
the ultimate limit on thermometry precision posed by the Quantum Fisher
Information (dashed black line) and the precision achieved by the optimal
Gaussian measurement (solid red), the worst Gaussian measurement (dotted red)
and single-mode Photo-Detection (blue) versus the temperature $T$. The shaded
area represents the precision achievable by the set of single-mode Gaussian
measurement. We also show the performance of Photo-Detection including
imperfect detectors, characterized by the parameter $e$ (dashed and dotted
blue lines).

In Figure 2 the Fisher information of Gaussian measurements is depicted versus
temperature. We compare this to the quantum Fisher information. Unlike in the
high-temperature regime, where Gaussian measurements (Heterodyne specifically)
are close to optimal, at low temperatures they perform rather poorly. In
Section III.3 we will see that photo-detection measurements are a viable
alternative in this latter case, and specifically that they perform close to
optimal at low temperatures.

####  III.2.6 Optimal Gaussian measurements on multimode systems

When dealing with a multimode scenario, the problem becomes harder, and we
were unable to derive an analytical solution. The main difficulty arises from
the fact that the objective function to be optimised, i.e., the Fisher
information, is a nonlinear function of the measurement covariance matrix.
Nonetheless, one can design an algorithm to find the optimal measurement
numerically. To begin with, one can use the fact that the covariance matrix of
the optimal measurement is always pure and therefore belongs to the family of
symmetric symplectic transformations, i.e., ${\bm{\sigma}}_{\rm
max}^{M}\in\\{\bar{\bm{S}}^{M}\coloneqq{\bm{S}}^{M}({\bm{S}}^{M})^{\intercal}|{\bm{S}}^{M}{\bm{\Omega}}({\bm{S}}^{M})^{\intercal}={\bm{\Omega}}\\}$,
where ${\bar{\bm{S}}^{M}}$ is a symmetric symplectic transformation, by
definition. We can then run an optimisation program, e.g., “fmincon” on MATLAB
programing language, to maximise ${\cal F}^{\rm
C}({\bm{\sigma}};{\bm{S}}^{M}({\bm{S}}^{M})^{\intercal})$—alternatively to
minimise -${\cal F}^{\rm
C}({\bm{\sigma}};{\bm{S}}^{M}({\bm{S}}^{M})^{\intercal})$—subject to the
constraint that
${\bm{S}}^{M}{\bm{\Omega}}({\bm{S}}^{M})^{\intercal}={\bm{\Omega}}$. See the
appendix for an example code.

####  III.2.7 Multimode Gaussian states at thermal equilibrium: Local versus
global Gaussian measurements

Figure 3: The role of global measurements when the system is in a product
state i.e., the covariance matrix is a direct sum. (a) The system undergoes an
arbitrary [physical] joint operation then each mode is locally measured. We
have no restriction in the performed local measurements. In this case one can
show that the optimal scenario is a local one. That is, we do not need to
perform any joint operation; optimisation over arbitrary local measurements is
sufficient. This is so because of the additivity of the quantum Fisher
information. (b) Same as (a), but when the joint operation and the local
measurements are all Gaussian. Unlike (a), if we are restricted to perform
only Gaussian measurements, it is not clear whether joint Gaussian operations
are beneficial or not. We conjecture that at thermal equilibrium they are not.

If the system is at thermal equilibrium, one can always diagonalise its
covariance matrix by means of some temperature-independent symplectic
transformation—i.e., to write it in the normal-mode basis. The symplectic
transformation is a Gaussian process, and can be considered as part of the
Gaussian measurement which we will optimise anyway. Therefore, we can focus
only on initial thermal states with an uncorrelated covariance matrix
${\bm{\sigma}}=\oplus_{k}{\nu_{k}}{\bm{I}}_{2}$. As for the measurement, one
can first perform a joint Gaussian symplectic transformation on all
modes—s.t.,
${\bm{\sigma}}\to{\bm{S}}_{s}^{M}{\bm{\sigma}}({\bm{S}}_{s}^{M})^{\intercal}$—and
then perform a local Gaussian measurement
${\bm{\sigma}}_{s}^{M}={\oplus}_{k}{\bm{\sigma}}_{s_{k}}^{M}$. The idea is to
find the optimal ${\bm{S}}_{s}^{M}$ and ${\bm{\sigma}}_{s_{k}}^{M}~{}\forall
k$.

Let us start by a scenario that only allows for local Gaussian measurements,
i.e., ${\bm{S}}^{M}_{s}={\bm{I}}_{2m}$. In this case, one can easily check
that the classical Fisher information is additive, that is

$\displaystyle{\cal F}^{\rm
C}(\oplus_{k}{\bm{\sigma}}_{k};\oplus_{k}{\bm{\sigma}}_{s_{k}}^{M})=\sum_{k}{\cal
F}^{\rm C}({\bm{\sigma}}_{k};{\bm{\sigma}}_{s_{k}}^{M}).$ (46)

Therefore, when we optimise over local Gaussian measurements, we have the
optimal Fisher information

${\cal F}^{\rm
C}(\oplus_{k}{\bm{\sigma}}_{k};\oplus_{k}{\bm{\sigma}}_{k,\max}^{M})=\sum_{k}{\cal
F}^{\rm C}({\bm{\sigma}}_{k};{\bm{\sigma}}_{k,\max}^{M}),$

that is the sum of optimal Fisher information of individual modes. The
question that remains to be answered is: can the performance of Gaussian
measurements get better if we allow for joint measurements, i.e., if
${\bm{S}}_{s}^{M}\neq{\bm{I}}_{2m}$?

As depicted in Fig. 3, the global Gaussian measurements can be seen as the
operation of some joint Gaussian unitary on all modes—namely an $m$-mode
symplectic transformation—followed by local measurements. The question is then
whether the joint Gaussian operation can be beneficial. Mathematically, we
know that ${\cal F}^{\rm
C}(\oplus_{k}{\bm{\sigma}}_{k};{\bm{\sigma}}^{M}_{\max})\geqslant{\cal F}^{\rm
C}(\oplus_{k}{\bm{\sigma}}_{k};\oplus_{k}{\bm{\sigma}}_{k,\max}^{M})$, since
the local measurements are a subset of global ones. We now wonder whether this
inequality is in fact an equality, i.e., if

$\displaystyle{\cal F}^{\rm
C}(\oplus_{k}{\bm{\sigma}}_{k};{\bm{\sigma}}^{M}_{\max})\overset{?}{=}{\cal
F}^{\rm C}(\oplus_{k}{\bm{\sigma}}_{k};\oplus_{k}{\bm{\sigma}}_{k,\max}^{M}).$
(47)

Note that, if we were not restricted to Gaussian measurements and could
perform arbitrary measurements, then the equality would always hold. This is
so because in this case the optimal Fisher information would reduce to the
quantum Fisher information, which is additive and can be achieved by locally
optimal measurements, that is

$\displaystyle{\cal F}^{\rm C}(\oplus_{k}{\bm{\sigma}}_{k};M_{\max})={\cal
F}^{\rm C}(\oplus_{k}{\bm{\sigma}}_{k};\oplus_{k}M_{k,\max})={\cal F}^{\rm
Q}(\oplus_{k}{\bm{\sigma}}_{k})=\sum_{k}{\cal F}^{\rm Q}({\bm{\sigma}}_{k}),$
(48)

with $M_{\max}$ being the optimal joint measurement—that is not necessarily
Gaussian—performed on all modes, while $M_{k,\max}$ is the optimal local
measurement—again not necessarily Gaussian—performed on mode $k$.

Unfortunately, we cannot prove the equivalence Eq. (47) in the general case.
However, in the appendix we prove it analytically when all modes have
identical covariance matrices—i.e., when
${\bm{\sigma}_{k}}={\bm{\sigma}}_{j}\forall j,k$. Furthermore, our numerical
simulations support the equivalence for arbitrary product states. Based on
these observations we conjecture this is generally true, however, a rigorous
proof is missing currently.

###  III.3 Photo-Detection

Another type of measurement that will be shown to be relevant for thermometry
and which can be of interest for experimental implementation is Photo-
Detection. Photo-Detection is a 2-outcome (detected or not detected)
measurement and therefore it is not Gaussian. However, a link can be made
between the Photo-Detection POVM elements and the POVM elements of a Gaussian
measurement, characterized in Section II. In what follows we consider
realistic photo-detection, in which we allow for a probability $P_{D}$ of
random counts—see e.g., Rahimi-Keshari  _et al._ (2016); Barnett  _et al._
(1998) for the error sources contributing to the random counts. For a single-
mode system the POVM elements can then be written as

$\displaystyle{M}_{0|\rm{rpd}}$ $\displaystyle=(1-P_{D})|0\rangle\langle 0|,$
$\displaystyle{M}_{{\bar{0}}|{\rm rpd}}$ $\displaystyle=I-{M}_{0|\rm{rpd}},$
(49)

where $|0\rangle\langle 0|$ is merely the vacuum state, and $I$ is the
identity operator in the Hilbert space of the quantum system (not to be
mistaken with ${\bm{I}}_{2}$). The vacuum itself is a pure Gaussian state with
a CM equal to identity ${\bm{\sigma}}_{\rm vac}^{M}={\bm{I}}_{2}$ and null
displacement vector. The probability of finding any Gaussian system—with
density matrix $\rho$, covariance matrix ${\bm{\sigma}}$ and a vanishing
displacement vector—in vacuum is then given by

$\displaystyle P(0|{\bm{\sigma}};{\rm rpd})=\operatorname{Tr}[\rho{M}_{0|{\rm
rpd}}]=(1-P_{D})\operatorname{Tr}[\rho|0\rangle\langle
0|]=\frac{1-P_{D}}{\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}(\bm{I}_{\bm{2}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}(\bm{I}_{\bm{2}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}(\bm{I}_{\bm{2}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}(\bm{I}_{\bm{2}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}},$ (50)

where we use Eq. (II.1.1). Of course we also find the probablity of detection
to be $P(1|{\rm rpd};{\bm{\sigma}})=1-P(0|{\rm rpd};{\bm{\sigma}})$. The
Classical Fisher Information associated to this probability distribution can
be simply found from (22) and reads

$\displaystyle{\cal F}^{\rm C}({\bm{\sigma}};{\rm rpd})$
$\displaystyle=\frac{[\partial_{T}P(0|{\bm{\sigma}};{\rm
rpd})]^{2}}{P(0|{\bm{\sigma}};{\rm rpd})P(1|{\bm{\sigma}};{\rm rpd})},$ (51)

As an example, we can use photo-detection for thermometry of a single mode at
thermal equilibrium, with ${\bm{\sigma}}=\coth(\omega/2T){\bm{I}}_{2}$. Using
(50) one finds $P(0|{\bm{\sigma}};{\rm rpd})=(1-P_{D})P_{0}$ with
$P_{0}\coloneqq(1-\exp[-\omega/T])^{-1}$. Thus the Fisher information is given
by

$\displaystyle{\cal F}^{\rm C}({\bm{\sigma}};{\rm
rpd})=\frac{(1-P_{D})[\partial_{T}P_{0}]^{2}}{P_{0}\left[1-(1-P_{D})P_{0}\right]}.$
(52)

In Fig. 2 we plot the classical Fisher information of realistic photo-
detection for various $P_{D}$, and benchmark it against the optimal
measurement and the optimal Gaussian measurement. Ideal photo-detection
performs close to optimal specifically at low temperatures. However, as soon
as detection error is introduced and $P_{D}>0$, the optimality at low
temperatures is lost. The intuition behind photo-detection being optimal is as
follows: Firstly, we note that at thermal equilibrium the optimal measurement
is in the basis of Hamiltonian, i.e., projection in the number basis.
Secondly, at very low $T$ the system is basically frozen in its lowest energy
eigenstates, namely it is in the vacuum state with some probability of
populating the first excited state. As such, the projection in the number
basis can be well approximated by ideal photo-detection. Since the probability
of observing the system in an excited state is already very low, then any
random count will hugely impact our statistics and therefore realistic photo-
detection will eventually fail to estimate low enough temperatures. This
argument can be made rigorous mathematically as well. At low temperature, one
can see

$\displaystyle{\cal F}^{\rm C}({\bm{\sigma}};{\rm rpd})$
$\displaystyle\approx\begin{cases}\frac{\omega^{2}e^{-\omega/T}}{T^{4}},&\hskip
14.22636ptP_{D}=0\\\
\frac{(1-P_{D})\omega^{2}e^{-2\omega/T}}{P_{D}T^{4}},&\hskip
14.22636ptP_{D}\neq 0\end{cases}$ (53)

Note that the exponential term in the case $P_{D}\neq 0$ vanishes more
rapidly, which is why the realistic photo-detection performs far from optimal
for small enough $T$. Moreover, by using Eq. (45) one can confirm that the
leading order of quantum Fisher information at $T\to 0$ has the same value as
the ideal photo-detection, i.e., to the leading order ${\cal F}^{\rm
Q}(\bm{\sigma})\approx\omega^{2}\exp[-\omega/T]/T^{4}$.

####  III.3.1 Multimode photo-detection

For general m-mode systems, we can perform photo-detection individually on
each mode. The $2^{\textit{m}}$ POVM elements can then be written as tensor
product $M_{{\bm{a}}|{\rm rpd}}={M}_{{\bm{a}}_{1}|{\rm
rpd}}^{1}\otimes{M}_{{\bm{a}}_{2}|{\rm
rpd}}^{2}\otimes...\otimes{M}_{{\bm{a}}_{\textit{M}}|{\rm rpd}}^{\textit{m}}$
where ${\bm{a}}_{i}\in\\{0,1\\}\forall i$. In what follows, we use ${\bm{0}}$
to denote the vector of all zeros, ${{\bm{0}}}_{k}$ to denote the vector of
all zeros except the $k$th element being ${\bar{0}}$, ${{\bm{0}}}_{k,l}$ to
denote the vector of all zeros except elements $k$ and $l$ which are
${\bar{0}}$ and so on. Given an $m$-mode Gaussian system $\rho$ with the
covariance matrix ${\bm{\sigma}}$, the probability of finding, for example,
all modes in the vacuum state is (recall that all first moments of the system
vanish)

$\displaystyle P({\bm{0}}|{\rm
rpd};{\bm{\sigma}})=\operatorname{Tr}[\rho{M}_{{\bm{0}}|{\rm rpd}}]$
$\displaystyle=\frac{\Pi_{i}(1-P_{D}^{i})}{\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}(\bm{I}_{{2m}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}(\bm{I}_{{2m}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}(\bm{I}_{{2m}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}(\bm{I}_{{2m}}+\bm{\sigma})\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}},$ (54)

where $P_{D}^{i}$ is the probability of random count for mode $i$, and
$\bm{I}_{{2m}}$ is the $2m\times 2m$ identity matrix. Similarly, by using
$M_{{\bm{0}}_{k}|{\rm rpd}}=M_{0|{\rm rpd}}^{1}\dots\otimes(I-M_{0|{\rm
rpd}}^{k})\dots\otimes M_{0|{\rm rpd}}^{m}$, the probability of finding the
$k$th mode in an excited state while all the other modes at vacuum reads

$\displaystyle P({{\bm{0}}_{k}}|{\rm rpd};{\bm{\sigma}})=\frac{\Pi_{i\neq
k}(1-P_{D}^{i})}{\mathchoice{{\hbox{$\displaystyle\sqrt{{\rm
det}(\operatorname{tr}_{k}[\bm{I}_{{2m}}+\bm{\sigma}])\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt{{\rm
det}(\operatorname{tr}_{k}[\bm{I}_{{2m}}+\bm{\sigma}])\,}$}\lower
0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt{{\rm
det}(\operatorname{tr}_{k}[\bm{I}_{{2m}}+\bm{\sigma}])\,}$}\lower
0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{\rm
det}(\operatorname{tr}_{k}[\bm{I}_{{2m}}+\bm{\sigma}])\,}$}\lower
0.4pt\hbox{\vrule height=3.75pt,depth=-3.00002pt}}}}-P({{\bm{0}}}|{\rm
rpd};{\bm{\sigma}}),$ (55)

where, for any covariance matrix ${\bm{\sigma}}$, the operation
$\operatorname{tr}_{k}[{\bm{\sigma}}]$ is understood as eliminating the $2$
rows and columns associated to the mode $k$. Following similar steps one can
find the probability of other outcomes as well. However, one can restrict to a
simplified photo-detection with $m+2$ POVM elements, $M_{{\rm I}|{\rm
rpd}}\coloneqq M_{{\bm{0}}|{\rm rpd}}$, $M_{{\rm II}_{k}|{\rm rpd}}\coloneqq
M_{{\bm{0}}_{k}|{\rm rpd}},\forall k$ and $M_{{\rm III}|{\rm rpd}}\coloneqq
I-M_{{\rm I}|{\rm rpd}}-\sum_{k}M_{{\rm II}_{k}|{\rm rpd}}$. One can consider
yet simpler version, in which one does not distinguish what mode is excited
therefore we have only three POVM elements with $M_{{\rm II}|{\rm
rpd}}=\sum_{k}M_{{\rm II}_{k}|{\rm rpd}}$ while the other two elements remain
the same. These simplified versions should be fairly comparable to the non-
simplified version at low temperatures, becasue at these temperatures the
other POVM elements should not play a significant role anyway.

MOHA: Can you see if this is analytically doable also out of equilibrium

##  IV Conclusions

The problem of thermometry of Gaussian states with experimentally feasible
measurements was addressed. The classical Fisher information of Gaussian
measurements was charachterized. Finding the best Gaussian measurement can be
expressed as an optimisation problem in which one looks for a pure covariance
matrix that maximises the classical Fisher information (36). This may be done
numerically in general. For thermal states one may find analytical solutions.
In particular, for single-mode states at thermal equilibrium the optimal
Gaussian measurement was found analytically to be either Homodyne or
Heterodyne detection depending on the temperature. This is unlike the optimal
non-Gaussian measurement, which is always a projection in the eigenbasis of
the Hamiltonian of the system, regardless of the temperature. For a multimode
scenario at thermal equilibrium we conjecture that the optimal Gaussian
measurement is composed of a Gaussian unitary that brings the system into its
normal modes basis, followed by a local optimal measurement on each of the
normal modes. The proof of this is equivalent to proving Eq. (47), which
remains an open question. Our numerical results support this conjecture, and
we have proven this analytically for a scenario with all normal modes being
equal.

We further compare the performance of the Gaussian measurements as well as
photo-detection-like measurements with the quantum Fisher information, for
systems at thermal equilibrium. We see that at low temperatures photo-
detection-like measurements perform almost optimally, whereas at high
temperatures Gaussian measurements are close to optimal and saturate the
quantum Cramér-Rao bound. We expect that for systems that are out of thermal
equilibrium these results do not hold. For example, it is shown in Correa  _et
al._ (2017); Khan  _et al._ (2021); Mehboudi  _et al._ (2019b) that at
ultracold temperatures and for non-equilibrium systems, sometimes position
measurement—which is a Gaussian measurement—is optimal. A deeper investigation
in this direction will be the be subject of future works.

Lastly, it would be interesting to consider the problem of thermometry with
Gaussian measurements in the context of the Bayesian formalism, which has been
the subject of few recent works Rubio  _et al._ (2020); Mehboudi  _et al._
(2021); Jorgensen  _et al._ (2021); Alves and Landi (2021). In particular,
since the optimal Gaussian measurement is temperature dependent, this adds an
extra challenge into designing optimal thermometry protocols when initial
uncertainty about the true temperature is non negligible.

## Acknowledgements

We are thankful to S. Rahimi-Keshari, F. Bakhshinezhad, J. Kolodynski, and P.
Sekatski for constructive discussions. This work was supported by Swiss
National Science Foundation NCCR SwissMAP, the Government of Spain
(FIS2020-TRANQI, ConTrAct FIS2017-83709-R, and Severo Ochoa CEX2019-000910-S),
Fundacio Cellex, Fundacio Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR
1381 and QuantumCAT), the ERC AdG CERQUTE, the AXA Chair in Quantum
Information Science. This project has received funding from the “Presidencia
de la Agencia Estatal de Investigación” within the “Convocatoria de
tramitación anticipada, correspondiente al año 2020, de las ayudas para
contratos predoctorales (Ref. PRE2020-094374 PhD Fellowship) para la formación
de doctores contemplada en el Subprograma Estatal de Formación del Programa
Estatal de Promoción del Talento y su Empleabilidad en I+D+i, en el marco del
Plan Estatal de Investigación Científica y Técnica y de Innovación 2017-2020,
cofinanciado por el Fondo Social Europeo”. LL acknowledges support from
Alexander von Humboldt Stiftung/Foundation.

##  Appendix A More about Gaussian measurements

###  A.1 General form of Gaussian measurements

In general, Gaussian measurements can be performed over an output state of any
Gaussian map on the system $\rho_{s}^{M}$, which is in general not a pure
state, and will then be represented by more general POVM elements that read

${M}_{\bm{a}|{s}}=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,D\big{(}\bm{y}\bm{T}_{s}\big{)}e^{{-\bm{y}\bm{N}_{s}\bm{y}^{\intercal}/2-i\bm{a}\bm{\Omega}\bm{y}^{\intercal}}},$
(56)

where the complete positivity condition
$\bm{N}_{s}-i\bm{T}_{s}\bm{\Omega}\bm{T}_{s}^{\intercal}/2\geqslant 0$ must
hold. The outcome probabilities of this Gaussian measurement are then given by

$\displaystyle p(\bm{a}|{s},\rho)$
$\displaystyle=\operatorname{tr}[\rho{M}_{\bm{a}|{s}}]=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,\operatorname{tr}[D\big{(}\bm{y}\bm{T}_{s}\big{)}\rho]e^{{-\bm{y}\bm{N}_{s}\bm{y}^{\intercal}/2-i\bm{a}\bm{\Omega}\bm{y}^{\intercal}}}$
$\displaystyle=\int\frac{d^{\,2m}\\!\bm{y}\,}{(2\pi)^{2m}}\,e^{{-\bm{y}\bm{T}_{s}\bm{\Omega}\bm{\sigma}\bm{\Omega}\bm{T}_{s}^{\intercal}\bm{y}^{\intercal}/2}+i\bm{d}\bm{\Omega}\bm{T}_{s}^{\intercal}\bm{y}^{\intercal}}e^{{-\bm{y}\bm{N}_{s}\bm{y}^{\intercal}/2-i\bm{a}\bm{\Omega}\bm{y}^{\intercal}}}$
(57) $\displaystyle=\frac{1}{(2\pi)^{2m}}\int
d^{2m}\bm{y}e^{-\frac{1}{2}\bm{y}{\tilde{\bm{\sigma}}}\bm{y}^{\intercal}+i{\tilde{\bm{d}}}\bm{y}^{\intercal}}.$
(58)

where we have defined
$\tilde{\bm{\sigma}}\coloneqq\bm{T}_{s}\bm{\Omega}\bm{\sigma}\bm{\Omega}\bm{T}_{s}^{\intercal}+\bm{N}_{s}$
and
$\tilde{\bm{d}}\coloneqq\bm{d}\bm{\Omega}\bm{T}_{s}^{\intercal}-\bm{a}\bm{\Omega}$.
Integration gives

$\displaystyle p(\bm{a}|{s},\rho)$
$\displaystyle=\frac{1}{(2\pi)^{m}\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\tilde{\bm{\sigma}}}\,}$}\lower
0.4pt\hbox{\vrule
height=7.22223pt,depth=-5.77782pt}}}{{\hbox{$\textstyle\sqrt{\det{\tilde{\bm{\sigma}}}\,}$}\lower
0.4pt\hbox{\vrule
height=7.22223pt,depth=-5.77782pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\tilde{\bm{\sigma}}}\,}$}\lower
0.4pt\hbox{\vrule
height=7.22223pt,depth=-5.77782pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\tilde{\bm{\sigma}}}\,}$}\lower
0.4pt\hbox{\vrule
height=7.22223pt,depth=-5.77782pt}}}}e^{-\frac{1}{2}\tilde{\bm{d}}{\tilde{\bm{\sigma}}}^{-1}\tilde{\bm{d}}^{\intercal}}.$
(59)

The previous scenario is revived by choosing
$\bm{N}_{s}=\bm{\Omega}{\bm{\sigma}_{s}^{M}}\bm{\Omega}$ and
$\bm{T}_{s}=\bm{I}_{2m}$.

###  A.2 Local Gaussian measurements, reduced states and covariance matrices

Suppose we have a multi-mode state and we perform a Gaussian measurement on a
subsystem. What is the post measurement state of the remaining modes? In
principle, after performing some POVM on the system and observing an outcome
associated with the element ${M}_{\bm{a}|{s}}$, the post measurement state
is—up to a unitary—described by the density matrix

$\displaystyle\rho({\bm{a}}|s)=\frac{\mathchoice{{\hbox{$\displaystyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}\rho\mathchoice{{\hbox{$\displaystyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}}{\operatorname{tr}[\mathchoice{{\hbox{$\displaystyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}\rho\mathchoice{{\hbox{$\displaystyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}]}=\frac{\mathchoice{{\hbox{$\displaystyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}\rho\mathchoice{{\hbox{$\displaystyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{{M}_{\bm{a}|{s}}\,}$}\lower
0.4pt\hbox{\vrule height=3.41666pt,depth=-2.73334pt}}}}{P(\bm{a}|{s},\rho)},$
(60)

Despite the unitary freedom in the post measurement state, the reduced state
of the modes over which we did not perform the measurement is independent of
this unitary.

Suppose we have a two-party density operator, $\rho_{AB}$, of $n+m$ modes,
where party $A$ has access to the first $m$ modes and party $B$ has the
remaining $n$ modes. A POVM measurement $\\{{M}_{\bm{b}|{s}}\\}$ is performed
on part $B$. That is, the overall POVM elements are given by
$\\{I_{A}\otimes{M}_{\bm{b}|{s}}\\}$. Given that we observe the result
${\bm{b}}$, the post-measurement [reduced] state of party $A$ reads

$\displaystyle\rho_{A}({\bm{b}},s)\coloneqq\operatorname{tr}_{B}[\rho_{AB}({\bm{b}}|s)]=\frac{\operatorname{tr}_{B}[I_{A}\otimes{M}_{\bm{b}|{s}}~{}\rho_{AB}]}{P(\bm{b}|{s},\rho_{AB})}.$
(61)

We are interested in knowing the covariance matrix of this reduced state. To
this aim, we need to find the characteristic function, from which one can
identify the covariance matrix.

### Characteristic function of the reduced density matrix

By definition, we have

$\displaystyle\chi(\bm{y}|\rho_{A}({\bm{b}},s))=\operatorname{tr}_{A}[D_{A}(\bm{y})\rho_{A}({\bm{b}},s)]=\frac{\operatorname{tr}[D_{A}(\bm{y})\otimes{M}_{\bm{b}|{s}}~{}\rho_{AB}]}{P(\bm{b}|{s},\rho)}.$
(62)

Let us focus on the nominator since the denominator is only a normalisation
factor independent of ${\bm{y}}$. By noting that $D_{B}({\bm{0}})=I_{B}$, this
can be written as

$\displaystyle\operatorname{tr}[D_{A}(\bm{y})\otimes{M}_{\bm{b}|{s}}~{}\rho_{AB}]$
$\displaystyle=\operatorname{tr}[D_{A}(\bm{y})\otimes
D_{B}({\bm{0}})~{}I_{A}\otimes{M}_{\bm{b}|{s}}~{}\rho_{AB}]\eqqcolon\operatorname{tr}[D_{AB}(\bm{y},{\bm{0}})~{}\varrho_{AB}~{}\rho_{AB}],$
(63)

where we denote $D_{AB}(\bm{y},{\bm{z}})=D_{A}(\bm{y})\otimes
D_{B}({\bm{z}})$, and $\varrho_{AB}\coloneqq I_{A}\otimes{M}_{\bm{b}|{s}}$. We
can expand both $\varrho_{AB}$ and $\rho_{AB}$ in terms of the displacement
operator by using Eq. (2), such that

$\displaystyle\varrho_{AB}\rho_{AB}$
$\displaystyle=\frac{1}{(2\pi)^{n+m}}\iint
d^{2(n+m)}\bm{x}d^{2(n+m)}\bm{z}~{}\operatorname{tr}[\varrho_{AB}D(\bm{x})]\operatorname{tr}[\rho_{AB}D(\bm{z})]D^{\dagger}(\bm{x})D^{\dagger}(\bm{z})$
$\displaystyle=\frac{1}{(2\pi)^{n+m}}\iint
d^{2(n+m)}\bm{x}d^{2(n+m)}\bm{z}~{}\chi(\bm{x}|\varrho_{AB})\chi(\bm{z}|\rho_{AB})D^{\dagger}(\bm{x})D^{\dagger}(\bm{z}),$
(64)

where $n$ is the number of modes over which we performed the measurement. Thus
we can find the characteristic function of $\varrho_{AB}\rho_{AB}$ using the
above expression

$\displaystyle\operatorname{tr}[D_{AB}(\bm{y})\varrho_{AB}\rho_{AB}]$
$\displaystyle=\frac{1}{(2\pi)^{n+m}}\iint
d^{2(n+m)}\bm{x}d^{2(n+m)}\bm{z}~{}\chi(\bm{x}|\varrho_{AB})\chi(\bm{z}|\rho_{AB})\operatorname{tr}[D(\bm{y})D^{\dagger}(\bm{x})D^{\dagger}(\bm{z})]$
$\displaystyle=\frac{1}{(2\pi)^{n+m}}\iint
d^{2(n+m)}\bm{x}d^{2(n+m)}\bm{z}~{}\chi(\bm{x}|\varrho_{AB})\chi(\bm{z}|\rho_{AB})$
$\displaystyle~{}~{}\times
e^{-i\bm{x}\bm{\Omega}\bm{z}^{\intercal}/2}e^{i\bm{y}\bm{\Omega}(\bm{z}+\bm{x})^{\intercal}/2}~{}\operatorname{tr}[D(\bm{y}-\bm{x}-\bm{z})]$
$\displaystyle=\iint
d^{2(n+m)}\bm{x}d^{2(n+m)}\bm{z}~{}\chi(\bm{x}|\varrho_{AB})\chi(\bm{z}|\rho_{AB})e^{-i\bm{x}\bm{\Omega}\bm{z}^{\intercal}/2}~{}\delta^{2(n+m)}(\bm{y}-\bm{x}-\bm{z})$
$\displaystyle=\int
d^{2(n+m)}\bm{x}~{}\chi(\bm{x}|\varrho_{AB})\chi(\bm{y}-\bm{x}|\rho_{AB})e^{-i\bm{x}\bm{\Omega}\bm{y}^{\intercal}/2},$
(65)

where we used that
$\operatorname{tr}[D(\bm{a})]=(2\pi)^{n+m}\delta^{2(n+m)}(\bm{a})$. Let us
denote with $\\{\bm{\sigma}_{AB},\bm{d}_{AB}\\}$ and
$\\{{\bm{\sigma}_{s}^{M}},\bm{d}_{s}^{M}\\}$ the covariance matrix and the
displacement vectors of $\rho_{AB}$ and $\varrho_{AB}$, respectively.
Therefore

$\displaystyle\operatorname{tr}[D_{AB}(\bm{y})\varrho_{AB}\rho_{AB}]$ (66)
$\displaystyle=\int
d^{2(n+m)}\bm{x}e^{-\frac{1}{2}\bm{x}\bm{\Omega}{\bm{\sigma}_{s}^{M}}\bm{\Omega}\bm{x}^{\intercal}+i\bm{d}_{s}^{M}\bm{\Omega}\bm{x}^{\intercal}}e^{-\frac{1}{2}(\bm{y}-\bm{x})\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}(\bm{y}-\bm{x})^{\intercal}+i\bm{d}_{AB}\bm{\Omega}(\bm{y}-\bm{x})^{\intercal}}e^{-i\bm{x}\bm{\Omega}\bm{y}^{\intercal}/2}$
$\displaystyle=e^{-\frac{1}{2}\bm{y}\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}\bm{y}^{\intercal}+i\bm{d}_{AB}\bm{\Omega}\bm{y}^{\intercal}}\int
d^{2(n+m)}\bm{x}e^{-\frac{1}{2}\bm{x}\bm{\Omega}(\bm{\sigma}_{AB}+{\bm{\sigma}}_{s}^{M})\bm{\Omega}\bm{x}^{\intercal}+(i\bm{d}_{s}^{M}-i\bm{d}_{AB}+i\bm{y}/2+\bm{y}\bm{\Omega}\bm{\sigma}_{AB})\bm{\Omega}\bm{x}^{\intercal}}$
$\displaystyle=e^{-\frac{1}{2}\bm{y}\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}\bm{y}^{\intercal}+i\bm{d}_{AB}\bm{\Omega}\bm{y}^{\intercal}}\int
d^{2(n+m)}\bm{x}e^{-\frac{1}{2}\bm{x}\bm{\tilde{\sigma}}\bm{x}^{\intercal}}e^{i\bm{\tilde{d}}\bm{x}^{\intercal}},$
(67)

with
$\bm{\tilde{\sigma}}\coloneqq\bm{\Omega}(\bm{\sigma}_{AB}+{{\bm{\sigma}_{s}^{M}}})\bm{\Omega}$,
and
$\bm{\tilde{d}}\coloneqq(\bm{d}_{s}^{M}-\bm{d}_{AB}+\bm{x}/2-i\bm{x}\bm{\Omega}\bm{\sigma}_{AB})\bm{\Omega}$.
The integral can also be evaluated using (II.1.1), which leads to

$\displaystyle\int
d^{2(n+m)}\bm{x}e^{-\frac{1}{2}\bm{x}\bm{\tilde{\sigma}}\bm{x}^{\intercal}}e^{i\bm{\tilde{d}}\bm{x}^{\intercal}}$
$\displaystyle=\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{(2\pi)^{2(n+m)}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=15.71112pt,depth=-12.56895pt}}}{{\hbox{$\textstyle\sqrt{\frac{(2\pi)^{2(n+m)}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=12.11111pt,depth=-9.68893pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{(2\pi)^{2(n+m)}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=10.11111pt,depth=-8.08893pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{(2\pi)^{2(n+m)}}{{\rm
det}{\tilde{\bm{\sigma}}}}\,}$}\lower 0.4pt\hbox{\vrule
height=10.11111pt,depth=-8.08893pt}}}e^{-\frac{1}{2}\bm{\tilde{d}}{\bm{\tilde{\sigma}}}^{-1}\bm{\tilde{d}}^{\intercal}}.$
(68)

For simplicity let us assume that $\bm{d}_{AB}=\bm{d}_{s}^{M}=0$. Then we have

$\displaystyle\operatorname{tr}[D(\bm{y})\varrho_{AB}\rho_{AB}]$
$\displaystyle\propto
e^{-\frac{1}{2}\bm{y}\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}\bm{y}^{\intercal}}e^{-\frac{1}{2}\bm{\tilde{d}}{\bm{\tilde{\sigma}}}^{-1}\bm{\tilde{d}}^{\intercal}}$
$\displaystyle=e^{-\frac{1}{2}\bm{y}\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}\bm{y}^{\intercal}}e^{\frac{1}{2}\bm{y}(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})(\bm{\sigma}_{AB}+{{\bm{\sigma}_{s}^{M}}})^{-1}(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})^{\intercal}\bm{y}^{\intercal}}$
$\displaystyle=\exp\left\\{{-\frac{1}{2}\bm{y}\left[\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}-(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})(\bm{\sigma}_{AB}+{\bm{\sigma}}_{s})^{-1}(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})^{\intercal}\right]\bm{y}^{\intercal}}\right\\}.$
(69)

We can now simply find the post-measurement covariance matrix of party $A$
after the Gaussian measurement ${M}_{\bm{a}|{s}}$ of party $B$ by setting
$D_{AB}(\bm{y})\to D_{AB}(\bm{y}_{A},{\bm{0}}_{B})$ in the exponent of the
above relation, and by noticing that the covariance matrix of the measurement
operator can be written as
${\bm{\sigma}_{s}^{M}}=\lim_{R\to\infty}R~{}{\bm{I}}_{A}\oplus\bm{\sigma}_{s,B}$.
Using the fact that

$\displaystyle(\bm{\sigma}_{AB}+{\bm{\sigma}_{s}^{M}})^{-1}=\lim_{R\to\infty}(\bm{\sigma}_{AB}+R~{}{\bm{I}}_{A}\oplus\bm{\sigma}_{s,B})^{-1}={\bm{0}}_{A}\oplus(\bm{\sigma}_{B}+\bm{\sigma}_{s,B})^{-1},$
(70)

we have

$\displaystyle\operatorname{tr}[D_{AB}(\bm{y}_{A},{\bm{0}}_{B})\varrho_{AB}\rho_{AB}]$
$\displaystyle\propto\exp\left\\{{-\frac{1}{2}(\bm{y}_{A},{\bm{0}}_{B})\left[\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}-(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})(\bm{\sigma}_{AB}+{\bm{\sigma}_{s}^{M}})^{-1}(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})^{\intercal}\right](\bm{y}_{A},{\bm{0}}_{B})^{\intercal}}\right\\}$
$\displaystyle=\exp\left\\{{-\frac{1}{2}(\bm{y}_{A},{\bm{0}}_{B})\left[\bm{\Omega}\bm{\sigma}_{AB}\bm{\Omega}-(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})({\bm{0}}_{A}\oplus(\bm{\sigma}_{B}+\bm{\sigma}_{s,B})^{-1})(\bm{I}/2-i\bm{\Omega}\bm{\sigma}_{AB})^{\intercal}\right](\bm{y}_{A},{\bm{0}}_{B})^{\intercal}}\right\\}$
$\displaystyle=\exp\left\\{-\frac{1}{2}\bm{y}_{A}\bm{\Omega}_{A}\left[\bm{\sigma}_{A}-\bm{\sigma}_{AB}^{\rm
Corr}(\bm{\sigma}_{B}+\bm{\sigma}_{s,B})^{-1}(\bm{\sigma}_{AB}^{\rm
Corr})^{\intercal}\right]\bm{\Omega}_{A}\bm{y}_{A}^{\intercal}\right\\}$ (71)

in which $\bm{\sigma}_{AB}^{\rm Corr}$ is the correlation block of the
covariance matrix $\bm{\sigma}_{AB}$, such that

$\displaystyle\bm{\sigma}_{AB}=\left[\begin{array}[]{cc}\bm{\sigma}_{A}&\bm{\sigma}_{AB}^{\rm
corr}\\\ (\bm{\sigma}_{AB}^{\rm corr})^{\intercal}&\bm{\sigma}_{B}\\\
\end{array}\right].$ (74)

Thus, the covariance matrix of $A$, after $B$’s measurement is independent of
the outcome and reads

$\displaystyle\bm{\sigma}_{A|s}=\bm{\sigma}_{A}-\bm{\sigma}_{AB}^{\rm
Corr}(\bm{\sigma}_{B}+\bm{\sigma}_{s,B})^{-1}(\bm{\sigma}_{AB}^{\rm
Corr})^{\intercal}.$ (75)

One can check that the same result holds for
$\bm{d}_{AB}\neq\bm{d}_{s}^{M}\neq 0$.

##  Appendix B The classical Fisher information of Gaussian probability
distributions

Given a Gaussian probability distribution $p({\bm{a}}|\theta)$ which can be
written as

$\displaystyle p(\bm{a}|\theta)$
$\displaystyle=\frac{e^{-\frac{1}{2}\bm{a}{\bm{\sigma}}^{-1}\bm{a}^{\intercal}}}{(2\pi)^{M}\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule height=3.47221pt,depth=-2.77779pt}}}},$ (76)

we want to find the classical Fisher information. Let us define
${\tilde{\bm{\sigma}}}\coloneqq{\bm{\sigma}}^{-1}$. Firstly, note that

$\displaystyle\partial_{\theta}\log
p(\bm{a}|\theta)=-\frac{\partial_{\theta}\det{\bm{\sigma}}}{2\det{\bm{\sigma}}}-\frac{1}{2}\bm{a}_{k}\bm{a}_{l}\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl},$
(77)

where summation over repeated indices is understood. Recall the Jacobi’s
formula for the derivative of determinant of [symmetric] matrices

$\displaystyle\partial_{\theta}\det{\bm{A}}$
$\displaystyle=\operatorname{Tr}[{\bm{A}}^{-1}\partial_{\theta}{\bm{A}}]\det{\bm{A}},$
(78) $\displaystyle\partial_{{\bm{A}}_{kl}}\det{\bm{A}}$
$\displaystyle={\bm{A}}^{-1}_{kl}\det{\bm{A}}.$ (79)

By substitution in the expression of classical Fisher information one gets

$\displaystyle F$ $\displaystyle=\int
d^{2M}{\bm{a}}p(\bm{a}|\theta)[\partial_{\theta}\log
p(\bm{a}|\theta)]^{2}=\int
d^{2M}{\bm{a}}p(\bm{a}|\theta)\left[-\frac{1}{2}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]-\frac{1}{2}{\bm{a}}_{k}{\bm{a}}_{l}\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\right]^{2}$
$\displaystyle=\frac{1}{4}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}+\frac{1}{2}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\int
d^{2M}{\bm{a}}~{}{\bm{a}}_{k}{\bm{a}}_{l}p(\bm{a}|\theta)+\frac{1}{4}\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\partial_{\theta}{\tilde{\bm{\sigma}}}_{mn}\int
d^{2M}{\bm{a}}~{}{\bm{a}}_{k}{\bm{a}}_{l}{\bm{a}}_{m}{\bm{a}}_{n}p(\bm{a}|\theta)$
$\displaystyle=\frac{1}{4}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}+\frac{\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}}{2(2\pi)^{M}\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}(-2\partial_{\tilde{\bm{\sigma}}_{kl}})\int
d^{2M}{\bm{a}}~{}e^{-\frac{1}{2}\bm{a}{\tilde{\bm{\sigma}}}\bm{a}^{\intercal}}$
$\displaystyle+\frac{\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\partial_{\theta}{\tilde{\bm{\sigma}}}_{mn}}{4(2\pi)^{M}\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}(4\partial_{\tilde{\bm{\sigma}}_{kl}}\partial_{\tilde{\bm{\sigma}}_{mn}})\int
d^{2M}{\bm{a}}~{}e^{-\frac{1}{2}\bm{a}{\tilde{\bm{\sigma}}}\bm{a}^{\intercal}}$
$\displaystyle=\frac{1}{4}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}-\frac{\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}}{\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}(\partial_{\tilde{\bm{\sigma}}_{kl}})\det{\tilde{\bm{\sigma}}}^{-\frac{1}{2}}+\frac{\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\partial_{\theta}{\tilde{\bm{\sigma}}}_{mn}}{\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}(\partial_{\tilde{\bm{\sigma}}_{kl}}\partial_{\tilde{\bm{\sigma}}_{mn}})\det{\tilde{\bm{\sigma}}}^{-\frac{1}{2}}$
$\displaystyle=\frac{1}{4}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}+\frac{1}{2}{\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}}[{\tilde{\bm{\sigma}}}^{-1}]_{kl}-\frac{\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\partial_{\theta}{\tilde{\bm{\sigma}}}_{mn}}{2\mathchoice{{\hbox{$\displaystyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\textstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=6.94444pt,depth=-5.55559pt}}}{{\hbox{$\scriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=4.8611pt,depth=-3.8889pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\det{\bm{\sigma}}\,}$}\lower
0.4pt\hbox{\vrule
height=3.47221pt,depth=-2.77779pt}}}}(\partial_{\tilde{\bm{\sigma}}_{mn}})([{\tilde{\bm{\sigma}}}^{-1}]_{kl}\det{\tilde{\bm{\sigma}}}^{-\frac{1}{2}})$
$\displaystyle=\frac{1}{4}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}+\frac{1}{2}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]\operatorname{Tr}[{\bm{\sigma}}\partial_{\theta}\tilde{\bm{\sigma}}]+\partial_{\theta}{\tilde{\bm{\sigma}}}_{kl}\partial_{\theta}{\tilde{\bm{\sigma}}}_{mn}\left(\frac{1}{4}[{\tilde{\bm{\sigma}}}^{-1}]_{kl}[{\tilde{\bm{\sigma}}}^{-1}]_{mn}+\frac{1}{2}[{\tilde{\bm{\sigma}}}^{-1}]_{mk}[{\tilde{\bm{\sigma}}}^{-1}]_{nl}\right)$
$\displaystyle=\frac{1}{4}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}-\frac{1}{2}\operatorname{Tr}[{\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}]^{2}+\frac{1}{4}\operatorname{Tr}[{\bm{\sigma}}\partial_{\theta}{\tilde{\bm{\sigma}}}]^{2}+\frac{1}{2}\operatorname{Tr}[(\partial_{\theta}{\tilde{\bm{\sigma}}}){\bm{\sigma}}(\partial_{\theta}{\tilde{\bm{\sigma}}}){\bm{\sigma}}]$
$\displaystyle=\frac{1}{2}\operatorname{Tr}[{\bm{\sigma}}^{-1}(\partial_{\theta}{\bm{\sigma}}){\bm{\sigma}}^{-1}(\partial_{\theta}{\bm{\sigma}})]=\frac{1}{2}\operatorname{Tr}[({\bm{\sigma}}^{-1}\partial_{\theta}{\bm{\sigma}})^{2}]$
(80)

where we benefited multiple uses of
${\tilde{\bm{\sigma}}}\partial_{\theta}{\bm{\sigma}}=-{\bm{\sigma}}\partial_{\theta}\tilde{\bm{\sigma}}$,
and also

$\displaystyle\partial_{\tilde{\bm{\sigma}}_{mn}}[{\tilde{\bm{\sigma}}}^{-1}]_{kl}$
$\displaystyle=-[{\tilde{\bm{\sigma}}}^{-1}(\partial_{\tilde{\bm{\sigma}}_{mn}}{\tilde{\bm{\sigma}}}){\tilde{\bm{\sigma}}}^{-1}]_{kl}=-[{\tilde{\bm{\sigma}}}^{-1}]_{kk^{\prime}}[\partial_{\tilde{\bm{\sigma}}_{mn}}{\tilde{\bm{\sigma}}}]_{k^{\prime}l^{\prime}}[{\tilde{\bm{\sigma}}}^{-1}]_{l^{\prime}l}=-[{\tilde{\bm{\sigma}}}^{-1}]_{kk^{\prime}}\delta_{mk^{\prime}}\delta_{nl^{\prime}}[{\tilde{\bm{\sigma}}}^{-1}]_{l^{\prime}l}$
$\displaystyle=-[{\tilde{\bm{\sigma}}}^{-1}]_{km}[{\tilde{\bm{\sigma}}}^{-1}]_{nl}.$
(81)

##  Appendix C Additivity of the Gaussian-measured Fisher information in a
special case

######  Proposition 1.

Let ${\bm{\sigma}}(T)$ be a single-mode covariance matrix, and assume that
${\bm{\sigma}}=\nu{\bm{I}}_{2}$ and
$\partial_{T}{\bm{\sigma}}=\mu{\bm{I}}_{2}$ are both proportional to the
identity. Then for all $n$ it holds that

$\max_{{\bm{\sigma}}^{M}_{s}}{\cal F}^{\rm C}\left({\bm{\sigma}}^{\oplus
n};{\bm{\sigma}}^{M}_{s}\right)=n\max_{{\bm{\gamma}}^{M}_{s}}{\cal F}^{\rm
C}\left({\bm{\sigma}};{\bm{\gamma}}^{M}_{s}\right)=n\mu^{2}\max\left\\{\frac{1}{\nu^{2}},\,\frac{2}{(1+\nu)^{2}}\right\\},$
(82)

where the optimization on the left-hand side is over all (multimode)
covariance matrices of the measurement ${\bm{\sigma}}^{M}_{s}$.

###### Proof.

Let ${\bm{\gamma}}^{M}_{s}$ be an arbitrary covariance matrix of an $n$-mode
Gaussian pure state. Then the spectrum of ${\bm{\gamma}}^{M}_{s}$ is well
known to be of the form
$\kappa_{1},\ldots,\kappa_{n},\kappa_{n}^{-1},\ldots,\kappa_{1}^{-1}$, where
$\kappa_{1},\ldots,\kappa_{n}\geq 1$. Thus,

$\displaystyle{\cal F}^{\rm
C}\left({\bm{\sigma}};{\bm{\sigma}}_{s}^{M}\right)$
$\displaystyle=\frac{1}{2}\operatorname{Tr}\left[\left(\left({\bm{\sigma}}^{\oplus
n}+{\bm{\gamma}}_{s}^{M}\right)^{-1}(\partial_{T}{\bm{\sigma}})^{\oplus
n}\right)^{2}\right]$
$\displaystyle=\frac{\mu^{2}}{2}\operatorname{Tr}\left[\left(\nu{\bm{I}}_{2n}+{\bm{\gamma}}_{s}^{M}\right)^{-2}\right]$
$\displaystyle=\frac{\mu^{2}}{2}\sum_{j=1}^{n}\left(\left(\nu+\kappa_{j}\right)^{-2}+\left(\nu+\kappa^{-1}\right)^{-2}\right).$

Hence, maximizing over ${\bm{\gamma}}^{M}_{s}$, that is to say, maximizing
over all choices of $\kappa_{j}$ subjected to the above constraints, we obtain
that

$\displaystyle\max_{{\bm{\sigma}}^{M}_{s}}{\cal F}^{\rm
C}\left({\bm{\sigma}}^{\oplus n};{\bm{\sigma}}^{M}_{s}\right)$
$\displaystyle=\max_{\kappa_{1},\ldots,\kappa_{n}\geq
1}\frac{\mu^{2}}{2}\sum_{j=1}^{n}\left(\left(\nu+\kappa_{j}\right)^{-2}+\left(\nu+\kappa^{-1}\right)^{-2}\right)$
$\displaystyle=\frac{\mu^{2}}{2}\sum_{j=1}^{n}\max_{\kappa_{j}\geq
1}\left(\left(\nu+\kappa_{j}\right)^{-2}+\left(\nu+\kappa^{-1}\right)^{-2}\right)$
$\displaystyle=\frac{n\mu^{2}}{2}\max\left\\{\frac{1}{\nu^{2}},\,\frac{2}{(1+\nu)^{2}}\right\\}$
$\displaystyle=n\max_{{\bm{\gamma}}^{M}_{s}}{\cal F}^{\rm
C}\left({\bm{\sigma}};{\bm{\gamma}}^{M}_{s}\right).$

This establishes the additivity of the Gaussian-measured Fisher information in
this special case. ∎

##  Appendix D Numerical optimisation, and the algorithm

We aim to identify the measurement CM ${\bm{\sigma}}_{s}^{M}$ that maximises
(36). We are dealing with a nonlinear optimisation problem subject to the
positive semi-definiteness criterion
${\bm{\sigma}}_{s}^{M}+i{\bm{\Omega}}\geqslant 0$. We already proved that we
should have ${\bm{\sigma}}_{s}^{M}={\bm{S}}^{M}({\bm{S}}^{M})^{T}$, with some
symplectic transformation satisfying
${\bm{S}}^{M}{\bm{\Omega}}({\bm{S}}^{M})^{T}={\bm{\Omega}}$. One can find the
optimal ${\bm{S}}^{M}$ in MATLAB by using the following code. Notice that this
is true for arbitrary parameter estimation problem and is not limited to
thermometry. One just needs to replace ${\rm dsigma}$ with the derivative of
the covariance matrix w.r.t. the parameter to be estimated.

⬇

1function [CFI,sigma_s] = optim_Gaussian(sigma_modes,dsigma_modes,num_modes)

2%——————-

3% sigma_modes is the Cov. Mat. of the system, which has

4% ‘num_modes’ modes

5% dsigma_modes is the derivative of sigma_modes

6% w.r.t the parameter

7%——————-

8%The Symplectic form

9Omega = [0 1;-1 0];Omega = kron(eye(num_modes),Omega);

10%Objective: (minus) the Fisher info.

11%Here, the measurement CM can be written as sigma^M=SS^T

12fisher = @(S) -1/2*trace( (dsigma_modes…

13 /(sigma_modes + S*transpose(S)))^2 );

14%The constraint on S, i.e., it should be symplectic

15nonlequ = @(S) S*Omega*transpose(S)-Omega;

16%the [] means we don’t have inequalities as constraints.

17%nonlequ(S) is an equality constraint

18nonlcon = @(S) deal([],nonlequ(S));

19%Initial guess for the optimisation variable;

20S0 = eye(2*num_modes);

21%fmincon options; change to suite your problem, if needed

22options = optimoptions(’fmincon’, ’StepTolerance’, 1e-30,…

23 ’FiniteDifferenceStepSize’, 1e-5,’MaxFunctionEvaluations’,10^7,…

24 MaxIterations,10^3,’Algorithm’,’interior-point’);

25%The actual minimisation

26[S,fval] = fmincon( fisher,S0,[],[],[],[],[],[],nonlcon,options);

27%The optimal measurement’s Cov. Mat.

28sigma_s = S*transpose(S);

29%The Fisher information corresponding to the optimal measurement

30CFI=-fval;

31end

## References

  * Mehboudi  _et al._ (2019a) M. Mehboudi, A. Sanpera, and L. A. Correa, Journal of Physics A: Mathematical and Theoretical 52, 303001 (2019a).
  * De Pasquale and Stace (2018) A. De Pasquale and T. M. Stace, in  _Thermodynamics in the Quantum Regime_ (Springer, 2018) pp. 503–527.
  * Giovannetti  _et al._ (2011) V. Giovannetti, S. Lloyd, and L. Maccone, Nature photonics 5, 222 (2011).
  * Giovannetti  _et al._ (2006) V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. Lett. 96, 010401 (2006).
  * Tóth and Apellaniz (2014) G. Tóth and I. Apellaniz, Journal of Physics A: Mathematical and Theoretical 47, 424006 (2014).
  * Correa  _et al._ (2015) L. A. Correa, M. Mehboudi, G. Adesso, and A. Sanpera, Phys. Rev. Lett. 114, 220405 (2015).
  * Jorgensen  _et al._ (2020) M. R. Jorgensen, P. P. Potts, M. G. A. Paris, and J. B. Brask, Phys. Rev. Research 2, 033394 (2020).
  * Hovhannisyan  _et al._ (2021) K. V. Hovhannisyan, M. R. J0rgensen, G. T. Landi, A. M. Alhambra, J. B. Brask, and M. Perarnau-Llobet, PRX Quantum 2, 020322 (2021).
  * 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, Rev. Mod. Phys. 84, 621 (2012).
  * Lewenstein  _et al._ (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Advances in Physics 56, 243 (2007), https://doi.org/10.1080/00018730701223200 .
  * Lewenstein  _et al._ (2012) M. Lewenstein, A. Sanpera, and V. Ahufinger,  _Ultracold Atoms in Optical Lattices: Simulating quantum many-body systems_ (Oxford University Press, 2012).
  * Serafini (2017a) A. Serafini,  _Quantum continuous variables: a primer of theoretical methods_ (CRC press, 2017).
  * Diamanti and Leverrier (2015) E. Diamanti and A. Leverrier, Entropy 17, 6072 (2015).
  * Navascués and Acín (2005) M. Navascués and A. Acín, Phys. Rev. Lett. 94, 020505 (2005).
  * Monras (2013) A. Monras, arXiv preprint arXiv:1303.3682 (2013).
  * Jiang (2014) Z. Jiang, Phys. Rev. A 89, 032128 (2014).
  * Nichols  _et al._ (2018) R. Nichols, P. Liuzzo-Scorpo, P. A. Knott, and G. Adesso, Phys. Rev. A 98, 012114 (2018).
  * Mehboudi  _et al._ (2019b) M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. A. García-March, and M. Lewenstein, Phys. Rev. Lett. 122, 030403 (2019b).
  * Planella  _et al._ (2020) G. Planella, M. F. B. Cenni, A. Acin, and M. Mehboudi, arXiv preprint arXiv:2001.11812 (2020).
  * Correa  _et al._ (2017) L. A. Correa, M. Perarnau-Llobet, K. V. Hovhannisyan, S. Hernández-Santana, M. Mehboudi, and A. Sanpera, Phys. Rev. A 96, 062103 (2017).
  * Hovhannisyan and Correa (2018) K. V. Hovhannisyan and L. A. Correa, Phys. Rev. B 98, 045101 (2018).
  * Mukherjee  _et al._ (2019) V. Mukherjee, A. Zwick, A. Ghosh, X. Chen, and G. Kurizki, Communications Physics 2, 1 (2019).
  * Martín-Martínez  _et al._ (2013) E. Martín-Martínez, A. Dragan, R. B. Mann, and I. Fuentes, New Journal of Physics 15, 053036 (2013).
  * Holevo (2011) A. S. Holevo,  _Probabilistic and Statistical Aspects of Quantum Theory_ , Publications of the Scuola Normale Superiore (Scuola Normale Superiore, Pisa, Italy, 2011).
  * Barnett and Radmore (2002) S. Barnett and P. M. Radmore,  _Methods in Theoretical Quantum Optics_ , Oxford Series in Optical and Imaging Sciences (Clarendon Press, 2002).
  * Serafini (2017b) A. Serafini,  _Quantum Continuous Variables: A Primer of Theoretical Methods_ (CRC Press, Taylor & Francis Group, Boca Raton, USA, 2017).
  * Newey and McFadden (1994) W. K. Newey and D. McFadden, Handbook of econometrics 4, 2111 (1994).
  * Oh  _et al._ (2019) C. Oh, C. Lee, C. Rockstuhl, H. Jeong, J. Kim, H. Nha, and S.-Y. Lee, npj Quantum Information 5, 1 (2019).
  * Morelli  _et al._ (2021) S. Morelli, A. Usui, E. Agudelo, and N. Friis, Quantum Science and Technology 6, 025018 (2021).
  * Kundu  _et al._ (2021) N. K. Kundu, M. R. McKay, and R. K. Mallik, “Machine learning based parameter estimation of gaussian quantum states,”  (2021), arXiv:2108.06061 [quant-ph] .
  * Giedke and Cirac (2002) G. Giedke and I. J. Cirac, Phys. Rev. A 66, 032316 (2002).
  * Williamson (1936) J. Williamson, Am. J. Math. 58, 141 (1936).
  * Lami  _et al._ (2017) L. Lami, C. Hirche, G. Adesso, and A. Winter, IEEE Trans. Inf. Theory 63, 7553 (2017).
  * Lampo  _et al._ (2017) A. Lampo, S. H. Lim, M. Á. García-March, and M. Lewenstein, Quantum 1, 30 (2017).
  * Khan  _et al._ (2021) M. M. Khan, M. Mehboudi, H. Tercas, M. Lewenstein, and M.-A. Garcia-March, arXiv preprint arXiv:2108.01754 (2021).
  * Rahimi-Keshari  _et al._ (2016) S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Phys. Rev. X 6, 021039 (2016).
  * Barnett  _et al._ (1998) S. M. Barnett, L. S. Phillips, and D. T. Pegg, Optics Communications 158, 45 (1998).
  * Rubio  _et al._ (2020) J. Rubio, J. Anders, and L. A. Correa, arXiv preprint arXiv:2011.13018 (2020).
  * Mehboudi  _et al._ (2021) M. Mehboudi, M. R. Jorgensen, S. Seah, J. B. Brask, J. Kolodynski, and M. Perarnau-Llobet,  (2021), arXiv:2108.05932 [quant-ph] .
  * Jorgensen  _et al._ (2021) M. R. Jorgensen, J. Kolodynski, M. Mehboudi, M. Perarnau-Llobet, and J. B. Brask,  (2021), arXiv:2108.05901 [quant-ph] .
  * Alves and Landi (2021) G. O. Alves and G. T. Landi, arXiv preprint arXiv:2106.12072 (2021).