\section{Relation of QFI to multipartite entanglement \label{sec:Entanglement}}

The QFI has a deep connection with multipartite entanglement~\cite{Hyllus2012,Toth2012}. 
Consider a system of $N$ spins with length $S=\frac 1 2$, and a linear observable $\hat{\cal O} = \frac{1}{2}\sum_l {\mathbf{n}_l}\cdot \boldsymbol\sigma_l$, where $\boldsymbol\sigma_l=\left(\sigma_l^x,\sigma_l^y,\sigma_l^z\right)$ is a vector containing the Pauli matrices $\sigma_l^\alpha$ associated to spin $l$, and $\mathbf{n}_l=(n_l^x,n_l^y,n_l^z)$ is a unit vector on the Bloch sphere.  
Then the system hosts at least $m+1$-partite entanglement if the QFI associated to $\hat{\cal O}$ fulfils \cite{Hyllus2012}
\begin{equation}
	\label{eq:QFIandEntanglement}
	F_Q > \left\lfloor \frac{N}{m} \right\rfloor m^2 + \left(N - \left\lfloor \frac{N}{m} \right\rfloor m\right)^2 \,,
\end{equation}
where $\left\lfloor X \right\rfloor$ is the largest integer smaller than or equal to $X$. 
For $m$ a divisor of $N$, the condition \eqref{eq:QFIandEntanglement} attains the simple form \eqref{eq:QFIDensityandEntanglement} when expressed through the quantum Fisher information density $f_Q\equiv F_Q/N$. 
Note that typically response functions that measure linear operators contain a factor of $1/N$ relative to our definition, Eq.~\eqref{eq:susceptIntro}, so that they will observe directly $f_Q$ rather than $F_Q$. 

The proof of Ref.~\cite{Hyllus2012} for bounding the multipartite entanglement with Eq.~(\ref{eq:QFIandEntanglement}) can be directly translated to degrees of freedom other than spins $1/2$, as long as $\hat{\cal O}$ represents a sum of local operators with bounded spectrum \cite{Pezze2014}. If $h_{\mathrm{max}}$ and $h_{\mathrm{min}}$ denote the largest, respectively smallest, eigenvalue of $\hat{\cal O}$,   then the right hand side of the condition \eqref{eq:QFIandEntanglement} acquires the prefactor $(h_{\mathrm{max}}-h_{\mathrm{min}})^2$. 
Therefore, the QFI can also witness spatial entanglement in systems other than spins-$1/2$, such as larger spins or fermions.  
Note that the relation of the QFI to response functions, Eq.~\eqref{eq:FQImChiIntro}, is independent of any such microscopic details of the underlying quantum many-body system or the hermitian operator $\hat{\cal O}$, which may even be non-local or unbounded.\\ 


 
\section{Solvability of the considered models \label{sec:solvability}}

All models used in this article for illustrating the main concepts are exactly solvable. 
The one-dimensional Ising chain in a transverse field [Eq.~\eqref{eq:Ising} with $\mathcal{J}_{lj}=\delta_{j,l+1}$] can be mapped to a free-fermion problem \cite{Sachdev1999}. 
Dynamical susceptibilities can then be calculated via Wick's decomposition of expectation values \cite{Derzhko1997}, and the Pfaffians appearing in the resulting expressions can be evaluated efficiently using the algorithm described in Ref.~\cite{Wimmer2012}.

The infinite-range Ising Hamiltonian [Eq.~\eqref{eq:Ising} with $\mathcal{J}_{lj}=1/N$] commutes with both $\mathbf{S}^2 = (S^x)^2 + (S^y)^2 + (S^z)^2$ and $S^z$ where $S^\alpha=\sum_{l} \sigma_l^\alpha/2$, $\alpha=x,y,z$. As a consequence, the Hamiltonian decomposes into disconnected blocks when represented in the common eigenbasis of $\mathbf{S}^2$ and $S^z$. Each block grows linearly with particle number $N$ and can be diagonalized efficiently. Here, we consider the largest of these blocks, with dimension $N+1$. 
The infinite-range Ising model has a thermal phase transition for $\theta<\pi/4$, with critical temperature \cite{Das2006}   
\begin{equation}
	\frac{T_c}{J}=\frac{\sin(\theta)}{\log[  (1+\tan(\theta)) / (1-\tan(\theta)) ]}\,.
\end{equation}
In all models, we set $\hbar=k_B=a=1$, where $a$ is the lattice spacing.  




 
 
 
 
 


 





 
 
\section{Scaling and universality \label{sec:Scaling}}


In the main text, we have demonstrated a close relationship of the QFI to thermodynamic response functions. 
These response functions underlie universal scaling behaviour close to quantum critical points, and consequently directly imprint a universal scaling behaviour onto the QFI. 
Aspects of critical scaling of the QFI and related quantum metrics have also been discussed, e.g., in 
Refs.~\cite{Gu2010,Ma2009,Liu2013QuantumFisher,Wang2014,Zheng2015}. 
The purpose of the present supplementary section is to present a general, unified scaling theory for the QFI. 

\subsection{Basic principles}

Before discussing the scaling of the QFI in detail, let us briefly recall the basic principles of the scaling Ansatz \cite{bookCardy}.  
Here, we focus on lattice models with lattice spacing $a$ that have a continuous phase transition at zero temperature, \emph{viz} a quantum phase transition (QPT). 
Generalisation to continuum models is straightforward.
The dimensionality of the system is $d$, its linear size $L$, and the total particle number $N=L^d$ (where we set $a=1$). 
To simplify the discussion, we focus on the simplest scenario, QPTs obtained by tuning a single relevant field, say the magnetic field $h$, to the critical point $h_c$. In this scenario, any small variation of $\tilde{h} = |h-h_c|$ drives the system in the thermodynamic limit to its stable fixed point with zero correlation length. Nevertheless, in the scaling regime, i.e., sufficiently close to the phase transition, all low-energy physics  becomes universal. As a remarkable consequence of this universality, very diverse microscopic models develop the same collective emerging behaviour at low energy, i.e., they display the same power-law decay of correlations, which is described by a small set of universal exponents $\set{\Delta_{\alpha}}$.  


In finite systems, or systems at finite temperature, the power-law of the correlations is cut off by system size or the thermal correlation length. 
In such cases, one can extract the universal properties of the QPT by performing an appropriate analysis of scale transformations, i.e., transformations that modify the lattice spacing $a\to \lambda a$, where $\lambda>1$ (this scale transformation can be envisioned as grouping subsets of $\lambda$ spins together). 
Close enough to the critical point one can expand all local operators in the basis of scaling operators, i.e., those that transform under rescaling  with a simple power of $\lambda$, $\hat{\cal O}^{\alpha}\to \lambda^{-\Delta_{\alpha}}{\hat{\cal O}}^{\alpha}$, with scaling exponent $\Delta_{\alpha} \geq 0$. 
The  two-point correlation function of these operators decays as 
\begin{equation}
\label{eq:scalingOfO}
\langle{{\hat{\cal O}}^{\alpha}(0){\hat{\cal O}}^{\alpha}(r)}\rangle-\langle{{\hat{\cal O}}^{\alpha}(0)}\rangle\langle{{\hat{\cal O}}^{\alpha}(r)}\rangle\propto r^{-2\Delta_{\alpha}}\,.
\end{equation}
In particular, the fixed-point Hamiltonian is a scaling operator with scaling exponent $\Delta_{H}$.

Scaling operators that transform with $\Delta_{\alpha} > \Delta_{H}$  are called \emph{irrelevant} operators, since, when added as a weak perturbation, coarse graining decreases their importance relative to the original Hamiltonian. Those operators with  $\Delta_{\alpha} <  \Delta_{H}$, on the other hand,  are called  \emph{relevant}, since they become larger and larger under coarse-graining transformations, eventually driving the system to a different fixed point. 
In the \emph{scaling regime},  the regime where relevant perturbations are still small compared to the original Hamiltonian, one can extract the universal low-energy physics by studying the response of the system to scaling transformations. 

Both the inverse system size $L^{-1}$ and the finite temperature $T$ are  relevant perturbations. They increase under a rescaling as $L^{-1} \to \lambda L^{-1}$ and $T \to \lambda^z T$, where $z$ is the dynamical critical exponent. The scaling regime for a finite and cold system is thus defined by $L^{-1} \ll 1$ and $T \ll 1$.
We will now discuss the scaling behaviour of the QFI with respect to these perturbations, as well as to a deviation $\tilde{h}$ from criticality. 


\subsection{Scaling of the quantum Fisher information}

To construct a scaling theory for the QFI close to a quantum critical point, we start from Eq.~\eqref{eq:QFIpure}, which defines the QFI at zero temperature, and consider a local generator $\hat{\cal O}=\sum_l {\hat{\cal O}}^\alpha_l$. 
By inserting Eq.~\eqref{eq:scalingOfO} into Eq.~\eqref{eq:QFIpure}, one finds that $f_Q$ transforms under rescaling by a factor $\lambda$ as 
\begin{equation}
 f_Q(T',L'^{-1},\tilde{h}')=\lambda^{d-2\Delta_{\alpha}}{\phi}_Q(T \lambda^z, L^{-1}\lambda,  \lambda^{1/\nu}\tilde{h} )+C\,.
 \label{eq:sca_fisher}
\end{equation}
Here, $C$ is a constant that is related to the non-universal short-distance correlations, and $\nu$ is the critical exponent associated to the correlation length. 
According to Eq.~\eqref{eq:sca_fisher}, the QFI has the scaling exponent $\Delta_Q = d-2\Delta_{\alpha}$. 
The scaling function ${\phi}_Q(T \lambda^z, L^{-1}\lambda,  \lambda^{1/\nu}\tilde{h} )$ encodes the large-distance part of the correlation functions, which dominates the scaling if the condition $\Delta_{Q}>0$ is fulfilled. In the following discussion, we assume that this is the case, allowing us to neglect the non-universal constant $C$, leading to Eq.~\eqref{eq:sca_fisher_maintext}. 
Below we will encounter, in the guise of the hard-core boson chain, an example where $\Delta_Q\leq 0$. In such cases, derivatives of $f_Q$ still allow one to extract the purely universal contribution. 
Note that when $\Delta_{Q}=0$ there can be logarithmic corrections to Eq.~\eqref{eq:sca_fisher}.  

We can use the scaling form Eq.~\eqref{eq:sca_fisher} as long as all relevant perturbations are small, that is as long as we stay in the scaling regime. 
The most relevant perturbation will constitute the most significant breaking of scale invariance, a fact that enables us to extract valuable information about the involved critical exponents. 
For example, when both $L^{-1}$ and $T$ are sufficiently small, the deviation $\tilde{h}$ from the critical point represents the most relevant perturbation: Scale invariance will be broken at the scale $\lambda\sim\tilde{h}^{-\nu}$. This implies that we can trade $\lambda$ by $\tilde{h}^{-\nu}$, and varying $\tilde{h}$ in this regime has the same effect as a rescaling transformation. In this regime, we expect $f_Q$ therefore to behave as 
\begin{equation}
f_Q(T,L^{-1},\tilde{h})= \tilde{h}^{-\nu \Delta_{Q}}{\phi}_{Q1}\left(T \tilde{h}^{-\nu z}, L^{-1}\tilde{h}^{-\nu}\right)\,. 
\label{eq:fQscalingwithh}
\end{equation}
For sufficiently relevant operators, for which $\Delta_Q >0$, $f_Q$ thus diverges  as $ \tilde{h}^{-\nu \Delta_Q}$ when approaching the critical point.
In the above expression, ${\phi}_{Q1}(T \tilde{h}^{-\nu z}, L^{-1}\tilde{h}^{-\nu})$ is a scaling function, embodying the fact that the normalised quantum Fisher information density, $f_Q/\tilde{h}^{-\nu \Delta_{Q}}$, is a well defined function of the scaling variables $ T \tilde{h}^{-\nu z}$ and $ L^{-1}\tilde{h}^{-\nu}$ rather than of $L$, $T$, and $\tilde{h}$  separately. 

Upon decreasing $\tilde{h}$, the system will reach a crossover with the next relevant perturbation. 
If $L^{-1} < T^{1/z}$, the crossover will occur around $T^{1/z} \sim \tilde{h}^{\nu}$ and carry over into a thermally dominated regime, for which the scaling of the QFI is dictated by 
\begin{equation}
f_Q(T,L^{-1}, \tilde{h})= T^{-\Delta_Q/z}{\phi}_{Q2}\left(L^{-1} T^{-1/z},\tilde{h}T^{-1/(z\nu)}\right)\,.
\label{eq:fQscalingwithT}
\end{equation} 
This scaling behaviour characterises the robustness of critical entanglement against non-zero temperatures (see, e.g., Ref.~\cite{Escher2011} for a discussion of the robustness in presence of noise).

In the opposite case $L^{-1} > T^{1/z}$, upon reducing $\tilde{h}$ one enters a regime dominated by finite-size effects. Here, $f_Q$ behaves as 
\begin{equation}
\label{eq:fQscalingwithsize}
f_Q(T,L^{-1}, \tilde{h})= L^{\Delta_Q}{\phi}_{Q3}\left(T L^{ z},L^{1/\nu}\tilde{h}\right)\,.
\end{equation} 
In this regime, $f_Q$ saturates to a finite, $L$ dependent value, which diverges with system size as $L^{\Delta_Q}$. The crossover is located around $L^{-1} \sim \tilde{h}^{\nu}$ and should thus approach the critical value $\tilde h=0$ as $L^{-1/\nu}$. 
Another interesting scenario appears when $T$ and $1/L$ are comparable, giving rise to further complex crossover phenomena. This can be seen in Fig.~\ref{fig:FQIsing} of the main text, where we study the universal scaling behaviour of the QFI in the paradigmatic model of the quantum Ising chain. 

This scaling analysis immediately permits us to identify strongly entangled quantum critical points. These are those critical points where $\Delta_Q>0$, i.e., where at low temperatures and close to criticality, $f_Q$ diverges with system size. 
Additionally, the scaling also allows us to extract a length scale of the entanglement present in the system \cite{Verstraete2004c}. 
Combining Eq.~\eqref{eq:sca_fisher} and main-text inequality~\eqref{eq:QFIDensityandEntanglement}, one sees that the number of entangled particles scales as $N_Q\geq f_Q \sim \lambda^{d-2\Delta_{\alpha}}$, which defines a multipartite entanglement length scale 
\begin{equation}
	l_{\mathrm{ent}}\sim N_Q^{1/d}\gtrsim \lambda^{1-2\Delta_{\alpha}/d}\,.
\end{equation}
Here, we assumed an isotropic distribution of entanglement; for anisotropic situations, the longest distance along which particles are witnessed to be entangled is actually larger. 
This entanglement length is related to the correlation length $\xi$ of the generator $\hat{\cal O}$. Since at criticality the correlation length scales as $\xi\sim \lambda$, one has $l_{\mathrm{ent}}\sim \xi^{1-2\Delta_\alpha/d}$. 
For $2\Delta_{\alpha}/d<1$, the length scale over which particles are entangled with each other thus diverges with system size. 


 
 
\subsection{Universal scaling of in the Mott-insulator--superfluid transition}

The scaling behaviour of $f_Q$, Eq.~\eqref{eq:sca_fisher}, appears in all relevant operators, not only the order parameter. 
As an illustration for this fact, we choose a system that is relevant for experiments on cold gases, a system of ultracold bosonic atoms confined in a one-dimensional optical lattice. 
In the limit of strong on-site interactions, this system can be modelled by the Hamiltonian~\cite{Paredes2004} 
\begin{equation}
	\label{eq:Hhardcorebosons}
	H_{\mathrm{hcb}} = - J \sum_{l=1}^{N-1} \left(\hb_l^\dagger \hb_{l+1} + \mathrm{h.c.} \right) - 2\mu \sum_{l=1}^N \hb_l^\dagger\hb_l \,,
\end{equation} 
where the bosonic creation (annihilation) operators $\hb_l^\dagger$ ($\hb_l$) obey the hard-core constraint $\hb_l^\dagger \hb_l^\dagger=0$. 
This model has a QPT at $\mu_c/J=1$ from a critical superfluid phase to a Mott-insulator at uniform filling, with critical exponents $z=2$ and $\nu=1/2$ \cite{Sachdev1999}. 

The transition has several peculiarities. 
First, the exact ground state is fully separable at the QPT \cite{Giampaolo2008}. 
Second, one of its order parameters, the particle density $\hat{n}= \sum_l \hb_l^\dagger \hb_l/N$, commutes with the Hamiltonian. 
We can, however, use a modified observable that is still a relevant operator, the staggered density $\hat{\cal O}= \sum_l (-1)^l \hb_l^\dagger \hb_l$. 
In optical-lattice experiments, the dynamic susceptibility of such density fluctuations may be probed by Bragg spectroscopy \cite{Stoferle2004,Ernst2010}. 

The Hamiltonian Eq.~\eqref{eq:Hhardcorebosons} is exactly solvable with the same methods as for the Ising chain in a transverse field \cite{Sachdev1999,Derzhko2000}. 
Moreover, under periodic boundary conditions, the dynamic structure factor $S(q,\omega)$ associated to density fluctuations $\hat{\cal O}=\sum_l \ensuremath{{\rm{e}}}^{i q r_l} \hb_l^\dagger \hb_l$ at wave vector $q$ can be evaluated analytically \cite{Roux2013}. 
It is gapless at ${q}=0$ and acquires a maximum at ${q}=\pi$, where
\begin{equation}
	S(\pi,\omega) = 2\pi \sum_k \delta(\omega -\epsilon_k + \epsilon_{k+\pi}) f(\epsilon_{k+\pi}) [1-f(\epsilon_k)]\,,
\end{equation}
with $\epsilon_k=-2[J+\mu \cos(k)]$ the single-particle energy of the diagonalised Hamiltonian and $f(\epsilon_k)=1/[\exp(\epsilon_k/ T)+1]$ the Fermi--Dirac occupation. 
As a consequence of the gaplessness, the QFI evaluated at ${q=0}$ vanishes, while it assumes its maximum value at ${q=\pi}$, i.e., for the staggered density. 
Using the fluctuation-dissipation theorem and evaluating the integral in Eq.~\eqref{eq:FQImChiIntro}, the corresponding quantum Fisher information density becomes 
\begin{equation}
	f_Q = \frac 4 N \sum_k \tanh^2\left(\frac{\epsilon_k - \epsilon_{k+\pi}}{2 T}\right) f(\epsilon_{k+\pi}) [1-f(\epsilon_k)]\,.
\end{equation}
This quantity is plotted in Fig.~\ref{fig:FQXY}(a) as a function of temperature and chemical potential, for $N=1000$. 
The entanglement witnessed by the QFI extends up to $ T\approx 0.91 J$, improving over other witnesses showing entanglement up to $ T\approx 0.7 J$ \cite{Guehne2006,Hide2007}. 
In contrast to Ref.~\cite{Guehne2006}, however, the QFI for the staggered magnetisation does not witness genuine multipartiteness of entanglement, but saturates at the value $f_Q=2$. 
(Notably, the QFI for other generators such as $b_l+b_l^\dagger$ does show divergent behaviour.) 
Towards the phase transition, the QFI for the staggered density decreases, consistent with the presence of the factorisation point. 
Despite this absence of divergences, the phase transition imprints its character onto the universal behaviour of the QFI and its derivatives, as explained in what follows.  


\begin{figure}
\centering
\includegraphics[width = 0.995\columnwidth]{figs/hardcore_bosons}
\caption{
(color online) Quantum Fisher information density, $f_Q$, for the staggered density of hard-core bosons in a one-dimensional optical lattice. 
{\bf{(a)}} 
In part of the superfluid phase, $f_Q$ witnesses two-site entanglement. Data for $N=1000$. 
{\bf{(b)}} 
Inset: $\partial f_Q/{\partial \mu}$ shows a peak that, with decreasing $T$, diverges and moves towards the critical point. 
Main panel: The predicted scaling exponents produce a perfect data collapse for $\partial f_Q/{\partial \mu}$. 
Data for $N=100$ and, from blue to red, 8 different temperatures distributed on a logarithmic scale between $T=0.01-0.2 J$ (inset: 4 temperatures in $0.01-0.13 J$).  
}
\label{fig:FQXY}
\end{figure}


At $T=0$ and $N=\infty$, one can use Eq.~\eqref{eq:QFIpure} to derive a particularly simple analytical formula for the QFI of the staggered density, 
\begin{equation}
f_Q = \left\{ 
	\begin{array}{ll} 
	4 n & \mu \leq 0 \\
	4(1-n)& \mu > 0
	\end{array}
	\right.\,.
\end{equation}
Here, we defined the mean particle density $n\equiv \braket{\hat{n}}$, which is $n = 1 - \mathrm{acos}(\mu/J)/\pi$ for $|\mu|<J$ and $0$ else. 
The quantum Fisher information density attains its maximal value of $f_Q= 2$ at $\mu=0$ and decreases towards the phase-transition point, where it vanishes. 
The above analytic expression allows us to obtain the scaling dimension of $f_Q$ by expanding it close to the critical point, $f_Q \propto \tilde{h}^{1/2}$, with $\tilde{h}=1-\mu/J$. 
Comparison with the scaling form Eq.~\eqref{eq:fQscalingwithh} gives the exact scaling exponent $\Delta_Q=-1$. 
Since $\Delta_Q<0$, $f_Q$ of the staggered magnetisation displays no divergent behaviour. 

Nevertheless, as seen in the inset of Fig.~\ref{fig:FQXY}(b), the non-analyticity of $f_Q$ at $\tilde{h}=0$ leads to a divergence in its derivatives. 
Considering the temperature as the most relevant perturbation, Eq.~\eqref{eq:fQscalingwithT}, we obtain the universal behaviour
\begin{equation}
\frac{\partial f_Q}{\partial {\tilde{h}}}= T^{-\Delta_Q/z+1/(z\nu)}{\phi}_{Q3}(\tilde{h}T^{-1/(z\nu)})\,.
\end{equation} 
That is, $\partial f_Q/\partial {\tilde{h}}$ has a peak that diverges as $T^{-\Delta_Q/z+1/(z\nu)}$ and the position of which lies at a line $\tilde h \sim T^{1/(z \nu)}$. 
Using the analytically known exponents $-\Delta_Q/z+1/(z\nu)=1/2$ and $1/(z \nu) =1$, an excellent scaling collapse is achieved up to rather large temperatures [main panel of Fig.~\ref{fig:FQXY}(b)]. 
We also performed fits to the numerical data for the peak height and peak position. These yield exponents of ${0.49}$ and ${1.02}$, respectively---even without high-precision data, we are able to reproduce the exact scaling exponents within two percent accuracy. 

 
 
 





\subsection{Universal scaling in the infinite-range quantum Ising model}

\begin{figure}
\centering
\includegraphics[width = 0.7\columnwidth]{figs/scaling_longrange_Ising}
\caption{(color online) Universal scaling of the quantum Fisher information density, $f_Q=F_Q/N$, calculated for the order parameter in the quantum Ising model with infinite-range interactions. 
Finite-temperature and finite-size scaling at criticality $\theta=\pi/4$. Bullets are from lighter to darker blue for $N=100,\dots,3000$. 
The universal scaling laws produce an excellent data collapse. The divergence $f_Q\sim N^{1/3}$ is much weaker than in the one-dimensional quantum Ising model, reflecting the more mean-field-like character of the infinite-range model. Nevertheless, it still implies a diverging multipartiteness of entanglement.}   
\label{fig:scaling_longrangeIsing}
\end{figure}

In the main text, we considered the infinite-range Ising model in a transverse field, main-text Eq.~\eqref{eq:Ising} with $\mathcal{J}_{lj}=1/N$, to illustrate the absence of features in the QFI at thermal transitions (see Fig.~\ref{fig:FQ_LRIsing}). 
Here, we provide its scaling analysis in the vicinity of its quantum phase transition.  

The infinite-range Ising model describes the universal properties above the upper critical dimension, with the longitudinal magnetization $\hat{\cal O} = \sum_{l=1}^N \sigma_l^x$ as the order parameter. 
In Fig.~\ref{fig:scaling_longrangeIsing}, we plot the QFI as a function of temperature for different system sizes at the system's quantum critical point, $\tilde h=0$. Using that $\Delta_Q=1/3$ and $z=1/3$~\cite{Botet1983} one obtains the scaling of the quantum Fisher information density $ f_Q \sim N^{1/3} \phi(T N^{1/3})$. From a scaling collapse, shown in Fig.~\ref{fig:scaling_longrangeIsing}, we determined numerically that $\Delta_Q = 0.34$ and $z=0.33$ which is in very good agreement with the exact solution.

 
 
 
 
 
 \section{Classical contribution to QFI and behaviour at thermal phase transitions}

In the main text, we study how well thermal states that lie close in a phase diagram are distinguishable, focusing entirely on quantum effects. 
More precisely, we ask how fast the probability distribution of a thermal state $\rho=\sum_\lambda p_\lambda \ket{\lambda}\bra{\lambda}$, with $p_\lambda=\exp(-E_\lambda/T)/Z$, changes due to the unitary transformation with $\ensuremath{{\rm{e}}}^{i\delta\vartheta \hat{\cal O}}$. 
It is instructive to extend this analysis to any source of distinguishability, be it quantum or classical. 

Consider a canonical ensemble $\rho$ at temperature $T$ for a Hamiltonian $H$. We wish to compare this state with the equilibrium ensemble $\rho'$ at the same temperature but for $H' = H+{\cal J}{\delta\vartheta} \hat{\cal O}$, where $\delta\vartheta\ll1$ and where ${\cal J}$ represents a suitable energy scale. 
There are two sources that render $\rho$ and $\rho'$ distinguishable. 
The first is the one considered in the main text. It can be rephrased as an adiabatic evolution of the state by slowly changing the Hamiltonian from $H$ to $H'$. 
If the evolution is adiabatic and unitary, the populations $p_\lambda$ remain unchanged. Any difference between initial and final density operator is then solely due to a modification of the energy eigenstates. Taking a unit time step $t=1/{\cal J}$, the change amounts to $\ket{\lambda}\to\ensuremath{{\rm{e}}}^{i\delta\vartheta \hat{\cal O}}\ket{\lambda}$. This is formally equivalent to the phase estimation scenario considered in this article. 

However, one may extend the scenario by including a sufficiently long waiting time after the transfer, such that the populations at the new Hamiltonian $H'$ can rethermalise with an external bath. This change of populations provides a second source of distinguishability between $\rho$ and $\rho'$, which to leading order remains disparate from the one described above.  
We can incorporate this effect by considering the full QFI \cite{BraunsteinCaves1994} 
\begin{eqnarray} 
\tilde{F}_Q&=&F_C+F_Q\\
&=&\sum_\lambda \frac{(\partial p_\lambda/\partial\vartheta)^2}{p_\lambda} + 2\sum_{\lambda,\lambda^\prime} \frac{(p_\lambda-p_{\lambda^\prime})^2}{p_\lambda+p_{\lambda^\prime}} \left|\bra{\lambda}\hat{\cal O}\ket{\lambda^\prime}\right|^2 \nonumber \,. 
\end{eqnarray}
The second term, $F_Q$, the subject of this work, is the term that is relevant at a QPT at $T=0$, where the first term, $F_C$, vanishes identically. 

Both terms are related to different aspects of the susceptibility of $\rho$ towards the generator $\hat{\cal O}$. 
As demonstrated in Eq.~(2) of the main text, $F_Q$ is in one-to-one correspondence to the Kubo dynamic susceptibility at non-zero frequencies. In Kubo's linear response theory, a system is perturbed weakly in a time-dependent way, generating a unitary dynamics, but---to leading order---without changing the populations of energy levels \cite{JensenBookChapter3}. 

The change of probabilities, instead, is related to the static, isothermal susceptibility, as has already been pointed out within purely classical statistical physics \cite{ProkopenkoObst2011}. 
Here, we demonstrate this relation in a quantum statistical calculation. 
Consider a static perturbation, $H \to H' = H + {\delta\vartheta'}\hat{\cal O}$, with $\delta\vartheta'=T\delta\vartheta$ (i.e., we choose the energy scale ${\cal J}=T$). 
The expectation value of the observable $\hat{\cal O}$ is given by a derivative of the free entropy $\ln Z$, with $Z=\sum_\lambda \exp(-E_\lambda/T)$ the partition function, 
\begin{equation}
\braket{\hat{\cal O}}=T \frac{\partial \ln Z}{\partial \vartheta'} = - \sum_\lambda p_\lambda \frac{\partial E_\lambda}{ \partial\vartheta'}\,,
\end{equation}
and the corresponding static, isothermal susceptibility is 
\begin{equation}
\chi^T \equiv \frac{\partial\braket{\hat{\cal O}}}{\partial \vartheta'} = - \sum_\lambda  \frac{\partial p_\lambda}{ \partial\vartheta'} \frac{\partial E_\lambda}{ \partial\vartheta'}  - \sum_{\lambda}   p_\lambda \frac{\partial^2 E_\lambda}{ (\partial\vartheta')^2} \,.
\end{equation} 
The derivative of $p_\lambda$ reads  
\begin{equation}
\label{eq:derivativeofplambda}
\frac{\partial p_\lambda}{ \partial\vartheta'}=-\frac {p_\lambda} T \left( \frac{\partial E_\lambda}{ \partial\vartheta'}-\sum_{\lambda'} \frac{\partial E_{\lambda'}}{ \partial\vartheta'}\right)
\end{equation}
One can use perturbation theory to evaluate the change of energy levels following the static perturbation. Working in a basis where degenerate eigenlevels $\ket{\lambda}$ are prediagonalised with respect to the perturbation operator $\hat{\cal O}$, we obtain    
\begin{eqnarray}
\chi^T &=&  \frac 1 T \left( \sum_\lambda p_\lambda \hat{\cal O}_{\lambda,\lambda}^2 - \braket{\hat{\cal O}}^2\right) + \sum_{\lambda,\lambda'}^{E_\lambda\neq E_\lambda'} \frac{|\hat{\cal O}_{\lambda,\lambda'}|^2 (p_\lambda-p_\lambda')}{E_{\lambda'}-E_\lambda} \nonumber \\
&=& \chi_{\mathrm{el}} + \chi_{\mathrm{vV}}\,.
\end{eqnarray}
Here, we defined the abbreviation $\hat{\cal O}_{\lambda,\lambda'}\equiv \bra{\lambda}\hat{\cal O}\ket{\lambda'}$. 
The first term is the elastic contribution, or classical Curie contribution, arising from a change of populations. It vanishes in a pure state, such as the quantum ground state at $T=0$. 
The second term is the quantum-mechanical van-Vleck correction, arising from matrix elements between different energy levels. It is connected smoothly to the real part of the Kubo dynamic susceptibility at non-zero frequencies and vanishes in a purely classical theory. 
Importantly, the elastic peak is a singular contribution at exactly $\omega=0$, and does thus not influence $F_Q$. Since it is this elastic peak that diverges at thermal phase transitions, these do not appear as particular features in $F_Q$. 

We can connect the elastic peak to the classical contribution to the Fisher information \cite{ProkopenkoObst2011}. Equation~\eqref{eq:derivativeofplambda} yields $F_C=\sum_\lambda \frac{(\partial p_\lambda/\partial\vartheta)^2}{p_\lambda}= \chi_{\mathrm{el}}/T$. 
Similarly to $\chi^T$, $\tilde{F}_Q$ thus has two parts: $F_Q$, the subject of this article, related to entanglement and divergent at QPTs; 
and $F_C$, of non-unitary, classical origin and divergent at thermal phase transitions. 
Both are related to the order parameter's susceptibility. Remarkably, by choosing the susceptibility at vanishing or non-vanishing frequencies, we can extract the purely classical or the purely quantum-mechanical contribution to the quantum Fisher information. 




 
 \section{Sum rule}
 
Our central relation Eq.~(2) allows us to derive a Thomas-Reiche-Kuhn-type \cite{Thomas1925,Reiche1925,Kuhn1925} sum rule that may be useful to bound the QFI. 
Since $\tanh\left(\frac{\omega}{2T}\right)\leq\frac{\omega}{2T}$, we have 
\begin{equation}
	F_Q(T) \leq \frac{1}{T}  \frac{2}{\pi}\int_{0}^{\infty} \ensuremath{{\rm{d}}} \omega \,\omega\, {\chi}''(\omega,T)\,.
\end{equation} 
Due to the Kramers--Kronig relations, we can now take advantage of a superconvergence theorem \cite{Lucarini2003}. 
This theorem states that two complex functions $g$ and $f$, for which 
$g(x) =  P \int_0^\infty \ensuremath{{\rm{d}}} \tilde{x}\, \frac{f(\tilde{x})}{x^2-\tilde{x}^2}$ holds, are related via 
$
\int_0^\infty \ensuremath{{\rm{d}}} x \, f(x) = \lim_{x\to\infty} \left(x^2 g(x)\right)\,,  
$
with $P$ the principal part of the integral. 
By setting $f(x)=-\frac 2 \pi x \chi''(x)$ and $g(x)=\chi'(x)\equiv\Re(\chi(x))$, one obtains
\begin{equation}
	\label{eq:FQsumRule}
	F_Q(T) \leq - \frac 1 T \lim_{\omega\to\infty}(\omega^2 \chi'(\omega,T))\,, 
\end{equation}
a Thomas-Reiche-Kuhn-type sum rule that upper bounds the QFI.  
Here, as in Eq.~\eqref{eq:FQImChiIntro} of the main text, the limit to infinity is to be understood as a frequency range much larger than all energy scales of the considered degrees of freedom (but smaller than energy scales at which resonances to additional degrees of freedom appear, the entanglement between which, however, could be analysed analogously).   
It is interesting to note that, if we take the classical limit by sending $T\to \infty$, the bound \eqref{eq:FQsumRule} becomes an equality relating the QFI to the high-frequency limit of $\chi'(\omega)$.\\ 








 


