# Bending the rules of low-temperature thermometry with periodic driving

Jonas Glatthard  Department of Physics and Astronomy, University of Exeter,
Exeter EX4 4QL, United Kingdom<EMAIL_ADDRESS>Luis A. Correa
Department of Physics and Astronomy, University of Exeter, Exeter EX4 4QL,
United Kingdom

###### Abstract

There exist severe limitations on the accuracy of low-temperature thermometry,
which poses a major challenge for future quantum-technological applications.
Low-temperature sensitivity might be manipulated by tailoring the interactions
between probe and sample. Unfortunately, the tunability of these interactions
is usually very restricted. Here, we focus on a more practical solution to
boost thermometric precision—driving the probe. Specifically, we solve for the
limit cycle of a periodically modulated linear probe in an equilibrium sample.
We treat the probe–sample interactions exactly and hence, our results are
valid for arbitrarily low temperatures $T$ and any spectral density. We find
that weak near-resonant modulation strongly enhances the signal-to-noise ratio
of low-temperature measurements, while causing minimal back action on the
sample. Furthermore, we show that near-resonant driving changes the power law
that governs thermal sensitivity over a broad range of temperatures, thus
‘bending’ the fundamental precision limits and enabling more sensitive low-
temperature thermometry. We then focus on a concrete example—impurity
thermometry in an atomic condensate. We demonstrate that periodic driving
allows for a sensitivity improvement of several orders of magnitude in sub-
nanokelvin temperature estimates drawn from the density profile of the
impurity atoms. We thus provide a feasible upgrade that can be easily
integrated into low-$T$ thermometry experiments.

##  1 Introduction

Recent experimental progress is making it possible to cool systems down to
unprecedented low temperatures. For instance, cold atomic gases have been
pushed into the sub-nanokelvin regime, and even down to few picokelvins [1, 2,
3, 4]. Yet measuring such record-breaking temperatures accurately remains
remarkably difficult [3]. This is not only due to technical issues, nor
specific to cold atomic gases. There are fundamental physical limitations that
render low-temperature thermometry inefficient [5]. Furthermore, at ultracold
temperatures, infinitesimal-coupling treatments that rely on local
thermalisation of probe systems become inadequate, since probe and sample
build quantum correlations which can push their marginals far from the Gibbs
state [6, 7, 8].

These shortcomings have motivated the stand-alone theoretical study of
thermometric precision in the quantum regime, termed quantum thermometry (see
[9, 10] and references therein). In particular, results on the scaling laws
that govern thermal sensitivity as the temperature $T\rightarrow 0$ [5, 11,
12] paint a grim picture. The signal-to-noise ratio of estimates on the
temperature of gapped samples is exponentially suppressed as $T$ decreases
[5], while in gapless systems it decays as a power law in $T$ instead111A
temperature-independent signal-to-noise ratio is possible in principle [5, 13,
14], but it requires engineering the probe–sample Hamiltonian, working with
limiting-case spectral densities, or using adaptive strategies. [5, 11, 12].
It is therefore essential to devise protocols that can enhance low-temperature
sensitivity as much as possible.

Engineering the probe–sample interactions [7, 15] can greatly improve
thermometric precision, but the degree of control over the spectral properties
of the sample is arguably limited. Alternatively, Mukherjee et al. [13]
proposed addressing the probe with a periodic driving field. Specifically, it
was shown that periodic driving allows to open new dissipation channels at
tuneable frequencies and to tailor the weight of their contribution to the
overall dissipative dynamics, ultimately yielding improved thermometric
bounds. Periodic driving has also proven useful in other metrological settings
[16, 17].

Technically, the analysis in Ref. [13] combines weak-coupling
Gorini–Kossakowski–Lindblad–Sudarshan (GKLS) quantum master equations [18, 19,
20] with Floquet theory [21, 22]. Imposing the GKLS form in a master equation
guarantees formal thermodynamic consistency and, for vanishingly weak
dissipation, they may yield an excellent approximation to the steady state of
a probe coupled to an ultracold sample. However, microscopically derived
Born–Markov secular GKLS master equations are known to break down at any
finite dissipation strength, as one moves to the $T\rightarrow 0$ limit [7,
15]. Namely, GKLS equations invariably predict local thermalisation of any
probe at the sample temperature in the absence of driving. At ultracold
temperatures, however, these do build strong correlations, which results in
non-thermal marginals. Furthermore, using GKLS equations on open systems with
tailored spectral densities may be questioned from a microscopic-derivation
viewpoint whenever the resulting correlation functions in the sample are long
lived [23].

Here we study the impact of periodic driving on low-temperature thermometric
precision, circumventing all these issues. We exploit the exact techniques
from [24, 25, 26] and solve for the non-equilibrium state of a periodically
driven Brownian particle in a harmonic potential and dissipatively coupled to
a linear equilibrium sample. Importantly, our results do not rely on any
assumptions on the strength of the dissipation, the spectral structure of the
sample, or its temperature. We provide explicit formulas for the limit-cycle
state of the probe, exact up to any order in the strength of the drive, and
use them to calculate the theoretical upper bound on the signal-to-noise ratio
of temperature estimates drawn from arbitrary measurements on the probe [10].
Ultimately, we analyse the low-$T$ scaling of the precision of such estimates.

Firstly, we find that our calculations qualitatively agree with the findings
in Ref. [13]; that is, the largest precision enhancements occur when the
external drive is weak and near-resonant. For our model, thermometric
precision without driving is known to obey a power law, with an exponent
dictated by the low-frequency behaviour (Ohmicity) of the spectral density of
the sample [5]. Interestingly, we find that near-resonant periodic driving
effectively modifies this scaling over a broad range of temperatures,
resulting in a large improvement of low-$T$ thermometry. We also demonstrate
that, in order to benefit from this precision boost, it is not necessary to
synchronise the measurements with the driving field—the advantage survives
time-averaging of the oscillating asymptotic probe state.

In order to illustrate the practical potential of this technique on a concrete
platform, we considered impurity thermometry in a Bose–Einstein condensate. We
exploit the fact that tightly confined atomic impurities immersed in a cold
atomic gas with a large condensed fraction can be modelled with our Brownian-
motion model [27, 28, 15]. Instead of evaluating the precision limits, we
adopt a more applied perspective. Working with the asymptotic time-averaged
state of the probe, we evaluate the responsiveness of the density profile of
the impurity cloud to small temperature fluctuations. For typical experimental
parameters, we show that near-resonant driving leads to a sensitivity
enhancement of several orders of magnitude at temperatures
$T\lesssim$1\text{\,}\mathrm{nK}$$.

Finally, we take into account a side-effect of addressing the probe with an
external field—the heating of the sample. Specifically, we calculate the
coarse grained long-time probe-sample heat flow. We find that, even when
driving near resonance, the back-action on the sample is a fourth-order effect
in the strength of the external drive.

This paper is structured as follows: The model is introduced in Sec. 2.1 and
an outline of the main steps in the calculation of its non-equilibrium limit-
cycle state is given in Sec. 2.2 (the full calculation is in Appendix A).
After briefly introducing the necessary estimation-theory tools in Sec. 3, we
present our main results in Sec. 4.1, and apply them to impurity thermometry
in Sec. 4.2. Sample heating is discussed in Sec. 4.3, with technical details
deferred to Appendices B and C. Finally, in Sec. 5 we summarise and draw our
conclusions.

##  2 The model and its limit cycle

###  2.1 Hamiltonian

Let the Hamiltonian of our probe be

$\boldsymbol{H}_{S}(t)=\frac{1}{2}\omega^{2}(t)\,\boldsymbol{x}^{2}+\frac{1}{2}\,\boldsymbol{p}^{2}.$
(1)

In all what follows, boldface symbols denote operators and we shall work in
units of $m=\hbar=k_{B}=1$ unless stated otherwise. We choose a time-dependent
frequency of the form

$\omega(t)^{2}=\omega_{0}^{2}+\upsilon\sin{\omega_{d}t}.$ (2)

The role of the sample will be played by a linear heat bath with Hamiltonian

$\boldsymbol{H}_{B}=\sum\nolimits_{\mu}\frac{1}{2}\omega_{\mu}^{2}m_{\mu}\boldsymbol{x}_{\mu}^{2}+\frac{1}{2m_{\mu}}\boldsymbol{p}_{\mu}^{2},$
(3)

while the dissipative coupling between probe and sample is

$\boldsymbol{H}_{I}=\boldsymbol{x}\,\sum\nolimits_{\mu}g_{\mu}\boldsymbol{x}_{\mu}.$
(4)

The interaction constants $g_{\mu}$ shape the spectral density

$J(\omega)\coloneqq\pi\sum_{\mu}\frac{g_{\mu}^{2}}{2m_{\mu}\omega_{\mu}}\delta(\omega-\omega_{\mu}),$

which we take to be Ohmic with algebraic cutoff

$J(\omega)=\frac{2\gamma\,\omega}{1+(\omega/\omega_{c})^{2}}.$ (5)

Here $\omega_{c}$ stands for the cutoff frequency of the bath. We refer such
spectrum as Ohmic since $J(\omega)\sim\omega^{s}$ at low frequency with $s=1$.
In contrast, $s>1$ and $s<1$ would correspond to super- and sub-Ohmic spectral
densities, respectively.

The exact Heisenberg equations of motion for the system can be written in
compact form as (cf. Appendix A.1)

$\ddot{\boldsymbol{x}}(t)+\big{(}\omega^{2}(t)+\omega_{R}^{2}\big{)}\,\boldsymbol{x}(t)\\\
-\int_{0}^{t}\chi(t-\tau)\,\boldsymbol{x}(\tau)d\tau=\boldsymbol{F}(t),$ (6)

where we assume that the bath’s marginal is initially a state of thermal
equilibrium. At this point, the initial state of the probe can be fully
general, though we will later restrict ourselves to Gaussian initial states.
The dissipation kernel $\chi(t)$ in Eq. (6) is

$\chi(t)=\frac{2}{\pi}\Theta(t)\int_{0}^{\infty}J(\omega)\sin{(\omega
t)}\,d\omega,$ (7)

and the stochastic force $\boldsymbol{F}(t)$ is defined in Appendix A.1.
$\Theta(t)$ stands for the Heaviside step function. The frequency shift
$\omega_{R}$ in Eq. (6) is

$\omega_{R}^{2}\coloneqq\frac{2}{\pi}\int_{0}^{\infty}\frac{J(\omega)}{\omega}d\omega.$
(8)

This is introduced in the Hamiltonian to compensate for the distortion on the
probe’s potential due to its coupling to the bath [29]. For our $J(\omega)$ in
Eq. (5), it evaluates to $\omega_{R}^{2}=2\,\gamma\,\omega_{c}$.

###  2.2 Limit cycle

Let $g(t,t^{\prime})$ be the Green’s function for Eq. (6), i.e.,

$\partial_{t}^{2}g(t,t^{\prime})+\big{(}\omega^{2}(t)+\omega_{R}^{2}\big{)}\,g(t,t^{\prime})\\\
-\int_{0}^{t}\chi(t-\tau)\,g(\tau,t^{\prime})\,d\tau=\delta(t-t^{\prime}).$
(9)

If a stable limit cycle exists (cf. Appendix C) $g(t,t^{\prime})$ can be cast
as

$g(t,t^{\prime})=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\int_{-\infty}^{\infty}a_{k}(\omega)\mathrm{e}^{i\omega(t-t^{\prime})}\mathrm{e}^{ik\omega_{d}t}d\omega.$
(10)

The detailed derivation of (10) is deferred to Appendix A.2 (see also Ref.
[26]). As we shall now see, knowledge of the amplitudes $a_{k}(\omega)$ allows
us to fully characterise the limit cycle of the probe.

To proceed further, let us insert (10) into (9), which leads to the following
relationship between the amplitudes $a_{k}(\omega)$

$\hat{g}_{0}[i(\omega+k\omega_{d})]^{-1}a_{k}(\omega)+\\\ \sum\nolimits_{l\neq
0}\,b_{l}\,a_{k-l}(\omega)=\delta_{k0},$ (11)

where $\delta_{ij}$ is a Kronecker delta and the periodic driving has been
Fourier expanded as

$\omega(t)^{2}=\omega_{0}^{2}+\sum\nolimits_{l\neq
0}b_{l}\,\mathrm{e}^{i\,l\,\omega_{d}\,t}.$ (12)

The notation $\hat{g}_{0}(s)$ in Eq. (11) stands for

$\hat{g}_{0}(s)\coloneqq\big{[}s^{2}+\omega_{0}^{2}+\omega_{R}^{2}-\hat{\chi}(s)\big{]}^{-1},$
(13)

where $\hat{\chi}(s)$ is the Laplace transform of the dissipation kernel,
i.e., $\hat{\chi}(s)=\int_{0}^{\infty}dt\,\mathrm{e}^{-st}\,\chi(t)$. We shall
implicitly work with $\hat{f}(i\omega)=\hat{f}(i\omega+0^{+})$. The choice of
notation in (13) is due to the fact that $\hat{g}_{0}(s)$ is the Laplace
transform of the Green’s function in Eq. (9) in the undriven case.

For our choice of spectral density we have

$\hat{\chi}(i\omega)=\frac{2\gamma\omega_{c}^{2}}{\omega_{c}+i\omega}.$ (14)

In turn, since our driving in Eq. (2) is sinusoidal the only non-vanishing
coefficients in the expansion (12) are $b_{\pm 1}=\mp\,i\,\upsilon/2$.

One can solve Eq. (11) for $a_{k}(\omega)$ self-consistently as

$a_{k}^{(n)}(\omega)=\hat{g}_{0}[i(\omega+k\omega_{d})]\\\
\times\big{(}\delta_{k0}-\sum\nolimits_{l\neq
0}b_{l}\,a_{k-l}^{(n-1)}(\omega)\big{)},$ (15)

where $a_{k}^{(n)}(\omega)$ is the solution after $n$ iterations, starting
from

$a_{k}^{(0)}(\omega)=\hat{g}_{0}[i(\omega+k\omega_{d})]\,\delta_{k0}.$

Each iteration of Eq. (15) adds only terms of higher order in $\upsilon$, and
$a_{k}=\lim_{n\rightarrow\infty}a_{k}^{(n)}$.

Substituting (12) and (15) into (10) yields the Green’s function
$g(t,t^{\prime})$ of the equation of motion of the system (6) in the limit
cycle. In fact, we can explicitly write the (time-dependent) asymptotic state
of the system in terms of the Fourier coefficients $a_{k}(\omega)$ and
$b_{k}$. To that end, we exploit the fact that the Hamiltonian is quadratic in
positions and momenta, thus generating a Gaussianity-preserving dynamics [30].
That is, given a Gaussian initial condition with vanishing first moments,
knowledge of the covariances $\sigma_{xx}=\langle\boldsymbol{x}^{2}\rangle$,
$\sigma_{xp}=\frac{1}{2}\langle\boldsymbol{x}\boldsymbol{p}+\boldsymbol{p}\boldsymbol{x}\rangle=\frac{1}{2}\langle\\{\boldsymbol{x},\boldsymbol{p}\\}\rangle$
and $\sigma_{pp}=\langle\boldsymbol{p}^{2}\rangle$ is enough to fully
characterise the state of the probe.

Specifically, the limit-cycle covariances of the probe in contact with an
equilibrium bath are

$\displaystyle\sigma_{\alpha\alpha}^{\text{lc}}(t)$
$\displaystyle=\mathrm{Re}\,\sum\nolimits_{j\,k}\varsigma_{jk}^{\alpha\alpha}\,\mathrm{e}^{i\,\omega_{d}(j-k)\,t},$
(16a) $\displaystyle\sigma_{xp}^{\text{lc}}(t)$
$\displaystyle=\mathrm{Im}\,\sum\nolimits_{j\,k}\varsigma_{jk}^{xp}\,\mathrm{e}^{i\,\omega_{d}(j-k)\,t},$
(16b)

where $\alpha=\\{x,p\\}$ and the super-index ‘lc’ indicates ‘limit cycle’. The
coefficients $\varsigma^{\alpha\beta}_{jk}$ are

$\displaystyle\varsigma_{jk}^{xx}$
$\displaystyle=\frac{1}{2}\int_{0}^{\infty}a_{j}(\omega)\,\hat{\mu}(\omega)\,a_{k}^{*}(\omega)\,d\omega,$
(17a) $\displaystyle\varsigma_{jk}^{pp}$
$\displaystyle=\frac{1}{2}\int_{0}^{\infty}(\omega+j\,\omega_{d})\,(\omega+k\,\omega_{d})$
$\displaystyle\qquad\qquad\qquad\times
a_{j}(\omega)\,\hat{\mu}(\omega)\,a_{k}^{*}(\omega)\,d\omega,$ (17b)
$\displaystyle\varsigma_{jk}^{xp}$
$\displaystyle=\frac{1}{2}\int_{0}^{\infty}(\omega+k\,\omega_{d})\,a_{j}(\omega)\hat{\mu}(\omega)\,a_{k}^{*}(\omega)\,d\omega.$
(17c)

Above, we have introduced the notation

$\hat{\mu}(\omega)\coloneqq\frac{2}{\pi}\,J(\omega)\,\coth{\left(\omega/2T\right)},$
(18)

and the asterisk denotes complex conjugation. The full derivation can be found
in Appendix A.2.

Note that, for our choice of driving, $a_{k}$ is $\pazocal{O}(\upsilon^{k})$.
Since $a_{k}(\omega)=a_{k}^{(2)}(\omega)+\pazocal{O}(\upsilon^{3})$, we may
stop after the second iteration in Eq. (15) for sufficiently small $\upsilon$
(i.e., weak periodic modulation). Explicitly, this gives us

$a_{0}(\omega)\simeq\hat{g}_{0}(i\omega)\\\ +\sum_{l\neq
0}\hat{g}_{0}(i\omega)\,b_{l}\,\hat{g}_{0}[i(\omega-l\omega_{d})]\,b_{-l}\,\hat{g}_{0}(i\omega).$
(19a) $a_{k\neq
0}(\omega)\simeq-\hat{g}_{0}[i(\omega+k\omega_{d})]\,b_{k}\,\hat{g}_{0}(i\omega)\\\
+\sum\nolimits_{l\neq\\{0,k\\}}\hat{g}_{0}[i(\omega+k\omega_{d})]b_{l}\\\
\times\hat{g}_{0}[i(\omega+(k-l)\omega_{d})]\,b_{k-l}\,\hat{g}_{0}(i\omega),$
(19b)

which facilitates calculations. For our simple driving, the only non-zero
amplitudes at second order are thus $a_{0}^{(2)}(\omega),a_{\pm
1}^{(2)}(\omega)$ and $a_{\pm 2}^{(2)}(\omega)$. We emphasise, however, that
Eqs. (16) and (17) are non-perturbative in the strength $\upsilon$ of the
external driving, so that aribitrarily high orders can be considered.

##  3 Ultimate thermometric precision

The Cramér–Rao inequality [31, 32] places a lower bound the statistical
uncertainty on any parameter inferred from a large number of measurements $N$
by means of an unbiased estimator. This bound has been extensively used in
high-precision sensing and, in particular, in low-temperature estimation (see
[10, 9] and references therein). Namely, the uncertainty in a temperature
estimate scales as

$\delta
T_{\mathrm{est}}\geq\frac{1}{\sqrt{N\pazocal{F}(\boldsymbol{\varrho}_{T})}},$
(20)

where $\pazocal{F}(\boldsymbol{\varrho})$ is the ‘quantum Fisher information’
calculated on the temperature-dependent state of the probe
$\boldsymbol{\varrho}_{T}$ as

$\pazocal{F}(\boldsymbol{\varrho}_{T})=-2\lim_{\tau\rightarrow
0}\frac{\partial^{2}\mathbb{F}(\boldsymbol{\varrho}_{T},\boldsymbol{\varrho}_{T+\tau})}{\partial\tau^{2}},$
(21)

and the Uhlmann fidelity
$\mathbb{F}(\boldsymbol{\varrho}_{1},\boldsymbol{\varrho}_{2})$ is defined as

$\mathbb{F}(\boldsymbol{\varrho}_{1},\boldsymbol{\varrho}_{2})=\mathrm{tr}\left(\sqrt{\sqrt{\boldsymbol{\varrho}_{1}}\,\boldsymbol{\varrho}_{2}\,\sqrt{\boldsymbol{\varrho}_{1}}}\right)^{2}.$

While the Cramér–Rao bound is inadequate when working with scarce data [33],
it has proven useful to understand the fundamental scaling laws of low-
temperature thermometry [7, 5, 11, 14].

The quantum Fisher information in Eq. (21) is particularly simple to evaluate
on (single-mode) Gaussian states. In that case, the fidelity between states
$\boldsymbol{\varrho}_{1}$ and $\boldsymbol{\varrho}_{2}$ with covariance
matrices $\mathsf{\Sigma}_{1}$ and $\mathsf{\Sigma}_{2}$, respectively, can be
cast as [34]

$\mathbb{F}(\mathsf{\Sigma}_{1},\mathsf{\Sigma}_{2})=2\left(\sqrt{\varkappa+1}-\sqrt{\lambda}\right)^{-1},$
(22)

where $\varkappa=4\det{\left(\Sigma_{1}+\Sigma_{2}\right)}$ and
$\lambda=(4\det\Sigma_{1}-1)(4\det\Sigma_{2}-1)$. In the formula above we have
taken into account that the first-order moments vanish in our case.

Figure 1: (a) Instantaneous (inverse) thermal sensitivity in the limit cycle.
Driving near resonance with $\omega_{d}=0.9\omega_{0}$ (dashed orange) and
$\omega_{d}=0.995\omega_{0}$ (dot-dashed green) improves the thermal
sensitivity of the probe. For comparison, the undriven case (solid blue) shows
the expected $1/T^{2}$ scaling as $T\rightarrow 0$ [5], while for a probe in a
Gibbs state with respect to its local Hamiltonian (dotted red), the relative
error grows exponentially. The limit cycle was sampled at 10 equispaced times
which span the shaded regions around the curves. (b) Lower bound on the
(inverse) thermal sensitivity from the time-averaged Fisher information
$\overline{\pazocal{F}(\boldsymbol{\varrho}_{T}(t))}$. In both plots
$\omega_{0}=1$, $\upsilon=0.1$, $\gamma=0.01$, and $\omega_{c}=100$.

##  4 Results and discussion

###  4.1 Enhancement of low-$T$ sensitivity

We are now in a position to assess the impact of periodic driving on the
thermal sensitivity of the probe at arbitrarily low temperatures. To do so, we
simply apply Eqs. (20–22) on the limit-cycle covariances from Eqs. (16) and
(17). It must be noted that, in contrast with master-equation methods, the
parameter ranges that we may study is not limited by time-scale separation
assumptions. That is, Eqs. (16) and (17) yield the exact stationary state,
only assuming Gaussianity for the initial condition, and that the sample is
initially in thermal equilibrium. For our calculations below we do approximate
the amplitudes $a_{k}(\omega)$ to second and fourth order in $\upsilon$, which
effectively limits our analysis to weak driving222More precisely, our
covariances $\sigma_{\alpha\beta}^{\text{lc}}(t)$ are accurate to second (or
fourth) order in $\upsilon$. Care must be taken, however, since the covariance
undergoes further manipulations in Eq. (22), which generates higher-order
terms which are not explicitly removed. Nonetheless, one can verify that the
thermal sensitivity calculated from the second-order covariances does agree
with the fourth-order calculations.—but not to weak dissipation or short-lived
bath correlations. Furthermore, as we argue below, weak driving is a practical
choice to minimise the heating of the sample (cf. B and C).

We take the signal-to-noise ratio $\frac{1}{\sqrt{N}}(T/\delta
T_{\text{est}})$ as the figure of merit, which we refer-to as thermal
sensitivity. The lower bound on the inverse of this quantity—as given by the
Cramér–Rao inequality—is plotted in Fig. 1(a) versus the temperature $T$, for
different modulations $\omega_{d}$. The calculation is repeated considering
various times within the limit cycle, since the asymptotic state of the probe
still depends (periodically) on time.

We observe that weak dissipation strength and near-resonant modulation
($\omega_{d}\approx\omega_{0}$) does boost sensitivity by several orders of
magnitude in the low-temperature regime. More interestingly, the power law
obeyed by the low-temperature sensitivity changes as a result of the driving.
It is known that the Fisher information on temperature of a Brownian particle
scales as $\pazocal{F}(\boldsymbol{\varrho}_{T})\sim T^{2\,s}$ at
$T\rightarrow 0$ [5] where $s$ is the Ohmicity parameter introduced in Sec.
2.1; in our case, $s=1$. In contrast, we find that for the power law is
effectively modified to $\pazocal{F}(\boldsymbol{\varrho}_{T})\sim T$ as
$\omega_{d}\rightarrow\omega_{0}$. Hence, near-resonant periodic driving
renders low-temperature thermometry much more efficient. This is our main
result.

By considering larger dissipation strengths we have also observed that
maintaining a substantial advantage over the undriven case requires stronger
driving $\upsilon$. This is due to the fact that strong dissipation is capable
of largely improving the low-temperature sensitivity [7] on its own, which
masks any improvements due to the drive.

As already noted, the thermal sensitivity of the driven probe oscillates in
time, spanning the shaded regions around the plotted curves in Fig. 1(a). In
order to get a clearer picture of the expected improvement in sensitivity, we
may prefer to work with the time-averaged limit-cycle

$\overline{\boldsymbol{\varrho}}_{T}\coloneqq\frac{\omega_{d}}{2\pi}\int_{0}^{2\pi/\omega_{d}}\boldsymbol{\varrho}_{T}(t)\,dt.$
(23)

Such average generally yields a non-Gaussian state lying on the convex hull of
Gaussian states. As a result, we can no longer use Eq. (22) to evaluate the
thermal sensitivity. Instead, we exploit of the extended convexity property of
the quantum Fisher information [35], which becomes

$\pazocal{F}_{T}\left[\int_{0}^{2\pi/\omega_{d}}\boldsymbol{\varrho}_{T}(t)\,dt\right]\leq\int_{0}^{2\pi/\omega_{d}}\pazocal{F}_{T}[\boldsymbol{\varrho}_{T}(t)]\,dt.$
(24)

Hence, by averaging the quantum Fisher information of the probe over the limit
cycle, we can bound the instantaneous thermal sensitivity from below. This is
illustrated in Fig. 1(b), which suggests that the scaling improvement for
near-resonant driving does survive time averaging. Indeed, in Sec. 4.2 we
give, in a concrete example, a lower bound to the left hand side of Ineq. (24)
which also exhibits an enhanced low-temperature scaling from the drive. This
has an important practical consequence: there would be no need to synchronise
the measurements on the probe with the drive in order to benefit from the
precision enhancement.

###  4.2 Application: Impurity thermometry in a Bose–Einstein condensate

Figure 2: Sensitivity improvement of the position quadrature of an impurity in
a cold atomic gas due to periodic modulation of the trapping potential. The
limit-cycle-averaged position dispersion of Yb impurities immersed in a BEC of
K is used to estimated the temperature of the latter. The trapping frequency
of for Yb is, on average, $\omega_{I}=2\pi\times$375\text{\,}\mathrm{Hz}$$,
and is sinusoidally modulated with amplitude $\upsilon=0.2\omega_{I}$ and
frequency $\omega_{d}=0.8\omega_{I}$. The resulting sensitivity (dashed
orange) can be greatly enhanced with respect to the unmodulated case (solid
blue). It is worth noting that in the unmodulated case, position dispersion is
a quasi-optimal temperature indicator [7]. The plotted temperature range of
$0.5\text{\,}\mathrm{nK}$–$2\text{\,}\mathrm{nK}$ is within experimental reach
[1, 3, 36]. The other parameters are $N_{B}=5000$,
$\omega_{B}=2\pi\times$750\text{\,}\mathrm{Hz}$$, $g_{IB}=0.55\times
10^{39}$\text{\,}\mathrm{J}\text{\,}\mathrm{m}$$ and $g_{B}=3\times
10^{39}$\text{\,}\mathrm{J}\text{\,}\mathrm{m}$$.

Thermometry in ultracold gases becomes particularly challenging as
temperatures decrease. The conventional time-of-flight absorption imaging
method can no longer be used on bosonic gases with a large condensed fraction,
and alternatives must be sought [3]. Specifically, the temperature of a
Bose–Einstein condensate (BEC) may be experimentally assessed using atomic
impurities as probes [37, 3, 38]. In the appropriate limits, a dilute gas of
impurity atoms co-trapped alongside a BEC can be described with our simple
quantum-Brownian-motion model [27, 28]. We may thus investigate whether
periodic modulation of the confining potential of the impurities can enhance
their thermal sensitivity at ultra-low temperatures. As we will see, it does.

Specifically, we consider an atomic impurity of mass $m_{I}$ in a 1D harmonic
trap with frequency $\omega_{I}$. In turn, the co-trapped BEC is comprised of
$N_{B}$ atoms of mass $m_{B}$ at temperature $T$, and is trapped in a harmonic
potential with frequency $\omega_{B}$. One may show that, if the temperature
is sufficiently low and the impurity is tightly confined around the centre of
the condensate trap, its dynamics is effectively described by the model in
Sec. 2.1, with spectral density [28]

$J(\omega)=2\,\gamma_{0}\,\omega^{4}\,\Theta(\omega_{B}-\omega),$ (25)

where the constant $\gamma_{0}$ is given by

$\gamma_{0}=\frac{\pi\,g_{B}}{\omega_{B}^{4}\,r^{3}}\left(\frac{g_{IB}\,\mu}{g_{B}\,\hbar\,\omega_{B}}\right)^{2},$

and $r$ and $\mu$ stand for the the Thomas–Fermi radius and the chemical
potential, respectively

$\displaystyle r$ $\displaystyle=\sqrt{2\mu/(m_{B}\omega_{B}^{2})},$
$\displaystyle\mu$
$\displaystyle=\left(\frac{3}{4\sqrt{2}}\,g_{B}\,N_{B}\,\omega_{B}\,\sqrt{m_{B}},\right)^{3/2}.$

Here, $g_{IB}$ is the interspecies (i.e., probe–sample) interactions, and
$g_{B}$ stands for the interatomic interactions within the sample. $J(\omega)$
also incorporates a hard cutoff through the Heaviside step function
$\Theta(\omega_{B}-\omega)$. Notice that, within this section, we work in SI
units.

Being able to tackle this problem without quantum master equations is
particularly important, since we are dealing with a super-Ohmic spectral
density. Environments with such power spectra have longer-lived correlation
functions than their Ohmic counterparts, which could easily compromise the
Born–Markov approximation underpinning most master-equation analyses.

Since the resulting spectral density differs from our Eq. (5), we must
recalculate the Laplace transform $\hat{\chi}$ of the corresponding
dissipation kernel. This reads [28]

$\hat{\chi}(i\omega)=\frac{\gamma_{0}\,\omega_{B}}{\pi}\\\
-\frac{2\,\gamma_{0}\,\omega^{2}}{\pi\,\omega_{B}^{3}}\left(\omega_{B}^{2}+\omega^{2}\,\log{\frac{\omega^{2}}{\omega^{2}+\omega_{B}^{2}}}\right).$
(26)

The only other difference with Sec. 2.1 is that the frequency shift
$\omega_{R}^{2}$ may be entirely neglected in the regime of validity of the
microscopic derivation of Eq. (25) [15]. Otherwise, Eqs. (16) and (17) are
directly applicable to this problem.

Rather than assuming a precise time control in the interrogation of the probe,
we work with the time-averaged state from Eq. (23), thus getting a better idea
of the sensitivity gain that may be expected in practice. Also, instead of
searching for ultimate precision bounds, we ask ourselves what is the
responsiveness of the position dispersion $\boldsymbol{x}^{2}$ of the
impurities to temperature fluctuations in the limit cycle. This is indeed an
experimentally accessible temperature estimator [39, 38], and its
responsiveness may quantified as [40, 41, 10]

$F_{T}(\boldsymbol{x}^{2})\coloneqq\frac{\big{|}\partial_{T}\langle\boldsymbol{x}^{2}\rangle_{\overline{\boldsymbol{\varrho}}_{T}}\big{|}^{2}}{(\Delta\boldsymbol{x}^{2})^{2}}\leq\pazocal{F}_{T},$
(27)

where the subindices emphasise that averages are taken over
$\overline{\boldsymbol{\varrho}}_{T}$ and thus,

$(\Delta\boldsymbol{x}^{2})^{2}=\langle\boldsymbol{x}^{4}\rangle_{\overline{\boldsymbol{\varrho}}_{T}}-(\langle\boldsymbol{x}^{2}\rangle_{\overline{\boldsymbol{\varrho}}_{T}})^{2}.$

Due to linearity, we can write
$\langle\boldsymbol{x}^{2}\rangle_{\overline{\boldsymbol{\varrho}}_{T}}$ as
$\overline{\sigma_{xx}^{\text{lc}}(t)}$, where the bar indicates limit-cycle
average. The temperature-dependence in Eqs. (17) is entirely contained in the
noise kernel $\hat{\mu}(\omega)$, which allows to readily calculate
$\partial_{T}\overline{\sigma_{xx}^{\text{lc}}}$. Furthermore, note that only
those terms with $j=k$ from Eq. (16) survive time averaging. On the other
hand, to evaluate the denominator of Eq. (27) we use the fact that the
instantaneous state of the probe is always Gaussian, so that the identity
$\langle\boldsymbol{x}^{4}(t)\rangle=3\langle\boldsymbol{x}(t)^{2}\rangle^{2}$
holds at all times. Exploiting linearity again, we may write

$\displaystyle\langle\boldsymbol{x}^{4}\rangle_{\overline{\boldsymbol{\varrho}}_{T}}-(\langle\boldsymbol{x}^{2}\rangle_{\overline{\boldsymbol{\varrho}}_{T}})^{2}=3\overline{\sigma_{xx}^{\text{lc}}(t)^{2}}-(\,\overline{\sigma_{xx}^{\text{lc}}(t)}\,)^{2}.$

The temperature-dependence of the responsiveness $F_{T}(\boldsymbol{x}^{2})$
is plotted in Fig. 2 with and without driving, for the cold atomic mixture
studied in [15]. As it can be seen, applying weak near-resonant driving to the
confining potential of the impurities yields sensitivity improvements of
several orders of magnitude at temperatures around
$k_{B}T\sim\hbar(\omega_{0}-\omega_{d})/4$. For typical parameters this
corresponds to the $0.1\text{\,}\mathrm{nK}$–$1\text{\,}\mathrm{nK}$ range,
which is of experimental interest [1, 3, 36], albeit extremely challenging for
conventional thermometry [3]. Closer to resonance, we again observe improved
scaling.

Although we have chosen cold atoms for illustration, the same basic principle
can be applied to other setups, like superconducting microresonators [42].
Switching to typical trapping frequencies in the $\text{\,}\mathrm{GHz}$ range
would result in improved temperature sensing at millikelvins, which is
directly relevant for this platform.

###  4.3 Heat dissipation

Figure 3: (top) Cycle-averaged heating rate of the sample as a function of the
driving frequency $\omega_{d}$. The heat current is calculated to second and
fourth order in $\upsilon$ (solid blue and dashed orange, respectively). As we
can see, the heating around resonance is a fourth-order effect. All parameters
are as in Fig. 1 and $T=0.025$. (bottom) Stability chart of a (classical)
damped Mathieu oscillator for varying driving frequency $\omega_{d}$ and
driving amplitude $\upsilon$, at fixed $\omega_{0}$. Unstable solutions of
exponentially growing amplitude appear in the shaded regions. In absence of
friction ($\gamma=0$; light grey) solutions are unstable at the resonance
frequencies $\omega_{d}=2\omega_{0}/n$ for integer $n$ (vertical dotted
lines), regardless of $\upsilon$. On the contrary, stable solutions may always
be found at sufficiently weak driving under non-zero friction. For instance,
unstable solutions for $\gamma=0.1$ are limited to the dark grey area.

So far we have seen how periodic driving can substantially boost the
sensitivity of low-temperature thermometry. This, however, comes at the cost
of heating the sample as the work pumped into the probe is dissipated away. In
order to assess the magnitude of this effect, we now focus on characterising
the probe–sample energy exchanges in the limit cycle.

The rate of change of the energy of the probe is simply
$\frac{d}{dt}\langle\boldsymbol{H}_{S}(t)\rangle=-i\langle[\boldsymbol{H}_{S}(t),\boldsymbol{H}_{I}]\rangle+\langle\partial_{t}\,\boldsymbol{H}_{S}(t)\rangle$,
where the brackets indicate averaging over the instantaneous state.
Specifically, the first term is the contribution due to the probe–sample
coupling. Let us denote it by

$\dot{Q}(t)=-i\langle[\boldsymbol{H}_{S}(t),\boldsymbol{H}_{I}]\rangle.$ (28)

In the limit cycle, this term may be cast as (see Appendix B)

$\dot{Q}(t)=\frac{1}{2}\frac{d}{dt}\sigma_{pp}(t)+\left(\omega^{2}(t)+\omega_{R}^{2}\right)\sigma_{xp}(t).$
(29)

Energy conservation and stationarity require that the asymptotic cycle average
of $\dot{Q}(t)$, i.e.,

$\overline{\dot{Q}}=\frac{\omega_{d}}{2\pi}\int_{0}^{2\pi/\omega_{d}}\dot{Q}^{\text{lc}}(t)\,dt,$

coincides with minus the average heating rate of the sample [26]. Explicitly,
$\overline{\dot{Q}}$ takes the form

$\overline{\dot{Q}}=-\frac{1}{\pi}\int_{0}^{\infty}\sum\nolimits_{k}\,k\,\omega_{d}\,\tilde{J}(\omega+k\,\omega_{d})\,a_{k}(\omega)\,a^{*}_{k}(\omega)\\\
\times J(\omega)\,\coth{(\omega/2T)}\,d\omega,$ (30)

where $\tilde{J}(\omega)$ is the spectral density extended to negative
frequencies as an odd function (cf. Appendix B).

Fig. 3 shows the limit-cycle-averaged heat current into the sample as a
function of the driving frequency. It has been calculated both to second and
fourth order in $\upsilon$, assuming again the Ohmic spectral density in Eq.
(5).

We see a large peak in the heat current around $\omega_{d}=2\omega_{0}$, and a
smaller one (a fourth-order effect) at $\omega_{d}=\omega_{0}$. Tuning the
drive close to resonance to enhance sensitivity thus comes at the cost of a
larger heating of the sample. The compromise between the two will be
ultimately dictated by the heat capacity of the sample. This is, however,
typically very large, so that fourth-order heating effects would be negligible
at small driving amplitude $\upsilon$. When it comes to the rest of parameters
in the problem, smaller dissipation strength and lower temperature both reduce
the heat current.

The dependence of the current on $\omega_{d}$ suggests a connection with
parametric resonance. Indeed, studying the closely related classical Mathieu
oscillator we find some striking similarities. The damped Mathieu equation
reads

$\ddot{x}+\gamma\,\dot{x}+\omega^{2}(t)\,x=0.$ (31)

In the limit of $\omega_{c}\rightarrow\infty$ the integral in Eq. (6)
generates a dissipation term of the form $\gamma\dot{x}$. The main difference
between the Mathieu equation and our model (besides the former being
classical) is the absence of thermal noise in Eq. (31).

The solutions of the Mathieu equation as time evolves are either decaying,
purely oscillatory, or have increasing amplitude (see Appendix C). The latter
are termed unstable. Fig. 3 shows the regions of instability of the solution
to Eq. (31) in parameter space. In the absence of damping (i.e., $\gamma=0$)
we recover the well-known instabilities at parametric resonance. Specifically,
at $\omega_{d}=2\omega_{0}/n$ for $n\in\\{1,2,...\\}$. Those are precisely the
points at which the heat current peaks in our model, thus suggesting a
connection between heating and instability. We do work, however, at finite
dissipation strength and small $\upsilon$, which would allow to avoid
instabilities even near resonance.

##  5 Conclusions

We have studied how the precision of low-temperature measurements can be
largely improved by periodic driving. More concretely, we have solved for the
long-time limit of a driven quantum Brownian temperature probe in an
equilibrium sample. We have found that the low-temperature sensitivity of the
Brownian probe may be substantially boosted by simply applying a weak
sinusoidal near-resonant modulation to its harmonic confining potential.

In order to quantify thermal sensitivity, we have calculated the upper limit
on the signal-to-noise ratio $\frac{1}{\sqrt{N}}\,(T/\delta T_{\text{est}})$.
For a large number of measurements $N$, this is given by the quantum Fisher
information of the limit-cycle probe state. Interestingly, we have shown that
the boost in precision does survive limit-cycle averaging. In other words,
enabling a precision improvement of several orders of magnitude in a
thermometry experiment would not require synchronisation between the
measurements and the drive.

Furthermore, we have illustrated how the exponent of the power-law-like
asymptotic behaviour of the signal-to-noise ratio at low $T$ [5] may be
manipulated by the driving. Effectively, this renders low-temperature
thermometry much more efficient than otherwise allowed by the spectral density
of the sample. Specifically, we show that the ultimate precision limit can
switch from a $\sim T^{2}$ scaling to a $\sim T$ behaviour as $T\rightarrow 0$
on a sample with Ohmic spectrum.

Finally, we have exploited the fact that an impurity thermometry setup in a
BEC can be microscopically modelled with our linear quantum Brownian motion
setup. We have shown that, for representative experimental parameters, weak
near-resonant driving of the impurity trap does lead to a very substantial
amplification of thermal sensitivity for $T\lesssim$1\text{\,}\mathrm{nK}$$.
In this case, we used a practically motivated proxy for thermal
sensitivity—the responsiveness of the time-averaged density profile of the
impurity gas to temperature fluctuations.

We have also investigated the impact of heat dissipation onto the sample as a
byproduct of the driving. We have found that, when driving near resonance,
heating is a fourth-order effect in the driving strength, thus producing
minimal back-action on the ultracold sample.

Crucially, the explicit formulae that we have derived for the asymptotic state
of the probe are exact to any order in the strength of the external drive.
Since we have avoided using Markovian quantum master equations, we have
bypassed the severely limiting assumptions of fast-decaying bath correlation
functions and weak system–bath coupling. Although the correct intuition can be
drawn from a simple Markovian calculation [13], our results extend to
arbitrary coupling and spectral density, and to arbitrarily low temperatures,
thus permitting the study of the $T\rightarrow 0$ scaling of thermometric
precision.

Ultimately, the physics behind the sensitivity improvement is the closing of
the Floquet quasi-energy gap, which has been exploited in other sensing
applications (see, e.g., [17]). We thus expect our results to qualitatively
hold beyond our particular model. Indeed, in the limit of weak dissipation,
periodic driving can be seen to have the same effect in non-harmonic probes
[13, 22]. Similarly, since the most crucial role for the question at hand is
played by the low-frequency tail of the spectral density, we also expect our
results to extend beyond the Ohmic and super-Ohmic spectral densities
considered here.

We have thus validated a robust solution to increase low-temperature
sensitivity in thermometry experiments, and quantified its impact on
thermometry in ultracold atomic gases. Studying the potential benefits on
other platforms and temperature ranges—e.g., $\text{\,}\mathrm{mK}$-cooled
superconducting circuits—remains an interesting open direction. Likewise, it
would be practically relevant to assess the impact of periodic driving when
working with small measurement records on both the accuracy and the speed of
convergence of a temperature estimate. This may be achieved with the newly
developed global (Bayesian) thermometry tools [33] and will be the subject of
future work.

## Acknowledgements

We thank M. Mehboudi, A. Riera-Campeny, N. Freitas, P. Potts, J. Anders, M.
Perarnau-Llobet and G. Barontini for useful discussions. JG is supported by a
scholarship from the College of Engineering, Mathematics and Physical Sciences
of the University of Exeter.

##  Appendix A Limit cycle of the probe

###  A.1 Equation of motion

The Heisenberg equations of motion from Eqs. (1–4) are

$\displaystyle\dot{\boldsymbol{x}}$ $\displaystyle=\boldsymbol{p},$ (32a)
$\displaystyle\dot{\boldsymbol{p}}$
$\displaystyle=-\big{(}\omega(t)^{2}+\omega_{R}^{2}\big{)}\boldsymbol{x}-\sum\nolimits_{\mu}g_{\mu}\,\boldsymbol{x}_{\mu},$
(32b) $\displaystyle\dot{\boldsymbol{x}}_{\mu}$
$\displaystyle=\boldsymbol{p}_{\mu}/m_{\mu},$ (32c)
$\displaystyle\dot{\boldsymbol{p}}_{\mu}$
$\displaystyle=-m_{\mu}\omega_{\mu}^{2}\,\boldsymbol{x}_{\mu}-g_{\mu}\,\boldsymbol{x}.$
(32d)

To make this paper self-contained, we shall now provide full detail on how to
obtain their limit cycle, following Ref. [26]. To simplify the notation we
define
$\boldsymbol{\mathsf{Z}}_{\mu}=(\boldsymbol{x}_{\mu},\boldsymbol{p}_{\mu})^{\mathsf{T}}$
and $\boldsymbol{\mathsf{Z}}=(\boldsymbol{x},\boldsymbol{p})^{\mathsf{T}}$.
Hence, dynamical equations for the sample degrees of freedom can be compacted
as

$\dot{\boldsymbol{\mathsf{Z}}}_{\mu}(t)+\begin{pmatrix}0&-m_{\mu}^{-1}\\\
m_{\mu}\omega_{\mu}^{2}&0\end{pmatrix}\,\boldsymbol{\mathsf{Z}}_{\mu}(t)=\begin{pmatrix}0&0\\\
-g_{\mu}&0\end{pmatrix}\,\boldsymbol{\mathsf{Z}}(t),$ (33)

where $\boldsymbol{\mathsf{Z}}(t)$ acts as a source term in the equations for
$\boldsymbol{\mathsf{Z}}_{\mu}(t)$. The corresponding Green’s function
$\mathsf{G}_{\mu}(t,t^{\prime})$ thus satisfies

$\partial_{t}\,\mathsf{G}_{\mu}(t,t^{\prime})+\begin{pmatrix}0&-m_{\mu}^{-1}\\\
m_{\mu}\omega_{\mu}^{2}&0\end{pmatrix}\mathsf{G}_{\mu}(t,t^{\prime})\\\
=\delta(t-t^{\prime})\,\mathbbm{1}_{2},$

where $\mathbbm{1}_{2}$ is the $2\times 2$ identity matrix and $\delta(x)$
stands for the Dirac delta. Hence,

$\mathsf{G}_{\mu}(t,t^{\prime})=\Theta(t-t^{\prime})\\\
\times\begin{pmatrix}\cos{(\omega_{\mu}(t-t^{\prime}))}&\frac{\sin{(\omega_{\mu}(t-t^{\prime}))}}{m_{\mu}\,\omega_{\mu}}\\\
-m_{\mu}\omega_{\mu}\sin{(\omega_{\mu}(t-t^{\prime}))}&\cos{(\omega_{\mu}(t-t^{\prime}))}\end{pmatrix}.$

Therefore, Eq. (33) is then solved by

$\boldsymbol{\mathsf{Z}}_{\mu}(t)=\mathsf{G}_{\mu}(t,0)\,\boldsymbol{\mathsf{Z}}_{\mu}(0)+\int_{0}^{t}\mathsf{G}_{\mu}(t,t^{\prime})\,\mathsf{C}_{\mu}\,\boldsymbol{Z}(t^{\prime})\,dt^{\prime},$
(34)

where $\mathsf{C}_{\mu}=\begin{pmatrix}0&0\\\ -g_{\mu}&0\end{pmatrix}$.

Similarly, we can now look at the equations of motion of the degrees of
freedom of the probe

$\dot{\boldsymbol{\mathsf{Z}}}(t)+\begin{pmatrix}0&-1\\\
\omega(t)^{2}+\omega_{R}^{2}&0\end{pmatrix}\boldsymbol{\mathsf{Z}}(t)\\\
-\int_{0}^{t}\begin{pmatrix}0&0\\\
\chi(t-t^{\prime})&0\end{pmatrix}\boldsymbol{\mathsf{Z}}(t^{\prime})dt^{\prime}=\boldsymbol{\mathsf{F}}(t),$
(35)

where the dissipation kernel $\chi(t-t^{\prime})$ was introduced in Eq. (7) in
the main text, and
$\boldsymbol{\mathsf{F}}(t)=\sum_{\mu}\mathsf{C}_{\mu}\,\mathsf{G}_{\mu}(t,0)\,\boldsymbol{\mathsf{Z}}_{\mu}(0)$.
Note that Eq. (35) can be recast as a second-order equation for
$\boldsymbol{x}$ of the form advanced in Eq. (6) of the main text, with
stochastic force $\boldsymbol{F}(t)$ given by

$\boldsymbol{F}(t)=-\sum_{\mu}g_{\mu}\,\\\
\times\left(\boldsymbol{x}_{\mu}(0)\cos{(\omega_{\mu}t)}+\frac{\boldsymbol{p}_{\mu}(0)}{m_{\mu}\omega_{\mu}}\sin{(\omega_{\mu}t)}\right).$

Back to Eq. (35), the corresponding Green’s function
$\mathsf{G}(t,t^{\prime})$ satisfies

$\partial t\,\mathsf{G}(t,t^{\prime})+\begin{pmatrix}0&-1\\\
\omega(t)^{2}+\omega_{R}^{2}&0\end{pmatrix}\mathsf{G}(t,t^{\prime})\\\
-\int_{0}^{t}\begin{pmatrix}0&0\\\
\chi(t-\tau)&0\end{pmatrix}\mathsf{G}(\tau,t^{\prime})\,d\tau=\delta(t-t^{\prime})\,\mathbbm{1}_{2}.$
(36)

The equation has the formal solution

$\boldsymbol{\mathsf{Z}}(t)=\mathsf{G}(t,0)\,\boldsymbol{\mathsf{Z}}(0)+\int_{0}^{t}\mathsf{G}(t,t^{\prime})\,\boldsymbol{\mathsf{F}}(t^{\prime})\,dt^{\prime},$
(37)

which will allow us to find $\mathsf{G}(t,t^{\prime})$ in the limit cycle.

###  A.2 Covariance matrix in the limit cycle

The covariance matrix takes the form

$\mathsf{\Sigma}(t)=\textrm{Re}\,\langle\boldsymbol{\mathsf{Z}}(t)\boldsymbol{\mathsf{Z}}(t)^{\mathsf{T}}\rangle-\langle\boldsymbol{\mathsf{Z}}(t)\rangle\langle\boldsymbol{\mathsf{Z}}(t)^{\mathsf{T}}\rangle\\\
=\begin{pmatrix}\sigma_{xx}(t)&\sigma_{xp}(t)\\\
\sigma_{xp}(t)&\sigma_{pp}(t)\end{pmatrix}.$

For initial states with
$\langle\boldsymbol{\mathsf{Z}}(0)\rangle=\langle\boldsymbol{\mathsf{Z}}_{\mu}(0)\rangle=0$,
the first moments vanish for all times, as can be seen from (37). Since
$\langle\boldsymbol{\mathsf{Z}}_{\mu}(0)\rangle=0$ is automatically fulfilled
for a sample in a thermal state, we only need to assume
$\langle\boldsymbol{\mathsf{Z}}(0)\rangle=0$.

From Eq. (37) the time evolution of the covariance matrix can be written as

$\mathsf{\Sigma}(t)=\mathsf{G}(t,0)\,\mathsf{\Sigma}(0)\,\mathsf{G}(t,0)^{\mathsf{T}}+\mathsf{G}(t,0)\,\textrm{Re}\,\langle\boldsymbol{\mathsf{Z}}(0)\boldsymbol{\mathsf{B}}(t)^{\mathsf{T}}\rangle+\textrm{Re}\,\langle\boldsymbol{\mathsf{B}}(t)\boldsymbol{\mathsf{Z}}(0)^{\mathsf{T}}\rangle\,\mathsf{G}(t,0)^{\mathsf{T}}+\textrm{Re}\,\langle\boldsymbol{\mathsf{B}}(t)\boldsymbol{\mathsf{B}}(t)^{\mathsf{T}}\rangle,$

with
$\boldsymbol{\mathsf{B}}(t)=\int_{0}^{\infty}\,\mathsf{G}(t,t^{\prime})\,\boldsymbol{\mathsf{F}}(t^{\prime})\,dt^{\prime}$.
The second and third term are propagations of initial probe–sample
correlations and therefore vanish in our case. In turn, the first term decays
for $t\rightarrow\infty$ for exponentially stable systems. Here, we assume
that this is indeed the case and later discuss when this assumption fails;
namely, at parametric resonance (see section 4.3 and Appendix C). We are thus
left with

$\displaystyle\mathsf{\Sigma}^{\text{lc}}(t)=\textrm{Re}\,\langle\boldsymbol{\mathsf{B}}(t)\boldsymbol{\mathsf{B}}(t)^{\mathsf{T}}\rangle=\int_{0}^{t}\int_{0}^{t}\mathsf{G}(t,t_{1})\sum\nolimits_{\mu,\nu}\mathsf{C}_{\mu}\mathsf{G}_{\mu}(t_{1},0)\,\textrm{Re}\,\langle\boldsymbol{\mathsf{Z}}_{\mu}(0)\boldsymbol{\mathsf{Z}}_{\nu}(0)^{\mathsf{T}}\rangle\mathsf{G}_{\nu}(t_{2},0)^{\mathsf{T}}\mathsf{C}^{\mathsf{T}}_{\nu}\mathsf{G}(t,t_{2})^{\mathsf{T}}\,dt_{1}\,dt_{2},$
(38)

where, as stated in the main text, the superscript ‘lc’ stands for ‘limit
cycle’. Since the sample is in thermal equilibrium at temperatures $T$

$\textrm{Re}\,\langle\boldsymbol{\mathsf{Z}}_{\mu}(0)\,\boldsymbol{\mathsf{Z}}_{\nu}(0)^{\mathsf{T}}\rangle\\\
=\delta_{\mu,\nu}\begin{pmatrix}\frac{1}{2\,m_{\mu}\omega_{\mu}}\coth\frac{\omega_{\mu}}{2T}&0\\\
0&\frac{m_{\mu}\omega_{\mu}}{2}\coth\frac{\omega_{\mu}}{2T}\end{pmatrix},$

where $\delta_{\mu,\nu}$ is a Kronecker delta. Eq. (38) thus simplifies to

$\mathsf{\Sigma}^{\text{lc}}(t)=\frac{1}{2}\,\int_{0}^{t}\int_{0}^{t}\mathsf{G}(t,t_{1})\\\
\times\begin{pmatrix}0&0\\\
0&\mu(t_{1}-t_{2})\end{pmatrix}\mathsf{G}(t,t_{2})^{\mathsf{T}}dt_{1}dt_{2},$
(39)

with the ‘noise kernel’ $\mu(t)$ given by

$\displaystyle\mu(t)=\frac{2}{\pi}\Theta(t)\int_{0}^{\infty}J(\omega)\,\cos{\omega
t}\,\coth{\frac{\omega}{2T}}\,d\omega.$ (40)

Eq. (39) reduces to the following expressions for the elements of the
covariance matrix

$\displaystyle\sigma_{xx}^{\text{lc}}(t)$
$\displaystyle=\frac{1}{2}\int_{0}^{t}\int_{0}^{t}g(t,t_{1})\mu(t_{1}-t_{2})$
$\displaystyle\quad\qquad\qquad\qquad\times g(t,t_{2})\,dt_{1}\,dt_{2},$ (41a)
$\displaystyle\sigma_{xp}^{\text{lc}}(t)$
$\displaystyle=\frac{1}{2}\int_{0}^{t}\int_{0}^{t}g(t,t_{1})\mu(t_{1}-t_{2})$
$\displaystyle\quad\qquad\qquad\qquad\times\partial_{t}g(t,t_{2})\,dt_{1}\,dt_{2},$
(41b) $\displaystyle\sigma_{pp}^{\text{lc}}(t)$
$\displaystyle=\frac{1}{2}\int_{0}^{t}\int_{0}^{t}\partial_{t}\,g(t,t_{1})\mu(t_{1}-t_{2})$
$\displaystyle\quad\qquad\qquad\qquad\times\partial_{t}g(t,t_{2})\,dt_{1}\,dt_{2},$
(41c)

where we have introduced the notation
$g(t,t^{\prime})\coloneqq[\mathsf{G}(t,t^{\prime})]_{12}$, and used
$[\mathsf{G}(t,t^{\prime})]_{22}=\partial_{t}[\mathsf{G}(t,t^{\prime})]_{12}$,
which follows from Eq. (36). Exploiting this identity, one further sees that

$\partial_{t}^{2}g(t,t^{\prime})+\big{(}\omega^{2}(t)+\omega_{R}^{2}\big{)}\,g(t,t^{\prime})\\\
-\int_{0}^{t}\chi(t-\tau)\,g(\tau,t^{\prime})\,d\tau=\delta(t-t^{\prime}),$
(42)

that is, $g(t,t^{\prime})$ is the Green’s function of the equation of motion
(6) for the system’s $\boldsymbol{x}(t)$.

Eqs. (41) can be manipulated further by casting the noise kernel (40) as

$\mu(t)=\textrm{Re}\,\int_{0}^{\infty}\hat{\mu}(\omega)\,\mathrm{e}^{i\omega
t}\,d\omega,$

with $\hat{\mu}(\omega)\coloneqq\frac{2}{\pi}\,J(\omega)\,\coth{(\omega/2T)}$,
as defined in the main text. This gives

$\displaystyle\sigma_{xx}^{\text{lc}}(t)$
$\displaystyle=\frac{1}{2}\,\textrm{Re}\,\int_{0}^{\infty}q(t,\omega)\,\hat{\mu}(\omega)\,q(t,\omega)^{*}d\omega,$
(43a) $\displaystyle\sigma_{xp}^{\text{lc}}(t)$
$\displaystyle=\frac{1}{2}\,\textrm{Re}\,\int_{0}^{\infty}q(t,\omega)\,\hat{\mu}(\omega)\,\partial_{t}q(t,\omega)^{*}d\omega,$
(43b) $\displaystyle\sigma_{pp}^{\text{lc}}(t)$
$\displaystyle=\frac{1}{2}\,\textrm{Re}\,\int_{0}^{\infty}\partial_{t}q(t,\omega)\,\hat{\mu}(\omega)\,\partial_{t}q(t,\omega)^{*}d\omega,$
(43c)

where $q(t,\omega)=\int_{0}^{t}g(t,t^{\prime})\,\mathrm{e}^{i\omega
t^{\prime}}\,dt^{\prime}$.

Recalling that our probe is periodically driven, i.e.,
$\boldsymbol{H}_{S}(t+\tau)=\boldsymbol{H}_{S}(t)$ with
$\tau=2\pi/\omega_{d}$, we must have
$g(t+\tau,t^{\prime}+\tau)=g(t,t^{\prime})$. Therefore,

$\displaystyle q(t+\tau,\omega)$
$\displaystyle=\int_{0}^{t+\tau}g(t+\tau,t^{\prime})\,\mathrm{e}^{i\omega
t^{\prime}}\,dt^{\prime}$
$\displaystyle=\int_{-\tau}^{t}g(t+\tau,t^{\prime}+\tau)\,\mathrm{e}^{i\omega(t^{\prime}+\tau)}\,dt^{\prime}$
$\displaystyle=\mathrm{e}^{i\omega\tau}\,\int_{-\tau}^{t}g(t,t^{\prime})\,\mathrm{e}^{i\omega
t^{\prime}}\,dt^{\prime}.$

The probe reaching a limit cycle requires $g(t,t^{\prime})$ going to zero
sufficiently fast for $\mathinner{\\!\left\lvert t-t^{\prime}\right\rvert}\gg
1$. In the limit cycle, i.e, $t\gg\tau$, we therefore have

$\int_{-\tau}^{t}g(t,t^{\prime})\,\mathrm{e}^{i\omega
t^{\prime}}\,dt^{\prime}\simeq\int_{0}^{t}g(t,t^{\prime})\,\mathrm{e}^{i\omega
t^{\prime}}\,dt^{\prime}=q(t,\omega).$

That is, $q(t+\tau,\omega)=\mathrm{e}^{i\omega\tau}\,q(t,\omega)$. Hence, the
function $p(t,\omega)\coloneqq\mathrm{e}^{-i\omega t}\,q(t,\omega)$ is
$\tau$-periodic in the limit cycle. As a result, we may write its Fourier
series as

$p(t,\omega)=\sum_{k=-\infty}^{\infty}a_{k}(\omega)\,\mathrm{e}^{ik\omega_{d}t}\\\
=\mathrm{e}^{-i\omega t}\int_{0}^{t}g(t,t^{\prime})\mathrm{e}^{i\omega
t^{\prime}}\,dt^{\prime}.$ (44)

Using $\delta(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{e}^{i\omega
t}d\omega$ it is easy to see that

$g(t,t^{\prime})=\frac{1}{2\pi}\int_{-\infty}^{\infty}p(t,\omega)\,\mathrm{e}^{i\omega(t-t^{\prime})}\,d\omega\\\
=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\int_{-\infty}^{\infty}a_{k}(\omega)\,\mathrm{e}^{i\omega(t-t^{\prime})}\,\mathrm{e}^{ik\omega_{d}t}\,d\omega,$
(45)

which is Eq. (10) from the main text.

Inserting now (45) into (42) gives

$\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega\,\mathrm{e}^{i\omega(t-t^{\prime})}\\\
\times\sum\nolimits_{k}\mathrm{e}^{ik\omega_{d}t}\big{[}a_{k}(\omega)\,\hat{g}_{0}[i(\omega+k\omega_{d})]^{-1}\\\
+\sum\nolimits_{l\neq 0}a_{k-l}(\omega)b_{l}\big{]}=\delta(t-t^{\prime}).$

Using again
$\delta(t-t^{\prime})=\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{e}^{i\omega(t-t^{\prime})}d\omega$
readily brings us to Eq. (11) in the main text, i.e.,

$\hat{g}_{0}[i(\omega+k\omega_{d})]^{-1}a_{k}(\omega)+\sum\nolimits_{l\neq
0}\,b_{l}\,a_{k-l}(\omega)=\delta_{k0},$ (46)

where the frequency has been expanded as

$\omega^{2}(t)=\omega_{0}^{2}+\sum\nolimits_{l\neq
0}b_{l}\,\mathrm{e}^{i\,l\,\omega_{d}\,t}.$

Finally, writing $q(t,\omega)=e^{i\omega t}p(t,\omega)$ and using the Fourier
series (44) allows to rewrite Eqs. (43) as Eqs. (16) and (17) in the main
text.

##  Appendix B Heating of the sample

The rate of change of the energy of the probe can be expressed as

$\displaystyle\frac{d}{dt}\langle\boldsymbol{H}_{s}(t)\rangle$
$\displaystyle=-i\langle\left[\boldsymbol{H}_{S}(t),\boldsymbol{H}(t)\right]\rangle+\partial_{t}\langle\boldsymbol{H}_{S}(t)\rangle$
$\displaystyle=-i\langle\left[\boldsymbol{H}_{S}(t),\boldsymbol{H}(t)\right]\rangle+\frac{1}{2}\partial_{t}\omega(t)^{2}\,\sigma_{xx}(t).$
(47)

Here, the last term is the input power from the external modulation of the
trapping frequency. Hence, we write

$\displaystyle\dot{W}(t)$
$\displaystyle\coloneqq\frac{1}{2}\partial_{t}\omega(t)^{2}\,\sigma_{xx}(t).$

In turn, as we saw in Sec. 4.3 above, the remaining term is the change of
energy due to the the probe–sample interaction, denoted $\dot{Q}(t)$.
Evaluating the commutator in (28) gives

$\displaystyle\dot{Q}(t)$
$\displaystyle=-\langle\textstyle\sum\nolimits_{\mu}g_{\mu}\,\boldsymbol{x}_{\mu}\,\boldsymbol{p}\rangle=-\frac{1}{2}\langle\\{\textstyle\sum_{\mu}g_{\mu}\boldsymbol{x}_{\mu},\boldsymbol{p}\\}\rangle.$

Using now the equation of motion (32b) allows us to write

$\displaystyle\dot{Q}(t)$
$\displaystyle=\frac{1}{2}\langle\\{\dot{\boldsymbol{p}},\boldsymbol{p}\\}\rangle+(\omega(t)^{2}+\omega_{R}^{2})\,\frac{1}{2}\langle\\{\boldsymbol{x},\boldsymbol{p}\\}\rangle$
$\displaystyle=\frac{1}{2}\frac{d}{dt}\sigma_{pp}(t)+(\omega(t)^{2}+\omega_{R}^{2})\,\sigma_{xp}(t).$
(48)

We now place ourselves at the limit cycle, thus working with
$\sigma_{pp}^{\text{lc}}(t)$ and $\sigma_{xp}^{\text{lc}}(t)$, and carry out
the cycle average

$\overline{\dot{Q}}\coloneqq\frac{2\pi}{\omega_{d}}\int_{0}^{2\pi/\omega_{d}}\dot{Q}^{\text{lc}}(t)\,dt.$

Note that the contribution to $\overline{\dot{Q}}$ from the first term on the
right-hand side of (B) would then vanish, as per Eq. (16a). In turn, using
Eqs. (12) and (16b) leads to

$\overline{\dot{Q}}=\frac{1}{2}\text{Im}\,\int_{0}^{\infty}\sum\nolimits_{j,k}b_{k-j}\,(\omega+k\omega_{d})\\\
\times a_{j}(\omega)\,\hat{\mu}(\omega)\,a_{k}^{*}(\omega)\,d\omega,$

where $b_{0}=\omega_{0}^{2}+\omega_{R}^{2}$. We now note that

$\text{Im}\,\sum\nolimits_{j,k}b_{k-j}\,a_{j}(\omega)\,a_{k}^{*}(\omega)=0,$

since $b_{j-k}={(b_{k-j})}^{*}$. Therefore, given that $\hat{\mu}(\omega)$ is
also real, we may write the simplified expression

$\overline{\dot{Q}}=\frac{\omega_{d}}{2}\,\text{Im}\,\int_{0}^{\infty}\sum\nolimits_{j,k}k\,b_{k-j}\,a_{j}(\omega)\,\hat{\mu}(\omega)\,a_{k}^{*}(\omega)\,d\omega.$

Using equation (46) we can write

$\sum\nolimits_{j}b_{k-j}\,a_{j}(\omega)=\delta_{k0}-(\hat{g}_{0}[i(\omega+k\omega_{d})]^{-1}-b_{0})\,a_{k}(\omega),$

so that the limit-cycle-averaged heat current becomes

$\overline{\dot{Q}}=-\frac{\omega_{d}}{2}\,\text{Im}\,\int_{0}^{\infty}\sum\nolimits_{k}k\,(\hat{g}_{0}[i(\omega+k\omega_{d})]^{-1}-b_{0})\\\
\times a_{k}(\omega)\,\hat{\mu}(\omega)\,a_{k}^{*}(\omega)\,d\omega.$

To complete the derivation we use the identity [7]

$\text{Im}\,(\hat{g}_{0}(i\omega)^{-1}-b_{0})=\text{Im}\,\hat{\chi}(i\omega)\\\
=J(\omega)\,\Theta(\omega)-J(-\omega)\,\Theta(-\omega)\coloneqq\tilde{J}(\omega).$

Writing $\hat{\mu}(\omega)$ explicitly as per Eq. (18) finally gives us Eq.
(30) from the main text; namely,

$\dot{\bar{Q}}=-\frac{1}{\pi}\int_{0}^{\infty}\sum_{k}k\omega_{d}J(\omega+k\omega_{d})a_{k}(\omega)\\\
\times J(\omega)a^{*}_{k}(\omega)\coth{\frac{\omega}{2T}}\,d\omega.$ (49)

##  Appendix C Mathieu equation and parametric resonance

As stated above, the damped Mathieu equation reads

$\ddot{x}+\gamma\,\dot{x}+\omega^{2}(t)\,x=0,$

with $\omega^{2}(t)=\omega_{0}^{2}+\upsilon\cos{\omega_{d}t}$. Introducing the
dimensionless parameters $\tilde{t}=\omega_{d}\,t/2$,
$\tilde{\omega}_{0}^{2}=4\omega_{0}^{2}/\omega_{d}^{2}$,
$\tilde{\upsilon}=2\upsilon/\omega_{d}^{2}$, and
$\tilde{\gamma}=2\gamma/\omega_{d}$, the equation can be recast as

$\ddot{x}+\tilde{\gamma}\,\dot{x}+\tilde{\omega}^{2}(\tilde{t})\,x=0,$

with
$\tilde{\omega}^{2}(t)=\tilde{\omega}_{0}^{2}+2\,\tilde{\upsilon}\,\cos{2\tilde{t}}$.
Applying now the transformation
$x=y\,\mathrm{e}^{-\tilde{\gamma}\,\tilde{t}/2}$ allows to absorb the damping
into a frequency shift, bringing the Mathieu equation into the form

$\ddot{y}+(\tilde{\omega}_{0}^{2}-\tilde{\gamma}^{2}/4+2\,\tilde{\upsilon}\,\cos{2\,\tilde{t}})\,y=0.$

This is solved by a function of the form [24, 43]

$y(\tilde{t})=c_{1}\,\mathrm{e}^{i\,\nu\,\tilde{t}}\,p(\tilde{t})+c_{2}\,\mathrm{e}^{-i\,\nu\,\tilde{t}}\,p(-\tilde{t}),$

where $\nu$ is the characteristic exponent, whose dependence on the parameters
$\tilde{\omega}_{0}$, $\tilde{\gamma}$, and $\tilde{v}$ can be readily
explored numerically. The functions $p(\tilde{t})$ here are periodic. Undoing
the previous change of variables back to $x$ then gives

$x(\tilde{t})=c_{1}\,\mathrm{e}^{(i\,\nu-\tilde{\gamma}/2)\,\tilde{t}}p(\tilde{t})+c_{2}\,\mathrm{e}^{-(i\,\nu-\tilde{\gamma}/2)\,\tilde{t}}p(-\tilde{t}).$

The asymptotic behaviour of this solution depends on $\text{Im}\,\nu$; namely,
if $\text{Im}\,\nu>\tilde{\gamma}/2$ the amplitude grows exponentially; if
$\text{Im}\,\nu<\tilde{\gamma}/2$ the amplitude decays to zero; and if
$\text{Im}\,\nu=\tilde{\gamma}/2$, the solution is purely oscillatory.

## References

  * Leanhardt et al. [2003] Aron E. Leanhardt, Thomas A. Pasquini, Michele Saba, André Schirotzek, Yong-il Shin, Dave Kielpinski, David E. Pritchard, and Wolfgang Ketterle.  Cooling bose-einstein condensates below 500 picokelvin.  _Science_ , 301(5639):1513–1515, 2003.  doi: 10.1126/science.1088827.
  * Bloch et al. [2012] Immanuel Bloch, Jean Dalibard, and Sylvain Nascimbene.  Quantum simulations with ultracold quantum gases.  _Nat. Phys._ , 8(4):267–276, 2012.  doi: 10.1038/nphys2259.
  * Olf et al. [2015] Ryan Olf, Fang Fang, G. Edward Marti, Andrew MacRae, and Dan M. Stamper-Kurn.  Thermometry and cooling of a bose gas to 0.02 times the condensation temperature.  _Nat. Phys._ , 11(9):720–723, 2015.  doi: 10.1038/nphys3408.
  * Deppner et al. [2021] Christian Deppner, Waldemar Herr, Merle Cornelius, Peter Stromberger, Tammo Sternke, Christoph Grzeschik, Alexander Grote, Jan Rudolph, Sven Herrmann, Markus Krutzik, André Wenzlawski, Robin Corgier, Eric Charron, David Guéry-Odelin, Naceur Gaaloul, Claus Lämmerzahl, Achim Peters, Patrick Windpassinger, and Ernst M. Rasel.  Collective-mode enhanced matter-wave optics.  _Phys. Rev. Lett._ , 127:100401, 2021.  doi: 10.1103/PhysRevLett.127.100401.
  * Hovhannisyan and Correa [2018] Karen V. Hovhannisyan and Luis A. Correa.  Measuring the temperature of cold many-body quantum systems.  _Phys. Rev. B_ , 98:045101, 2018.  doi: 10.1103/PhysRevB.98.045101.
  * De Pasquale et al. [2016] Antonella De Pasquale, Davide Rossini, Rosario Fazio, and Vittorio Giovannetti.  Local quantum thermal susceptibility.  _Nat. Commun._ , 7(1):1–8, 2016.  doi: 10.1038/ncomms12782.
  * Correa et al. [2017] Luis A. Correa, Martí Perarnau-Llobet, Karen V. Hovhannisyan, Senaida Hernández-Santana, Mohammad Mehboudi, and Anna Sanpera.  Enhancement of low-temperature thermometry by strong coupling.  _Phys. Rev. A_ , 96:062103, 2017.  doi: 10.1103/PhysRevA.96.062103.
  * Miller and Anders [2018] Harry J. D. Miller and Janet Anders.  Energy-temperature uncertainty relation in quantum thermodynamics.  _Nat. Commun._ , 9(1):1–8, 2018.  doi: 10.1038/s41467-018-04536-7.
  * de Pascuale and Stace [2019] Antonella de Pascuale and Tom Stace.  _Thermodynamics in the Quantum Regime. Fundamental Theories of Physics, vol. 195_ , chapter Quantum thermometry.  Springer, Cham, 2019.  ISBN 978-3-319-99046-0.  doi: 10.1007/978-3-319-99046-0_21.
  * Mehboudi et al. [2019a] Mohammad Mehboudi, Anna Sanpera, and Luis A. Correa.  Thermometry in the quantum regime: recent theoretical progress.  _J. Phys. A Math. Theor._ , 52:011611, 2019a.  doi: 10.1088/1751-8121/ab2828.
  * Potts et al. [2019] Patrick P. Potts, Jonatan Bohr Brask, and Nicolas Brunner.  Fundamental limits on low-temperature quantum thermometry with finite resolution.  _Quantum_ , 3:161, 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-07-09-161.
  * Jørgensen et al. [2020] Mathias R. Jørgensen, Patrick P. Potts, Matteo G. A. Paris, and Jonatan B. Brask.  Tight bound on finite-resolution quantum thermometry at low temperatures.  _Phys. Rev. Res._ , 2:033394, 2020.  doi: 10.1103/PhysRevResearch.2.033394.
  * Mukherjee et al. [2019] Victor Mukherjee, Analia Zwick, Arnab Ghosh, Xi Chen, and Gershon Kurizki.  Enhanced precision bound of low-temperature quantum thermometry via dynamical control.  _Commun. Phys._ , 2:162, 2019.  doi: 10.1038/s42005-019-0265-y.
  * Henao et al. [2021] Ivan Henao, Karen V Hovhannisyan, and Raam Uzdin.  Thermometric machine for ultraprecise thermometry of low temperatures.  _arXiv preprint arXiv:2108.10469_ , 2021.  doi: 10.48550/arXiv.2108.10469.
  * Mehboudi et al. [2019b] Mohammad Mehboudi, Aniello Lampo, Christos Charalambous, Luis A. Correa, Miguel Ángel García-March, and Maciej Lewenstein.  Using polarons for sub-nk quantum nondemolition thermometry in a bose-einstein condensate.  _Phys. Rev. Lett._ , 122:030403, 2019b.  doi: 10.1103/PhysRevLett.122.030403.
  * Mishra and Bayat [2021a] Utkarsh Mishra and Abolfazl Bayat.  Integrable quantum many-body sensors for ac field sensing.  _arXiv preprint arXiv:2105.13507_ , 2021a.  doi: 10.48550/arXiv.2105.13507.
  * Mishra and Bayat [2021b] Utkarsh Mishra and Abolfazl Bayat.  Driving enhanced quantum sensing in partially accessible many-body systems.  _Phys. Rev. Lett._ , 127:080504, 2021b.  doi: 10.1103/PhysRevLett.127.080504.
  * Davies [1974] E. Brian Davies.  Markovian master equations.  _Commun. Math. Phys._ , 39(2):91 – 110, 1974\.  doi: 10.1007/BF01608389.
  * Gorini et al. [1976] Vittorio Gorini, Andrzej Kossakowski, and Ennackal Chandy George Sudarshan.  Completely positive dynamical semigroups of n-level systems.  _J. Math. Phys._ , 17(5):821–825, 1976.  doi: 10.1063/1.522979.
  * Lindblad [1976] Goran Lindblad.  On the generators of quantum dynamical semigroups.  _Commun. Math. Phys._ , 48(2):119–130, 1976.  doi: 10.1007/BF01608499.
  * Alicki et al. [2006] Robert Alicki, Daniel A. Lidar, and Paolo Zanardi.  Internal consistency of fault-tolerant quantum error correction in light of rigorous derivations of the quantum markovian limit.  _Phys. Rev. A_ , 73:052311, 2006.  doi: 10.1103/PhysRevA.73.052311.
  * Alicki et al. [2012] Robert Alicki, David Gelbwaser-Klimovsky, and Gershon Kurizki.  Periodically driven quantum open systems: Tutorial.  _arXiv preprint arXiv:1205.4552_ , 2012.  doi: 10.48550/arXiv.1205.4552.
  * Breuer et al. [2002] Heinz-Peter Breuer, Francesco Petruccione, et al.  _The theory of open quantum systems_.  Oxford University Press on Demand, 2002.  ISBN 9780199213900.  doi: 10.1093/acprof:oso/9780199213900.001.0001.
  * Zerbe and Hänggi [1995] Christine Zerbe and Peter Hänggi.  Brownian parametric quantum oscillator with dissipation.  _Phys. Rev. E_ , 52:1533–1543, 1995.  doi: 10.1103/PhysRevE.52.1533.
  * Freitas and Paz [2014] Nahuel Freitas and Juan Pablo Paz.  Analytic solution for heat flow through a general harmonic network.  _Phys. Rev. E_ , 90:042128, 2014.  doi: 10.1103/PhysRevE.90.042128.
  * Freitas and Paz [2017] Nahuel Freitas and Juan Pablo Paz.  Fundamental limits for cooling of linear quantum refrigerators.  _Phys. Rev. E_ , 95:012146, 2017.  doi: 10.1103/PhysRevE.95.012146.
  * Lampo et al. [2017] Aniello Lampo, Soon Hoe Lim, Miguel Ángel García-March, and Maciej Lewenstein.  Bose polaron as an instance of quantum Brownian motion.  _Quantum_ , 1:30, 2017.  ISSN 2521-327X.  doi: 10.22331/q-2017-09-27-30.
  * Lampo et al. [2018] Aniello Lampo, Christos Charalambous, Miguel Ángel García-March, and Maciej Lewenstein.  Non-markovian polaron dynamics in a trapped bose-einstein condensate.  _Phys. Rev. A_ , 98:063630, 2018.  doi: 10.1103/PhysRevA.98.063630.
  * Weiss [1999] Ulrich Weiss.  _Quantum dissipative systems_ , volume 13.  World Scientific, 2nd edition, 1999.  ISBN 978-981-4374-91-0.  doi: 10.1142/8334.
  * Ferraro et al. [2005] Alessandro Ferraro, Stefano Olivares, and Matteo G. A. Paris.  _Gaussian States in Quantum Information_.  Napoli Series on physics and Astrophysics. Bibliopolis, 2005.  ISBN 88-7088-483-X.  doi: 10.48550/arXiv.quant-ph/0503237.
  * Cramér [2016] Harald Cramér.  _Mathematical Methods of Statistics (PMS-9)_.  Princeton University Press, 2016.  ISBN 9781400883868.  doi: doi:10.1515/9781400883868.
  * Rao [1945] Calyampudi Radhakrishna Rao.  Information and the accuracy attainable in the estimation of statistical parameters.  _Reson. J. Sci. Educ_ , 20:78–90, 1945.  doi: 10.1007/978-1-4612-0919-5_16.
  * Rubio et al. [2021] Jesús Rubio, Janet Anders, and Luis A. Correa.  Global quantum thermometry.  _Phys. Rev. Lett._ , 127:190402, 2021.  doi: 10.1103/PhysRevLett.127.190402.
  * Scutaru [1998] Horia Scutaru.  Fidelity for displaced squeezed thermal states and the oscillator semigroup.  _J. Phys. A Math. Gen._ , 31(15):3659–3663, 1998\.  doi: 10.1088/0305-4470/31/15/025.
  * Alipour and Rezakhani [2015] Sahar Alipour and Ali T. Rezakhani.  Extended convexity of quantum fisher information in quantum metrology.  _Phys. Rev. A_ , 91:042104, 2015.  doi: 10.1103/PhysRevA.91.042104.
  * Aveline et al. [2020] David C. Aveline, Jason R. Williams, Ethan R. Elliott, Chelsea Dutenhoffer, James R. Kellogg, James M. Kohel, Norman E. Lay, Kamal Oudrhiri, Robert F. Shotwell, Nan Yu, and Robert J. Thompson.  Observation of Bose–Einstein condensates in an Earth-orbiting research lab.  _Nature_ , 582(7811):193–197, 2020.  doi: 10.1038/s41586-020-2346-1.
  * Spethmann et al. [2012] Nicolas Spethmann, Farina Kindermann, Shincy John, Claudia Weber, Dieter Meschede, and Artur Widera.  Dynamics of single neutral impurity atoms immersed in an ultracold gas.  _Phys. Rev. Lett._ , 109:235301, 2012.  doi: 10.1103/PhysRevLett.109.235301.
  * Hohmann et al. [2016] Michael Hohmann, Farina Kindermann, Tobias Lausch, Daniel Mayer, Felix Schmidt, and Artur Widera.  Single-atom thermometer for ultracold gases.  _Phys. Rev. A_ , 93:043607, 2016.  doi: 10.1103/PhysRevA.93.043607.
  * Catani et al. [2012] Jacopo Catani, Giacomo Lamporesi, Devang Naik, Michael Gring, Massimo Inguscio, Fransesco Minardi, Aadrian Kantian, and Thierry Giamarchi.  Quantum dynamics of impurities in a one-dimensional bose gas.  _Phys. Rev. A_ , 85:023623, 2012.  doi: 10.1103/PhysRevA.85.023623.
  * 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.  doi: 10.1103/PhysRevLett.72.3439.
  * Tóth and Apellaniz [2014] Géza Tóth and Iagoba Apellaniz.  Quantum metrology from a quantum information science perspective.  _J. Phys. A Math. Theor._ , 47:424006, 2014.  doi: 10.1088/1751-8113/47/42/424006.
  * Zmuidzinas [2012] Jonas Zmuidzinas.  Superconducting microresonators: Physics and applications.  _Annu. Rev. Condens. Matter Phys._ , 3(1):169–214, 2012.  doi: 10.1146/annurev-conmatphys-020911-125022.
  * Nayfeh and Mook [1995] Ali Hasan Nayfeh and Dean T. Mook.  _Nonlinear Oscillations_.  John Wiley and Sons, Ltd, 1995.  ISBN 9780471121428.  doi: 10.1002/9783527617586.