# Taking the temperature of a pure quantum state

Mark T. Mitchison<EMAIL_ADDRESS>School of Physics, Trinity College
Dublin, College Green, Dublin 2, Ireland  Archak Purkayastha  School of
Physics, Trinity College Dublin, College Green, Dublin 2, Ireland  Marlon
Brenes  School of Physics, Trinity College Dublin, College Green, Dublin 2,
Ireland  Alessandro Silva  SISSA, Via Bonomea 265, I-34135 Trieste, Italy
John Goold<EMAIL_ADDRESS>School of Physics, Trinity College Dublin, College
Green, Dublin 2, Ireland

###### Abstract

Temperature is a deceptively simple concept that still raises deep questions
at the forefront of quantum physics research. The observation of
thermalisation in completely isolated quantum systems, such as cold-atom
quantum simulators, implies that a temperature can be assigned even to
individual, pure quantum states. Here, we propose a scheme to measure the
temperature of such pure states through quantum interference. Our proposal
involves interferometry of an auxiliary qubit probe, which is prepared in a
superposition state and subsequently undergoes decoherence due to weak
coupling with a closed, thermalised many-body system. Using only a few basic
assumptions about chaotic quantum systems — namely, the eigenstate
thermalisation hypothesis and the emergence of hydrodynamics at long times —
we show that the qubit undergoes pure exponential decoherence at a rate that
depends on the temperature of its surroundings. We verify our predictions by
numerical experiments on a quantum spin chain that thermalises after absorbing
energy from a periodic drive. Our work provides a general method to measure
the temperature of isolated, strongly interacting systems under minimal
assumptions.

Advances in our understanding of thermodynamic concepts have always been
inspired by the technologies of the time, from steam engines in the nineteenth
century to ultra-cold atom simulators in the twenty-first. Irrespective of the
historical era, the importance of measuring temperature cannot be overstated.
In 1798, the American military man and scientist, Count Rumford, noticed that
he could generate heat from friction while boring cannons in the arsenal of
the Bavarian army he was tasked with reorganising. Rumford reported the
systematic temperature increase of the water in which the cannon barrels were
immersed Thompson (1798), challenging the prevailing caloric theory of heat
and inspiring James Joule to perform the decisive experiments that established
energy conservation as the first law of a new thermodynamic theory. In his
famous paddle-bucket experiment, Joule measured the mechanical equivalent of
heat by observing the temperature change induced by stirring fluid in a
thermally isolated container Joule (1850). Here, we show that recasting
Joule’s experiment as a fully quantum-mechanical process leads to a
fundamentally new scheme to measure the temperature of an isolated quantum
many-body system. Our proposal relies on entangling the system with an
auxiliary qubit that undergoes decoherence with a temperature-dependent rate.
This thermometer scale is defined entirely through quantum interference and
allows the measurement of temperature for generic systems in pure quantum
states.

In the last two decades, experimental progress in cold-atom physics has
enabled coherent quantum dynamics to persist over extraordinary timescales:
long enough to observe isolated many-body systems thermalise without coupling
to any external bath Trotzky  _et al._ (2012); Clos  _et al._ (2016); Kaufman
_et al._ (2016); Bordia  _et al._ (2017); Tang  _et al._ (2018). The emergence
of thermodynamics in this context is elegantly explained by the eigenstate
thermalisation hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol  _et
al._ (2008). The ETH posits that, in a sufficiently complex and chaotic
system, each energy eigenstate encodes the properties of the equilibrium
ensemble. As a result, local observables in a far-from-equilibrium scenario
eventually thermalise under unitary evolution D’Alessio  _et al._ (2016). The
final temperature is set by the energy density of the initial condition, which
may be effectively a pure quantum state. Thermal fluctuations thus arise
locally because of quantum entanglement between different parts of the system
Goldstein  _et al._ (2006); Popescu  _et al._ (2006) rather than by any
classical statistical mixing. This begs the question: can the temperature of a
pure state also be measured in a completely quantum-mechanical way?

Figure 1: Illustration of an experiment where work is performed on a thermally
isolated system — such as (a) a bucket of water or (b) an ultracold atomic gas
— thus driving it into a non-equilibrium state. After the external force is
removed, collisions between particles lead to irreversible thermalisation at a
temperature determined by the energy density of the initial state, even though
in (b) the global evolution is unitary and the system is described by a pure
quantum state. The final temperature can be inferred by entangling the system
to a qubit probe and measuring the resulting decoherence rate. Figure 2:
Unitary heating of a quantum spin-$\tfrac{1}{2}$ chain. The Hamiltonian reads
$\hat{H}=J\sum_{j=1}^{L}\left(\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}+\Delta\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}\right)+h\sum_{j\,\rm
odd}\hat{\sigma}_{j}^{z}$, where $\hat{\sigma}^{x,y,z}_{j}$ are Pauli
operators pertaining to lattice site $j$. The exchange coupling $J$ and
anisotropy $J\Delta$ respectively describe the kinetic and interaction energy
of conserved spin excitations, while $h$ is a staggered field acting on odd
sites of the chain. We take periodic boundary conditions, half filling
($N=L/2$ excitations), and choose parameters $\Delta=0.55$ and $h=J$ for which
the model is known to be robustly non-integrable Brenes  _et al._ (2020a).
Results are for $L=26$ sites unless otherwise indicated. Numerical methods are
described in the Supplemental Material SM . (a) Mean energy
$\bar{E}=\braket{\psi(t_{\rm prep})}{\hat{H}}{\psi(t_{\rm prep})}$ of the
chain as a function of the preparation time $t_{\rm prep}$ under local
periodic driving,
$\hat{H}(t)=\hat{H}+a\sin(\omega_{0}t)\hat{\sigma}_{j_{0}}^{z}$, with
amplitude $a=2J$ and frequency $\omega_{0}=8J$ applied to an odd-numbered
site, $j_{0}=L/2$ (or $j_{0}=L/2+1$ if $L/2$ is even). The inset shows the
energy fluctuations, $\Delta E^{2}=\braket{\psi(t_{\rm
prep})}{(\hat{H}-\bar{E})^{2}}{\psi(t_{\rm prep})}$, relative to
$\bar{E}_{*}=\bar{E}-E_{0}$ with $E_{0}$ the ground-state energy, as a
function of system size at fixed microcanonical temperature $T=10J$. (b)
Energy distribution of the prepared state,
$|\psi(E)|^{2}=\sum_{n}|\braket{E_{n}}{\psi(t_{\rm
prep})}|^{2}\delta(E-E_{n})$, where $\ket{E_{n}}$ is an eigenstate of
$\hat{H}$ with energy $E_{n}$. (c) Equilibration of the local magnetisation on
site $j_{0}$ after the drive is switched off. Dashed lines show the
corresponding microcanonical average. (d) Time-averaged local magnetisation
after equilibration (black dots, obtained by averaging over a time interval
$\delta t\geq 20J^{-1}$) compared with the microcanonical average (blue line).
(e) Auto-correlation function of the local operator
$\hat{A}=\sum_{j}u_{j}\hat{\sigma}^{z}_{j}$, where
$u_{j}\propto\mathrm{e}^{-(j-j_{0})^{2}}$ is a Gaussian profile
($\sum_{j}u_{j}=1$). Lines show the result for $t-t_{\rm prep}=100J^{-1}$,
after the system has relaxed to equilibrium, while squares indicate near-
identical values for $t-t_{\rm prep}=110J^{-1}$. (f) Inverse temperature
estimated by fitting the low-frequency noise and response functions to the FDT
$\tilde{\chi}^{\prime\prime}(\omega)/\tilde{S}(\omega)=\tanh(\beta\omega/2)$
(black dots) and the corresponding microcanonical prediction (blue line).

Our pure-state thermometry scheme, depicted in Fig. 1, draws inspiration from
Joule’s pioneering experiment, for which thermal isolation was vital. We
consider the extreme case of an isolated quantum system such as an ultra-cold
atomic gas. Work is performed by changing some external constraint, thus
driving the system out of equilibrium in analogy to Joule’s paddles. The
driving force is then removed and the system relaxes under unitary evolution.
Local observables thermalise to a temperature governed by the work performed,
i.e. the mechanical equivalent of heat. Joule’s apparatus included an in situ
thermometer to measure the temperature change of the insulated fluid. In our
setup, this role is played by an auxiliary qubit that becomes entangled with
the many-body system. Assuming only the ETH and the equations of diffusive
hydrodynamics, we show that the qubit undergoes pure exponential decoherence
at a rate with a universal temperature dependence which can be
interferometrically measured Cetina  _et al._ (2015, 2016); Skou  _et al._
(2021), providing a uniquely quantum thermometer for pure states. Our
inherently non-equilibrium protocol offers a general solution to the problem
of thermometry for isolated quantum systems, while avoiding the limitations of
small thermal probes that equilibrate with the system Mehboudi  _et al._
(2019).

The quantum equivalent of Joule’s paddle bucket is best illustrated by a
specific example, although our scheme is general. Fig. 2 details an in silico
experiment where a thermally isolated many-body system is heated by periodic
driving Bunin  _et al._ (2011); D’Alessio and Rigol (2014); Lazarides  _et
al._ (2014). We simulate an archetypal model of a quantum chaotic system: a
Heisenberg spin-$\tfrac{1}{2}$ chain Jepsen  _et al._ (2020); Scheie  _et al._
(2021) with a staggered magnetic field that breaks integrability Brenes  _et
al._ (2020a). Units with $\hbar=k_{B}=1$ are adopted throughout. At time
$t=0$, the chain is prepared in its ground state with energy $E_{0}$. An
oscillatory field is then applied locally, pumping energy steadily into the
system until the drive is switched off at time $t_{\rm prep}$ [Fig. 2(a)].
This procedure generates a class of non-equilibrium pure states whose average
energy $\bar{E}$ can be selected by tuning the preparation time. These states
have a structured energy distribution featuring sharp peaks spaced by the
drive frequency [Fig. 2(b)]. Importantly, the corresponding energy
fluctuations $\Delta E$ are sub-extensive, meaning that $\Delta
E/(\bar{E}-E_{0})\equiv\Delta E/\bar{E}_{*}$ decreases with system size [Fig.
2(a) inset].

After the drive is switched off, the system evolves autonomously and local
observables relax to equilibrium [Fig. 2(c)], exhibiting small fluctuations
around a value that is close to the prediction of the microcanonical ensemble
[Fig. 2(d)]. This ensemble is characterised by a single parameter: the average
energy, $\bar{E}$, with the corresponding inverse temperature
$T^{-1}\equiv\beta=\beta(\bar{E})$ given by the fundamental definition
$\beta(E)=\mathrm{d}\mathcal{S}/\mathrm{d}E$, where $\mathcal{S}(E)$ is the
microcanonical entropy. Similar thermal behaviour is observed in correlation
functions like
$C(t^{\prime},t)=\braket{\hat{A}(t^{\prime})\hat{A}(t)}-\braket{\hat{A}(t^{\prime})}\braket{\hat{A}(t)}$,
with $\hat{A}$ a local observable, which become approximately stationary at
long times, i.e. $C(t+\tau,t)\approx C(\tau)$ [Fig. 2(e)]. Conventionally, one
writes $C(\tau)$ in terms of the symmetrised noise function
$S(\tau)=\,\mathrm{Re}\,[C(\tau)]$ and the dissipative response function
$\chi^{\prime\prime}(\tau)=\mathrm{i}\,\mathrm{Im}\,[C(\tau)]$. After
relaxation, their Fourier transforms are related by the fluctation-dissipation
theorem (FDT),
$\tilde{S}(\omega)=\coth(\beta\omega/2)\tilde{\chi}^{\prime\prime}(\omega)$,
as expected in thermal equilibrium [Fig. 2(f)].

While the thermalisation of these “paddle-bucket” preparations is striking in
light of the highly non-equilibrium energy distribution displayed in Fig.
2(b), it is also completely generic. Any other initial condition with sub-
extensive energy fluctuations — e.g. as created by a sudden parameter quench —
will similarly thermalise in a non-integrable system at energies well above
the ground state D’Alessio  _et al._ (2016). This behaviour is fully explained
by the ETH, which can be formulated as an ansatz for the matrix elements of an
arbitrary local observable, $\hat{A}$, in the energy eigenbasis Srednicki
(1999), i.e. $A_{mn}=\langle E_{m}|\hat{A}|E_{n}\rangle$, where
$\hat{H}\ket{E_{n}}=E_{n}\ket{E_{n}}$. The ansatz reads as

$A_{mn}=\begin{cases}\hskip
5.69054ptA(E_{n})+\mathcal{O}(\mathcal{D}^{-1/2}),&m=n,\\\ \hskip
5.69054pt\mathrm{e}^{-\mathcal{S}(E_{mn})/2}f(E_{mn},\omega_{mn})R_{mn}+\mathcal{O}(\mathcal{D}^{-1}),&m\neq
n,\end{cases}$ (1)

where the meaning of each term on the right-hand side is illustrated in Fig. 3
and explained as follows. The diagonal matrix elements ($m=n$) are given by a
smooth function of energy, $A(E)$, up to fluctuations that scale inversely
with the square root of the Hilbert-space dimension, $\mathcal{D}$ [Fig.
3(a)]. Therefore, all energy eigenstates near a given energy $E$ yield the
same expectation value $\langle\hat{A}\rangle=A(E)$ in the thermodynamic
limit. This identifies $A(E)$ as the microcanonical average of $\hat{A}$ at
inverse temperature $\beta(E)$. Meanwhile, the off-diagonal matrix elements
($m\neq n$) are exponentially small and erratically distributed [Fig. 3(b)
inset], as described in Eq. (1) by a Hermitian matrix $R_{mn}$ of random
numbers with zero mean and unit variance. Underlying this distribution is a
smooth spectral function $f(E,\omega)$ of the mean energy,
$E_{mn}=\tfrac{1}{2}(E_{m}+E_{n})$ and transition frequency,
$\omega_{mn}=E_{m}-E_{n}$, which is revealed in the variance of the matrix
elements within small energy and frequency windows [Fig. 3(b) main panel].

Figure 3: Eigenstate thermalisation in the staggered-field Heisenberg spin
chain [see Fig. 2 caption for details]. (a) The diagonal elements of the local
operator $\hat{A}$ concentrate around a smooth function of the energy density,
$\epsilon_{n}=(E_{n}-E_{\rm min})/(E_{\rm max}-E_{\rm min})$. The black line
shows a running average for $L=20$ within microcanonical windows of width
$\delta\epsilon=0.02$. The inset of (a) shows the variance of matrix elements
evaluated within the central 10% of the spectrum for different system sizes,
plotted against the Hilbert-space dimension, $\mathcal{D}$. The scaling
$\mathcal{D}^{-1}$ is indicated by the dashed red line. (b) The low-frequency
spectral function for $L=18$ and three different temperatures. The data are
generated by finding all off-diagonal matrix elements at energy $E_{mn}$
consistent with a given temperature, $T=[\beta(E_{mn})]^{-1}$, then computing
the variance, ${\rm
var}[A_{mn}]=\mathrm{e}^{-\mathcal{S}(E_{mn})}|f(E_{mn},\omega_{mn})|^{2}$,
within small frequency windows $\delta\omega\sim 0.2J$. The inset shows the
magnitude of the matrix elements near $T=5J$ (grey points; only one in every
hundred is shown) and a running average of $|A_{mn}|$ (black line).

Aside from energy eigenstates, the ETH also describes the ergodic dynamics of
non-equilibrium pure states seen in Fig. 2. Using Eq. (1), the time average of
an observable is found to be
$\overline{\braket{\hat{A}}}=\sum_{n}|\braket{E_{n}}{\psi}|^{2}A_{nn}=A(\bar{E})+\mathcal{O}(\Delta
E^{2}/\bar{E}^{2}_{*})$, while temporal fluctuations away from this value are
proportional to the off-diagonal elements $A_{mn}$ and thus are exponentially
suppressed D’Alessio  _et al._ (2016). Moreover, thermal statistics arise for
correlation functions because both noise and response are determined (up to
sub-extensive corrections) by the spectral function as D’Alessio  _et al._
(2016); Brenes  _et al._ (2020a)

$\displaystyle\tilde{S}(\omega)=2\pi\cosh(\beta\omega/2)|f(\bar{E},\omega)|^{2},$
(2)
$\displaystyle\tilde{\chi}^{\prime\prime}(\omega)=2\pi\sinh(\beta\omega/2)|f(\bar{E},\omega)|^{2},$
(3)

immediately implying the FDT. While these features of the ETH have long been
understood, the low-frequency behaviour of the spectral function has only
recently been identified as a sensitive indicator of quantum chaos in many-
body systems Brenes  _et al._ (2020b); Pandey  _et al._ (2020). For a generic
observable in a non-integrable system, $f(E,0)$ is non-zero and may vary
significantly with temperature [Fig. 3(b)]. This observation forms the basis
of our thermometry scheme.

A traditional strategy is to let the thermometer exchange energy with its
surroundings and simply wait for equilibration. This working principle was
used in Joule’s instruments Joule (1850), for example, but becomes challenging
to implement at low temperature, where a precise thermometer needs small
energy scales and correspondingly long thermalisation times Correa  _et al._
(2015). We thus take a fundamentally different approach based on pure
decoherence Johnson  _et al._ (2016); Razavian  _et al._ (2019); Mitchison
_et al._ (2020), which is independent of the probe’s energy. Our thermometer
comprises a qubit with energy eigenstates $\ket{\uparrow}$ and
$\ket{\downarrow}$, coupled to the system by an interaction of the form
$\hat{H}_{\rm int}=\ket{\uparrow}\bra{\uparrow}\otimes g\hat{A}$ for some
local observable $\hat{A}$ and coupling constant $g$. This kind of interaction
— which can be engineered, for example, using Feshbach resonances in ultracold
gases Cetina  _et al._ (2016) — conserves the qubit’s energy and ensures that
it does not participate in the dynamics while in its ground state
$\ket{\downarrow}$.

Suppose that at time $t_{0}$, the thermal system of interest is described by
the pure state $\ket{\psi(t_{0})}=\ket{\psi_{0}}$. The protocol begins by
exciting the qubit into a superposition
$\ket{+}=\tfrac{1}{\sqrt{2}}(\ket{\uparrow}+\ket{\downarrow})$ with a
$\pi/2-$pulse, which prepares the joint product state
$\ket{\Psi(t_{0})}=\ket{+}\ket{\psi_{0}}$. In a frame rotating at the qubit
precession frequency, the Schrödinger evolution is then
$\ket{\Psi(t)}=\frac{1}{\sqrt{2}}\left(\mathrm{e}^{-\mathrm{i}\hat{H}(t-t_{0})}\ket{\downarrow}\ket{\psi_{0}}+\mathrm{e}^{-\mathrm{i}(\hat{H}+g\hat{A})(t-t_{0})}\ket{\uparrow}\ket{\psi_{0}}\right)$.
Entanglement develops between the probe and the system, leading to a loss of
distinguishability quantified by the fidelity

$|\varv(t)|^{2}=|\braket{\psi_{0}}{\mathrm{e}^{\mathrm{i}\hat{H}(t-t_{0})}\mathrm{e}^{-\mathrm{i}(\hat{H}+g\hat{A})(t-t_{0})}}{\psi_{0}}|^{2}.$
(4)

The resulting decrease in interference contrast is reflected in the off-
diagonal elements of the qubit density matrix, $\hat{\rho}_{\rm
q}(t)=\mathrm{Tr}_{\rm sys}\,\ket{\Psi(t)}\bra{\Psi(t)}$, which decay in time
according to $\braket{\downarrow\rvert\hat{\rho}_{\rm
q}(t)\lvert\uparrow}=\tfrac{1}{2}\varv(t)$. This decoherence is finally probed
by applying a second $\pi/2$-pulse with a phase $\theta$ relative to the first
one, then measuring the excited-state probability of the qubit,
$P_{\uparrow}=\tfrac{1}{2}(1+\,\mathrm{Re}\,[\mathrm{e}^{\mathrm{i}\theta}\varv(t)]).$
The time-dependent overlap $\varv(t)$ is thus reconstructed by varying
$\theta$.

To assess the temperature dependence of the interference contrast, we focus on
the weak-coupling regime and approximate the fidelity (4) by a cumulant
expansion to second order in $g$ SM . We obtain
$|\varv(t)|^{2}=\mathrm{e}^{-\Gamma(t)}$, where

$\displaystyle\Gamma(t)=4g^{2}\int\frac{\mathrm{d}\omega}{2\pi}\,\tilde{S}(\omega)\frac{\sin^{2}[\omega(t-t_{0})/2]}{\omega^{2}}.$
(5)

At weak coupling, the largest effects are seen for $t-t_{0}\gg\tau_{c}$, where
$\tau_{c}$ is the characteristic timescale for the correlation function
$C(\tau)$ to decay to zero. The integral in Eq. (5) is then dominated by the
contribution near $\omega=0$, which implies pure exponential decoherence,
$|\varv(t)|^{2}\sim\mathrm{e}^{-\gamma(t-t_{0})}$, with an asymptotic decay
rate $\gamma=g^{2}\tilde{S}(0)\propto|f(\bar{E},0)|^{2}$. We numerically
confirm this behaviour in Fig. 4(a), which shows the fidelity for a probe
coupled to a spin chain heated by the procedure of Fig. 2. Even for moderate
coupling strengths, we observe near-perfect exponential decay with a
temperature-dependent rate in close agreement with the weak-coupling
prediction. The decoherence is associated with a growth in the entanglement
entropy $\mathsf{S}[\hat{\rho}_{\rm q}]=-\mathrm{Tr}[\hat{\rho}_{\rm
q}\ln\hat{\rho}_{\rm q}]$, which saturates to the temperature-independent
value $\mathsf{S}[\hat{\rho}_{\rm q}]\to\ln 2$ characterising a maximally
entangled state [Fig. 4(b)]. This distinguishes our non-equilibrium protocol
from a thermalisation process. In Fig. 4(c), the temperature dependence of the
decoherence rate is analysed in more detail. We find that $\gamma$ depends
almost linearly on energy density [Fig. 4(c) inset], which translates into a
non-linear variation with temperature [Fig. 4(c) main panel] that is greatest
at low temperatures.

Figure 4: Decoherence of a qubit with coupling $g=0.2J$ to a spin-chain
environment prepared in a pure thermal state as shown in Fig. 2. (a) The
fidelity and (b) the entanglement entropy for three different temperatures.
Solid lines show an exact calculation of Eq. (4), while the dashed lines show
the weak-coupling approximation
$|\varv(t)|^{2}=\mathrm{e}^{-\gamma(t-t_{0})}$, with
$\gamma=g^{2}\tilde{S}(0)$ extracted from dynamical calculations of $C(\tau)$.
We take $t_{0}-t_{\rm prep}=100J^{-1}$ for $\bar{E}=-4J,-8J$ and $t_{0}-t_{\rm
prep}=200J^{-1}$ for $\bar{E}=-16J$. (c) Temperature dependence of the
asymptotic decoherence rate, $\gamma=g^{2}\tilde{S}(0)$, with the energy
dependence as an inset. (d) QFI (black dots) and Fisher information for a
measurement in the qubit eigenbasis (red triangles), computed within the weak-
coupling approximation as a function of temperature. Both quantities are
evaluated at the time $t^{*}$ where the QFI is maximised, i.e.
$\mathcal{F}^{Q}_{T}=\mathcal{F}^{Q}_{T}(t^{*})\equiv\max_{t}\mathcal{F}^{Q}_{T}(t)$
and $\mathcal{F}^{\parallel}_{T}=\mathcal{F}^{\parallel}_{T}(t^{*})$.

We quantify the temperature information that can be extracted from our
protocol using the quantum Fisher information (QFI). Consider a temperature
estimate constructed from $M$ independent measurements in a given basis,
$\mu$, on identical qubit preparations. For large $M$, the statistical error
of any unbiased estimate is asymptotically bounded by $\Delta T^{2}\geq
1/M\mathcal{F}^{\mu}_{T}\geq 1/M\mathcal{F}_{T}^{Q}$. Here,
$\mathcal{F}^{\mu}_{T}$ is the Fisher information for the chosen basis while
the QFI, $\mathcal{F}_{T}^{Q}=\max_{\mu}\mathcal{F}_{T}^{\mu}$, is the maximum
over all measurements and thus describes the ultimate uncertainty limit
imposed by quantum mechanics Braunstein and Caves (1994). The temperature can
be inferred from the exponential decay of $|\varv(t)|$ by measuring in the
eigenbasis of $\hat{\rho}_{\rm q}(t)$, i.e. by applying a final $\pi/2$-pulse
with phase $\theta=-\arg\varv(t)$ SM . Fig. 4(d) shows the corresponding
Fisher information, $\mathcal{F}^{\parallel}_{T}$, in the weak-coupling limit.
Since $\mathcal{F}_{T}^{\parallel}\approx\mathcal{F}^{Q}_{T}$, we conclude
that the decoherence rate captures almost all temperature information
available from the probe in this example. For instance, we obtain the value
$T^{2}\mathcal{F}^{\parallel}_{T}\approx 0.2$ at temperature $T=5J$, implying
that $M=500$ measurements could suffice to achieve a precision of $\Delta
T/T\gtrsim 10\%$. Note that a single ultra-cold gas sample may host thousands
of independent impurities Cetina  _et al._ (2016).

Our results show that the temperature of an isolated system can be measured
using the most primitive features of quantum dynamics: namely, unitarily
evolving wavefunctions and entanglement between subsystems. The scale of our
thermometer is defined not through the energetic fluctuations of some
statistical mixture, but by the rate of entanglement growth in a quantum
decoherence process. While this rate should generally increase with
temperature, the precise dependence is system- and observable-specific.
However, since a generic system is expected to display hydrodynamic behaviour
at long timescales, we can predict a general form for $\gamma(T)$ assuming
that the probe couples to diffusive modes of a conserved density. In three
spatial dimensions (3D), we obtain SM

$\gamma(T)=\frac{2\bar{g}^{2}\chi_{0}T}{D},$ (6)

where $D$ is the diffusion coefficient, $\chi_{0}$ is the thermodynamic
susceptibility to long-wavelength density perturbations and $\bar{g}$ is a
renormalised coupling that depends only on the spatial profile of the probe.
According to Eq. (6), the qubit’s decoherence rate provides an ideal, linear
thermometer scale within any temperature range where $D$ and $\chi_{0}$ are
approximately constant, and allows for accurate thermometry in general
whenever $D$ and $\chi_{0}$ are known as a function of temperature.

Due to the complexity of simulating 3D quantum-chaotic dynamics, we have
instead presented a numerical example in 1D. Here, as in the 3D case,
exponential decoherence arises at long times due to the finite low-frequency
spectral function, $|f(\bar{E},0)|^{2}$, a generic feature of chaotic systems
with finite volume Brenes  _et al._ (2020b); Pandey  _et al._ (2020). In the
thermodynamic limit, however, diffusion in 1D implies that
$|f(\bar{E},0)|^{2}$ is proportional to the system size, with the low-
frequency noise diverging as $\tilde{S}(\omega)\sim\omega^{-1/2}$ SM . Eq. (5)
then predicts a crossover from non-exponential decoherence, $\Gamma(t)\sim
t^{3/2}$, to exponential decay, $\Gamma(t)\sim t$, after a timescale
commensurate with the Thouless time, which grows with volume (a similar
conclusion holds in 2D with different scalings) SM . Although the crossover
occurs too rapidly to clearly distinguish at system sizes accessible to our
numerics, the temperature of very large, low-dimensional, diffusive systems
must therefore be inferred from the dominant non-exponential decoherence at
intermediate times. This interesting competition of timescales calls for
further research to characterise how Markovian dynamics Srednicki (1999);
Nation and Porras (2019); Parra-Murillo  _et al._ (2021) and thermodynamics
Iyoda  _et al._ (2017); Riera-Campeny  _et al._ (2021) emerges for open
quantum systems in chaotic environments.

Accurate, in situ thermometry of isolated quantum systems is an outstanding
problem in cold-atom physics, where strong, short-ranged correlations confound
destructive global measurement techniques such as time-of-flight imaging.
Conversely, a small quantum probe facilitates local, minimally destructive
temperature measurements, in principle Hohmann  _et al._ (2016); Bouton  _et
al._ (2020). Our proposal to infer temperature from decoherence dynamics does
not require thermalisation of the qubit nor fine-tuning of its energy levels,
and is applicable to generic many-body systems in arbitrary non-equilibrium
preparations (including mixed states) with sub-extensive energy fluctuations.
This opens a pathway for the toolbox of quantum-enhanced thermometry Mehboudi
_et al._ (2019) to probe the ultimate limit of an isolated system in a pure
quantum state.

Acknowledgements. We thank S. R. Clark, C. Jarzynski, A. Polkovnikov, and J.
Richter for useful feedback on the manuscript. M. B. and J. G. thank M. Rigol
for stimulating their interest in the ETH. A. P. acknowledges funding from the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No. 890884. Some calculations were
performed on the Lonsdale cluster maintained by the Trinity Centre for High
Performance Computing. This cluster was funded through grants from Science
Foundation Ireland (SFI). We acknowledge the DJEI/DES/SFI/HEA Irish Centre for
High-End Computing (ICHEC) for the provision of computational facilities,
Project No. TCPHY138A. This work was supported by a SFI-Royal Society
University Research Fellowship (J. G.) and the Royal Society (M. B.). J. G.
and M. T. M. acknowledge funding from the European Research Council Starting
Grant ODYSSEY (Grant Agreement No. 758403) and the EPSRC-SFI joint project
QuamNESS.

## References

  * Thompson (1798) B. Thompson, Philos. Trans. R. Soc. London 88, 80 (1798).
  * Joule (1850) J. P. Joule, Philos. Trans. R. Soc. London 140, 61 (1850).
  * Trotzky  _et al._ (2012) S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Nature Phys. 8, 325 (2012).
  * Clos  _et al._ (2016) G. Clos, D. Porras, U. Warring, and T. Schaetz, Phys. Rev. Lett. 117, 170401 (2016).
  * Kaufman  _et al._ (2016) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
  * Bordia  _et al._ (2017) P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Nature Phys. 13, 460 (2017).
  * Tang  _et al._ (2018) Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Phys. Rev. X 8, 021030 (2018).
  * Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Rigol  _et al._ (2008) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * Goldstein  _et al._ (2006) S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì, Phys. Rev. Lett. 96, 050403 (2006).
  * Popescu  _et al._ (2006) S. Popescu, A. J. Short, and A. Winter, Nature Phys. 2, 754 (2006).
  * Brenes  _et al._ (2020a) M. Brenes, S. Pappalardi, J. Goold, and A. Silva, Phys. Rev. Lett. 124, 040605 (2020a).
  * (15) See the Supplemental Material for further details on the numerical examples, the weak-coupling expansion of the fidelity, the quantum Fisher information, and the predictions of diffusive hydrodynamics for different spatial dimensionalities.
  * Cetina  _et al._ (2015) M. Cetina, M. Jag, R. S. Lous, J. T. M. Walraven, R. Grimm, R. S. Christensen, and G. M. Bruun, Phys. Rev. Lett. 115, 135302 (2015).
  * Cetina  _et al._ (2016) M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. M. Walraven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt, M. Knap, and E. Demler, Science 354, 96 (2016).
  * Skou  _et al._ (2021) M. G. Skou, T. G. Skov, N. B. Jørgensen, K. K. Nielsen, A. Camacho-Guardian, T. Pohl, G. M. Bruun, and J. J. Arlt, Nature Phys. (2021), 10.1038/s41567-021-01184-5.
  * Mehboudi  _et al._ (2019) M. Mehboudi, A. Sanpera, and L. A. Correa, J. Phys. A: Math. Theor. 52, 303001 (2019).
  * Bunin  _et al._ (2011) G. Bunin, L. D’Alessio, Y. Kafri, and A. Polkovnikov, Nature Phys. 7, 913 (2011).
  * D’Alessio and Rigol (2014) L. D’Alessio and M. Rigol, Phys. Rev. X 4, 041048 (2014).
  * Lazarides  _et al._ (2014) A. Lazarides, A. Das, and R. Moessner, Phys. Rev. E 90, 012110 (2014).
  * Jepsen  _et al._ (2020) P. N. Jepsen, J. Amato-Grill, I. Dimitrova, W. W. Ho, E. Demler, and W. Ketterle, Nature 588, 403 (2020).
  * Scheie  _et al._ (2021) A. Scheie, N. E. Sherman, M. Dupont, S. E. Nagler, M. B. Stone, G. E. Granroth, J. E. Moore, and D. A. Tennant, Nature Phys. (2021), 10.1038/s41567-021-01191-6.
  * Srednicki (1999) M. Srednicki, J. Phys. A 32, 1163 (1999).
  * Brenes  _et al._ (2020b) M. Brenes, J. Goold, and M. Rigol, Phys. Rev. B 102, 075127 (2020b).
  * Pandey  _et al._ (2020) M. Pandey, P. W. Claeys, D. K. Campbell, A. Polkovnikov, and D. Sels, Phys. Rev. X 10, 041017 (2020).
  * Correa  _et al._ (2015) L. A. Correa, M. Mehboudi, G. Adesso, and A. Sanpera, Phys. Rev. Lett. 114, 220405 (2015).
  * Johnson  _et al._ (2016) T. H. Johnson, F. Cosco, M. T. Mitchison, D. Jaksch, and S. R. Clark, Phys. Rev. A 93, 053619 (2016).
  * Razavian  _et al._ (2019) S. Razavian, C. Benedetti, M. Bina, Y. Akbari-Kourbolagh, and M. G. A. Paris, Eur. Phys. J. Plus 134, 284 (2019).
  * Mitchison  _et al._ (2020) M. T. Mitchison, T. Fogarty, G. Guarnieri, S. Campbell, T. Busch, and J. Goold, Phys. Rev. Lett. 125, 080402 (2020).
  * Braunstein and Caves (1994) S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 3439 (1994).
  * Nation and Porras (2019) C. Nation and D. Porras, Phys. Rev. E 99, 052139 (2019).
  * Parra-Murillo  _et al._ (2021) C. A. Parra-Murillo, M. Bramberger, C. Hubig, and I. De Vega, Phys. Rev. A 103, 032204 (2021).
  * Iyoda  _et al._ (2017) E. Iyoda, K. Kaneko, and T. Sagawa, Phys. Rev. Lett. 119, 100601 (2017).
  * Riera-Campeny  _et al._ (2021) A. Riera-Campeny, A. Sanpera, and P. Strasberg, PRX Quantum 2, 010340 (2021).
  * Hohmann  _et al._ (2016) M. Hohmann, F. Kindermann, T. Lausch, D. Mayer, F. Schmidt, and A. Widera, Phys. Rev. A 93, 043607 (2016).
  * Bouton  _et al._ (2020) Q. Bouton, J. Nettersheim, D. Adam, F. Schmidt, D. Mayer, T. Lausch, E. Tiemann, and A. Widera, Phys. Rev. X 10, 011018 (2020).
  * Mondaini and Rigol (2017) R. Mondaini and M. Rigol, Phys. Rev. E 96, 012157 (2017).
  * Khatami  _et al._ (2013) E. Khatami, G. Pupillo, M. Srednicki, and M. Rigol, Phys. Rev. Lett. 111, 050403 (2013).
  * Brenes  _et al._ (2020c) M. Brenes, T. LeBlond, J. Goold, and M. Rigol, Phys. Rev. Lett. 125, 070605 (2020c).
  * LeBlond and Rigol (2020) T. LeBlond and M. Rigol, Phys. Rev. E 102, 062113 (2020).
  * Steinigeweg  _et al._ (2015) R. Steinigeweg, J. Gemmer, and W. Brenig, Phys. Rev. B 91, 104404 (2015).
  * Weiße  _et al._ (2006) A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, Rev. Mod. Phys. 78, 275 (2006).
  * Yang  _et al._ (2020) Y. Yang, S. Iblisdir, J. I. Cirac, and M. C. Bañuls, Phys. Rev. Lett. 124, 100602 (2020).
  * Kubo (1962) R. Kubo, J. Phys. Soc. Japan 17, 1100 (1962).
  * Jaynes (2003) E. T. Jaynes,  _Probability Theory_ (Cambridge University Press, 2003).
  * Kadanoff and Martin (1963) L. P. Kadanoff and P. C. Martin, Ann. Phys. 24, 419 (1963).
  * Forster (1990) D. Forster,  _Hydrodynamic fluctuations, broken symmetry, and correlation functions_ (Perseus Books, 1990).

## Supplemental Material

###  S1 Numerical methods

In this section we provide further details on the model and the methods used
to analyse the quantitative examples discussed in the main text.

####  S1.1 Hamiltonian and observables

The examples in the main text are based on the Hamiltonian

$\hat{H}=J\sum_{j=1}^{L}\left(\hat{\sigma}_{j}^{x}\hat{\sigma}_{j+1}^{x}+\hat{\sigma}_{j}^{y}\hat{\sigma}_{j+1}^{y}+\Delta\hat{\sigma}_{j}^{z}\hat{\sigma}_{j+1}^{z}\right)+h\sum_{j\,\rm
odd}\hat{\sigma}_{j}^{z},$ (S1)

with periodic boundary conditions. This Hamiltonian conserves the number of
spin excitations, $\hat{N}=\tfrac{1}{2}\sum_{j}(1+\hat{\sigma}^{z}_{j})$, and
we work in the half-filled symmetry sector of states with eigenvalue $N=L/2$.
The bulk parameters are chosen to be $h=J$ and $\Delta=0.55$, following Ref.
Brenes  _et al._ (2020a). Some of our examples focus on the local operator
$\hat{A}=\sum_{j}u_{j}\hat{\sigma}_{j}^{z}$, where
$u_{j}\propto\mathrm{e}^{-(j-j_{0})^{2}}$ is a Gaussian profile centred on
site $j_{0}$, where $j_{0}=L/2$ if $L/2$ is odd and $j_{0}=L/2+1$ if $L/2$ is
even. To improve numerical efficiency, we set $u_{j}=0$ on all sites where
$\mathrm{e}^{-(j-j_{0})^{2}}<10^{-3}$, and then normalise as
$\sum_{j}u_{j}=1$. This generates an observable with support restricted to
five sites of the lattice.

It is important to remark that, in order to obtain the function $f(E,\omega)$
from the coarse-grained average of the off-diagonal matrix elements of a local
operator Mondaini and Rigol (2017); Khatami  _et al._ (2013); Brenes  _et al._
(2020a, c, b), all symmetries of the model should be resolved. Resolving these
symmetries amounts to restricting the block-diagonal Hamiltonian to a single
symmetry sub-sector of states corresponding to a given eigenvalue of the
symmetry generator (performing a separate calculation for each block, if more
than one is required). If the operators admix the symmetry sub-sectors LeBlond
and Rigol (2020), resolving the corresponding symmetry is not required to
obtain $f(E,\omega)$ from the coarse-grained average of the off-diagonal
matrix elements.

The function $f(E,\omega)$ obtained from this procedure is shown in Fig. 3(b).
In contrast with the open-boundary chain, the model with periodic boundary
conditions is translation-invariant. Instead of resolving this symmetry, we
break it by augmenting the magnetic field acting on site $j=1$ by a small
amount $\delta h=0.1h$. Even with the addition of this small perturbation, in
the zero-magnetisation sector, an underlying spatial reflection symmetry
remains. This symmetry is broken by the operator $\hat{A}$ when $L/2$ is odd,
as in Fig. 3(b).

####  S1.2 Dynamical evolution

To compute dynamical quantities we solve the Schrödinger equation

$\mathrm{i}\partial_{t}\ket{\psi(t)}=\hat{H}\ket{\psi(t)},$ (S2)

to obtain the state vector $\ket{\psi(t)}$, from which any observable
$\braket{\hat{A}(t)}=\braket{\psi(t)}{\hat{A}}{\psi(t)}$ can be calculated.
The decoherence function itself is obtained from the overlap
$\varv(t)=\braket{\psi(t)}{\psi^{\prime}(t)}$, where the states
$\ket{\psi(t)}$ and $\ket{\psi^{\prime}(t)}$ are propagated under Eq. (S2)
with Hamiltonians $\hat{H}$ and $\hat{H}^{\prime}=\hat{H}+g\hat{A}$,
respectively. We integrate Eq. (S2) using a standard fourth-order Runge-Kutta
Steinigeweg  _et al._ (2015) algorithm with time step $J\delta t=0.01$. This
choice yields an excellent approximation to unitary evolution, e.g. after an
evolution time $Jt=50$ the normalisation $\braket{\psi(t)}{\psi(t)}$ drops by
less than 0.1% at temperature $T=20J$ and by less than 1% at $T=5J$. For very
long evolutions, such as those required to compute $\varv(t)$ at weak
coupling, we normalise the state after each time step to enforce unitarity.

To compute two-point correlation functions, we follow the procedure described,
for example, in Ref. Steinigeweg  _et al._ (2015). For $t^{\prime}>t$ we can
write

$\displaystyle\braket{\hat{A}(t^{\prime})\hat{A}(t)}$
$\displaystyle=\braket{\psi(t_{0})}{\mathrm{e}^{\mathrm{i}\hat{H}(t^{\prime}-t_{0})}\hat{A}\mathrm{e}^{-\mathrm{i}\hat{H}(t^{\prime}-t)}\hat{A}\mathrm{e}^{-\mathrm{i}\hat{H}(t-t_{0})}}{\psi(t_{0})}$
$\displaystyle=\braket{\psi(t^{\prime})}{\hat{A}}{\phi(t^{\prime})},$ (S3)

where $\ket{\psi(t_{0})}$ and $t_{0}$ are the (arbitrary) initial state and
time, and
$\ket{\phi(t^{\prime})}=\mathrm{e}^{-\mathrm{i}\hat{H}(t^{\prime}-t)}\hat{A}\ket{\psi(t)}$
is obtained by propagating the Schrödinger equation (S2) for a time interval
$t^{\prime}-t$ starting from the initial condition $\hat{A}\ket{\psi(t)}$. The
case $t^{\prime}<t$ is obtained by complex conjugation.

Figure S1: (a) Real part (solid line) and imaginary part (dotted line) of the
correlation function for a non-equilibrium pure state at energy
$\bar{E}=-16J$, after thermalisation as described in Fig. 2 of the main text.
The inset zooms in on the the real part at longer times for four different
system sizes at the same temperature, $T=5J$. (b) Noise and response functions
for $L=26$, obtained from the correlation function by numerical Fourier
transform up to a cutoff $J\tau_{*}=10$. The black solid line shows their
ratio, which approximately obeys the fluctuation-dissipation relation within
the range of frequencies where $\tilde{\chi}^{\prime\prime}(\omega)$ is non-
zero. The dashed line shows $\tanh(\beta\omega/2)$ for comparison, with
$\beta(\bar{E})=0.20J^{-1}$ obtained independently from a KPM calculation [see
Sec. S1.3].

Fig. S1(a) shows an example of a correlation function computed in this way.
Both real and imaginary parts of the correlation function decay over time,
after which they execute small fluctuations near zero. These long-time
fluctuations are a finite-size effect: they disappear rapidly with increasing
system size [Fig. S1(a) inset]. To avoid this effect when constructing the
noise and response functions by numerical Fourier transform, we sample the
correlation functions only up to a time $\tau_{*}$, chosen to be
$J\tau_{*}=10$. This is sufficiently large to capture all non-trivial
dynamical features of $C(\tau)$, while generating a relatively smooth
approximation to the noise and response function as shown in Fig. S1(b).
Larger values of $\tau_{*}$ tend to generate spurious features in the
frequency domain due to the long-time fluctuations of $C(\tau)$, which are
most prevalent at lower temperatures. To extract a temperature from the
Fourier data, as in Fig. 2(f), we average the slope of
$\tilde{\chi}^{\prime\prime}(\omega)/\tilde{S}(\omega)$ over the frequency
range $|\omega|\leq 2J$, within which
$\tanh(\beta\omega/2)\approx\beta\omega/2$ to an excellent approximation for
the temperatures of interest.

####  S1.3 Microcanonical predictions

Predictions of the microcanonical ensemble are evaluated using the kernel
polynomial method (KPM) Weiße  _et al._ (2006); Yang  _et al._ (2020). The
three quantities that we compute in this way are the density of states,
$\Omega(E)$, the microcanonical expectation value, $A(E)$, and the local
density of states, $|\psi(E)|^{2}$, given respectively by

$\displaystyle\Omega(E)=\sum_{n}\delta(E-E_{n}),$ (S4) $\displaystyle
A(E)=\frac{1}{\Omega(E)}\sum_{n}A_{nn}\delta(E-E_{n}),$ (S5)
$\displaystyle|\psi(E)|^{2}=\sum_{n}|\braket{E_{n}}{\psi}|^{2}\delta(E-E_{n}).$
(S6)

The microcanonical entropy and temperature are extracted from the density of
states via Boltzmann’s relation $\mathcal{S}(E)=\ln W(E)$, where
$W(E)=\Omega(E)\mathrm{d}E$ corresponds to the number of microstates in a
small energy interval $\mathrm{d}E$.

The kernel polynomial method works by expanding the above functions in the
basis of orthogonal Chebyshev polynomials, $T_{n}(E)$. Since these polynomials
are defined only on the interval $E\in[-1,1]$, we first rescale the
Hamiltonian spectrum to lie within this interval. Then, following the standard
procedure detailed in Ref. Weiße  _et al._ (2006), we approximate a continuous
function of energy $\Xi(E)$ by evaluating a finite number of Chebyshev
moments, $\mu_{m}=\int\mathrm{d}E\,\Xi(E)T_{m}(E)$. The function is then
reconstructed as the truncated basis expansion

$\Xi(E)\approx\frac{1}{\pi\sqrt{1-E^{2}}}\left[g_{0}\mu_{0}+2\sum_{m=1}^{M_{\rm
Cheb}}g_{m}\mu_{m}T_{m}(E)\right],$ (S7)

where $g_{m}$ are coefficients that decay with increasing $m$, which smooth
the high-frequency oscillations (Gibbs phenomenon) that would otherwise result
from truncating the expansion at finite order, $M_{\rm Cheb}$. We use the
values of $g_{m}$ corresponding to the Jackson kernel, which is known to be
optimal and effectively broadens the delta functions entering Eqs. (S4)–(S6)
to Gaussians of width $\lesssim\pi/M_{\rm Cheb}$; see Ref. Weiße  _et al._
(2006) for details.

For example, the density of states is evaluated from the moments

$\mu_{m}=\mathrm{Tr}[T_{m}(\hat{H})]\approx\frac{1}{R}\sum_{r=1}^{R}\braket{r}{T_{m}(\hat{H})}{r}.$
(S8)

The last step approximates the trace using a set of $R$ random vectors, where
each component of the vector $\ket{r}$ is chosen independently from a Gaussian
distribution. In a large Hilbert space only a few random vectors, say
$R\lesssim 10$ for spin chain of length $L=26$, are typically needed to
achieve excellent convergence within the relevant energy range for
thermalisation. In the form (S8), the moments can be computed recursively
using the defining three-term relation of the Chebyshev polynomials,

$\displaystyle T_{m+1}(\hat{H})=2\hat{H}T_{m}(\hat{H})-T_{m-1}(\hat{H}),$ (S9)

for $m\geq 1$, with $T_{0}(\hat{H})=1$ and $T_{1}(\hat{H})=\hat{H}$. Acting
this expression on a vector as $\ket{v_{m}}=T_{m}(\hat{H})\ket{r}$ generates
the recursion relation $\ket{v_{m+1}}=2\hat{H}\ket{v_{m}}-\ket{v_{m-1}}$. Each
iteration thus requires a single matrix-vector multiplication, making the KPM
a very efficient method for large, sparse matrices such as Eq. (S1).

Figure S2: Density of states for a spin chain of $L=20$ sites. The bars show a
normalised histogram of the energy levels obtained from exact diagonalisation.
The blue solid line is the density of states obtained from KPM with $M_{\rm
Cheb}=100$ moments averaged over $R=150$ random vectors. Dividing $\Omega(E)$
by the Hilbert-space dimension ensures normalisation:
$\int\mathrm{d}E\,\Omega(E)/\mathcal{D}=1$.

In all our calculations of the inverse temperature and microcanonical
averages, we take $M_{\rm Cheb}=100$ moments in total. This is sufficient to
obtain an accurate, smooth approximation, as shown in Fig. S2 for the density
of states, for example. Larger values of $M_{\rm Cheb}$ can lead to spurious
oscillations due to the underlying discreteness of the finite-dimensional
Hilbert space, especially at low temperatures. For the local density of states
shown in Fig. 2(b), we take $M_{\rm Cheb}=250$ in order to better resolve the
sharp features.

###  S2 Weak-coupling and long-time limit of the decoherence function

In this section we derive the weak-coupling expansion of the decoherence
function, $\varv(t)$, and discuss its convergence to the long-time limit. To
simplify the notation we set the initial time of the protocol to be $t_{0}=0$
in this section. Our starting point is the general expression for the
decoherence function [c.f. Eq. (4)]

$\varv(t)=\Braket{\mathrm{e}^{\mathrm{i}\hat{H}t}\mathrm{e}^{-\mathrm{i}(\hat{H}+g\hat{A})t}}=\Braket{{\rm
T}\exp\left[-\mathrm{i}g\int_{0}^{t}\mathrm{d}t^{\prime}\hat{A}(t^{\prime})\right]},$
(S10)

where we recognise the definition of the time-ordered exponential, with
$\hat{A}(t)=\mathrm{e}^{\mathrm{i}\hat{H}t}\hat{A}\mathrm{e}^{-\mathrm{i}\hat{H}t}$
the Heisenberg-picture operator with respect to the Hamiltonian $\hat{H}$.
Here, and in the following, expectation values are taken with respect to the
initial state of the many-body system,
$\braket{\bullet}=\mathrm{Tr}[\bullet\hat{\rho}]$, e.g.
$\hat{\rho}=\ket{\psi_{0}}\bra{\psi_{0}}$ in the case of a pure initial state
as considered in the main text. Note, however, that our framework can be
equally applied to any mixed state for which the one- and two-point functions
of the operator $\hat{A}$ are approximately thermal.

We proceed by expanding Eq. (S10) into time-ordered cumulants Kubo (1962) and
neglecting terms of order $\mathcal{O}(g^{3})$, which yields

$\displaystyle-\ln\varv(t)$
$\displaystyle\approx\mathrm{i}g\int_{0}^{t}\mathrm{d}t^{\prime}\,\braket{\hat{A}(t^{\prime})}+g^{2}\int_{0}^{t}\mathrm{d}t^{\prime}\int_{0}^{t^{\prime}}\mathrm{d}t^{\prime\prime}\,C(t^{\prime},t^{\prime\prime})$
(S11) $\displaystyle\approx\mathrm{i}gt\overline{\braket{\hat{A}}}$
$\displaystyle\quad+g^{2}\int_{0}^{t}\mathrm{d}t^{\prime}\int_{0}^{t^{\prime}}\mathrm{d}t^{\prime\prime}\,\left[S(t^{\prime}-t^{\prime\prime})+\chi^{\prime\prime}(t^{\prime}-t^{\prime\prime})\right],$
$\displaystyle=\mathrm{i}\Phi(t)+\frac{1}{2}\Gamma(t)$ (S12)

where the first line is the second-order cumulant expansion, while the second
line follows by assuming the noise statistics are approximately stationary. On
the final line, we moved to the Fourier domain, e.g.
$S(\tau)=\int\mathrm{d}\omega\,\mathrm{e}^{-\mathrm{i}\omega\tau}\tilde{S}(\omega)/2\pi$,
and defined

$\displaystyle\Gamma(t)=4g^{2}\int\frac{\mathrm{d}\omega}{2\pi}\,\tilde{S}(\omega)\frac{\sin^{2}(\omega
t/2)}{\omega^{2}},$ (S13)
$\displaystyle\Phi(t)=gt\overline{\braket{\hat{A}}}+g^{2}\int\frac{\mathrm{d}\omega}{2\pi}\,\tilde{\chi}^{\prime\prime}(\omega)\frac{\sin(\omega
t)-\omega t}{\omega^{2}}.$ (S14)

In order to analyse convergence to the long-time limit, it is convenient to
remain in the time domain. Consider the second-order contribution to Eq. (S11)
in the stationary approximation

$\displaystyle\int_{0}^{t}\mathrm{d}t^{\prime}\int_{0}^{t^{\prime}}\mathrm{d}t^{\prime\prime}\,C(t^{\prime}-t^{\prime\prime})=\int_{0}^{t}\mathrm{d}\tau\,(t-\tau)C(\tau),$
(S15)

where we introduced the variables $\tau=t^{\prime}-t^{\prime\prime}$ and
$\bar{t}=(t^{\prime}+t^{\prime\prime})/2$ and performed the trivial integral
over $\bar{t}\in[\tau/2,t-\tau/2]$. Let $\tau_{c}$ denote the correlation time
after which $C(\tau)$ has decayed to zero. Assuming that the improper integral
$\int_{0}^{\infty}\mathrm{d}\tau\,\tau C(\tau)$ exists, it gives a sub-leading
(i.e. constant) correction for large $t$. For $t\gg\tau_{c}$ we can therefore
neglect this term and take the upper integration limit to infinity, obtaining

$\int_{0}^{t}\mathrm{d}\tau\,(t-\tau)C(\tau)\approx
t\int_{0}^{\infty}\mathrm{d}\tau\,C(\tau)=\frac{1}{2}\left[\tilde{S}(0)-\mathrm{i}\chi_{\hat{A}}\right],$
(S16)

where
$\chi_{\hat{A}}=\int\mathrm{d}\omega\,\tilde{\chi}^{\prime\prime}(\omega)/\pi\omega$
is the thermodynamic susceptibility corresponding to the observable $\hat{A}$.
This result also follows from taking $t\to\infty$ directly in Eqs. (S13) and
(S14). Eq. (S16) implies the emergence of pure exponential decay for times
$t\gg\tau_{c}$, with the rate
$\gamma=\lim_{t\to\infty}\mathrm{d}\Gamma/\mathrm{d}t=g^{2}\tilde{S}(0)$.
Therefore, self-consistency of the exponential approximation requires that the
decoherence function evolves slowly so that $\gamma\tau_{c}\ll 1$. This
condition is well satisfied by our examples, where the correlation function
fully decays after a time $J\tau_{c}\lesssim 10$ [Fig. S1(a)], while the
characteristic timescale for $\varv(t)$ is an order of magnitude longer [Figs.
4(a)].

Note that the above arguments break down whenever $C(\tau)\sim\tau^{-p}$ with
$p\leq 1$ for large $\tau$. In particular, diffusion in the thermodynamic
limit implies that $S(\tau)=\,\mathrm{Re}\,[C(\tau)]\sim(D\tau)^{-d/2}$ in $d$
spatial dimensions with diffusion coefficient $D$, as discussed in Sec. S4.
For $d=1$, therefore, Eq. (S15) is dominated by the second term in parentheses
on the right-hand side, which implies the long-time behaviour $\Gamma(t)\sim
t^{3/2}$. For $d=2$, we obtain the asymptotic scaling $\Gamma(t)\sim
t\ln(Dt/\ell^{2})$, where $\ell$ is the length scale characterising the probe
[see Eq. (S4)]. In a finite system, however, diffusive dynamics persists up to
the Thouless time, $t_{\rm T}\sim L^{2}/D$ ($L$ is the linear dimension of the
system, defined in Sec. S4), after which the correlation function $C(\tau)$
drops to zero (up to small fluctuations, see Fig. S1). One therefore expects
Eq. (S16) to hold for low-dimensional diffusive systems after a correlation
time scaling as $\tau_{c}\sim t_{\rm T}$.

###  S3 Quantum Fisher information

In this section we discuss the quantum Fisher information and its
contributions from the norm and phase of the decoherence function. Let us
first briefly recap the meaning of the Fisher information in the context of
parameter estimation. Suppose that the qubit probe is in the state
$\hat{\rho}_{\rm q}(T)$, which depends on the temperature $T$. We consider a
measurement described by a set of positive operators, $\\{\hat{\Pi}(\xi)\\}$,
such that $\int\mathrm{d}\xi\,\hat{\Pi}(\xi)=1$, where the possible
measurement outcomes are labelled by the continuous index $\xi$ without loss
of generality. A temperature estimate $T_{\rm est}(\bm{\xi})$ is constructed
from the outcomes $\bm{\xi}=\\{\xi_{1},\ldots,\xi_{M}\\}$ of a large number,
$M$, of identical measurements on independent preparations of the state
$\hat{\rho}_{\rm q}(T)$. We consider (asymptotically) unbiased estimators such
as maximum likelihood estimation, which satisfy $\mathbb{E}[T_{\rm est}]=T$,
where

$\mathbb{E}[T_{\rm
est}]=\int\mathrm{d}\xi_{1}\cdots\int\mathrm{d}\xi_{M}p(\xi_{1}|T)\cdots
p(\xi_{M}|T)T_{\rm est}(\bm{\xi}),$ (S17)

while $p(\xi|T)=\mathrm{Tr}[\hat{\Pi}(\xi)\hat{\rho}_{\rm q}(T)]$ denotes the
probability of obtaining outcome $\xi$ in a single measurement. The expected
error in the temperature estimate is thus $\Delta T^{2}=\mathbb{E}[(T_{\rm
est}-T)^{2}]$. This obeys the Cramér-Rao bound Jaynes (2003), $\Delta
T^{2}\geq 1/M\mathcal{F}_{T}$, where the Fisher information is given by

$\mathcal{F}_{T}=\int\mathrm{d}\xi\,p(\xi|T)\left(\frac{\partial\ln
p(\xi|T)}{\partial T}\right)^{2},$ (S18)

which measures the sensitivity of the distribution to changes in the parameter
$T$. The Fisher information depends on the choice of measurement basis, and is
upper-bounded by the quantum Fisher information Braunstein and Caves (1994)
(QFI), $\mathcal{F}_{T}\leq\mathcal{F}_{T}^{Q}$. The bound is saturated by the
measurement of a specific observable: the symmetric logarithmic derivative
(SLD), $\hat{\Lambda}_{T}$.

For a pure dephasing evolution, the qubit state is of the form
$\hat{\rho}_{\rm q}=\tfrac{1}{2}(1+\bm{v}\cdot\bm{\hat{\sigma}}),$ where
$\bm{\hat{\sigma}}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})^{\sf
T}$ is a vector of Pauli operators and
$\bm{v}=(\,\mathrm{Re}\,[\varv],-\,\mathrm{Im}\,[\varv],0)^{\sf T}$ is the
Bloch vector. Parameterising the decoherence function as
$\varv=|\varv|\mathrm{e}^{-\mathrm{i}\phi}$, the QFI takes the form Mitchison
_et al._ (2020)

$\mathcal{F}_{T}^{Q}=\frac{1}{1-|\varv|^{2}}\left(\frac{\partial|\varv|}{\partial
T}\right)^{2}+|\varv|^{2}\left(\frac{\partial\phi}{\partial
T}\right)^{2}=\mathcal{F}^{\parallel}_{T}+\mathcal{F}^{\perp}_{T}.$ (S19)

These two terms respectively correspond to the Fisher information for
measurements of
$\hat{\sigma}_{\parallel}=\cos(\phi)\hat{\sigma}_{x}-\sin(\phi)\hat{\sigma}_{y}$
and
$\hat{\sigma}_{\perp}=\cos(\phi)\hat{\sigma}_{y}+\sin(\phi)\hat{\sigma}_{x}$,
i.e. the bases parallel and perpendicular to $\bm{v}$ in the equatorial plane
of the Bloch sphere [see Fig. 1 of the main text]. Up to irrelevant additive
and multiplicative factors, the SLD is given by

$\displaystyle\hat{\Lambda}_{T}\propto\cos(\varphi)\hat{\sigma}_{\parallel}+\sin(\varphi)\hat{\sigma}_{\perp},$
(S20)
$\displaystyle\tan(\varphi)=\frac{|\varv|(1-|\varv|)^{2}\partial_{T}\phi}{\partial_{T}|\varv|}.$

Figure S3: (a) Phase of the decoherence function with coupling strength
$g=0.2J$ at three different temperatures, as in Fig. 4(a) of the main text.
Solid lines show the exact calculation and dashed lines show the corresponding
asymptotic approximation,
$\phi=\Phi(t)\to(g\overline{\braket{\hat{A}}}-g^{2}\chi_{\hat{A}}/2)t$ [see
Eq. (S14)], obtained from a dynamical calculation of $C(\tau)$ [see Sec. S1.2]
(b) Temperature dependence of the asymptotic phase accumulation rate,
$\dot{\Phi}=\mathrm{d}\Phi/\mathrm{d}t$, in the weak-coupling limit, as
extracted from a dynamical calculation of $\chi^{\prime\prime}(\tau)$.

The main text considers thermometry using measurements of the norm
$|\varv(t)|$, neglecting the phase $\phi$. In Fig. S3(a) we plot the phase at
the same coupling $g=0.2J$ as considered in the main text. The results are
compared to the asymptotic weak-coupling approximation,
$\phi\approx\Phi(t)\to(g\overline{\braket{\hat{A}}}-g^{2}\chi_{\hat{A}}/2)t$
in Fig. S3. The exact results noticeably diverge from the weak-coupling
prediction at longer times, presumably because of higher-order effects not
captured by the perturbative expansion to order $\mathcal{O}(g^{2})$. The
difficulty of predicting the phase accurately with perturbative approximations
is well known, e.g. from studies of the Fermi edge singularity where the
cumulant expansion for $\phi$ diverges already at second order Mitchison  _et
al._ (2020). In Fig. 4(b) we plot the temperature dependence of phase
contributions in the weak-coupling limit, showing that they depend more weakly
on temperature than the corresponding decoherence rate in this example. As a
consequence, the QFI (S19) is dominated by the first contribution,
$\mathcal{F}^{\parallel}_{T}$, as shown in Fig. 4(d) of the main text. The
corresponding measurement of $\hat{\sigma}_{\parallel}$ could be enacted by
tuning the phase of the final $\pi/2$-pulse, $\theta$, to equal the phase of
the decoherence function, i.e. $\theta=\phi$, as discussed in the main text.

To obtain a smooth prediction for the Fisher information in the weak-coupling
limit, as shown in Fig. 4(d), we construct a spline interpolation of the data
shown in Figs. 4(c) and S3(b). The derivatives in Eq. (S19) are then
approximated by a first-order finite difference with $\delta T\approx 0.2J$.

###  S4 Hydrodynamic response function

In this section, we discuss the form of the response function that follows
from diffusive hydrodynamics, and the resulting temperature dependence of the
decoherence rate for large systems. We closely follow the classic analysis of
Kadanoff & Martin Kadanoff and Martin (1963); Forster (1990), generalising it
to a large but finite system in $d$ spatial dimensions. We consider a
translation-invariant system of volume $L^{d}$ under periodic boundary
conditions, described by the Hamiltonian $\hat{H}$. Let
$\hat{A}=\int\mathrm{d}\mathbf{r}\,u(\mathbf{r})\hat{n}(\mathbf{r})$ be a
local observable written in terms of the density, $\hat{n}(\mathbf{r})$, of a
globally conserved charge,
$\hat{N}=\int\mathrm{d}\mathbf{r}\,\hat{n}(\mathbf{r})$, such that
$[\hat{H},\hat{N}]=0$. The slowly varying function $u(\mathbf{r})$ represents
the spatial profile of the probe, which we normalise as
$\int\mathrm{d}\mathbf{r}\,u(\mathbf{r})=1$ without loss of generality.

Translation invariance allows us to separately analyse each Fourier component
of the density,
$\hat{n}_{\mathbf{k}}=\int\mathrm{d}\mathbf{r}\,\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\hat{n}(\mathbf{r})$,
where the discrete wavevector $\mathbf{k}$ has components $k_{j}=2\pi n_{j}/L$
with $n_{j}\in\mathbb{Z}$ and $j=1,\ldots d$. For any state that is invariant
under space and time translations we can then define the density response
function via

$\frac{1}{2L^{d}}\braket{[\hat{n}_{\mathbf{k}}(t+\tau),\hat{n}_{-\mathbf{k}^{\prime}}(t)]}=\delta_{\mathbf{k}\mathbf{k}^{\prime}}\chi^{\prime\prime}_{\mathbf{k}}(\tau),$
(S21)

with $\tilde{\chi}^{\prime\prime}_{\mathbf{k}}(\omega)$ the Fourier transform
of $\chi^{\prime\prime}_{\mathbf{k}}(\tau)$. Note that
$\chi^{\prime\prime}_{\mathbf{k}=0}(\tau)=0$ identically as a consequence of
the conservation of total charge, $\hat{N}=\hat{n}_{\mathbf{k}=0}$. The probe
observable reads
$\hat{A}=L^{-d}\sum_{\mathbf{k}}c_{-\mathbf{k}}\hat{n}_{\mathbf{k}}$, with the
corresponding response function
$\tilde{\chi}^{\prime\prime}(\omega)=L^{-d}\sum_{\mathbf{k}\neq
0}|u_{\mathbf{k}}|^{2}\tilde{\chi}_{\mathbf{k}}^{\prime\prime}(\omega)$, where
$u_{\mathbf{k}}$ is the Fourier transform of $u(\mathbf{r})$. For example, if
$u(\mathbf{r})\propto e^{-r^{2}/2\ell^{2}}$ is a Gaussian of width $\ell$,
$u_{\mathbf{k}}=\mathrm{e}^{-\ell^{2}k^{2}/2}$ cuts off wavevectors
$k\gg\ell^{-1}$; other smooth profiles show similar behaviour. Therefore, so
long as $u(\mathbf{r})$ is slowly varying (i.e. $\ell$ is large), only long-
wavelength diffusive modes contribute significantly to
$\tilde{\chi}^{\prime\prime}(\omega)$.

To find the response function associated with diffusion, a small density
modulation is introduced by adiabatically switching on a weak force,
$F(\mathbf{r})$, and then suddenly removing it at $t=0$. This is modelled by
the potential

$\hat{V}(t)=-\Theta(-t)\mathrm{e}^{\varepsilon
t}\int\mathrm{d}\mathbf{r}\,F(\mathbf{r})\hat{n}(\mathbf{r}),$ (S22)

where $\Theta(t)$ is the unit step function and $\varepsilon\to 0$ is a small,
non-negative convergence parameter. The resulting density deviation is then
computed from linear-response theory, starting from an equilibrium state at
the initial time $t_{0}\to-\infty$. For $t<0$, the system adiabatically
follows the slow perturbation, thus remaining in thermal equilibrium,
$\langle\hat{n}_{\mathbf{k}}(t\leq
0)\rangle=\langle\hat{n}_{\mathbf{k}}\rangle_{\rm eq}$. Meanwhile, linear-
response theory yields $\langle\hat{n}_{\mathbf{k}}(t\leq
0)\rangle=\chi_{\mathbf{k}}F_{\mathbf{k}}$, where

$\chi_{\mathbf{k}}=\int\mathrm{d}\omega\,\frac{\tilde{\chi}^{\prime\prime}_{\mathbf{k}}(\omega)}{\pi\omega}=\left.\frac{\partial\langle\hat{n}_{\mathbf{k}}\rangle_{\rm
eq}}{\partial F_{\mathbf{k}}}\right\rvert_{F=0},$ (S23)

which is identified as the thermodynamic susceptibility. For $t>0$, we obtain

$\braket{\hat{n}_{\mathbf{k}}(t)}=\int\mathrm{d}\omega\,\frac{\tilde{\chi}^{\prime\prime}_{\mathbf{k}}(\omega)F_{\mathbf{k}}}{\pi\omega}\mathrm{e}^{-\mathrm{i}\omega
t},$ (S24)

which, assuming a slowly varying $F(\mathbf{r})$, should evolve according to
the diffusion equation
$(\partial_{t}+Dk^{2})\braket{\hat{n}_{\mathbf{k}}(t)}=0$. This is readily
solved by
$\braket{\hat{n}_{\mathbf{k}}(t)}=\chi_{\mathbf{k}}F_{\mathbf{k}}\mathrm{e}^{-Dk^{2}t}$,
given the initial condition at $t=0$. Comparing these two solutions for
$\braket{\hat{n}_{\mathbf{k}}(t)}$ in the Laplace domain, we get

$\int\frac{\mathrm{d}\omega^{\prime}}{\mathrm{i}\pi}\frac{\tilde{\chi}^{\prime\prime}_{\mathbf{k}}(\omega^{\prime})}{\omega^{\prime}(\omega^{\prime}-z)}=\frac{\chi_{\mathbf{k}}}{Dk^{2}-\mathrm{i}z},$
(S25)

where the Laplace variable obeys $\,\mathrm{Im}\,z>0$. Finally, by continuing
$z\to\omega+\mathrm{i}0$ to the real axis one deduces the density response
function $\tilde{\chi}^{\prime\prime}_{\mathbf{k}}(\omega)$. The response
function for $\hat{A}$ then follows as

$\tilde{\chi}^{\prime\prime}(\omega)=\frac{1}{L^{d}}\sum_{\mathbf{k}\neq
0}\frac{\chi_{\mathbf{k}}|u_{\mathbf{k}}|^{2}Dk^{2}\omega}{\omega^{2}+(Dk^{2})^{2}}.$
(S26)

In a sufficiently large system, the summation is well approximated by an
integral using the standard prescription
$L^{-d}\sum_{\mathbf{k}}\to(2\pi)^{-d}\int\mathrm{d}\mathbf{k}$. If we also
assume that $u_{\mathbf{k}}$ samples only small wavevectors such that
$\chi_{\mathbf{k}}\approx\lim_{\mathbf{k}\to
0}\chi_{\mathbf{k}}\equiv\chi_{0}$, then the integral can be carried out
explicitly. At small positive frequencies, we obtain the limiting behaviour

$\tilde{\chi}^{\prime\prime}(\omega)\sim\begin{cases}\sqrt{\omega}&(d=1)\\\
-\omega\ln(\ell^{2}\omega/D)&(d=2)\\\ \omega&(d\geq 3).\end{cases}$ (S27)

For $d=2$ only, the low-frequency response function depends sensitively on the
ultraviolet (UV) cutoff scale, $D/\ell^{2}$, where $\ell$ is the width of the
probe function $u(\mathbf{r})$ discussed above. For $d\geq 3$, the response is
Ohmic and the zero-frequency limit of the symmetrised noise
$\tilde{S}(\omega)=\coth(\beta\omega/2)\tilde{\chi}^{\prime\prime}(\omega)$ is
well behaved. The long-time decoherence rate thus follows from the limit
$\gamma=\lim_{\omega\to 0}2\tilde{\chi}^{\prime\prime}(\omega)/\beta\omega$,
which can be taken directly in Eq. (S26) to obtain

$\gamma=\frac{2g^{2}T}{D}\int\frac{\mathrm{d}\mathbf{k}}{(2\pi)^{3}}\frac{\chi_{\mathbf{k}}|u_{\mathbf{k}}|^{2}}{k^{2}}\approx\frac{2\bar{g}^{2}T\chi_{0}}{D},$
(S28)

where the final approximation follows from the assumption that
$u_{\mathbf{k}}$ samples only long wavelengths, and we defined the
renormalised coupling
$\bar{g}^{2}=g^{2}\int\mathrm{d}\mathbf{k}|u_{\mathbf{k}}|^{2}/(8\pi^{3}k^{2}),$
which depends only on properties of the probe.

For $d<3$, the low-frequency limit of Eq. (S26) is ill-behaved in an infinite
system. For any finite system, however, the lower $k$-integration limit is cut
off by the exclusion of $\mathbf{k}=0$. In 1D this yields

$\gamma=\frac{2g^{2}T}{\pi
D}\int_{2\pi/L}^{\infty}\mathrm{d}k\,\frac{\chi_{k}|u_{k}|^{2}}{k^{2}}\sim L,$
(S29)

since for large $L$ the integral is dominated by its lower limit. A similar
argument in 2D yields $\gamma\sim\ln(L/\ell)$, where again the UV cutoff
appears explicitly for dimensional reasons. This diverging zero-frequency
noise in low dimensions originates from the long diffusive tail of the
correlation function in time. To see this explicitly, we take the Fourier
transform of Eq. (S26) to obtain, for $\tau>0$,

$\displaystyle\mathrm{i}\chi^{\prime\prime}(\tau)$
$\displaystyle=\frac{1}{2L^{d}}\sum_{\mathbf{k}\neq
0}\chi_{\mathbf{k}}|u_{\mathbf{k}}|^{2}Dk^{2}\mathrm{e}^{-Dk^{2}\tau}$ (S30)

Assuming as above that
$\chi_{\mathbf{k}}|u_{\mathbf{k}}|^{2}\approx\chi_{0}\mathrm{e}^{-\ell^{2}k^{2}/2}$,
we obtain $\chi(\tau)\propto(\ell^{2}+2D\tau)^{-(d/2+1)}$ in the thermodynamic
limit. The leading-order contribution to $C(\tau)$ arises from the
corresponding symmetrised noise, which for $\tau\to\infty$ can be approximated
from the low-frequency response as

$\displaystyle S(\tau)$
$\displaystyle=\int\frac{\mathrm{d}\omega}{2\pi}\coth(\beta\omega/2)\tilde{\chi}^{\prime\prime}(\omega)\mathrm{e}^{-\mathrm{i}\omega\tau}$
$\displaystyle\approx\int\frac{\mathrm{d}\omega}{2\pi}\frac{2\tilde{\chi}^{\prime\prime}(\omega)}{\beta\omega}\mathrm{e}^{-\mathrm{i}\omega\tau}$
$\displaystyle=\frac{1}{\beta L^{d}}\sum_{\mathbf{k}\neq
0}\chi_{\mathbf{k}}|u_{\mathbf{k}}|^{2}\mathrm{e}^{-Dk^{2}\tau},$ (S31)

from which we deduce $S(\tau)\sim\left(\ell^{2}+2D\tau\right)^{-d/2}.$