# Operational significance of nonclassicality in nonequilibrium Gaussian
quantum thermometry

Safoura Mirkhalaf<EMAIL_ADDRESS>SSM and MM contributed equally.
Department of Physics, University of Tehran, P.O. Box 14395-547, Tehran, Iran
School of Nano Science, Institute for Research in Fundamental Sciences (IPM),
P.O. Box 19395-5531, Tehran, Iran  Centre for Quantum Optical Technologies,
Centre of New Technologies, University of Warsaw, Banacha 2c, 02-097 Warszawa,
Poland  Mohammad Mehboudi<EMAIL_ADDRESS>Technische
Universität Wien, 1020 Vienna, Austria  Département de Physique Appliquée,
Université de Genève, 1211 Genève, Switzerland  Zohre Nafari Qaleh  School of
Physics, Institute for Research in Fundamental Sciences (IPM), P.O. Box
19395-5531, Tehran, Iran  Saleh Rahimi-Keshari<EMAIL_ADDRESS>Department of
Physics, University of Tehran, P.O. Box 14395-547, Tehran, Iran  School of
Physics, Institute for Research in Fundamental Sciences (IPM), P.O. Box
19395-5531, Tehran, Iran

###### Abstract

We provide new operational significance of nonclassicality in nonequilibrium
temperature estimation of bosonic baths with Gaussian probe states and
Gaussian dynamics. We find a bound on the thermometry performance using
classical probe states. Then we show that by using nonclassical probe states,
single-mode and two-mode squeezed vacuum states, one can profoundly improve
the classical limit. Interestingly, we observe that this improvement can also
be achieved by using Gaussian measurements. Hence, we propose a fully Gaussian
protocol for enhanced thermometry, which can simply be realized and used in
quantum optics platforms.

##  I Introduction

Temperature estimation is essential in characterizing and engineering quantum
systems for any technological applications. Given finite experimental
resources, quantum thermometry designs probes with maximum information gain
Mehboudi  _et al._ (2019); De Pasquale and Stace (2018). Along with the number
of probes, time is an important resource. Without any time restriction, it is
often best to let a temperature probe thermalise with the system before
reading it out. In this procedure, known as equilibrium thermometry, the
initial state of the probe does not play a role and therefore the potential of
quantum resources in the state preparation is not used. Nonetheless, the
impact of initial probe states in a nonequilibrium (dynamical) scenario can be
important and has been the subject of previous studies Brunelli  _et al._
(2011, 2012); Correa  _et al._ (2015); Jevtic  _et al._ (2015); Razavian  _et
al._ (2019); De Pasquale  _et al._ (2017); Mitchison  _et al._ (2022); Cavina
_et al._ (2018); Guo  _et al._ (2015); Mitchison  _et al._ (2020); Seah  _et
al._ (2019); Shu  _et al._ (2020); Mancino  _et al._ (2020); Henao and Uzdin
(2021); Hovhannisyan  _et al._ (2021); Sekatski and Perarnau-Llobet (2021)

Continuous-variable systems in Gaussian states are the relevant model in many
physical platforms for thermometry such as quantum optics, quantum gases,
Josephson junctions, and mechanical resonators Weedbrook  _et al._ (2012). In
this context, a framework for using Gaussian measurements was recently
proposed with some preliminary results at thermal equilibrium Cenni  _et al._
(2022). On the other hand, in a recent paper Sekatski and Perarnau-Llobet
(2021), the limits of nonequilibrium thermometry in the Markovian environments
were established. The results of Sekatski and Perarnau-Llobet (2021) are based
on bounds in terms of operator norm that are only applicable to probe systems
with a finite-dimensional Hilbert space. Therefore, of particular interest is
to understand the limits and advantages of continuous-variable probe systems
with infinite-dimensional Hilbert spaces—that are readily available in the
laboratory in Gaussian states—in nonequilibrium thermometry.

In this paper, we propose a fully Gaussian setup with simple physical
realizations that exploits the power of nonclassicality for nonequilibrium
thermometry. Specifically, in our setup, the interaction between a continuous-
variable probe system and a thermal bath is described by a class of Gaussian
dynamics that can be described as a Brownian motion Mancino  _et al._ (2020);
Correa  _et al._ (2015); Schaller (2014). We first formalize the limit of
classical probes for thermometry and show that this limit can be profoundly
improved with the use of nonclassical Gaussian states and measurements. Our
formalism provides a new insight into the operational significance of
nonclassical states in nonequilibrium thermometry, and can readily be employed
e.g., in quantum optical platforms.

The structure of the paper is as follows. In Section II we formalise the
problem, the setup, and the figure of merit. In Section III, we review the
Gaussian formalism, including the Gaussian dynamics. We then proceed with our
main results in Section IV. Finally, in Section V, we close the paper with
some remarks. Derivations and details of some results, as well as more
simulations, are presented in the appendices in order to keep the paper’s main
part coherent.

##  II The setup

As shown in Fig. 1, the probe mode, which can be in an entangled state with an
auxiliary mode, interacts with a bath at temperature $T$. After some
interaction time $t$, the probe system is measured to infer the temperature of
the bath. We assume that the total time of running the experiment is fixed to
$\tau$. During this time, one can reset and repeat the whole process
$M=\tau/t$ times. In what follows we assume that the time required for
preparation and measurement protocols is not relevant compared to the time
required for the parametrization (the probe’s dynamics). Given that the
density matrix of the system right after parametrization reads ${\cal
E}_{t}(\rho)$, one performs a POVM measurement $\Pi$ with elements
$\\{\Pi_{k}\\}$. After $M$ rounds, we can collect the outcomes to the dataset
${\bf x}=\\{x_{1},\dots,x_{M}\\}$. An estimator function ${\tilde{T}}({\bf
x})$ maps the dataset into an estimate of the true parameter. From the Cramér-
Rao bound we know that the mean square error of any unbiased estimator is
bounded from below Cramér (1999); Newey and McFadden (1994)

$\displaystyle{\rm
MSE}({\tilde{T}};\rho;\Pi;t)\\!\coloneqq\\!\left<\\!\big{(}\\!{\tilde{T}}({\bf
x})-T_{0}\big{)}^{\\!2}\right>_{\\!\bf
x}\\!\\!\\!\geqslant\\!\frac{t}{\tau{\cal F}^{\rm C}(\rho;\Pi;t)},$ (1)

where the bound is saturable by choosing the maximum likelihood estimator.
Here, $T_{0}$ is the mean value of the estimator over all measurement
outcomes, say $T_{0}=\langle\tilde{T}(x)\rangle_{x}$; ${\cal F}^{\rm
C}(\rho;\Pi;t)$ is the classical Fisher information (CFI) associated with the
measurement $\\{\Pi\\}$ and the initial state $\rho$. The CFI is upper bounded
by the quantum Fisher information (QFI) defined as ${\cal F}^{\rm
Q}(\rho;t)\coloneqq\max_{\Pi}{\cal F}^{\rm C}(\rho;\Pi;t)$ Braunstein  _et
al._ (1996); Braunstein and Caves (1994); Paris (2009); Giovannetti  _et al._
(2011); Hayashi (2006). Motivated by (1) and following Chin  _et al._ (2012);
Correa  _et al._ (2015); Sekatski and Perarnau-Llobet (2021); del Campo  _et
al._ (2013) we set the rate of CFI, ${\tilde{\cal F}^{\rm
C}}(\rho;\Pi;t)\coloneqq{\cal F}^{\rm C}(\rho;\Pi;t)/t$, and the rate of QFI,
$\tilde{\cal F}^{\rm Q}(\rho;t)\coloneqq{\cal F}^{\rm Q}(\rho;t)/t$, as our
figures of merit. There are three main variables to optimize over, the initial
state $\rho$, the interrogation time $t$, and the measurement $\Pi$.

Figure 1: Our setting for estimating the temperature of a system: the probe is
prepared in a state, which may or may not be correlated to auxiliary units,
and then interacts with the sample for some time $t$. A measurement on the
probe (and possibly the auxiliary units) will reveal some information about
the temperature. One then resets the probe’s state and repeats this process
$M=\tau/t$ times, aiming to maximize the acquired information.

Our main aims are two folded (I) to understand the impact of Gaussian,
nonclassical probe states and (II) to study the performance of Gaussian
measurements. These are motivated by experimental feasibility; in quantum
optics, preparation of Gaussian states is always possible via single-mode
displacement and squeezing operations together with linear-optical networks,
and Gaussian measurements are readily available via homodyne detection and a
linear-optical network as well. We show that nonclassicality in Gaussian
states, namely single-mode squeezing and two-mode squeezing can significantly
improve the classical bound on thermometry precision. Moreover, Gaussian
measurements often perform as well as the optimal measurement, especially in
short time scales.

##  III Gaussian formalism

We assume that probes are initially prepared in Gaussian states and their
marginal states remain Gaussian through the interaction. Let us denote the
vector of canonical operators by $R=(x_{1},p_{1},\dots,x_{m},p_{m})^{T}$,
which also defines the matrix elements of the symplectic form
$[R_{j},R_{k}]\eqqcolon i\hbar\Omega_{jk}$. A Gaussian state is fully
described by the vector of the mean values $d\coloneqq{\rm Tr}[\rho R]$ and
the covariance matrix ($\hbar=1$)

$\displaystyle\sigma\coloneqq{\rm Tr}\\!\left[\rho R\circ
R^{T}\right]-2dd^{T},$ (2)

where we use the notation $R\circ R^{T}\equiv RR^{T}+(RR^{T})^{T}$ with
$RR^{T}$ being a matrix of operators that is not symmetric.

In the most general case, a Lindbladian Gaussian dynamics ${\cal E}_{t}$
evolves any probe state as Demoen  _et al._ (1977); Lindblad (2000); Cerf  _et
al._ (2007); Holevo and Werner (2001)

$\displaystyle d_{t}=X_{t}d,\hskip 28.45274pt\sigma_{t}=X_{t}\sigma
X^{T}_{t}+Y_{t}.$ (3)

Thus, the double
$(X_{t}\in{\mathbb{R}}^{2m}\times{\mathbb{R}}^{2m},Y_{t}\in{\mathbb{R}}^{2m}\times{\mathbb{R}}^{2m})$
fully characterise the dynamics. They satisfy $Y_{t}^{T}=Y_{t}$ and
$\Omega+Y+X_{t}\Omega X_{t}^{T}\geqslant 0$ to guarantee complete positivity
of the dynamics.

For the standard Brownian motion the temperature information enters the
dynamics only through the noise matrix $Y_{t}$ while the drift matrix $X_{t}$
is independent of temperature. This is the case in some situations such as the
damped harmonic oscillator or loss of an optical cavity and has been
previously used to study various problems in quantum thermometry Mancino  _et
al._ (2020); Correa  _et al._ (2015); Schaller (2014). In Appendix A, we
provide a microscopic derivation of this master equation where we prove the
term $X_{t}$ is indeed temperature independent. More specifically, such
Gaussian dynamical evolution can be described by $X_{t}=\exp(-\gamma
t/2)O_{t}$, and $Y_{t}=(1-\exp(-\gamma t))\sigma_{T}$, where $O_{t}$ is an
orthogonal matrix and $\sigma_{T}=\coth(\omega/2T)I_{2}$ is the covariance
matrix of a bosonic mode with frequency $\omega$ at thermal equilibrium with
temperature $T$. Here, $I_{2}$ is the $2\times 2$ identity matrix. The
parameter $\gamma=J(\omega)$ is related to the spectral density of the
environment $J(\omega)$ and generally depends on the system’s bare frequency
$\omega$; however, it is independent of the environment temperature. Also, in
our model, the orthogonal transformation $O_{t}$ corresponds to the coherent
dynamics of the probe system that is independent of temperature. Hence, we can
always work in the interaction picture with $O_{t}=I_{2}$ and incorporate the
phase rotation in time into the measurement. This also implies that, without
loss of generality, we can assume that the initial covariance matrix is
diagonal, as ${\cal E}_{t}(U(\theta)\rho U^{\dagger}(\theta))=U(\theta){\cal
E}_{t}(\rho)U^{\dagger}(\theta)$ and the phase rotation $U(\theta)$ can
further be absorbed into the measurement. In this case, the covariance matrix
$\sigma_{t}$ remains diagonal and is given by

$\sigma_{t}=e^{-\gamma t}\sigma+\big{(}1-e^{-\gamma t}\big{)}\nu I_{2}$ (4)

where $\nu=\coth(\omega/2T)$.

Having the Gaussian state of the probe at the interrogation time $t$, one can
calculate the QFI Monras (2013); Pinel  _et al._ (2013); Jiang (2014); Nichols
_et al._ (2018)

$\displaystyle\begin{split}{\cal F}^{\rm
Q}(\sigma,d;t)&\\!=2\partial_{T}d_{t}^{T}\sigma_{t}^{-1}\partial_{T}d_{t}+\frac{1}{2}\langle\partial_{T}\sigma_{t}|(\sigma_{t}\otimes\sigma_{t}-\Omega\otimes\Omega)^{-1}|\partial_{T}\sigma_{t}\rangle,\end{split}$
(5)

where $\partial_{T}$ denotes the partial derivative with respect to
temperature, and we used the vectorization notation
$|\partial_{T}\sigma_{t}\rangle\coloneqq([\partial_{T}\sigma_{t}]_{11},\dots,[\partial_{T}\sigma_{t}]_{1m},\dots,[\partial_{T}\sigma_{t}]_{mm})^{T}$.
Furthermore, if we perform a Gaussian measurement described by the covariance
matrix $\sigma^{M}$, then the corresponding CFI reads Cenni  _et al._ (2022);
Monras (2013); Malagò and Pistone (2015)

$\displaystyle{\cal F}^{\rm
C}(\sigma,d;\sigma^{M};t)=2\partial_{T}d_{t}^{T}(\sigma_{t}+\sigma^{M})^{-1}\partial_{T}d_{t}+\frac{1}{2}{\rm
Tr}\\!\left[\left(({\sigma}_{t}+{\sigma}^{M})^{-1}\partial_{T}{\sigma}_{t}\right)^{2}\right].$
(6)

However, as $X_{t}$ is temperature independent, no information about the
temperature is imprinted on $d_{t}$; see Fig. 2. Therefore, the first term in
the above expressions for QFI and CFI vanishes.

Figure 2: (a) Schematic representation of the probe’s evolution in phase
space. In our model, the quadrature mean values tend to zero (dashed black
arrow) in a temperature-independent manner. The temperature rather determines
how much noise is being added to the state (blue arrows). (b) If the probe is
in a highly squeezed-vacuum state, the squeezed quadrature (narrow dark blue
region) will be highly sensitive to temperature, while the temperature
information imprinted in the conjugate quadrature is almost negligible
compared with the initial large uncertainty. Therefore, homodyne measurement
on the squeezed quadrature—which is a Gaussian measurement—is the optimal
measurement.

##  IV Main results

Our main results are as follows. (i) The precision of using classical probe
states is upper bounded by that of using the vacuum probe state. (ii) This
bound can be overcome by using nonclassical states, namely single and two
mode-squeezed states. In particular, we prove that single-mode squeezed states
perform significantly better than the vacuum state, which previously was
thought to be the best reparation Correa  _et al._ (2015). Moreover,
entanglement between the probe and an auxiliary mode (two-mode squeezed
states) further improves the thermometry precision. (iii) We prove the (extent
of) optimality of Gaussian measurements, specifically homodyne detection for
non-equilibrium thermometry.

###  IV.1 Result (i): An upper bound for classical states

To begin with, we characterize the performance of classical states. A quantum
state is classical if and only if its Glauber-Sudarshan $P$-function is a
probability density distribution, and the state can be written as a mixture of
coherent states, $\rho^{\rm cl}\coloneqq\int d^{2}\alpha
P(\alpha)|\alpha\rangle\langle\alpha|$ Titulaer and Glauber (1965); Mandel
(1986). As discussed, in our model, no information about the temperature is
imprinted on the first-order moments of the probe state. Therefore, the
estimation precision of using all coherent probe states (including the vacuum
probe state), which have the same covariance matrix $\sigma_{1}=I_{2}$, are
the same. Also, as the Fisher information is convex—i.e., statistical mixture
reduces the Fisher information—we conclude that the performance of all
classical probe states is upper bounded by that of the vacuum probe state,
that is,

$\displaystyle{\cal F}^{\rm Q}(\rho^{\rm cl};t)\leqslant{\cal F}^{\rm
Q}(|\alpha\rangle\langle\alpha|;t)={\cal F}^{\rm Q}(|0\rangle\langle
0|;t),\hskip 5.69046pt\forall~{}t.$ (7)

Note that this result holds for nonGaussian classical states as well.

###  IV.2 Result (ii): Beating the classical bound with nonclassical states

As the first-order moments and a phase rotation do not matter, we consider the
squeezed-vacuum state with the covariance matrix $\sigma_{r}={\rm
diag}(r,1/r)$, where the squeezing parameter is between $r=0$ for infinitely
$x$-quadrature squeezed state and $r=1$ for the vacuum state, as potentially
optimal pure single-mode Gaussian probe state. Note that $p$-quadrature
squeezed states can be obtained by a phase rotation from the $x$-quadrature
squeezed states, so without loss of generality we assume $0\leqslant
r\leqslant 1$. By using Equations (4) and (5), we can find an analytical
expression for the QFI rate

$\displaystyle\tilde{\cal F}^{\rm Q}(\sigma_{r};t)=\frac{(1-e^{-\gamma
t})^{2}(\partial_{T}\nu)^{2}\left([\sigma_{t}]_{11}^{2}+[\sigma_{t}]_{22}^{2}+2\right)}{2t([\sigma_{t}]_{11}^{2}[\sigma_{t}]_{22}^{2}-1)}.$
(8)

In Figure. 3 we depict the (rate of the) QFI against time for the vacuum state
and for a squeezed state. It appears that the QFI rate is optimal when we
perform the measurement as quickly as possible. In the limit of short times
and high squeezing ($r\ll 1$) and by keeping the most relevant terms in $t$
and $r$, we have

$\displaystyle\tilde{\cal F}^{\rm Q}(\sigma_{r\ll 1};\gamma t\ll
1)=\frac{\gamma(\partial_{T}\nu)^{2}}{2(\gamma t\nu^{2}+2\nu r)}.$ (9)

It is clear from this expression that the optimal measurement time approaches
zero. Also, note that the rate can get arbitrarily large if we choose $\gamma
t$ and $r$ sufficiently small. This implies that nonclassicality in terms of
squeezing can significantly enhance precision. While the behavior of Fig. 3
depends on the temperature and the frequency (in particular the ratio
$\omega/T$ is determinant through the parameter $\nu=\coth(\omega/2T)$) our
numerics in Appendix B show that that nonclassicality in the form of squeezing
improves the QFI rate for other temperatures as well.

Figure 3: (a) The quantum Fisher information and (b) its rate are depicted in
solid. The probe is initially prepared in the vacuum or the squeezed state.
The (rate of the) classical Fisher information of optimal Gaussian
measurements are also plotted in dashed. All classical probe states have a
precision (or rate) that lies within the shaded green area, whose border is
determined by the vacuum state. One can see that squeezed states can overcome
this classical bound at short times. Interestingly, this improvement can still
be attained by using Gaussian measurements, which are readily available in the
laboratory. The optimal measurement should be performed as quickly as
possible. As for the optimal Gaussian measurement, however, one may wait until
the rate reaches its maximum, before measuring the system—better seen from
Fig. 4 that depicts short-time behavior. Here we set $T=1$, $\omega=1$,
$\gamma=0.2$, $r=10^{-3}$.

Let us now explore the role of entanglement. To this aim, we assume possession
of a secondary mode that does not undergo the dissipative dynamics but is
rather initially entangled with the probe system—see Fig. 1. The dynamics can
be simply extended to include the second mode by taking $X_{t}\to X_{t}\oplus
I_{2}$, and $Y_{t}\to Y_{t}\oplus 0_{2}$, where $0_{2}$ is the $2\times 2$
null matrix. Equation (4) is generalized as

$\displaystyle\sigma_{t}=\left(\begin{array}[]{cc}e^{-\gamma t}A+(1-e^{-\gamma
t})\nu I_{2}&e^{-\frac{\gamma t}{2}}C\\\ e^{-\frac{\gamma
t}{2}}C^{T}&B\end{array}\right),$ (12)

where $A$, $B$ and $C$ are $2\times 2$ matrices making up the initial
covariance matrix at $t=0$. As we can see, at long times the correlations
vanish. However, one can harness the correlations for enhanced thermometry at
short evolution times. As the initial input, we choose the two-mode squeezed
vacuum state given by

$\displaystyle\sigma_{2,r}\coloneqq\left(\begin{array}[]{ccc}r^{-1}\,{I}&\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}}\,{Z}\\\
\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower 0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}\,{Z}&r^{-1}\,{I}\end{array}\right),$ (15)

with $0<r\leqslant 1$ being the squeezing parameter and $Z={\rm diag}(1,-1)$.
This covariance matrix represents a pure entangled state iff $r<1$. Notice
that this scenario is at least as good as the one with a single mode squeezed
state with the same squeezing parameter $\sigma_{r}$. If one performs a local
Gaussian measurement described by $\sigma^{\rm M}$ on the auxiliary mode of
(15), the probe state collapses into Gaussian states with the covariance
matrix $\sigma_{{\rm PS}}=r^{-1}I-(r^{-2}-1)Z(r^{-1}I+\sigma^{\rm M})^{-1}Z$.
Setting $\sigma^{\rm M}=\lim_{s\to 0}{\rm diag}(s,1/s)$ (homodyne), the
covariance matrix of the probe states becomes $\sigma_{{\rm PS}}=\sigma_{r}$,
that is a single-mode squeezed state. However, by performing joint
measurements more information may be extracted.

By using Eq. (5), the QFI rate for the two-mode probe can be calculated and
compared with the result of the single-mode probe. Our simulations depicted in
Fig. 4 confirm that indeed at short times the QFI rate can gain up to a two-
fold improvement over the single mode squeezed state. Note, however, that this
improvement requires joint measurements, which may not be practical in
general. In the following, we show that by restricting to Gaussian
measurements, which are readily available, one can still exploit the
nonclassicality in probe states for thermometry precision.

Figure 4: Same as Fig. 3, but comparing the single-mode squeezed states with
the two-mode squeezed states for the probe. By using entanglement the QFI
saturates to its steady state value faster—at $\gamma t\gg 1$ the solid red
curve will catch up with the solid blue one—making them even better than
single mode squeezed states. Here, we examine the performance of a specific
joint Gaussian measurement, consisting of a $50:50$ beamsplitter followed by
two homodyne measurements. Evidently, even this simple Gaussian measurement
can exploit the improvement offered by the entangled Gaussian state, making
them comparable to the single mode squeezing interrogated with optimal
nonGaussian measurement. The rest of the parameters are taken similarly to
Fig. 3. $T=1$, $\gamma=0.2$, $r=10^{-3}$, $\omega=1$.

###  IV.3 Result (iii): Optimality of Gaussian measurements

We now restrict to a fully Gaussian scenario where the measurements are
Gaussian as well. In Appendix C, we analytically find the optimal Gaussian
measurement for single-mode Gaussian probe states, extending some of the
results of Cenni  _et al._ (2022) to dynamical scenarios. We prove that the
optimal Gaussian measurement has a diagonal covariance matrix $\sigma^{\rm
M}_{L}={\rm diag}(L,1/L)$, which becomes homodyne for short times. By using
Eqs. (4) and (6) for the single-mode squeezed vacuum probe state with
$\sigma_{r}$ and homodyne measurement on the $x$-quadrature ($L=0$), the CFI
rate is given by

$\displaystyle\tilde{\cal F}^{\rm C}(\sigma_{r};\sigma^{\rm
M}_{0};t)=\frac{(1-e^{-\gamma t})^{2}({\partial_{T}\nu})^{2}}{2t(re^{-\gamma
t}+(1-e^{-\gamma t})\nu)^{2}}.$ (16)

We can check that $\tilde{\cal F}^{\rm C}(\sigma_{r\ll 1};\sigma^{\rm
M}_{0};t\ll 1)\approx\gamma^{2}t(\partial_{T}\nu)^{2}/[2(r+\gamma t\nu)^{2}]$.
By comparing this with (9) and for $r\to 0$ one sees that $\tilde{\cal F}^{\rm
Q}(\sigma_{r\ll 1};t\ll 1)\approx\tilde{\cal F}^{\rm C}(\sigma_{r\ll
1};\sigma^{\rm M}_{0};t\ll 1)\approx(\partial_{T}\nu)^{2}/(2t\nu^{2})$ which
suggests that homodyne measurement is the optimal measurement in the regime of
large squeezing and at short times. For a finite squeezing parameter, unlike
the QFI rate, the best interrogation time is not zero and roughly satisfies
$t^{*}={r}/(\gamma\nu)$, which decreases by increasing the squeezing and the
dissipation rate $\gamma$. We have $\tilde{\cal F}^{\rm
C}(\sigma_{r};\sigma^{\rm
M}_{0};t^{*})\approx\gamma(\partial_{T}\nu)^{2}/(8r\nu)$.

Let us now discuss the role of entanglement in two-mode squeezed vacuum states
using joint Gaussian measurements. In general, joint Gaussian measurements can
be modeled in terms of two local phase rotation operations, a beamsplitter and
two single-mode Gaussian measurements Rahimi-Keshari  _et al._ (2015).
However, instead of optimizing over all joint Gaussian measurements, we
consider a measurement that boosts the performance with respect to the single-
mode scenario. Specifically, we take $\sigma^{\rm M}_{2,0}\coloneqq S_{\rm
BS}\sigma_{xp}S_{BS}^{T}$ where the probe mode after the interaction overlaps
with the auxiliary mode on a balanced beamsplitter with the symplectic
transformation $S_{\rm BS}$. Then at the output two homodyne measurements, one
on $x$-quadrature and one on $p$-quadrature with the CM
$\sigma_{xp}=\lim_{L\to 0}{\rm diag}(L,1/L,1/L,L)$, are performed. In this
case, the CFI rate reads (see Appendix D for details)

$\displaystyle\tilde{\cal F}^{\rm C}(\sigma_{2,r};\sigma^{\rm
M}_{2,0};t)=\frac{(\partial_{T}\lambda)^{2}}{t\lambda^{2}},$ (17)

where $\lambda\coloneqq(1+e^{-\gamma t})r^{-1}+\nu(1-e^{-\gamma
t})-2e^{\frac{-\gamma
t}{2}}\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}}$. As depicted in Fig. 4
this fully Gaussian setting can exploit entanglement to overperform the best
(nonclassical) single mode setting, including nonGaussian measurements. At
short times and for small values of $r$, we have

$\displaystyle\tilde{\cal F}^{\rm C}(\sigma_{2,r};\sigma^{\rm M}_{2,0};t\ll
1)\approx\frac{16r^{2}\gamma^{2}t(\partial_{T}\nu)^{2}}{[(\gamma
t)^{2}+4r^{2}+4r\gamma t\nu]^{2}},$ (18)

which shows the optimal time is $t^{*}\approx{r}/({\gamma\nu})$. We then have
$\tilde{\cal F}^{\rm C}(\sigma_{2,r};\sigma^{\rm M}_{2,0};t^{*})\approx
16\gamma\nu^{3}(\partial_{T}\nu)^{2}/(r(1+8\nu^{2})^{2})$, which roughly shows
two times improvement in the CFI rate compared to the single-mode squeezed
probe state and homodyne measurement.

##  V Discussion

Our formalism for nonequilibrium thermometry of bosonic baths by means of
continuous-variable systems can be used to describe many platforms such as
impurities in cold gases, mechanical resonators, Josephson junctions, and
quantum-optical systems. We set a bound on the thermometric performance of all
classical states (i.e., those with non-negative $P$-function), based on
showing that all coherent states are equally useful. We then showed that
nonclassical squeezed-vacuum states can beat this bound. An interesting lesson
is that, increasing the energy of probes in terms of the displacement in phase
space (or the amplitude of oscillation) is not beneficial for thermometry. But
rather the energy should be put into increasing the degree of squeezing. We
have shown, in particular, that fully Gaussian scenarios, which can simply be
implemented in the laboratory in terms of single-mode and two-mode squeezed
vacuum probe states and homodyne measurements, are extremely useful in
boosting thermometry precision when time is a resource.

Our results also complement recent developments in nonequilibrium thermometry
Sekatski and Perarnau-Llobet (2021), which describes a fully non-restricted
scenario. There are three major differences between the two scenarios. (1) In
Sekatski and Perarnau-Llobet (2021), the probe system has a finite
dimensional-Hilbert space and the optimal preparation is always the eigenstate
of the Hamiltonian with the maximum energy. Extending this to the infinite-
dimensional Hilbert space with unbounded operators would be non-trivial.
However, in our method, we consider the Gaussian states of a continuous-
variable probe system with an infinite-dimensional Hilbert space. Gaussian
states can be easily prepared using linear optics, and we find the optimal
probe state within the Gaussian family. An interesting future direction is to
generalize the upper bound of Sekatski and Perarnau-Llobet (2021) to
continuous variable probes by considering a constraint on the average energy
of the probe states. (2) When there is no restriction on the measurement, it
is best to perform the measurement as fast as possible. However, when
restricting to the Gaussian measurements, there exists an optimal measurement
time—which as we have shown scales with the initial squeezing in the state of
the probe. This is because while the QFI for the optimal input state grows
linearly with time, the CFI of Gaussian measurements grows quadratically with
time (except in the unphysical case where the probe is initially infinitely
squeezed and the CFI scales linearly with time). (iii) unlike the findings of
Sekatski and Perarnau-Llobet (2021), there is no Lamb-shift effect in our
model since the Lamb-shift is temperature independent—See Appendix A. We
expect that our formalism finds applications in estimating the temperature in
quantum optical platforms or within cold Bosonic gases where state-of-the-art
techniques allow for the implementation of Gaussian measurements such as
homodyne/heterodyne measurements or time-of-flight measurements.

## Acknowledgements

We are thankful to P. Sekatski and M. Perarnau-Llobet for constructive
discussions. MM was supported the DFG/FWF Research Unit FOR 2724 ‘Thermal
machines in the quantum world’ and by Swiss National Science Foundation NCCR
SwissMAP. SM is supported by a grant from Basic Sciences Research, Iran, Fund
(No. BSRF-phys-399-05).

##  Appendix A Case study: The quantum Brownian motion

Consider a Harmonic oscillator—playing the role of a Brownian particle, namely
an impurity atoms or a light mode—that is embedded inside a bosonic
bath—representing, e.g., cold bosonic gases or radiation fields. The Brownian
particle is used as our probe to take the temperature of the bath. The total
Hamiltonian reads

$\displaystyle H=H_{p}+H_{B}+H_{I}.$ (19)

Here, $H_{p}$ is the free Hamiltonian of the oscillator

$\displaystyle H_{p}$ $\displaystyle=\omega_{0}a^{\dagger}a,$ (20)

with $a$ and $a^{\dagger}$ being the Bosonic creation and annihilation
operators and $\omega_{0}$ is the bare frequency of the probe. Moreover,
$H_{B}$ is the Hamiltonian of the bath given by

$\displaystyle H_{B}=\sum_{\mu}\omega_{\mu}b^{\dagger}_{\mu}b_{\mu},$ (21)

with $a_{\mu}$, $b^{\dagger}_{\mu}$, and $\omega_{\mu}$ being the creation,
annihilation,and the bare frequency of the $\mu$th oscillator in the bath.
Note that the creation and annihilation operators satisfy the commutation
relation $[a,a^{\dagger}]=1$, $[b_{\mu},b^{\dagger}_{\nu}]=\delta_{\mu,\nu}$.
Finally, $H_{I}$ is the linear interaction between the probe and the bath that
we assume to read

$\displaystyle H_{I}$
$\displaystyle=a^{\dagger}\otimes\sum_{\mu}g_{\mu}b_{\mu}+\text{ h.c.},$ (22)

where $g_{\mu}$ is the coupling strength and h.c. refers to Hermitian
conjugate. Note that by defining the quadrature operators in the usual way,
one can rewrite the interaction Hamiltonian in terms of Hermitian canonical
operators

$H_{I}=\sum_{k=1}^{2}A_{k}\otimes B_{k}$ (23)

where

$\displaystyle A_{1}$
$\displaystyle=\frac{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}}}}(a^{\dagger}+a)=x,\qquad
B_{1}=\frac{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}}}}\sum_{\mu}g_{\mu}(b_{\mu}^{\dagger}+b_{\mu})=\sum_{\mu}g_{\mu}x_{\mu},$
$\displaystyle\hskip 5.69046ptA_{2}$
$\displaystyle=\frac{i}{\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}}}}(a^{\dagger}-a)=p,\qquad
B_{2}=\frac{i}{\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}}}}\sum_{\mu}g_{\mu}(b^{\dagger}_{\mu}-b_{\mu})=\sum_{\mu}g_{\mu}p_{\mu}.$

In the limit of continuous bath modes, we will use the following general
description for the spectral density

$\displaystyle J(\omega)=2\pi\sum_{k}g^{2}_{k}\delta(\omega-\omega_{k}).$ (24)

For now, we do not assume any specific shape for the spectral density and keep
the analysis as general as possible. As we prove below, the
Gorini–Kossakowski–Sudarshan–Lindblad (GKLS) master equation for $\rho_{p}$
the state of probe mode reads

$\displaystyle\frac{d\rho_{p}}{dt}=-i\big{[}H_{p}~{}+\Delta
H_{p}\,,\rho_{p}\big{]}+\Gamma_{\rm
out}\bigg{(}\\!a\,\rho_{p}\,a^{\dagger}-\frac{1}{2}\big{\\{}a^{\dagger}a,\rho_{p}\big{\\}}\\!\bigg{)}+\Gamma_{\rm
in}\bigg{(}\\!a^{\dagger}\rho_{p}\,a-\frac{1}{2}\big{\\{}aa^{\dagger},\rho_{p}\big{\\}}\\!\bigg{)},$
(25)

with $\Gamma_{\rm out}\coloneqq J(\omega_{0})(N(\omega_{0},T)+1)$,
$\Gamma_{\rm in}\coloneqq J(\omega_{0})N(\omega_{0},T)$, and
$N(\omega_{0},T)=[\exp(\omega_{0}/T)-1]^{-1}$. The term $\Delta H_{P}$
represents a [Lamb] shift, which as we show below is temperature independent
and thus the coherent part has no impact on the thermometry performance. In
the following subsection, we derive the master equation governing the damped
harmonic oscillator dynamics in two ways, and particularly prove that the
shift term is independent of temperature.

###  A.1 Derivation of the GKLS master equation (25)

In the first approach, one can start from the von-Neumann equation Rivas and
Huelga (2011); Breuer and Petruccione (2002). We start by writing the
evolution of the total system in the interaction picture

$\displaystyle\dfrac{d\bm{\rho}(t)}{dt}=-i\big{[}\bm{H}_{I}(t),\bm{\rho}(t)\big{]},$
(26)

where we use boldface letters for the operators in the interaction picture

$\displaystyle\bm{\rho}(t)$ $\displaystyle=U^{\dagger}_{p}(t,0)\otimes
U^{\dagger}_{B}(t,0)~{}\rho(t)~{}U_{p}(t,0)\otimes U_{B}(t,0),$ (27)
$\displaystyle\bm{H}_{I}(t)$ $\displaystyle=U^{\dagger}_{p}(t,0)\otimes
U^{\dagger}_{B}(t,0)~{}H_{I}~{}U_{p}(t,0)\otimes U_{B}(t,0).$ (28)

with $U_{p}(t,0)=e^{-itH_{p}}$ and $U_{B}(t,0)=e^{-itH_{B}}$, and $\rho(t)$ is
the joint state of the probe and bath systems. By integrating Eq. (26) we
obtain

$\displaystyle\bm{\rho}(t)=\bm{\rho}(0)-i\int_{0}^{t}ds\big{[}\bm{H}_{I}(s),\bm{\rho}(s)\big{]},$
(29)

which by recursive replacement reads

$\displaystyle\bm{\rho}(t)=\bm{\rho}(0)-i\int_{0}^{t}ds\big{[}\bm{H}_{I}(s),\bm{\rho}(0)\big{]}-\int_{0}^{t}ds\int_{0}^{s}dr\big{[}\bm{H}_{I}(s),\left[\bm{H}_{I}(r),\bm{\rho}(r)\right]\\!\big{]}.$
(30)

If we differentiate again, we have

$\displaystyle\dfrac{d\bm{\rho}(t)}{dt}=-i\big{[}\bm{H}_{I}(t),\bm{\rho}(0)\big{]}-\int^{t}_{0}ds\big{[}\bm{H}_{I}(t),\left[\bm{H}_{I}(s),\bm{\rho}(s)\right]\\!\big{]}.$
(31)

Now, by taking partial trace over the degrees of freedom of bath we have

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}=-i\,\mathrm{Tr}_{B}\\!\big{[}\bm{H}_{I}(t),\bm{\rho}(0)\big{]}-\int^{t}_{0}ds\,\mathrm{Tr}_{B}\\!\big{[}\bm{H}_{I}(t),\left[\bm{H}_{I}(s),\bm{\rho}(s)\right]\\!\big{]}.$
(32)

Considering that the initial state of the total system is a product state
$\bm{\rho}(0)=\bm{\rho}_{p}(0)\otimes\bm{\rho}_{B}$, and that the bath is
initially in a state which commutes with $H_{B}$ such as the thermal state,
i.e. $\rho_{B}=e^{-H_{B}/T}/\mathrm{Tr}_{B}(e^{-H_{B}/T})$, the first term of
Eq. (32) vanishes Breuer and Petruccione (2002). In this step, we employ the
Born approximation, where the coupling between the probe and the bath is
considered weak such that the effect of the system on the bath can be
neglected; this implies that for the sufficiently large bath one can assume
the bath state is time-independent and the probe-bath remain uncorrelated at
all times, $\bm{\rho}(s)\approx\bm{\rho}_{p}(s)\otimes\bm{\rho}_{B}$ Breuer
and Petruccione (2002); Lidar (2020). Thus, Eq. (32) becomes

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}$
$\displaystyle=-\int^{t}_{0}ds\,\mathrm{Tr}_{B}\\!\big{[}\bm{H}_{I}(t),\left[\bm{H}_{I}(s),\bm{\rho}_{p}(s)\otimes\bm{\rho}_{B}\right]\\!\big{]}$
$\displaystyle=-\int^{t}_{0}ds\,\mathrm{Tr}_{B}\\!\Big{[}\bm{H}_{I}(t)\bm{H}_{I}(t-s)\bm{\rho}_{p}(t-s)\otimes\bm{\rho}_{B}-\bm{H}_{I}(t-s)\bm{\rho}_{p}(t-s)\otimes\bm{\rho}_{B}\bm{H}_{I}(t)+\text{h.c.}\Big{]}$
$\displaystyle=-\sum_{kk^{\prime}}\int_{0}^{t}ds\Big{(}\\!{\bm{A}}_{k}(t){\bm{A}}_{k^{\prime}}(t-s)\bm{\rho}_{p}(t-s)\\!-\\!{\bm{A}}_{k^{\prime}}(t-s)\bm{\rho}_{p}(t-s){\bm{A}}_{k}(t)\\!\Big{)}\mathrm{Tr}_{B}\\!\big{[}{\bm{B}}_{k}(t){{\bm{B}}_{k^{\prime}}(t-s)\bm{\rho}_{B}}\big{]}+\text{h.c.},$
$\displaystyle=-\sum_{kk^{\prime}}\int_{0}^{t}ds\Big{(}\\!{\bm{A}}_{k}(t){\bm{A}}_{k^{\prime}}(t-s)\bm{\rho}_{p}(t)\\!-\\!{\bm{A}}_{k^{\prime}}(t-s)\bm{\rho}_{p}(t){\bm{A}}_{k}(t)\\!\Big{)}\mathrm{Tr}_{B}\\!\big{[}{\bm{B}}_{k}(t){{\bm{B}}_{k^{\prime}}(t-s)\bm{\rho}_{B}}\big{]}+\text{h.c.},$
(33)

where we have applied $s\rightarrow t-s$ in the second line; in the third line
we used the general form of the interaction Hamiltonian in the Schrödinger
picture as $H_{I}=\sum_{k}A_{k}\otimes B_{k}$, with $A_{k}$ and $B_{k}$
defined in Eq. (23), and thus $\bm{A}_{k}(t)$ and $\bm{B}_{k}(t)$ are
corresponding operators in the interaction picture

$\displaystyle\begin{split}\bm{A}_{1}(t)&=\frac{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}}}}\Big{(}e^{i\omega_{0}t}a^{\dagger}+e^{-i\omega_{0}t}a\Big{)},\qquad\bm{B}_{1}(t)=\dfrac{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}}}}\sum_{\mu}g_{\mu}\Big{(}e^{it\omega_{\mu}}b_{\mu}^{\dagger}+e^{-it\omega_{\mu}}b_{\mu}\Big{)},\\\
\bm{A}_{2}(t)&=\frac{i}{\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}}}}\Big{(}e^{i\omega_{0}t}a^{\dagger}-e^{-i\omega_{0}t}a\Big{)},\qquad\bm{B}_{2}(t)=\dfrac{i}{\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}}}}\sum_{\mu}g_{\mu}\Big{(}e^{it\omega_{\mu}}b_{\mu}^{\dagger}-e^{-it\omega_{\mu}}b_{\mu}\Big{)};\end{split}$
(34)

finally, in the last line, we use the Markov approximation by replacing
$\bm{\rho}_{p}(t-s)\to\bm{\rho}_{p}(t)$.

To proceed further, we assume that $t\gg\tau_{B}$, where the bath correlation
time scale $\tau_{B}$ is defined as the timescale beyond which the bath two-
time correlation functions
$\mathrm{Tr}_{B}[{{\bm{B}}_{k^{\prime}}(t-s){\bm{B}}_{k}(t)\bm{\rho}_{B}}]$
decay rapidly Lidar (2020); Albash  _et al._ (2012). As a result, we can take
the upper limit of the integral to infinity. We have

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}=\sum_{kk^{\prime}}\int_{0}^{\infty}\\!ds\Big{(}\\!{\bm{A}}_{k^{\prime}}(t-s)\bm{\rho}_{p}(t){\bm{A}}_{k}(t)\\!-\\!{\bm{A}}_{k}(t){\bm{A}}_{k^{\prime}}(t-s)\bm{\rho}_{p}(t)\\!\Big{)}\mathrm{Tr}_{B}\\!\big{[}{\bm{B}}_{k}(t){{\bm{B}}_{k^{\prime}}(t-s)\bm{\rho}_{B}}\big{]}+\text{h.c.}.$
(35)

By using the Bosonic relations

$\displaystyle\big{\langle}b^{\dagger}_{k}b_{l}\big{\rangle}$
$\displaystyle\equiv\mathrm{Tr}_{B}\\!\big{[}\rho_{B}b^{\dagger}_{k}b_{l}\big{]}=\delta_{kl}N(\omega_{k},T),$
(36) $\displaystyle\big{\langle}b_{k}b^{\dagger}_{l}\big{\rangle}$
$\displaystyle=\delta_{kl}\big{(}N(\omega_{k},T)+1\big{)},$ (37)
$\displaystyle\big{\langle}b^{\dagger}_{k}b^{\dagger}_{l}\big{\rangle}$
$\displaystyle=\big{\langle}b_{k}b_{l}\big{\rangle}=0,$ (38)

where, $N(\omega_{k},T)$ is the average number of the $k$-th mode, we put
together Eqs. (34) and (35) and obtain

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}=\sum_{k}g_{k}^{2}\int_{0}^{\infty}ds\Big{(}$
$\displaystyle
a^{\dagger}\bm{\rho}_{p}(t)aN(\omega_{k},T)e^{-is(\omega_{0}-\omega_{k})}+a\bm{\rho}_{p}(t)a^{\dagger}\big{(}N(\omega_{k},T)+1\big{)}e^{is(\omega_{0}-\omega_{k})}$
$\displaystyle-$ $\displaystyle
a^{\dagger}a\bm{\rho}_{p}(t)\big{(}N(\omega_{k},T)+1\big{)}e^{is(\omega_{0}-\omega_{k})}-aa^{\dagger}\bm{\rho}_{p}(t)N(\omega_{k},T)e^{-is(\omega_{0}-\omega_{k})}\Big{)}+\text{h.c.}.$
(39)

Now by using the following formula

$\displaystyle\int_{0}^{\infty}dse^{-is(\omega_{0}-\omega_{k})}=\pi\delta(\omega_{0}-\omega_{k})-i\mathbb{P}\dfrac{1}{\omega_{0}-\omega_{k}},$
(40)

where $\mathbb{P}$ denotes the Cauchy principal value, and also using the
spectral density function (24), equation (39) simplifies to

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}$
$\displaystyle=\frac{1}{2}a^{\dagger}\bm{\rho}_{p}(t)a\left(\\!J(\omega_{0})N(\omega_{0},T)-\frac{i}{\pi}\mathbb{P}\int_{-\infty}^{\infty}\\!d\omega\,\dfrac{J(\omega)N(\omega,T)}{\omega_{0}-\omega}\right)$
$\displaystyle~{}+\frac{1}{2}a\bm{\rho}_{p}(t)a^{\dagger}\left(\\!J(\omega_{0})\left[N(\omega_{0},T)+1\right]+\frac{i}{\pi}\mathbb{P}\int_{-\infty}^{\infty}d\omega\,\dfrac{J(\omega)\left[N(\omega,T)+1\right]}{\omega_{0}-\omega}\right)$
$\displaystyle~{}-\frac{1}{2}a^{\dagger}a\bm{\rho}_{p}(t)\left(\\!J(\omega_{0})\left[N(\omega_{0},T)+1\right]+\frac{i}{\pi}\mathbb{P}\int_{-\infty}^{\infty}\\!d\omega\,\dfrac{J(\omega)\left[N(\omega,T)+1\right]}{\omega_{0}-\omega}\right)$
$\displaystyle~{}-\frac{1}{2}aa^{\dagger}\bm{\rho}_{p}(t)\left(\\!J(\omega_{0})N(\omega_{0},T)-\frac{i}{\pi}\mathbb{P}\int_{-\infty}^{\infty}\\!d\omega\,\dfrac{J(\omega)N(\omega,T)}{\omega_{0}-\omega}\right)+\text{h.c.}~{}.$
(41)

Finally, the dynamics of the damped harmonic oscillator in the interaction
picture reads

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}=$
$\displaystyle-i\Delta_{1}\big{[}a^{\dagger}a,\bm{\rho}_{p}(t)\big{]}+i\Delta_{2}\left[aa^{\dagger},\bm{\rho}_{p}(t)\right]+J(\omega_{0})\big{(}N(\omega_{0},T)+1\big{)}\Big{(}a\bm{\rho}_{p}(t)a^{\dagger}-\dfrac{1}{2}\big{\\{}a^{\dagger}a,\bm{\rho}_{p}(t)\big{\\}}\\!\Big{)}$
$\displaystyle+J(\omega_{0})N(\omega_{0},T)\Big{(}a^{\dagger}\bm{\rho}_{p}(t)a-\dfrac{1}{2}\big{\\{}aa^{\dagger},\bm{\rho}_{p}(t)\big{\\}}\\!\Big{)},$
(42)

where

$\displaystyle\Delta_{1}$
$\displaystyle=\frac{1}{2\pi}\mathbb{P}\int_{-\infty}^{\infty}\\!d\omega\,\dfrac{J(\omega)\left[N(\omega,T)+1\right]}{\omega_{0}-\omega},$
(43) $\displaystyle\Delta_{2}$
$\displaystyle=\frac{1}{2\pi}\mathbb{P}\int_{-\infty}^{\infty}\\!d\omega\,\dfrac{J(\omega)N(\omega,T)}{\omega_{0}-\omega}.$
(44)

Since $[a,a^{\dagger}]=\mathbb{I}$, then:

$\displaystyle\dfrac{d\bm{\rho}_{p}(t)}{dt}=$
$\displaystyle-i\big{[}\left(\Delta_{1}-\Delta_{2}\right)a^{\dagger}a\,,\bm{\rho}_{p}(t)\big{]}+J(\omega_{0})\big{(}N(\omega_{0},T)+1\big{)}\Big{(}a\bm{\rho}_{p}(t)a^{\dagger}-\dfrac{1}{2}\big{\\{}a^{\dagger}a,\bm{\rho}_{p}(t)\big{\\}}\\!\Big{)}$
$\displaystyle+J(\omega_{0})N(\omega_{0},T)\Big{(}a^{\dagger}\bm{\rho}_{p}(t)a-\dfrac{1}{2}\big{\\{}aa^{\dagger},\bm{\rho}_{p}(t)\big{\\}}\\!\Big{)}.$
(45)

By using Eq. (27) the master equation in the Schrödinger picture can also be
obtained

$\displaystyle\dfrac{d\rho_{p}(t)}{dt}=-i\big{[}H_{\rm
LS}+\omega_{0}a^{\dagger}a,\rho_{p}(t)\big{]}+\mathcal{D}[\rho_{p}(t)],$ (46)

where

$\displaystyle\mathcal{D}[\rho_{p}(t)]=J(\omega_{0})\big{(}N(\omega_{0},T)+1\big{)}\Big{(}a\rho_{p}(t)a^{\dagger}-\dfrac{1}{2}\big{\\{}a^{\dagger}a,\rho_{p}(t)\big{\\}}\\!\Big{)}+J(\omega_{0})N(\omega_{0},T)\Big{(}a^{\dagger}\rho_{p}(t)a-\dfrac{1}{2}\big{\\{}aa^{\dagger},\rho_{p}(t)\big{\\}}\\!\Big{)},$
(47)

is the dissipator of the master equation. and

$\displaystyle H_{\rm LS}=(\Delta_{1}-\Delta_{2})a^{\dagger}a,$ (48)

is the Lamb-shift Hamiltonian which leads to a renormalization of bare
particle Hamiltonian due to the particle-bath interaction. In general, it has
been shown that the Lamb-shift terms have a significant impact on the dynamics
Nafari Qaleh and Rezakhani (2022), and even on thermometry Sekatski and
Perarnau-Llobet (2021). However, in our case, since

$\displaystyle\Delta_{1}-\Delta_{2}=\mathbb{P}\int_{-\infty}^{\infty}d\omega\dfrac{J(\omega)}{\omega_{0}-\omega},$
(49)

it is clear that the Lamb-shift term is independent of $N(\omega_{0},T)$ and
therefore independent of the bath temperature. Thus, the presence of the Lamb
shift does not affect the thermometry performance, in particular, the first
term in Eqs. (5) and (6) disappear and our no-go result holds.

Note that a similar result is expected if the interaction Hamiltonian was
considered to be of $H_{I}=(a+a^{\dagger})\sum_{\mu}g_{\mu}(b+b^{\dagger})$,
which can be effectively reduced to our model by means of the secular
approximation. See also Ref. Correa and Glatthard (2023) for a thorough
discussion on Lamb-Shift and its impact on the dynamics for the same model
with an Ohmic spectral density with Lorentz-Drude cutoff.

In the second approach for deriving the master equation (25), one can use the
standard closed form of the Lindblad equation to describe the time evolution
of the damped harmonic oscillator Breuer and Petruccione (2002); Lidar (2020).
For this purpose, we start from the following equation in the Schrödinger
picture

$\displaystyle\dfrac{d\rho_{p}(t)}{dt}=-i\big{[}H_{p},\rho_{p}(t)\big{]}+\sum_{\omega}\sum_{k,k^{\prime}}\Gamma_{kk^{\prime}}(\omega)\Big{(}\\!A_{k^{\prime}}(\omega)\rho_{p}(t)A_{k}^{\dagger}(\omega)-A_{k}^{\dagger}(\omega)A_{k^{\prime}}(\omega)\rho_{p}(t)\\!\Big{)}+\text{h.c.},$
(50)

where we assumed the interaction Hamiltonian (23) and define the coefficients
$\Gamma_{kk^{\prime}}(\omega)$ and operators $A_{k}(\omega)$ in the following.
In the above equation, we introduce the one-sided Fourier transformation of
the bath correlation functions as

$\displaystyle\Gamma_{kk^{\prime}}(\omega)=\int_{0}^{\infty}\\!ds\,e^{i\omega
s}\left\langle{\bm{B}}_{k}(s){\bm{B}}_{k^{\prime}}(0)\right\rangle,$ (51)

where
$\left\langle{\bm{B}}_{k}(s){\bm{B}}_{k^{\prime}}(0)\right\rangle=\mathrm{Tr}_{B}[{\bm{B}}_{k}(s){\bm{B}}_{k^{\prime}}(0)\rho_{B}]$
with ${\bm{B}}_{k}(t)$ being the bath operators in the interaction picture
(34). Substituting ${\bm{B}}_{1}(t)$ and ${\bm{B}}_{2}(t)$ into Eq. (51) and
using Eqs. (36)-(38) we find

$\displaystyle\Gamma_{11}(\omega)=\Gamma_{22}(\omega)=\sum_{k}\dfrac{g^{2}_{k}}{2}\bigg{(}\\!N(\omega_{k},T)\\!\int_{0}^{\infty}dse^{i(\omega+\omega_{k})s}+\big{(}N(\omega_{k},T)+1\big{)}\\!\int_{0}^{\infty}dse^{i(\omega-\omega_{k})s}\\!\bigg{)},$
(52)
$\displaystyle\Gamma_{12}(\omega)=-\Gamma_{21}(\omega)=\sum_{k}-i\dfrac{g^{2}_{k}}{2}\bigg{(}\\!N(\omega_{k},T)\\!\int_{0}^{\infty}dse^{i(\omega+\omega_{k})s}-\big{(}N(\omega_{k},T)+1\big{)}\\!\int_{0}^{\infty}dse^{i(\omega-\omega_{k})s}\\!\bigg{)}.$
(53)

By considering a continuous spectrum for the bath and using spectral density
function (24) and Eq. (40), we get

$\displaystyle\Gamma_{11}(\omega)=\Gamma_{22}(\omega)=\frac{1}{4}\big{(}N(\omega,T)+1\big{)}\big{(}J(\omega)-J(-\omega)\big{)}+\dfrac{i}{4\pi}\mathbb{P}\\!\int_{-\infty}^{\infty}\\!d\omega^{\prime}J(\omega^{\prime})\bigg{(}\dfrac{N(\omega^{\prime},T)}{\omega^{\prime}+\omega}+\dfrac{N(\omega^{\prime},T)+1}{\omega^{\prime}-\omega}\bigg{)},$
(54)
$\displaystyle\Gamma_{12}(\omega)=-\Gamma_{21}(\omega)=\frac{i}{4}\big{(}N(\omega,T)+1\big{)}\big{(}J(\omega)+J(-\omega)\big{)}+\dfrac{1}{4\pi}\mathbb{P}\\!\int_{-\infty}^{\infty}\\!d\omega^{\prime}J(\omega^{\prime})\bigg{(}\dfrac{N(\omega^{\prime},T)}{\omega^{\prime}+\omega}-\dfrac{N(\omega^{\prime},T)+1}{\omega^{\prime}-\omega}\bigg{)}.$
(55)

By using the eigenvectors and eigenvalues of the probe Hamiltonian,
$H_{p}|n\rangle=(n+\dfrac{1}{2})\omega_{0}|n\rangle$, the operators
$A_{k}(\omega)$ in (50) are also defined as

$\displaystyle
A_{k}(\omega)\coloneqq\sum_{(n^{\prime}-n)\omega_{0}=\omega}\\!\\!|n\rangle\langle
n|A_{k}|n^{\prime}\rangle\langle n^{\prime}|.$ (56)

For $A_{1}=(a^{\dagger}+a)/\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}}}$, we have

$\displaystyle A_{1}(\omega)$
$\displaystyle=\dfrac{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}}}}\sum_{(n^{\prime}-n)\omega_{0}=\omega}\\!\\!|n\rangle\langle
n|(a^{\dagger}+a)|n^{\prime}\rangle\langle
n^{\prime}|=\dfrac{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}}}}\\!\sum_{(n^{\prime}-n)\omega_{0}=\omega}\\!\\!\\!\Big{(}\\!\mathchoice{{\hbox{$\displaystyle\sqrt{n^{\prime}+1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{n^{\prime}+1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{n^{\prime}+1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{n^{\prime}+1\,}$}\lower
0.4pt\hbox{\vrule height=3.26389pt,depth=-2.61113pt}}}|n\rangle\langle
n^{\prime}|\delta_{n,n^{\prime}+1}+\mathchoice{{\hbox{$\displaystyle\sqrt{n^{\prime}\,}$}\lower
0.4pt\hbox{\vrule
height=5.8611pt,depth=-4.6889pt}}}{{\hbox{$\textstyle\sqrt{n^{\prime}\,}$}\lower
0.4pt\hbox{\vrule
height=5.8611pt,depth=-4.6889pt}}}{{\hbox{$\scriptstyle\sqrt{n^{\prime}\,}$}\lower
0.4pt\hbox{\vrule
height=4.125pt,depth=-3.30002pt}}}{{\hbox{$\scriptscriptstyle\sqrt{n^{\prime}\,}$}\lower
0.4pt\hbox{\vrule height=3.26389pt,depth=-2.61113pt}}}|n\rangle\langle
n^{\prime}|\delta_{n,n^{\prime}-1}\Big{)},$ (57)

which leads to two Lindblad operators with corresponding frequencies
$+\omega_{0}$ and $-\omega_{0}$ as

$\displaystyle A_{1}(+\omega_{0})$
$\displaystyle=\dfrac{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}}}}\sum_{n}\mathchoice{{\hbox{$\displaystyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule height=2.15277pt,depth=-1.72223pt}}}|n-1\rangle\langle
n|=\dfrac{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}}}}a,$ (58) $\displaystyle
A_{1}(-\omega_{0})$
$\displaystyle=\dfrac{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}}}}\sum_{n}\mathchoice{{\hbox{$\displaystyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}|n+1\rangle\langle
n|=\dfrac{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}}}}a^{\dagger}$ (59)

Similarly, for
$A_{2}=i(a^{\dagger}-a)/\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}}}$, we get other Lindblad
operators

$\displaystyle A_{2}(+\omega_{0})$
$\displaystyle=-\dfrac{i}{\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}}}}\sum_{n}\mathchoice{{\hbox{$\displaystyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\textstyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule
height=4.30554pt,depth=-3.44446pt}}}{{\hbox{$\scriptstyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule
height=3.01389pt,depth=-2.41113pt}}}{{\hbox{$\scriptscriptstyle\sqrt{n\,}$}\lower
0.4pt\hbox{\vrule height=2.15277pt,depth=-1.72223pt}}}|n-1\rangle\langle
n|=-\dfrac{i}{\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}}}}a,$ (60) $\displaystyle
A_{2}(-\omega_{0})$
$\displaystyle=\dfrac{i}{\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}}}}\sum_{n}\mathchoice{{\hbox{$\displaystyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{n+1\,}$}\lower
0.4pt\hbox{\vrule height=3.22221pt,depth=-2.57779pt}}}|n+1\rangle\langle
n|=\dfrac{i}{\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}}}}a^{\dagger}.$ (61)

It is clear that for the damped harmonic oscillator we have
$k,k^{\prime}\in\\{1,2\\}$ and $\omega\in\\{\omega_{0},-\omega_{0}\\}$, thus
Eq. (50) becomes

$\displaystyle\dfrac{d\rho_{p}(t)}{dt}=$
$\displaystyle-i\left[H_{p},\rho_{p}(t)\right]+\Gamma_{11}(\omega_{0})\Big{(}\\!A_{1}(\omega_{0})\rho_{p}(t)A_{1}^{\dagger}(\omega_{0})-A_{1}^{\dagger}(\omega_{0})A_{1}(\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{12}(\omega_{0})\Big{(}\\!A_{2}(\omega_{0})\rho_{p}(t)A_{1}^{\dagger}(\omega_{0})-A_{1}^{\dagger}(\omega_{0})A_{2}(\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{21}(\omega_{0})\Big{(}\\!A_{1}(\omega_{0})\rho_{p}(t)A_{2}^{\dagger}(\omega_{0})-A_{2}^{\dagger}(\omega_{0})A_{1}(\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{22}(\omega_{0})\Big{(}\\!A_{2}(\omega_{0})\rho_{p}(t)A_{2}^{\dagger}(\omega_{0})-A_{2}^{\dagger}(\omega_{0})A_{2}(\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{11}(-\omega_{0})\Big{(}\\!A_{1}(-\omega_{0})\rho_{p}(t)A_{1}^{\dagger}(-\omega_{0})-A_{1}^{\dagger}(-\omega_{0})A_{1}(-\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{12}(-\omega_{0})\Big{(}\\!A_{2}(-\omega_{0})\rho_{p}(t)A_{1}^{\dagger}(-\omega_{0})-A_{1}^{\dagger}(-\omega_{0})A_{2}(-\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{21}(-\omega_{0})\Big{(}\\!A_{1}(-\omega_{0})\rho_{p}(t)A_{2}^{\dagger}(-\omega_{0})-A_{2}^{\dagger}(-\omega_{0})A_{1}(-\omega_{0})\rho_{p}(t)\Big{)}$
$\displaystyle+\Gamma_{22}(-\omega_{0})\Big{(}\\!A_{2}(-\omega_{0})\rho_{p}(t)A_{2}^{\dagger}(-\omega_{0})-A_{2}^{\dagger}(-\omega_{0})A_{2}(-\omega_{0})\rho_{p}(t)\Big{)}+\text{h.c.}.$
(62)

By substituting Lindblad operators and considering
$\Gamma_{kk^{\prime}}(\omega)=\gamma_{kk^{\prime}}(\omega)+iS_{kk^{\prime}}(\omega)$,
Eq. (62) becomes

$\displaystyle\dfrac{d\rho_{p}(t)}{dt}=-i\left[H_{p},\rho_{p}(t)\right]$
$\displaystyle+\dfrac{1}{2}\big{(}\gamma_{11}(\omega_{0})+iS_{11}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-a^{\dagger}a\rho_{p}(t)\big{)}+\dfrac{1}{2}\big{(}\gamma_{11}(\omega_{0})-iS_{11}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-\rho_{p}(t)a^{\dagger}a\big{)}$
$\displaystyle-\dfrac{i}{2}\big{(}\gamma_{12}(\omega_{0})+iS_{12}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-a^{\dagger}a\rho_{p}(t)\big{)}+\dfrac{i}{2}\big{(}\gamma_{12}(\omega_{0})-iS_{12}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-\rho_{p}(t)a^{\dagger}a\big{]})$
$\displaystyle-\dfrac{i}{2}\big{(}\gamma_{12}(\omega_{0})+iS_{12}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-a^{\dagger}a\rho_{p}(t)\big{)}+\dfrac{i}{2}\big{(}\gamma_{12}(\omega_{0})-iS_{12}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-\rho_{p}(t)a^{\dagger}a\big{)}$
$\displaystyle+\dfrac{1}{2}\big{(}\gamma_{11}(\omega_{0})+iS_{11}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-a^{\dagger}a\rho_{p}(t)\big{)}+\dfrac{1}{2}\big{(}\gamma_{11}(\omega_{0})-iS_{11}(\omega_{0})\big{)}\big{(}a\rho_{p}(t)a^{\dagger}-\rho_{p}(t)a^{\dagger}a\big{)}$
$\displaystyle+\dfrac{1}{2}\big{(}\gamma_{11}(-\omega_{0})+iS_{11}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-aa^{\dagger}\rho_{p}(t)\big{)}+\dfrac{1}{2}\big{(}\gamma_{11}(-\omega_{0})-iS_{11}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-\rho_{p}(t)aa^{\dagger}\big{)}$
$\displaystyle+\dfrac{i}{2}\big{(}\gamma_{12}(-\omega_{0})+iS_{12}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-aa^{\dagger}\rho_{p}(t)\big{)}-\dfrac{i}{2}\big{(}\gamma_{12}(-\omega_{0})-iS_{12}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-\rho_{p}(t)aa^{\dagger}\big{)}$
$\displaystyle+\dfrac{i}{2}\big{(}\gamma_{12}(-\omega_{0})+iS_{12}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-aa^{\dagger}\rho_{p}(t)\big{)}-\dfrac{i}{2}\big{(}\gamma_{12}(-\omega_{0})-iS_{12}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-\rho_{p}(t)aa^{\dagger}\big{)}$
$\displaystyle+\dfrac{1}{2}\big{(}\gamma_{11}(-\omega_{0})+iS_{11}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-aa^{\dagger}\rho_{p}(t)\big{)}+\dfrac{1}{2}\big{(}\gamma_{11}(-\omega_{0})-iS_{11}(-\omega_{0})\big{)}\big{(}a^{\dagger}\rho_{p}(t)a-\rho_{p}(t)aa^{\dagger}\big{)}.$
(63)

The above equation then simplifies to

$\displaystyle\dfrac{d\rho_{p}(t)}{dt}=-i\left[H_{p}\,,\rho_{p}(t)\right]$
$\displaystyle+2\big{(}\gamma_{11}(\omega_{0})+S_{12}(\omega_{0})\big{)}\Big{(}a\rho_{p}(t)a^{\dagger}-\dfrac{1}{2}\big{\\{}a^{\dagger}a,\rho_{p}(t)\big{\\}}\\!\Big{)}-i\big{(}S_{11}(\omega_{0})-\gamma_{12}(\omega_{0})\big{)}\big{[}a^{\dagger}a,\rho_{p}(t)\big{]}$
$\displaystyle+2\big{(}\gamma_{11}(-\omega_{0})-S_{12}(-\omega_{0})\big{)}\Big{(}a^{\dagger}\rho_{p}(t)a-\dfrac{1}{2}\big{\\{}aa^{\dagger},\rho_{p}(t)\big{\\}}\\!\Big{)}-i\big{(}S_{11}(-\omega_{0})+\gamma_{12}(-\omega_{0})\big{)}\big{[}aa^{\dagger},\rho_{p}(t)\big{]}.$
(64)

The coefficients in the above equation are easily calculated to be

$\displaystyle\gamma_{11}(\omega_{0})+S_{12}(\omega_{0})=\frac{1}{2}J(\omega_{0})\big{(}N(\omega_{0},T)+1\big{)},$
(65) $\displaystyle S_{11}(\omega_{0})-\gamma_{12}(\omega_{0})=\Delta_{1},$
(66)
$\displaystyle\gamma_{11}(-\omega_{0})-S_{12}(-\omega_{0})=-\frac{1}{2}J(\omega_{0})\big{(}N(-\omega_{0},T)+1\big{)}=\frac{1}{2}J(\omega_{0})N(\omega_{0},T),$
(67) $\displaystyle S_{11}(-\omega_{0})+\gamma_{12}(-\omega_{0})=-\Delta_{2},$
(68)

where we have used $N(-\omega_{0},T)=-N(\omega_{0},T)-1$. Finally by using
$[a,a^{\dagger}]=\mathbb{I}$, the quantum master equation for the damped
harmonic oscillator in the Schrödinger picture is obtained

$\displaystyle\dfrac{d\rho_{p}(t)}{dt}=$
$\displaystyle-i\left[(\Delta_{1}-\Delta_{2}+\omega_{0})a^{\dagger}a,\rho_{p}(t)\right]+J(\omega_{0})\big{(}N(\omega_{0},T)+1\big{)}\Big{(}a\rho_{p}(t)a^{\dagger}-\dfrac{1}{2}\big{\\{}a^{\dagger}a,\rho_{p}(t)\big{\\}}\\!\Big{)}$
$\displaystyle+J(\omega_{0})N(\omega_{0},T)\Big{(}a^{\dagger}\rho_{p}(t)a-\dfrac{1}{2}\big{\\{}aa^{\dagger},\rho_{p}(t)\big{\\}}\\!\Big{)},$
(69)

which has the same form as in Eq. (46).

###  A.2 The master equation for the moments

To write down the master equation for the first and the second moments we need
to bring the master equation in the standard form Nicacio  _et al._ (2016). We
can define new jump operators to absorb the decay rates, i.e.,

$\displaystyle L^{+}$
$\displaystyle=\mathchoice{{\hbox{$\displaystyle\sqrt{\Gamma_{\rm
in}\,}$}\lower 0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{\Gamma_{\rm
in}\,}$}\lower 0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{\Gamma_{\rm
in}\,}$}\lower 0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\Gamma_{\rm
in}\,}$}\lower 0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}a=\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{\Gamma_{\rm
in}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=11.92552pt,depth=-9.54047pt}}}{{\hbox{$\textstyle\sqrt{\frac{\Gamma_{\rm
in}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=8.3746pt,depth=-6.69972pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{\Gamma_{\rm
in}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.36348pt,depth=-5.0908pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{\Gamma_{\rm
in}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.36348pt,depth=-5.0908pt}}}(1,~{}i)R\eqqcolon(c^{+})^{T}R,$ (70)
$\displaystyle L^{-}$
$\displaystyle=\mathchoice{{\hbox{$\displaystyle\sqrt{\Gamma_{\rm
out}\,}$}\lower 0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\textstyle\sqrt{\Gamma_{\rm
out}\,}$}\lower 0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.46667pt}}}{{\hbox{$\scriptstyle\sqrt{\Gamma_{\rm
out}\,}$}\lower 0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.82668pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\Gamma_{\rm
out}\,}$}\lower 0.4pt\hbox{\vrule
height=3.41666pt,depth=-2.73334pt}}}a^{\dagger}=\mathchoice{{\hbox{$\displaystyle\sqrt{\frac{\Gamma_{\rm
out}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=11.77776pt,depth=-9.42224pt}}}{{\hbox{$\textstyle\sqrt{\frac{\Gamma_{\rm
out}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=8.26904pt,depth=-6.61526pt}}}{{\hbox{$\scriptstyle\sqrt{\frac{\Gamma_{\rm
out}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.25792pt,depth=-5.00636pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\frac{\Gamma_{\rm
out}}{2}\,}$}\lower 0.4pt\hbox{\vrule
height=6.25792pt,depth=-5.00636pt}}}(1,~{}-i)R\eqqcolon(c^{-})^{T}R,$ (71)

where the decay rates are defined as $\Gamma_{\rm in}\coloneqq
J(\omega_{0})N(\omega_{0},T)$, and $\Gamma_{\rm out}\coloneqq
J(\omega_{0})[N(\omega_{0},T)+1]$, and $R=(x,~{}p)^{T}$ is the vector of the
canonical operators. Then, the Lindbladian master equation (25) can be
equivalently written for the displacement vector and covariance matrix:

$\displaystyle\frac{d}{dt}d$ $\displaystyle=A~{}d,$
$\displaystyle\frac{d}{dt}\sigma$
$\displaystyle=A~{}\sigma~{}+~{}\sigma~{}A^{T}~{}+~{}D,$ (72)

with $A=\Omega~{}[G-{\rm Im}(CC^{\dagger})]$, $D=-2\Omega~{}{\rm
Re}(CC^{\dagger})~{}\Omega$, and $C=\left(c^{+},c^{-}\right)$. Substituting
for $C$, we find that

$\displaystyle{\rm Re}~{}CC^{\dagger}$
$\displaystyle=\frac{J(\omega_{0})}{2}\big{(}2N(\omega_{0},T)+1\big{)}\\!\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)=\frac{J(\omega_{0})}{2}\coth\\!\Big{(}\frac{\omega_{0}}{2T}\Big{)}I_{2}$
(75) $\displaystyle{\rm Im}~{}CC^{\dagger}$
$\displaystyle=\frac{J(\omega_{0})}{2}\left(\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right).$ (78)

Putting everything together, we have

$\displaystyle
D=J(\omega_{0})\coth\\!\Big{(}\frac{\omega_{0}}{2T}\Big{)}I_{2},$ (79)

and

$\displaystyle
A=\left(\begin{array}[]{cc}-\frac{J(\omega_{0})}{2}&\omega_{0}\\\
-\omega_{0}&-\frac{J(\omega_{0})}{2}\end{array}\right).$ (82)

The time dependent solution to the “Lyapunov differential equation” (A.2) is
given by

$\displaystyle d(t)$ $\displaystyle=e^{tA}d,$ $\displaystyle\sigma(t)$
$\displaystyle=e^{tA}\sigma(0)e^{tA^{T}}+\int_{0}^{t}dse^{(t-s)A}De^{(t-s)A^{T}}.$
(83)

Note that since the real part of the eigenvalues of $A$ are always negative,
we will have a steady state. One can find the exponential of $A$, e.g., by
using its Jordan form. It reads

$\displaystyle{\rm e}^{tA}$ $\displaystyle={\rm
e}^{-J(\omega_{0})t}~{}\left(\begin{array}[]{cc}\cos(\omega_{0}~{}t)&\sin(\omega_{0}~{}t)\\\
-~{}\sin(\omega_{0}~{}t)&\cos(\omega_{0}~{}t)\end{array}\right)\eqqcolon
e^{-J(\omega_{0})t/2}O_{t},$ (86)

and thus

$\displaystyle d(t)$ $\displaystyle=e^{-J(\omega_{0})t/2}O_{t}d,$
$\displaystyle\sigma(t)$ $\displaystyle=e^{-J(\omega_{0})t}O_{t}\sigma
O^{T}_{t}+\big{(}1-e^{-J(\omega_{0})t}\big{)}\nu I_{2}.$ (87)

with $\nu=\coth{({\omega_{0}}/{2T})}$. Again, note that the rotation $O_{t}$
is temperature independent and has no thermometric value; thus we work in the
interaction picture. The steady state moments can be easily found to be the
thermal ones, i.e., $d(\infty)=0$, and
$\sigma(\infty)=\sigma_{T}=\coth{({\omega_{0}}/{2T})}I_{2}$. One sees that the
speed of converging to the thermal state is given by $J(\omega_{0})$. That is,
the spectral density at the frequency of the oscillator plays the role of the
decay rate $\gamma$ in Eq. (4).

##  Appendix B More simulations for the QFI rate

As shown in Fig. (3), the nonclassicality of the initial probe state, in the
form of squeezing, improves the QFI rate for $T=\omega$. However, this
improvement is not restricted to this specific temperate. In Fig. 5, we
observe a similar behavior at different values of $T/\omega$. Moreover, as we
see in Fig. 6, the two-mode squeezed vacuum states can also improve the
precision at different temperatures. In particular, at lower temperatures the
improvement is even more significant. Nonetheless, our ansatz Gaussian
measurement does not fully exploit this improvement.

Figure 5: The QFI rate for various temperatures, (a) $T=0.2\omega$
corresponding to the low temperature regime, and (b) $T=5\omega$ corresponding
to the high temperature regime. Similar to the main text, we observe an
improvement in using squeezed probe states. Although, for smaller $T/\omega$
one observes the enhancement at smaller times. Here, we set the rest of the
parameters to $\omega=1$, $\gamma=0.1$, and $r=10^{-3}$.

Figure 6: The QFI rate for various temperatures; comparison between two-mode
and single-mode squeezed probe staets. (a) $T=0.2\omega$ corresponding to the
low temperature regieme, and (b) $T=5\omega$ corresponding to the high
temperature regieme. Similar to the main text, we observe an enhancement in
using entangled two-mode squeezed vacuum states. For smaller values of
$T/\omega$, the enhancement can be even larger. Here, we set the rest of the
parameters to $\omega=1$, $\gamma=0.1$, and $r=10^{-3}$.

##  Appendix C The optimal classical Fisher information for single-mode
squeezed probe states

###  C.1 The covariance matrix of the optimal Gaussian measurement is always
diagonal

As discussed in the main text, by working in the interaction picture
$O_{t}=I_{2}$, the covariance matrix of the probe state, Eq. (4), becomes

$\displaystyle\sigma_{t}=e^{-\gamma t}\sigma_{r}+(1-e^{-\gamma t})\nu I_{2}.$
(88)

Here, $\sigma_{r}={\rm diag}(r,1/r)$, as proven in the main text, is the
optimal form of the initial preparation of the probe state. Thus, the
covariance matrix $\sigma_{t}={\rm
diag}\left([\sigma_{t}]_{11},[\sigma_{t}]_{22}\right)$ remains diagonal. Note
that, since we have symmetry in $x$\- and $p$-quadrature, without loss of
generality, we can assume that $0<r\leqslant 1$, which implies
$[\sigma_{t}]_{11}\leqslant[\sigma_{t}]_{22}$.

To obtain the optimal single-mode Gaussian measurements, we only need to
consider rank-1 measurements whose covariance matrix, in general, is described
by

$\displaystyle\sigma^{M}=\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{array}\right)\\!\left(\begin{array}[]{cc}L&0\\\
0&L^{-1}\end{array}\right)\\!\left(\begin{array}[]{cc}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{array}\right)=\left(\begin{array}[]{cc}a&c\\\
c&b\end{array}\right)$ (97)

where $0<L\leqslant 1$ and $0\leqslant\theta<\pi$. Thus, two parameters
$\theta$ and $r$ characterize all rank-1 Gaussian measurements. In particular,
homodyne measurement on $x$ ($p$) quadrature is characterized by $L\to 0$ and
$\theta=0$ ($L\to 0$ and $\theta=\pi/2$). heterodyne, the measurement in
coherent state basis, is identified with $L=1$.

The classical Fisher information (CFI) for a single-mode Gaussian probe state
and Gaussian measurement, by using Eqs. (88) and (97), is given by

$\displaystyle{\cal F}^{\rm C}(\sigma_{r};\sigma^{M};t)$
$\displaystyle=\partial_{T}d_{t}^{T}\sigma_{t}^{-1}\partial_{T}d_{t}+\frac{1}{2}{\rm
Tr}\\!\left[\left(({\sigma}_{t}+{\sigma}^{M})^{-1}\partial_{T}{\sigma}_{t}\right)^{2}\right]$
(98) $\displaystyle=\frac{1}{2}(1-e^{-\gamma
t})^{2}({\partial_{T}\nu})^{2}\,{\rm
Tr}\\!\left[\left(\begin{array}[]{cc}[\sigma_{t}]_{11}+a&c\\\
c&[\sigma_{t}]_{22}+b\end{array}\right)^{\\!\\!-2}\right].$ (101)

where we used $\partial_{T}d_{t}=0$. By diagonalizing the argument,

$\displaystyle\left(\begin{array}[]{cc}[\sigma_{t}]_{11}+a&c\\\
c&[\sigma_{t}]_{22}+b\end{array}\right)^{-1}=U\left(\begin{array}[]{cc}1/\lambda_{+}&0\\\
0&1/\lambda_{-}\end{array}\right)U^{\dagger},$ (106)

and using the eigenvalues

$\displaystyle\lambda_{\pm}=\frac{1}{2}\bigg{(}[\sigma_{t}]_{11}+[\sigma_{t}]_{22}+a+b\pm\mathchoice{{\hbox{$\displaystyle\sqrt{4c^{2}+\left([\sigma_{t}]_{11}-[\sigma_{t}]_{22}+a-b\right)^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.30444pt,depth=-7.44359pt}}}{{\hbox{$\textstyle\sqrt{4c^{2}+\left([\sigma_{t}]_{11}-[\sigma_{t}]_{22}+a-b\right)^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.30444pt,depth=-7.44359pt}}}{{\hbox{$\scriptstyle\sqrt{4c^{2}+\left([\sigma_{t}]_{11}-[\sigma_{t}]_{22}+a-b\right)^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.53888pt,depth=-5.23112pt}}}{{\hbox{$\scriptscriptstyle\sqrt{4c^{2}+\left([\sigma_{t}]_{11}-[\sigma_{t}]_{22}+a-b\right)^{2}\,}$}\lower
0.4pt\hbox{\vrule height=5.03888pt,depth=-4.03113pt}}}\bigg{)}$ (107)

and Eq. (97), we get the general form of the CFI as follows

$\displaystyle{\cal F}^{\rm C}(\sigma_{r};\sigma^{M};t)$
$\displaystyle=\frac{1}{2}(1-e^{-\gamma
t})^{2}({\partial_{T}\nu})^{2}(\lambda_{+}^{-2}+\lambda_{-}^{-2})$ (108)
$\displaystyle=2(1-e^{-\gamma t})^{2}({\partial_{T}\nu})^{2}$
$\displaystyle\times\frac{1+L\left({\rm
Tr}[\sigma_{t}]+L(L^{2}+[\sigma_{t}]_{22}^{2}+[\sigma_{t}]_{11}^{2}+L{\rm
Tr}[\sigma_{t}])\right)-L(L^{2}-1)\\!\left([\sigma_{t}]_{22}-[\sigma_{t}]_{11}\right)\cos(2\theta)}{\left({\rm
Tr}[\sigma_{t}]+L\left(2+2[\sigma_{t}]_{22}[\sigma_{t}]_{11}+L{\rm
Tr}[\sigma_{t}]\right)+\left(L^{2}-1\right)\\!\left([\sigma_{t}]_{22}-[\sigma_{t}]_{11}\right)\cos(2\theta)\right)^{2}}.$
(109)

There are two cases. First, for $r=1$, when the initial probe state is vacuum
or coherent state, by using Eq. (88) we have
$[\sigma_{t}]_{22}=[\sigma_{t}]_{11}$. Hence, the CFI (109) is independent of
$\theta$, and it can be set to zero to have a diagonal $\sigma^{M}$ for the
Gaussian measurement; in this case, Eq. (109) becomes

${\cal F}^{\rm C}(\sigma_{1};\sigma^{M};t)=\frac{1}{2}(1-e^{-\gamma
t})^{2}({\partial_{T}\nu})^{2}\bigg{(}\\!\frac{1}{([\sigma_{t}]_{11}+L)^{2}}+\frac{1}{([\sigma_{t}]_{11}+L^{-1})^{2}}\\!\bigg{)}.$
(110)

Second, for $0<r<1$, corresponding to squeezed-vacuum probe states, we have
$[\sigma_{t}]_{11}\leqslant[\sigma_{t}]_{22}$. In this case for any
$0<L\leqslant 1$, the CFI (109) is maximized when
$(L^{2}-1)([\sigma_{t}]_{22}-[\sigma_{t}]_{11})\cos(2\theta)<0$ and
$|(L^{2}-1)([\sigma_{t}]_{22}-[\sigma_{t}]_{11})\cos(2\theta)|$ is maximum.
Hence, the optimal value is $\theta=0$, and again the covariance matrix
$\sigma^{M}$ is diagonal; in this case, Eq. (109) becomes

${\cal F}^{\rm C}(\sigma_{r};\sigma^{M};t)=\frac{1}{2}(1-e^{-\gamma
t})^{2}({\partial_{T}\nu})^{2}\bigg{(}\\!\frac{1}{([\sigma_{t}]_{11}+L)^{2}}+\frac{1}{([\sigma_{t}]_{22}+L^{-1})^{2}}\\!\bigg{)}.$
(111)

Therefore, in both cases, if the covariance matrix of the initial state is
diagonal, the optimal rank-1 Gaussian measurement has a diagonal covariance
matrix as well.

###  C.2 The optimal Gaussian measurement at short times is homodyne detection

We start by using the CFI—Eq. (111)—for the state given by Eq. (88), and a
Gaussian measurement with diagonal covariance matrix $\sigma^{M}_{L}={\rm
diag}(L,1/L)$. The CFI reads

$\displaystyle{\cal F}_{\rm exact}^{\rm C}(\sigma_{r};\sigma^{M}_{L};t)$
$\displaystyle\coloneqq{\cal F}^{\rm
C}(\sigma_{r};\sigma^{M}_{L};t)=\frac{1}{2}\big{(}1-e^{-\gamma
t}\big{)}^{2}(\partial_{T}\nu)^{2}\\!\left(\\!\frac{{\mathrm{e}}^{2\,\gamma\,t}}{{{\left(r-\nu+L\,{\mathrm{e}}^{\gamma\,t}+\nu\,{\mathrm{e}}^{\gamma\,t}\right)}}^{2}}+\frac{L^{2}\,r^{2}\,{\mathrm{e}}^{2\,\gamma\,t}}{{{\left(L+r\,{\mathrm{e}}^{\gamma\,t}-L\,\nu\,r+L\,\nu\,r\,{\mathrm{e}}^{\gamma\,t}\right)}}^{2}}\\!\right)$
$\displaystyle=\frac{1}{2}(1-e^{-\gamma
t})^{2}(\partial_{T}\nu)^{2}\left(\\!\frac{{\mathrm{e}}^{2\,\gamma\,t}\,{{\left(L+r\,{\mathrm{e}}^{\gamma\,t}-L\,\nu\,r+L\,\nu\,r\,{\mathrm{e}}^{\gamma\,t}\right)}}^{2}+L^{2}\,r^{2}\,{\mathrm{e}}^{2\,\gamma\,t}\,{{\left(r-\nu+L\,{\mathrm{e}}^{\gamma\,t}+\nu\,{\mathrm{e}}^{\gamma\,t}\right)}}^{2}}{{{\left(r-\nu+L\,{\mathrm{e}}^{\gamma\,t}+\nu\,{\mathrm{e}}^{\gamma\,t}\right)}}^{2}\,{{\left(L+r\,{\mathrm{e}}^{\gamma\,t}-L\,\nu\,r+L\,\nu\,r\,{\mathrm{e}}^{\gamma\,t}\right)}}^{2}}\right).$
(112)

We can now expand both the numerator and the denominator for short times, and
keep the leading order terms in $\gamma t$. Note that in doing so we are not
allowed to ignore $\gamma t$ compared to $L$ or $r$, since they can be equally
small. However, we can ignore terms like $r\gamma t$ compared to $r$, or
ignore $L\gamma t$ compared to $L$. We have

$\displaystyle{\cal F}^{\rm C}(\sigma_{r};\sigma^{M}_{L};t\ll 1)$
$\displaystyle\approx\frac{(\gamma
t\partial_{T}\nu)^{2}}{2}\frac{(L+r+r(1+L\nu)\gamma
t)^{2}+L^{2}r^{2}(r+L+(L+\nu)\gamma t)^{2}}{(r+L+(L+\nu)\gamma
t)^{2}(L+r+r(1+L\nu)\gamma t)^{2}}$ $\displaystyle\approx\frac{(\gamma
t\partial_{T}\nu)^{2}}{2}\frac{(L+r)^{2}+L^{2}r^{2}(r+L+\nu\gamma
t)^{2}}{(r+L+\nu\gamma t)^{2}(L+r)^{2}}.$ (113)

To proceed further, we consider two specific cases. (i) If ${\cal
O}(L+r)\lesssim{\cal O}(\nu\gamma t)$ then the second term in the numerator is
irrelevant compared to the first term and therefore we have

$\displaystyle{\cal F}^{\rm C}(\sigma_{r};\sigma^{M}_{L};t\ll 1)$
$\displaystyle\approx\frac{(\gamma t\partial_{T}\nu)^{2}}{2(r+L+\nu\gamma
t)^{2}},$ (114)

which clearly reaches its maximum if we chose $L=0$. For this choice, the
Fisher information is maximised at $t^{*}\propto r/(\gamma\nu)$. (ii) If
${\cal O}(L+r)>{\cal O}(\nu\gamma t)$, then we can ignore the $\nu\gamma t$
terms in the parenthesis to get

$\displaystyle{\cal F}^{\rm C}(\sigma_{r};\sigma^{M}_{L};t\ll 1)$
$\displaystyle\approx\frac{(\gamma
t\partial_{T}\nu)^{2}(1+L^{2}r^{2})}{2(r+L)^{2}}.$ (115)

By taking derivative of this expression with respect to $L$, we find no
solutions, i.e., the function is monotonic. A simple comparison then shows
that ${\cal F}^{\rm C}(\sigma_{r};\sigma^{M}_{0};t\ll 1)>{\cal F}^{\rm
C}(\sigma_{r};\sigma^{M}_{1};t\ll 1)$, that is, homodyne detection is the best
Gaussian measurement in this case too. For this measurement, the Fisher
information is zero initially and grows quadratically with time—until ${\cal
O}(r)\lesssim{\cal O}(\nu\gamma t)$ when we have to use the expression (114).
Finally, note that when $L=0$, one can use equation (114) which covers both
cases

$\displaystyle{\cal F}^{\rm C}_{\rm approx}(\sigma_{r};\sigma^{M}_{0};t\ll
1)\coloneqq\frac{(\gamma t\partial_{T}\nu)^{2}}{2(r+\nu\gamma t)^{2}}.$ (116)

As depicted in Fig. 7, this approximation is in a very good agreement with the
exact value of the CFI given by Eq. (C.2).

Figure 7: The comparison between the exact and the approximation of the CFI,
respectively, Eqs. (C.2) and (116), shows a good agreement. In the left panel
we set $r=10^{-3}$ corresponding to a highly squeezed initial preparation,
while in the right panel, we set $r=1$ corresponding to vacuum. The rest of
the parameters are set to $T=1$, $\gamma=0.2$, $r=10^{-3}$, and $\omega=1$.

##  Appendix D Two-mode squeezed states and a feasible joint Gaussian
measurement

Figure 8: The schematic of a fully Gaussian setup that uses an auxiliary mode
in an entangled state with probe mode and a joint measurement for further
improving the thermometry precision. The dashed blue box on the left describes
the preparation procedure: A two-mode squeezed vacuum state is prepared by
injecting two orthogonally squeezed-vacuum states into a 50:50 beamsplitter.
The probe (upper) mode interacts with the bath and then both modes are jointly
measured. The joint measurement (the dashed red box) consists of a 50:50
beamsplitter and two homodyne measurements, one on the $x$-quadrature and the
other on the $p$-quadrature.

We describe a fully Gaussian setup that can exploit entanglement for enhanced
precision in temperature estimation—see Fig. 8. The joint system of probe and
an auxiliary bosonic mode is prepared in a two-mode squeezed vacuum state,
which is an entangled state. Such state can be generated by a spontaneous
parametric down-conversion source or equivalently, as shown in Fig. 8, by
overlapping two single-mode squeezed vacuum states, which are orthogonally
squeezed, on a 50:50 beamsplitter. By using the symplectic transformation
describing the 50:50 beamsplitter which reads

$S_{\rm BS}=\frac{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}}}}\left(\begin{array}[]{cccc}1&0&1&0\\\
0&1&0&1\\\ -1&0&1&0\\\ 0&-1&0&1\end{array}\right),$ (117)

the covariance matrix of two-mode squeezed vacuum states can be obtained as

$\displaystyle\sigma_{2,r}=\frac{1}{2}S_{\rm
BS}\left(\begin{array}[]{cccc}u&0&0&0\\\ 0&u^{-1}&0&0\\\ 0&0&u^{-1}&0\\\
0&0&0&u\end{array}\right)S_{\rm
BS}^{T}=\left(\begin{array}[]{cc}r^{-1}I_{2}&\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}}Z\\\
\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower 0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}&r^{-1}I_{2}\end{array}\right),$ (124)

where $Z={\rm diag}(1,-1)$ and we defined $r\coloneqq 2(u+1/u)^{-1}$. After
preparing the probe and the auxiliary systems in an entangled state, the probe
mode interacts with the bath and, as discussed in the main text, this
transforms the state into a two-mode Gaussian state that is described by the
covariance matrix

$\displaystyle\sigma_{t}=\left(\begin{array}[]{cccc}[(1-e^{-\gamma
t})\nu+e^{-\gamma t}r^{-1}]I_{2}&e^{-\gamma
t/2}\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}}Z\\\ e^{-\gamma
t/2}\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=3.44165pt,depth=-2.75334pt}}}Z&r^{-1}I_{2},\end{array}\right).$ (127)

Now, to exploit the power of entanglement, a joint measurement is required.
Our joint measurement consists of a 50:50 beamsplitter followed by two
homodyne measurements, as shown in Fig. 8. However, to obtain the Fisher
information, we apply the beamsplitter on the state, which gives the
covariance matrix $S_{\rm BS}\sigma_{t}S_{\rm BS}=\sigma^{\prime}_{t}$ with
the following elements,

$\displaystyle[\sigma_{t}^{\prime}]_{11}$
$\displaystyle=[\sigma_{t}^{\prime}]_{44}=\frac{1-e^{-\gamma
t}}{2}\nu+\frac{1+e^{-\gamma t}}{2}r^{-1}-e^{-\gamma
t/2}\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}},$ (128)
$\displaystyle[\sigma_{t}^{\prime}]_{13}$
$\displaystyle=[\sigma_{t}^{\prime}]_{24}=e^{-\gamma
t/2}\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}},$
$\displaystyle[\sigma_{t}^{\prime}]_{22}$
$\displaystyle=[\sigma_{t}^{\prime}]_{33}\frac{1-e^{-\gamma
t}}{2}\nu+\frac{1+e^{-\gamma t}}{2}r^{-1}+e^{-\gamma
t/2}\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}},$
$\displaystyle[\sigma_{t}^{\prime}]_{31}$
$\displaystyle=[\sigma_{t}^{\prime}]_{42}=\frac{1-e^{-\gamma
t}}{2}(\nu-r^{-1})\mathchoice{{\hbox{$\displaystyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\textstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.15558pt}}}{{\hbox{$\scriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.6089pt}}}{{\hbox{$\scriptscriptstyle\sqrt{r^{-2}-1\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.75334pt}}},$
$\displaystyle[\sigma_{t}^{\prime}]_{jk}$ $\displaystyle=0,\text{ otherwise}.$
(129)

Then, using two homodyne measurements, we measure the $x$-quadrature of the
first output mode and the $p$-quadrature of the second output mode. Note that
by performing this joint measurement (50:50 beamsplitter followed by two
homodyne measurements), we effectively sample from a two-dimensional Gaussian
probability distribution with the covariance matrix $\tilde{\sigma}={\rm
diag}([\sigma_{t}^{\prime}]_{11},[\sigma_{t}^{\prime}]_{44})=[\sigma_{t}^{\prime}]_{11}I_{2}$.
Therefore, the CFI is given by

$\displaystyle{\cal F}^{\rm C}(\sigma_{2,r};\sigma^{M}_{2,0};t)$
$\displaystyle=\frac{1}{2}{\rm
Tr}\\!\left[\left(\tilde{\sigma}^{-1}\partial_{T}\tilde{\sigma}\right)^{2}\right]=\frac{(1-e^{-\gamma
t})^{2}({\partial_{T}\nu})^{2}}{4[\sigma_{t}^{\prime}]_{11}^{2}},$ (130)

which is Eq. (13) of the main text.

## References

  * Mehboudi  _et al._ (2019) Mohammad Mehboudi, Anna Sanpera, and Luis A Correa, “Thermometry in the quantum regime: recent theoretical progress,” Journal of Physics A: Mathematical and Theoretical 52, 303001 (2019).
  * De Pasquale and Stace (2018) Antonella De Pasquale and Thomas M Stace, “Quantum thermometry,” in  _Thermodynamics in the Quantum Regime_ (Springer, 2018) pp. 503–527.
  * Brunelli  _et al._ (2011) Matteo Brunelli, Stefano Olivares, and Matteo G. A. Paris, “Qubit thermometry for micromechanical resonators,” Phys. Rev. A 84, 032105 (2011).
  * Brunelli  _et al._ (2012) Matteo Brunelli, Stefano Olivares, Mauro Paternostro, and Matteo G. A. Paris, “Qubit-assisted thermometry of a quantum harmonic oscillator,” Phys. Rev. A 86, 012125 (2012).
  * Correa  _et al._ (2015) Luis A. Correa, Mohammad Mehboudi, Gerardo Adesso, and Anna Sanpera, “Individual quantum probes for optimal thermometry,” Phys. Rev. Lett. 114, 220405 (2015).
  * Jevtic  _et al._ (2015) Sania Jevtic, David Newman, Terry Rudolph, and T. M. Stace, “Single-qubit thermometry,” Phys. Rev. A 91, 012331 (2015).
  * Razavian  _et al._ (2019) Sholeh Razavian, Claudia Benedetti, Matteo Bina, Yahya Akbari-Kourbolagh, and Matteo GA Paris, “Quantum thermometry by single-qubit dephasing,” The European Physical Journal Plus 134, 284 (2019).
  * De Pasquale  _et al._ (2017) Antonella De Pasquale, Kazuya Yuasa, and Vittorio Giovannetti, “Estimating temperature via sequential measurements,” Phys. Rev. A 96, 012316 (2017).
  * Mitchison  _et al._ (2022) Mark T. Mitchison, Archak Purkayastha, Marlon Brenes, Alessandro Silva, and John Goold, “Taking the temperature of a pure quantum state,” Phys. Rev. A 105, L030201 (2022).
  * Cavina  _et al._ (2018) Vasco Cavina, Luca Mancino, Antonella De Pasquale, Ilaria Gianani, Marco Sbroscia, Robert I. Booth, Emanuele Roccia, Roberto Raimondi, Vittorio Giovannetti, and Marco Barbieri, “Bridging thermodynamics and metrology in nonequilibrium quantum thermometry,” Phys. Rev. A 98, 050101 (2018).
  * Guo  _et al._ (2015) Li-Sha Guo, Bao-Ming Xu, Jian Zou, and Bin Shao, “Improved thermometry of low-temperature quantum systems by a ring-structure probe,” Phys. Rev. A 92, 052112 (2015).
  * Mitchison  _et al._ (2020) Mark T. Mitchison, Thomás Fogarty, Giacomo Guarnieri, Steve Campbell, Thomas Busch, and John Goold, “In situ thermometry of a cold fermi gas via dephasing impurities,” Phys. Rev. Lett. 125, 080402 (2020).
  * Seah  _et al._ (2019) Stella Seah, Stefan Nimmrichter, Daniel Grimmer, Jader P. Santos, Valerio Scarani, and Gabriel T. Landi, “Collisional quantum thermometry,” Phys. Rev. Lett. 123, 180602 (2019).
  * Shu  _et al._ (2020) Angeline Shu, Stella Seah, and Valerio Scarani, “Surpassing the thermal cramér-rao bound with collisional thermometry,” Phys. Rev. A 102, 042417 (2020).
  * Mancino  _et al._ (2020) Luca Mancino, Marco G. Genoni, Marco Barbieri, and Mauro Paternostro, “Nonequilibrium readiness and precision of gaussian quantum thermometers,” Phys. Rev. Research 2, 033498 (2020).
  * Henao and Uzdin (2021) Ivan Henao and Raam Uzdin, “Catalytic transformations with finite-size environments: applications to cooling and thermometry,” Quantum 5, 547 (2021).
  * Hovhannisyan  _et al._ (2021) Karen V. Hovhannisyan, Mathias R. Jorgensen, Gabriel T. Landi, Álvaro M. Alhambra, Jonatan B. Brask, and Martí Perarnau-Llobet, “Optimal quantum thermometry with coarse-grained measurements,” PRX Quantum 2, 020322 (2021).
  * Sekatski and Perarnau-Llobet (2021) Pavel Sekatski and Martí Perarnau-Llobet, “Optimal nonequilibrium thermometry in finite time,” arXiv preprint arXiv:2107.04425 (2021), https://doi.org/10.48550/arXiv.2107.04425.
  * Weedbrook  _et al._ (2012) Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J Cerf, Timothy C Ralph, Jeffrey H Shapiro, and Seth Lloyd, “Gaussian quantum information,” Reviews of Modern Physics 84, 621 (2012).
  * Cenni  _et al._ (2022) Marina F.B. Cenni, Ludovico Lami, Antonio Acín, and Mohammad Mehboudi, “Thermometry of Gaussian quantum systems using Gaussian measurements,” Quantum 6, 743 (2022).
  * Schaller (2014) Gernot Schaller,  _Open quantum systems far from equilibrium_, Vol. 881 (Springer, 2014).
  * Cramér (1999) Harald Cramér,  _Mathematical methods of statistics_ , Vol. 43 (Princeton university press, 1999).
  * Newey and McFadden (1994) Whitney K. Newey and Daniel McFadden, “Chapter 36 large sample estimation and hypothesis testing,”  Handbook of Econometrics, 4, 2111–2245 (1994).
  * Braunstein  _et al._ (1996) Samuel L. Braunstein, Carlton M. Caves, and G.J. Milburn, “Generalized uncertainty relations: Theory, examples, and lorentz invariance,” Annals of Physics 247, 135–173 (1996).
  * Braunstein and Caves (1994) Samuel L. Braunstein and Carlton M. Caves, “Statistical distance and the geometry of quantum states,” Phys. Rev. Lett. 72, 3439–3443 (1994).
  * Paris (2009) Matteo GA Paris, “Quantum estimation for quantum technology,” International Journal of Quantum Information 7, 125–137 (2009).
  * Giovannetti  _et al._ (2011) Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone, “Advances in quantum metrology,” Nature photonics 5, 222–229 (2011).
  * Hayashi (2006) Masahito Hayashi,  _Quantum information_ (Springer, 2006).
  * Chin  _et al._ (2012) Alex W. Chin, Susana F. Huelga, and Martin B. Plenio, “Quantum metrology in non-markovian environments,” Phys. Rev. Lett. 109, 233601 (2012).
  * del Campo  _et al._ (2013) A. del Campo, I. L. Egusquiza, M. B. Plenio, and S. F. Huelga, “Quantum speed limits in open system dynamics,” Phys. Rev. Lett. 110, 050403 (2013).
  * Demoen  _et al._ (1977) Bart Demoen, P Vanheuverzwijn, and André Verbeure, “Completely positive maps on ccr-algebra,” Letters in mathematical physics 2, 161–166 (1977).
  * Lindblad (2000) Göran Lindblad, “Cloning the quantum oscillator,” Journal of Physics A: Mathematical and General 33, 5059–5076 (2000).
  * Cerf  _et al._ (2007) N J Cerf, G Leuchs, and E S Polzik,  _Quantum Information with Continuous Variables of Atoms and Light_ (PUBLISHED BY IMPERIAL COLLEGE PRESS AND DISTRIBUTED BY WORLD SCIENTIFIC PUBLISHING CO., 2007).
  * Holevo and Werner (2001) A. S. Holevo and R. F. Werner, “Evaluating capacities of bosonic gaussian channels,” Phys. Rev. A 63, 032312 (2001).
  * Monras (2013) Alex Monras, “Phase space formalism for quantum estimation of gaussian states,” arXiv:1303.3682 (2013).
  * Pinel  _et al._ (2013) O. Pinel, P. Jian, N. Treps, C. Fabre, and D. Braun, “Quantum parameter estimation using general single-mode gaussian states,” Phys. Rev. A 88, 040102 (2013).
  * Jiang (2014) Zhang Jiang, “Quantum fisher information for states in exponential form,” Phys. Rev. A 89, 032128 (2014).
  * Nichols  _et al._ (2018) Rosanna Nichols, Pietro Liuzzo-Scorpo, Paul A. Knott, and Gerardo Adesso, “Multiparameter gaussian quantum metrology,” Phys. Rev. A 98, 012114 (2018).
  * Malagò and Pistone (2015) Luigi Malagò and Giovanni Pistone, “Information geometry of the gaussian distribution in view of stochastic optimization,” Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII , 150–162 (2015).
  * Titulaer and Glauber (1965) U. M. Titulaer and R. J. Glauber, “Correlation functions for coherent fields,” Phys. Rev. 140, B676–B682 (1965).
  * Mandel (1986) L Mandel, “Non-classical states of the electromagnetic field,” Physica Scripta T12, 34–42 (1986).
  * Rahimi-Keshari  _et al._ (2015) Saleh Rahimi-Keshari, Timothy C Ralph, and Carlton M Caves, “Operational discord measure for gaussian states with gaussian measurements,” New Journal of Physics 17, 063037 (2015).
  * Rivas and Huelga (2011) Á. Rivas and S.F. Huelga,  _Open Quantum Systems: An Introduction_ , SpringerBriefs in Physics (Springer Berlin Heidelberg, 2011).
  * Breuer and Petruccione (2002) Heinz-Peter Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).
  * Lidar (2020) Daniel A. Lidar, “Lecture notes on the theory of open quantum systems,” (2020), arXiv:1902.00967 [quant-ph] .
  * Albash  _et al._ (2012) Tameem Albash, Sergio Boixo, Daniel A Lidar, and Paolo Zanardi, “Quantum adiabatic markovian master equations,” New Journal of Physics 14, 123016 (2012).
  * Nafari Qaleh and Rezakhani (2022) Z. Nafari Qaleh and A. T. Rezakhani, “Enhancing energy transfer in quantum systems via periodic driving: Floquet master equations,” Phys. Rev. A 105, 012208 (2022).
  * Correa and Glatthard (2023) Luis A Correa and Jonas Glatthard, “Potential renormalisation, lamb shift and mean-force gibbs state–to shift or not to shift?” arXiv preprint arXiv:2305.08941 (2023), https://doi.org/10.48550/arXiv.2305.08941.
  * Nicacio  _et al._ (2016) F. Nicacio, M. Paternostro, and A. Ferraro, “Determining stationary-state quantum properties directly from system-environment interactions,” Phys. Rev. A 94, 052129 (2016).