#  Local dynamics and thermal activation in the transverse-field Ising chain

Jiahao Yang  Tsung-Dao Lee Institute & School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China  Weishi Yuan  Department
of Physics and Astronomy, McMaster University, Hamilton, Ontario, L8S 4M1,
Canada  Takashi Imai  Department of Physics and Astronomy, McMaster
University, Hamilton, Ontario, L8S 4M1, Canada  Qimiao Si<EMAIL_ADDRESS>Department of Physics & Astronomy, Rice University, Houston, Texas 77005, USA
Jianda Wu<EMAIL_ADDRESS>Tsung-Dao Lee Institute & School of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Márton
Kormos<EMAIL_ADDRESS>MTA-BME Quantum-Dynamics and Correlations
Research Group, Eötvös Loránd Research Network (ELKH), Budapest University of
Technology and Economics, 1111 Budapest, Budafoki út 8, Hungary

###### Abstract

There has been considerable recent progress in identifying candidate materials
for the transverse-field Ising chain (TFIC), a paradigmatic model for quantum
criticality. Here, we study the local spin dynamical structure factor of
different spin components in the quantum disordered region of the TFIC. We
show that the low-frequency local dynamics of the spins in the Ising- and
transverse-field directions have strikingly distinctive temperature
dependencies. This leads to the thermal-activation gap for the secular term of
the NMR $1/T_{2}^{\prime}$ relaxation rate to be half of that for the
$1/T_{1}$ relaxation rate. Our findings reveal a new surprise in the nonzero-
temperature dynamics of the venerable TFIC model and uncover a means to evince
the material realization of the TFIC universality.

Introduction.  While classical matter freezes at zero temperature, quantum
many-body systems often display multiple ground states due to the competition
between different couplings. Upon a continuous transformation from one ground
state to another, quantum criticality develops. It is anchored by a quantum
critical point (QCP) at zero temperature, in contrast to a classical critical
point that appears at a thermally-induced phase transition. Quantum
criticality has emerged as a general organizing principle to understand many
of the richest phenomena that have been observed in quantum materials Special
Issue: Quantum Criticality and Novel Phases (2013); Paschen and Si (2021).
These include the cuprates Keimer  _et al._ (2015); Ramshaw  _et al._ (2015),
heavy fermion metals Kirchner  _et al._ (2020); Coleman and Schofield (2005);
Löhneysen  _et al._ (2007) and iron pnictides Dai  _et al._ (2009);
Licciardello  _et al._ (2019). One of the prominent features of quantum
criticality is that it mixes spatial and temporal fluctuations. While this
intermixing complicates the description of quantum criticality, it also
implies that dynamical properties can be used to characterize the nature of
quantum criticality. As another outstanding feature of quantum criticality,
approaching the QCP by a non-thermal control parameter ($g$) and by
temperature ($T$) represent two independent ways to examine its universal
behavior. As such, it is instructive to probe quantum criticality by analyzing
dynamical properties as a function of temperature, which can be conveniently
studied experimentally.

A paradigmatic model for continuous quantum phase transitions is the
transverse field Ising model in one spatial dimension Niemeijer (1967); Pfeuty
(1970); Sachdev (2011); Barouch and McCoy (1971); Suzuki (1971, 1976); Kopp
and Chakravarty (2005). It represents a prototype setting to explore the
properties of quantum criticality. Yet, in spite of its venerable status,
there is much about it that remains to be understood. Suitable materials to
study this model have only been emerging recently Coldea  _et al._ (2010);
Wang  _et al._ (2018); Cui  _et al._ (2019); Zhang  _et al._ (2020); Zou  _et
al._ (2021). They have allowed new experiments that are providing puzzling
experimental results, especially on dynamical properties at nonzero
temperatures. At the same time, calculating dynamical quantities near a QCP
are always challenging. That is also the case for the transverse field Ising
chain, notwithstanding the considerable efforts Niemeijer (1967); Pfeuty
(1970); Wu  _et al._ (2018). One of the particularly interesting quantities is
the local dynamical structural factor Wu  _et al._ (2014); unp (2013);
Steinberg  _et al._ (2019), which can be measured by the NMR longitudinal
relaxation rate $1/T_{1}$ Kinross  _et al._ (2014).

In this work, we study the local dynamics that can be measured by the NMR
transverse relaxation rate $1/T_{2}$. We consider the quantum disordered
region at $g>g_{c}$, and calculate $S^{xx}$ and $S^{zz}$ vs. $T$ at
$\omega\ll{k_{\text{B}}}T/\hbar$ [for notations about $x$ and $z$, see Eq. (1)
below.] We find that they have thermally-activated behavior with activation
gaps that differ by a factor of $2$. This leads to the surprising and
experimentally testable finding that the activation gap for the secular term
of $1/T_{2}$, named $1/T_{2}^{\prime}$ (Ref. [Jaccarino, 1965]), will be half
of that for $1/T_{1}$.

The model.  The Hamiltonian of the transverse field Ising chain is given by
Pfeuty (1970)

$H_{I}=-J\sum_{i=1}^{N}\left({{\sigma_{i}^{z}\sigma_{i+1}^{z}}+g{\sigma_{i}^{x}}}\right)\,,$
(1)

where $\sigma_{i}^{x}$ and $\sigma_{i}^{z}$ are Pauli matrices associated with
the spin components $S^{\mu}=\sigma^{\mu}/2,({\mu=x,y,z})\;$on site $i$, and
$g$ is the coupling with the transverse field. Below we shall refer to the $z$
(Ising) direction as longitudinal and to the $x$ direction as transverse. At
zero temperature, the system undergoes a quantum phase transition when the
transverse field is tuned across its QCP $g=g_{c}=1$. The Hamiltonian can be
conveniently converted to fermionic operators, $c_{i}$ and $c_{i}^{\dagger}$
through Jordan-Wigner transformation. After introducing a Bogoliubov rotation,
as described in Appendix A, the Hamiltonian takes the canonical form
$H_{\text{I}}=\sum_{k}\epsilon_{k}(\gamma_{k}^{\dagger}\gamma_{k}-\frac{1}{2})$
with single-particle energy dispersion

$\epsilon_{k}=2J\sqrt{1+g^{2}-2g\cos k}\,,$ (2)

where momentum $k$ is dimensionless in all calculations. At zero momentum, the
gap $\Delta=2J|g-g_{c}|$ which vanishes at $g_{\text{c}}.$

In the vicinity of the QCP, the low energy effective description of the system
is given by an Ising field theory obtained as the scaling limit of the lattice
Hamiltonian Eq. (18). In this limit the lattice spacing goes to zero, $a\to
0,$ while $J\to\infty$ and $g\to 1$ such that the energy gap and the “speed of
light” are kept fixed, $2J(1-g)=\Delta,$ $2Ja/\hbar=c.$ The resulting
Hamiltonian describes a relativistic field theory of free Majorana fermions
with mass $m=\Delta/c^{2}$

${\mathcal{H}}_{\text{I}}=\int\mathrm{d}x\left[\hbar
c\frac{i}{2}(\psi\frac{\partial\psi}{\partial
x}-\bar{\psi}\frac{\partial\bar{\psi}}{\partial
x})\pm\Delta\frac{i}{2}(\bar{\psi}\psi-\psi\bar{\psi})\right]\,,$ (3)

where the sign of second term is $+$ ($-$) for the paramagnetic
(ferromagnetic) regime, corresponding to $g>g_{c}$ ($g<g_{c}$) in the lattice
model. The field operators are related to the lattice operators as
$\stackrel{{\scriptstyle\scriptscriptstyle{\,(-)}}}{{\psi}}\\!\\!\\!\\!\\!(ja)=\frac{1}{\sqrt{2a}}\left(e^{\mp
i\pi/4}c_{j}+e^{\pm i\pi/4}c_{j}^{\dagger}\right).$ The spectrum of the Ising
field theory is built from multi-particle free fermion states with
relativistic dispersion for a single particle

$\varepsilon(p)=\sqrt{\Delta^{2}+p^{2}c^{2}}=\Delta\cosh\theta\,,$ (4)

where $p=\hbar k/a$ and $\theta$ is the relativistic rapidity parameter. The
scaling limit of the transverse magnetization operator can be written as

$\sigma_{j}^{x}=-2a\,\varepsilon(x=ja)\equiv-2a\,i\,\bar{\psi}(x)\psi(x)\,.$
(5)

The field theory operator $\sigma(x)$ corresponding to the order parameter is
non-local with respect to the fermions and cannot be expressed in a simple way
in terms of the latter. Assuming the standard field theory normalization, it
is related to the lattice magnetization as ($\hbar=c=1$) Pfeuty (1970)

$\sigma(x=ja)=2^{-1/24}e^{-1/8}\mathcal{A}^{3/2}a^{-1/8}\sigma^{z}_{j}=\bar{s}J^{1/8}\sigma^{z}_{j}\,,$
(6)

where $\mathcal{A}=1.2824271291\dots$ is Glaisher’s constant.

Dynamic structure factor.  We compute the local spin dynamical structure
factor (DSF) for frequencies $\omega\ll k_{\text{B}}T$ which are relevant to
the NMR relaxation rates. The DSF of the spin component $\alpha$ is given by

$S^{\alpha\alpha}(\omega,q)=\frac{-2}{1-e^{-\beta\omega}}\mathrm{Im}\,\chi^{\alpha\alpha}(\omega,q)=\\\
\sum_{l}\int_{-\infty}^{\infty}\mathrm{d}t\,e^{i\omega
t-iqla}\langle\sigma^{\alpha}_{l+1}(t)\sigma^{\alpha}_{1}(0)\rangle_{T}\,,$
(7)

where $\beta=(k_{\text{B}}T)^{-1}$, and $\chi^{\alpha\alpha}(\omega,q)$ is the
dynamical spin susceptibility at the transferred energy $\omega$ and momentum
$q$. In the field theory, we consider the continuum operators Eqs. (5), (6)
and the summation over lattice sites is replaced by an integral over $x$. The
local DSF is
$S^{\alpha\alpha}(\omega)=\int\frac{\mathrm{d}q}{2\pi}S^{\alpha\alpha}(\omega,q)$.
Both on the lattice and in the field theory, we proceed with a Lehmann
spectral representation. Using the field theory language, it reads

$S^{\alpha\alpha}(\omega,q)=\frac{1}{\mathcal{Z}}\sum_{r,s=0}^{\infty}C^{\alpha\alpha}_{r,s}(\omega,q)$
(8)

with $\alpha=x,y,z$ sig  and

$C^{\alpha\alpha}_{r,s}(\omega,q)=\int\frac{\mathrm{d}\theta_{1}\dots\mathrm{d}\theta_{r}}{(2\pi)^{r}r!}\int\frac{\mathrm{d}\theta^{\prime}_{1}\dots\mathrm{d}\theta^{\prime}_{s}}{(2\pi)^{s}s!}e^{-\beta
E_{r}}(2\pi)^{2}\\\
\delta(\omega+E_{r}-E_{s})\delta(q+P_{r}-P_{s})\,|\langle\theta_{1}\dots\theta_{r}|\sigma^{\alpha}(0,0)|\theta^{\prime}_{1}\dots\theta^{\prime}_{s}\rangle|^{2}\,,$
(9)

where the energy and momentum eigenvalues of the multiparticle states
$|\theta_{1},\dots,\theta_{n}\rangle$ are
$E_{n}=\Delta\sum_{i}^{n}\cosh\theta_{i}$ and
$P_{n}=(\Delta/c)\sum_{i}^{n}\sinh\theta_{i}$ ($c$ is re-introduced here and
following). As described in Appendix B, this series is a low-temperature
expansion. It has term-by-term divergences that can be regularized in a linked
cluster expansion Essler and Konik (2009); Pozsgay and Takács (2010),
${S}^{\alpha\alpha}(t,x)=\sum_{r=0,s=0}^{\infty}D^{\alpha\alpha}_{r,s}(t,x),$
where the finite terms $D^{\alpha\alpha}_{r,s}$ are certain linear
combinations of $C^{\alpha\alpha}_{r,s}$ and terms appearing in the expansion
of the partition function.

Local transverse DSF in the scaling limit.  We now turn to truncated form
factor series results of $S^{xx}(\omega)$ in the quantum disordered region
with $\omega\ll k_{\text{B}}T\ll\Delta$. The field theory operator
corresponding to the transverse magnetization,
$\varepsilon=i\,\bar{\psi}(x)\psi(x)$ [cf. Eq. (5)], is quadratic in the
fermionic operators, so it only has nonzero form factors between states that
either have an equal number of particles or the particle number difference is
2. In particular,

$F_{2}^{\varepsilon}(\theta|\theta^{\prime})\equiv\langle\theta|\varepsilon(0)|\theta^{\prime}\rangle=-i\frac{\Delta}{c}\cosh\left(\frac{\theta-\theta^{\prime}}{2}\right)\,,$
(10)

as can be obtained from the plane wave expansion of the fields [cf. Eq. (19)]
in a straightforward way.

The first term in the series is the vacuum contribution $C_{00}^{xx}=|\langle
0|\varepsilon|0\rangle|^{2}\delta(\omega)\delta(k)$, i.e. 0 particle - 0
particle (0p - 0p) contribution. In the field theory this is a divergent
contribution and requires renormalization (e.g. by normal ordering). However,
since we are interested in the finite (but small) $\omega$ domain, we ignore
this term. The terms $D_{0,s}^{xx}$ and $D_{s,0}^{xx}$ contribute to
frequencies $\omega\geq s\Delta$ outside of our domain of interest,
$0<\omega\ll\Delta.$ The first contributing term is thus

$C_{11}^{xx}(\omega,q)=\int\int\frac{\mathrm{d}\theta}{2\pi}\frac{\mathrm{d}\theta^{\prime}}{2\pi}|F_{2}^{\varepsilon}(\theta|\theta^{\prime})|^{2}e^{-\beta\Delta\cosh\theta}(2\pi)^{2}\\\
\delta\left[q+\frac{\Delta}{c}(\sinh\theta-\sinh\theta^{\prime})\right]\delta[\omega+\Delta(\cosh\theta-\cosh\theta^{\prime})]\,.$
(11)

Then after performing two integrals, the leading behavior of local transverse
DSF at $\omega\ll k_{\text{B}}T\ll\Delta$ region is (cf. Appendix C)

$C_{11}^{xx}(\omega)\approx-\frac{\Delta}{\pi
c^{2}}e^{-\frac{\Delta}{k_{\text{B}}T}}\left[\ln\left(\frac{\omega}{4k_{\text{B}}T}\right)-\frac{k_{\text{B}}T}{2\Delta}+\gamma_{\text{E}}\right]\,,$
(12)

where $\gamma_{E}$ is Euler’s constant. The main features of the result are
the $\sim e^{-\Delta/k_{\text{B}}T}$ temperature dependence and the
$\sim\ln(\omega/k_{\text{B}}T)$ logarithmic frequency dependence.

We have computed $C_{11}^{xx}$, but since in
$D_{11}^{xx}=C_{11}^{xx}-\mathcal{Z}_{1}C_{00}^{xx}$ the second term
contributes to the static $\sim\delta(\omega)$ response only, the obtained
result is equal to $D_{11}^{xx}$ for any nonzero $\omega$. The higher particle
number contributions $D_{r,s}^{xx}$ with $\mathrm{max}(r,s)\geq 2$ contain the
factor $e^{-\beta E_{r}}\delta(E_{r}-E_{s}+\omega)$, so they are exponentially
suppressed at low temperature $\beta\Delta\gg 1$ for frequencies
$\omega\ll\Delta$, in particular, $D_{r,s}^{xx}\sim
e^{-\mathrm{max}(r,s)\beta\Delta}$. So we find that the leading order
contribution to the local transverse DSF $S^{xx}(\omega)$ in the quantum
disordered region with $\omega\ll k_{\text{B}}T\ll\Delta$ is given by Eq.
(12). An alternative derivation of the same result is given in Appendix C.1.

Local transverse DSF in the spin chain.  The truncated form factor expansion
can also be applied in this case, with the obvious modifications: the
relativistic dispersion relation has to be replaced by the lattice one in Eq.
(2), and integrations over the lattice momenta runs from $-\pi/a$ to $\pi/a.$
As the calculation is analogous to the field theory calculation, we report it
in Appendix D. Exploiting the fact that the transverse magnetization has a
local expression in terms of the free fermionic variables, we can go beyond
the truncated form factor expansion and obtain exact transverse DSF
$S^{xx}(\omega,q)$, Eq. (67), at any temperature (cf. Appendix E for details).

In the quantum disordered region with $\omega\ll k_{\text{B}}T\ll\Delta$, the
asymptotic expression for $S^{xx}(\omega)$ reads (cf. Appendix E.1)

$S^{xx}(\omega)\approx\frac{\Delta}{\pi
J^{2}}e^{-\frac{\Delta}{k_{\text{B}}T}}\left[-\ln\left(\frac{\omega}{4k_{\text{B}}T}\right)+\frac{k_{\text{B}}T}{2\Delta}-\gamma_{E}\right]\,.$
(13)

Using $J=c/(2a)$ and recalling the rescaling factor $2a$ between the
$\sigma^{x}$ and the field theory operator $\varepsilon,$ we find perfect
agreement with the result Eq. (12).

Figure 1:  The local transverse DSF vs. frequency (blue dots) at fixed
$\Delta=0.1J$ and fixed $k_{\text{B}}T=0.01J$. The frequency dependence is
fitted by $10^{6}\times S^{xx}(\omega)=-2.08-1.39\ln(\omega/\Delta)$, shown by
the red solid line.

Fig. (1) shows the frequency dependence of $S^{xx}(\omega)$, evaluated from
the lattice model, and its fit in the $\omega\ll k_{\text{B}}T\ll\Delta$
region. With parameters $\Delta=0.1J$ and $k_{\text{B}}T=0.01J$, the exact
lattice result Eq. (67) gives the expected logarithmic divergence in $\omega$:
$10^{6}\times S^{xx}(\omega)=-2.08-1.39\ln(\omega/\Delta).$ This agrees well
with the asymptotic result Eq. (13), $10^{6}\times
S^{xx}(\omega)=-1.99-1.38\ln(\omega/\Delta)$.

Similarly, we also study the temperature dependence of $S^{xx}$ with
$\omega=10^{-4}J$ and $\Delta=0.1J$. The exact lattice result, shown in Fig.
(2), exhibits an exponential decay in temperature together with a logarithmic
correction in the prefactor. A fit to the lattice result yields
$S^{xx}(T)=e^{-\Delta/(k_{\text{B}}T)}\left[0.24+0.03\ln(k_{\text{B}}T/\Delta)\right]$
which conforms with the result obtained from the asymptotic expression Eq.
(13),
$S^{xx}(T)=e^{-\Delta/(k_{\text{B}}T)}\left[0.23+0.02(k_{\text{B}}T/\Delta)+0.03\ln(k_{\text{B}}T/\Delta)\right]$.
In the fitting, the $k_{\text{B}}T/\Delta$ term is not taken into account
since it is negligible compared with other terms in the considered region.

Local longitudinal DSF at low frequencies.  We next turn to the DSF of the
order parameter field, i.e. $S^{zz}(\omega)$. This operator is highly nonlocal
in terms of the Jordan–Wigner fermions prohibiting an exact calculation based
on free fermion techniques. However, one can still use the truncated form
factor series approach. The calculation of the form factors of $\sigma^{z}$
are far from being trivial, but are known exactly even on the finite spin
chain Iorgov  _et al._ (2011). We consider the paramagnetic phase in the
scaling limit, focusing on the DSF of the continuum spin operator $\sigma(x)$
in Eq. (6). We obtain, for $k_{\text{B}}{T}\ll\Delta$ (cf. Appendix F for
details),

$S^{zz}(\omega)\approx\frac{\bar{\sigma}^{2}}{\Delta}\frac{3\sqrt{3}}{2\pi}\left(\frac{k_{\text{B}}T}{\Delta}\right)^{2}e^{-\frac{2\Delta}{k_{\text{B}}T}}\,,$
(14)

where $\bar{\sigma}=\bar{s}(\Delta/c^{2})^{1/8}$ (cf. Eq.(6)). This result
agrees with the scaling limit of the corresponding result found in Refs.
Steinberg  _et al._ (2019); unp (2013). In Fig. (2) we show the temperature
dependence of $S^{zz}(\omega)$ Eq. (88) by numerical evaluation with fixed
frequency $\omega=0.1\Delta$ on a logarithmic scale. The fitted line,
$-0.13571-1.99478\Delta/(k_{\text{B}}T)-2.06216\ln[\Delta/(k_{\text{B}}T)]$,
is consistent with the prediction of Eq. (14),
$-0.189959-2\Delta/(k_{\text{B}}T)-2\ln[\Delta/(k_{\text{B}}T)]$.

Figure 2:  The contrast of the thermal-activation gaps between $S^{xx}$ (blue
dots) and $S^{zz}$ (black dots); the gap for $S^{xx}$ is half of that for
$S^{zz}$.

Contrasting thermal activations of $S^{xx}$ and $S^{zz}$.  Our key result is
the different exponential temperature dependence,
$e^{-\Delta/(k_{\text{B}}T)}$ vs. $e^{-2\Delta/(k_{\text{B}}T)}$, for $S^{xx}$
and $S^{zz}$, respectively. This is eventually a consequence of the
$\mathbb{Z}_{2}:\sigma^{z}\to-\sigma^{z}$ symmetry of the model. For large
transverse fields $g\gg 1,$ the ground state corresponds to spins pointing in
the $x$ direction, so intuitively, excitations in the paramagnetic phase can
be thought of as spin flips in the transverse direction generated by
$\sigma^{z}$. This means that $\sigma^{z}$ creates and destroys particles that
carry the same quantum numbers as the operator itself, i.e. they are odd under
spin reversal in the $z$ direction. As a consequence, the nonzero matrix
elements of $\sigma^{z}$ in the paramagnetic phase are between states with
different particle number parity Berg  _et al._ (1979); Iorgov  _et al._
(2011). For the transverse magnetization $\sigma^{x}$ it is the opposite: its
only nonvanishing matrix elements are between states of the same parity
because it is quadratic in the fermionic operators.

As discussed in the introduction, at low frequencies $\omega\ll\Delta$ only
those matrix elements can contribute in the Lehmann-representation where the
energies of the two states are close to each other due to the Dirac-delta
expressing energy conservation. This implies that matrix elements between the
ground state and the excited states do not contribute. Moreover, independently
of which state carries the Boltzmann factor, the contribution of the matrix
element will be $\sim e^{-n\Delta/(k_{\text{B}}T)}$ where $n$ is the larger of
the particle numbers in the two states. Together with the parity property of
$\sigma^{z}$ this implies that the leading temperature dependence of the
longitudinal DSF is $e^{-2\Delta/(k_{\text{B}}T)}$ coming from 1p - 2p
contributions while that of the transverse DSF is
$e^{-\Delta/(k_{\text{B}}T)}$ coming from 1p - 1p matrix elements.

We note in passing that in the quantum $E_{8}$ integrable field theory, the
scaling limit of the Ising spin chain perturbed by a longitudinal magnetic
field at its QCP ($g=g_{c}=1$), even the longitudinal DSF features a
$e^{-\Delta/(k_{\text{B}}T)}$ decay Wu  _et al._ (2014). This is a consequence
of the fact that $\sigma^{z}$ has nonzero 1p - 1p matrix element Delfino and
Mussardo (1995) as the $\mathbb{Z}_{2}$ symmetry is explicitly broken by the
perturbing field.

NMR relaxation rates.  The transverse field applied along the $x$-axis in the
TFIC serves as the applied static magnetic field in an NMR setup Kinross  _et
al._ (2014). The longitudinal NMR relaxation rate in this field geometry is
given by Moriya (1956); Jaccarino (1965)

$\frac{1}{T_{1}}\sim|A_{y}|^{2}S^{yy}(\omega_{n})+|A_{z}|^{2}S^{zz}(\omega_{n})\,.$
(15)

Here, $A_{j}$ ($j$ = $x$, $y$ and $z$) is the scalar hyperfine coupling, which
we take as constants for simplicity, and $\omega_{n}$ is the resonant
frequency ($\sim$MHz) of NMR measurements. Therefore, $1/T_{1}$ probes local
spin dynamics through $S^{zz}$ and $S^{yy}$ along the two orientations
orthogonal to the transverse-field direction $x$. In the TFIC,
$S^{yy}(\omega)=S^{zz}(\omega)\omega^{2}/[4(gJ)^{2}]$ Wu  _et al._ (2014) so
the contribution from $S^{yy}$ can be ignored in an NMR setting (which is at a
very low frequency). As such, we expect $1/T_{1}\sim
S^{zz}(\omega=\omega_{n})\sim e^{-2\Delta/(k_{\text{B}}T)}$, for $\omega\ll
k_{\text{B}}T\ll\Delta$ Kinross  _et al._ (2014).

To probe $S^{xx}(\omega)$, we consider the measurement of $1/T_{2}$ with the
same field geometry as $1/T_{1}$. In general Moriya (1956); Jaccarino (1965)

$\displaystyle\frac{1}{T_{2}}=$
$\displaystyle\frac{1}{T_{1}^{\prime}}+\frac{1}{T_{2}^{\prime}}\,,$
$\displaystyle\frac{1}{T_{1}^{\prime}}=$ $\displaystyle A\frac{1}{T_{1}}\,,$
(16) $\displaystyle\frac{1}{T_{2}^{\prime}}=$
$\displaystyle|A_{x}|^{2}S^{xx}(\omega\rightarrow 0)\,.$

Here, the non-secular contribution $1/T_{1}^{\prime}$ can be estimated from
the result of $1/T_{1}$ measurement using the prefactor A calculated
theoretically based on Bloch-Wangsmann-Redfield theorey, as outlined in
Appendix I. The secular term $1/T_{2}^{\prime}$ can then be determined from
Eq. (16) (first line); in practice, this separation is the easiest at
relatively low temperatures, where the non-secular term $1/T_{1}^{\prime}$ has
been suppressed by the larger gap $2\Delta$ of $S^{zz}$ but
$1/T_{2}^{\prime}\sim S^{xx}(\omega\rightarrow 0)\sim
e^{-\Delta/(k_{\text{B}}T)}$ still remains sizable.

Conclusion.  To conclude, we determined the leading behavior of the local
transverse DSF in the quantum disordered region of the TFIC at small transfer
energy with temperature much smaller than the gap. It is shown that when the
transfer energy is much smaller than the temperature the local transverse DSF
exhibits a logarithmic singularity. We found that the extracted thermal
activation gap from the local transverse DSF is half of that for the
longitudinal one, which can be attributed to the different parities of
$\sigma^{x}$ and $\sigma^{z}$ in the quantum disordered region of the TFIC.
This sharp contrast can be directly tested in a proper NMR setup through
$1/T_{1}$ and $1/T_{2}$ relaxation rate measurements. It is known that the
thermodynamical quantity, the Grüneisen ratio of the TFIC exhibits unique
quantum critical behavior Wu  _et al._ (2018). Our work reveals a sharp
contrast of the temperature dependence of the transverse and longitudinal
local DSFs of the TFIC, which represents a unique and surprising dynamical
feature for the TFIC. Accessing the universality of the TFIC is a crucial step
toward a realization of the exotic quantum phenomena such as the new particles
in the quantum $E_{8}$ integrable model Wu  _et al._ (2014); Zhang  _et al._
(2020); Wang  _et al._ (2021); Zou  _et al._ (2021). Our work implies that
combined measurements of the dynamical and thermodynamic quantities provide
telltale experimental signs for the class of TFIC universality.

## Acknowledgments

The work at Shanghai Jiao Tong University was partially supported by Natural
Science Foundation of Shanghai with Grant No. 20ZR1428400 and Shanghai Pujiang
Program with Grant No. 20PJ1408100 (J.Y, J.W), and by National Research,
Development and Innovation Office (NKFIH) under the research grant K-16 No.
119204 and by the Fund TKP2020 IES (Grant No. BME-IE-NAT), under the auspices
of the Ministry for Innovation and Technology. J.W. acknowledges additional
support from a Shanghai talent program. The work at McMaster was supported by
NSERC. The work at Rice has been supported in part by the NSF Grant No.
DMR-1920740 and the Robert A. Welch Foundation Grant No. C-1411. Q.S.
acknowledges the hospitality of the Aspen Center for Physics, which is
supported by the NSF grant No.PHY-1607611. M.K. was supported by the ÚNKP-20-5
new National Excellence Program of the Ministry for Innovation and Technology
from the source of the National Research, Development and Innovation Fund, and
acknowledges support by a Bolyai János grant of the HAS.

##  Appendix A Diagonalization of the lattice and field theory Hamiltonians

After Jordan–Wigner transformation, $\sigma^{x}_{i}=1-2c^{\dagger}_{i}c_{i}$
and
$\sigma_{i}^{z}=-\prod\limits_{j<i}{\left({1-2c_{j}^{\dagger}c_{j}}\right)}\left({c_{i}+c_{i}^{\dagger}}\right)$,
the TFIC Hamiltonian Eq. (1) becomes

$H_{\text{I}}=-J\sum_{i=1}^{N}\left[(c_{i}^{\dagger}c_{i+1}+c_{i}^{\dagger}c_{i+1}^{\dagger}+h.c.)+g(1-2c_{i}^{\dagger}c_{i})\right]$
(17)

in terms of the fermionic operators $c_{i},c_{i}^{\dagger}.$ After Fourier
transformation, $c_{j}=\frac{1}{\sqrt{N}}\sum_{k}c_{k}e^{ikj},$ the
Hamiltonian is diagonalized by a Bogoliubov rotation
$\gamma_{k}=u_{k}c_{k}-iv_{k}c_{-k}^{\dagger}\,,$ where
$u_{k}=\cos(\theta_{k}/2)\,,v_{k}=\sin(\theta_{k}/2)$ with the Bogoliubov
angle $\tan(\theta_{k})=\frac{\sin k}{g-\cos k}$. After these steps, we arrive
at

$H_{\text{I}}=\sum_{k}\epsilon_{k}\left(\gamma_{k}^{\dagger}\gamma_{k}-\frac{1}{2}\right)$
(18)

with single-particle energy dispersion $\epsilon_{k}=2J\sqrt{1+g^{2}-2g\cos
k}.$

The field theory Hamiltonian Eq. (3) can be diagonalized by the plane wave
expansion of the Majorana fields which in the paramagnetic phase reads
(setting $\hbar=1$)

$\displaystyle\psi(x,t)$
$\displaystyle=\phantom{i}\sqrt{\frac{mc}{2}}\int\frac{\mathrm{d}\theta}{2\pi}e^{-\theta/2}\left[\alpha
a(\theta)e^{ip_{\theta}x-i\varepsilon_{\theta}t}+\alpha^{*}a^{\dagger}(\theta)e^{-ip_{\theta}x+i\varepsilon_{\theta}t}\right]\,,$
(19a) $\displaystyle\bar{\psi}(x,t)$
$\displaystyle=i\sqrt{\frac{mc}{2}}\int\frac{\mathrm{d}\theta}{2\pi}\,e^{\theta/2}\;\,\left[\alpha
a(\theta)e^{ip_{\theta}x-i\varepsilon_{\theta}t}-\alpha^{*}a^{\dagger}(\theta)e^{-ip_{\theta}x+i\varepsilon_{\theta}t}\right]\,,$
(19b)

where $\alpha=e^{-i\pi/4},$ and $p_{\theta}=mc\sinh(\theta)$ and
$\varepsilon_{\theta}=mc^{2}\cosh(\theta)$ are the momentum and energy in
terms of the rapidity variable $\theta.$ The creation/annihilation operators
obey the algebra

$\\{a(\theta),a^{\dagger}(\theta^{\prime})\\}=2\pi\delta(\theta-\theta^{\prime})$
(20)

and diagonalize the Hamiltonian which becomes

$\mathcal{H}_{\text{I}}=\int\frac{\mathrm{d}\theta}{2\pi}\,a^{\dagger}(\theta)a(\theta)\,mc^{2}\cosh\theta\,.$
(21)

##  Appendix B Form-factor method and linked cluster expansion

Exploiting the local nature of the transverse magnetization in terms of the
fermions, the transverse DSF $S^{xx}(\omega,q)$ can be obtained exactly. This
is however not true for the longitudinal DSF $S^{zz}(\omega,q).$ Still, in
both cases one can give a systematic low temperature expansion Essler and
Konik (2009); Pozsgay and Takács (2010). As we are mainly interested in the
low-temperature NMR relaxation rates, we first discuss this more general
approach, applicable both in the spin chain and in the field theory. We shall
use the field theory notations but everything can be translated to the spin
chain in a straightforward manner.

Our starting point is the Lehmann spectral representation,

$S^{\alpha\alpha}(\omega,q)=\sum_{l}\int_{-\infty}^{\infty}\mathrm{d}t\,e^{i\omega
t-iql}\langle\sigma_{l+1}^{\alpha}(t)\sigma^{\alpha}_{1}(0)\rangle_{T}=\frac{1}{\mathcal{Z}}\sum_{n,m}e^{-\beta
E_{n}}(2\pi)^{2}\delta(\hbar\omega+E_{n}-E_{m})\delta(q+P_{n}-P_{m})|\langle
n|\sigma^{\alpha}_{1}(0)|m\rangle|^{2}\,.$ (22)

Using the multiparticle energy eigenstates111More generally, in interacting
field theories the basis of asymptotic scattering states is used.
$|\theta_{1},\dots,\theta_{n}\rangle$ in the Lehmann representation (22) leads
to

$S^{\alpha\alpha}(\omega,q)=\frac{1}{\mathcal{Z}}\sum_{r,s=0}^{\infty}C^{\alpha\alpha}_{r,s}(\omega,q)$
(23)

with $\alpha=x,y,z$ and

$C^{\alpha\alpha}_{r,s}(\omega,q)=\int\frac{\mathrm{d}\theta_{1}\dots\mathrm{d}\theta_{r}}{(2\pi)^{r}r!}\int\frac{\mathrm{d}\theta^{\prime}_{1}\dots\mathrm{d}\theta^{\prime}_{s}}{(2\pi)^{s}s!}e^{-\beta
E_{r}}\\\
(2\pi)^{2}\delta(\hbar\omega+E_{r}-E_{s})\delta(q+P_{r}-P_{s})\,|\langle\theta_{1}\dots\theta_{r}|\sigma^{\alpha}(0,0)|\theta^{\prime}_{1}\dots\theta^{\prime}_{s}\rangle|^{2}\,,$
(24)

where the energy and momentum eigenvalues are
$E_{n}=mc^{2}\sum_{i}^{n}\cosh\theta_{i}$ and
$P_{n}=mc\sum_{i}^{n}\sinh\theta_{i}.$

This series is a low-temperature expansion in the following sense. The term
$C_{r,s}$ (we omit the αα superscript) carries a factor $e^{-\beta
E_{r}}<e^{-r\beta m}$, thus the small expansion parameter is $e^{-\beta m}.$
Its dependence on $s$ is less obvious but note that thanks to the energy
conserving Dirac-delta, the energies of the two states with $r$ and $s$
particles are related. For a fixed frequency $\omega$ and at low temperature
one can truncate the double sum in both $r$ and $s$. The partition function
can also be written as $\mathcal{Z}=\sum_{n=0}^{\infty}\mathcal{Z}_{n}$ where
$\mathcal{Z}_{n}$ has a factor of $e^{-\beta E_{n}}$.

Now the main problem to be solved is the regularization of the singularities
present in the partition function and in the matrix elements (form factors) of
the operators in question. In infinite volume all $\mathcal{Z}_{n}$ are
singular due to the fact that they contain a scalar product of two momentum
eigenstates. Similarly, $C_{r,s}$ inherits the kinematical poles of the form
factors whenever two rapidities in the two sets coincide, i.e.
$\theta_{i}=\theta^{\prime}_{j}$ for some $i,j$. Since the structure factor is
a well-defined physical quantity, these singularities must cancel each other.
In order to make this manifest we reshuffle the infinite series in a linked
cluster expansion Essler and Konik (2009); Pozsgay and Takács (2010)

$\chi^{\alpha\alpha}(t,x)=\sum_{r=0,s=0}^{\infty}D^{\alpha\alpha}_{r,s}(t,x)\,,$
(25)

where the terms

$\displaystyle D_{0,s}$ $\displaystyle=C_{0,s}\,,$ (26) $\displaystyle
D_{1,s}$ $\displaystyle=C_{1,s}-\mathcal{Z}_{1}C_{0,s-1}\,,$ (27)
$\displaystyle D_{2,s}$
$\displaystyle=C_{2,s}-\mathcal{Z}_{1}C_{1,s-1}+(\mathcal{Z}_{1}^{2}-\mathcal{Z}_{2})C_{0,s-2}+\dots$
(28)

are supposed to be finite, and equivalent relations hold with the indices
interchanged. In order to obtain a finite result one needs to regularize the
divergencies either in a continuum scheme by adding infinitesimal imaginary
parts to the rapidities Essler and Konik (2009), or by going to a large but
finite volume $L$ that satisfies $1\ll mL\ll e^{m\beta}$, so that the density
of thermally excited particles is small Essler and Konik (2009); Pozsgay and
Takács (2010). The singularities manifest themselves as positive powers of
$L$, while the final result for the $D_{r,s}$ should be
$\sim\mathcal{O}(L^{0})$. All the $D_{1n}(t,x)$ terms are given for any
massive relativistic diagonal scattering theory in Ref. Pozsgay and Takács
(2010), while the general expression for $D_{22}(t,x)$ can be found in Ref.
Szécsényi and Takács (2012).

The resulting series still can have diverging terms, signalling that the zero
temperature quantity is already singular. This happens around the single
particle dispersion relation $\hbar\omega\sim\epsilon(k)$, where the zero
temperature DSF is proportional to an on-shell Dirac-delta which broadens at
non-zero temperatures Essler and Konik (2009). In these cases a resummation of
infinitely many terms is necessary in order to obtain the finite result.
However, if we are interested in the small-$\omega$ behavior in the disordered
phase of the Ising model then due to $\hbar\omega\ll mc^{2}$ we are far from
the mass shell. In this case the individual terms are not singular and the
truncated series should give a good approximation.

##  Appendix C Detailed field theory calculation of local transverse DSF

In this appendix we provide the details of the calculations leading to Eq.
(12). Our starting expression is Eq. (11) from the main text,

$C_{11}(\omega,q)=\int\int\frac{\mathrm{d}\theta}{2\pi}\frac{\mathrm{d}\theta^{\prime}}{2\pi}|F_{2}^{\varepsilon}(\theta|\theta^{\prime})|^{2}e^{-\beta
m\cosh\theta}(2\pi)^{2}\delta[q+m(\sinh\theta-\sinh\theta^{\prime})]\delta[\omega+m(\cosh\theta-\cosh\theta^{\prime})]\,.$
(29)

Both integrals can be performed using the Dirac-delta constraints. The system
of equations can be brought to the following form in terms of
$x=e^{\theta^{\prime}}$ and $y=e^{\theta}$

$\displaystyle x-y$ $\displaystyle=\tilde{\omega}+\tilde{p}\equiv A\,,$ (30a)
$\displaystyle\frac{1}{x}-\frac{1}{y}$
$\displaystyle=\tilde{\omega}-\tilde{p}\equiv B\,,$ (30b)

where $\tilde{\omega}=\omega/m,\tilde{q}=q/m.$ This leads to the two solutions
$\\{x_{+},y_{+}\\}$ and $\\{x_{-},y_{-}\\}$

$\displaystyle x_{\pm}$ $\displaystyle=\frac{AB\pm\sqrt{AB(AB-4)}}{2B}\,,$
(31a) $\displaystyle y_{\pm}$
$\displaystyle=\frac{-AB\pm\sqrt{AB(AB-4)}}{2B}\,.$ (31b)

These roots must be real and positive, so their product must be positive,
implying

$x_{+}y_{+}=x_{-}y_{-}=-x_{+}x_{-}=-\frac{A}{B}=\frac{\tilde{q}+\tilde{w}}{\tilde{q}-\tilde{w}}>0\,,$
(32)

so $\tilde{q}^{2}-\tilde{w}^{2}>0.$ Their sum is also positive, leading to

$x_{\pm}+y_{\pm}=\frac{\pm\sqrt{AB(AB-4)}}{B}\,,$ (33)

so only the $``+"$ solution is valid for $B=\tilde{\omega}-\tilde{q}>0$ and
only the $``-"$ solution is valid for $B=\tilde{\omega}-\tilde{q}<0.$ Note
that since $AB<0,$ the expression under the square root is always positive so
the reality condition is automatically satisfied. Summarizing,

$\displaystyle e^{\theta_{\pm}}$
$\displaystyle=\frac{q^{2}-\omega^{2}\pm\sqrt{(q^{2}-\omega^{2})(q^{2}-\omega^{2}+4m^{2})}}{2m(\omega-q)}\,,$
(34) $\displaystyle e^{\theta^{\prime}_{\pm}}$
$\displaystyle=\frac{-q^{2}+\omega^{2}\pm\sqrt{(q^{2}-\omega^{2})(q^{2}-\omega^{2}+4m^{2})}}{2m(\omega-q)}\,,$
(35)

where the $"+"$ roots are valid for $q<-|\omega|$ and the $"-"$ roots are
valid for $q>|\omega|.$ The Jacobian of the change of variables

$\\{m(\cosh\theta^{\prime}-\cosh\theta),m(\sinh\theta^{\prime}-\sinh\theta)\\}\to\\{\theta,\theta^{\prime}\\}$
(36)

is $m^{2}\sinh(\theta-\theta^{\prime}).$ Conveniently,

$\cosh(\theta-\theta^{\prime})=1+(q^{2}-\omega^{2})/(2m^{2})\,,\qquad\cosh\theta_{\pm}=\left(-\omega\mp
q\sqrt{1+4m^{2}/(q^{2}-\omega^{2})}\right)/(2m)\,.$ (37)

Using
$|F^{(\varepsilon)}(\theta|\theta^{\prime})|^{2}=m^{2}[1+\cosh(\theta-\theta^{\prime})]/2$
we obtain

$C_{11}(\omega,q)=e^{-m\beta\cosh\theta_{\pm}}\frac{q^{2}-\omega^{2}+4m^{2}}{2\sqrt{(q^{2}-\omega^{2})(q^{2}-\omega^{2}+4m^{2})}}=\frac{1}{2}e^{\beta\omega/2}\sqrt{1+\frac{4m^{2}}{q^{2}-\omega^{2}}}e^{-\frac{\beta}{2}|q|\sqrt{1+\frac{4m^{2}}{q^{2}-\omega^{2}}}}\,,$
(38)

where we used $\pm q=-|q|.$ The local dynamic structure factor is given by

$C_{11}(\omega)=\int\frac{\mathrm{d}q}{2\pi}C_{11}(\omega,q)=2\int_{|\omega|}^{\infty}\frac{\mathrm{d}q}{2\pi}C_{11}(\omega,q)\,.$
(39)

Note that as long as $\omega$ is non-zero, the exponential cuts off the
diverging prefactor and the integrand remains finite. This is no longer true
for $\omega=0$ which signals a logarithmic singularity. We can extract the
leading behavior in the small-$\omega$ (and low-$T$) limit as

$\begin{split}&C_{11}(\omega)\approx\frac{1}{2}e^{\beta\omega/2}\int_{|\omega|}^{\infty}\frac{\mathrm{d}q}{\pi}\left(\frac{2m}{q}+\frac{q}{4m}\right)e^{-\frac{\beta
m}{2}\left(\frac{q^{2}}{4m^{2}}+2+\frac{\omega^{2}}{q^{2}}\right)}\\\
&\approx\frac{1}{2}e^{\beta\omega/2}\int_{\omega}^{p_{m}}\frac{\mathrm{d}q}{\pi}\left(\frac{2m}{q}+\frac{q}{4m}\right)e^{-\frac{\beta
m}{2}\left(2+\frac{\omega^{2}}{q^{2}}\right)}+\frac{1}{2}e^{\beta\omega/2}\int_{p_{m}}^{\infty}\frac{\mathrm{d}q}{\pi}\left(\frac{2m}{q}+\frac{q}{4m}\right)e^{-\frac{\beta
m}{2}\left(\frac{q^{2}}{4m^{2}}+2\right)}\\\
&=\frac{e^{\beta\omega/2}e^{-m\beta}}{32\pi}\left\\{\frac{2\omega}{m}e^{-\frac{\beta
m(m+\omega)}{2m+\omega}}\left(e^{\frac{\beta m}{2}}(2m+\omega)-\omega
e^{\frac{\beta
m\omega}{4m+2\omega}}\right)+\left(16m-\beta\omega^{2}\right)\left[\Gamma\left(0,\frac{m\beta\omega}{4m+2\omega}\right)-\Gamma\left(0,\frac{m\beta}{2}\right)\right]\right\\}\\\
&\approx-\frac{m}{\pi}e^{\beta\omega/2}e^{-m\beta}\left[\ln\left(\frac{\beta\omega}{4}\right)-\frac{1}{2m\beta}+\gamma_{\text{E}}\right]\,,\end{split}$
(40)

where the incomplete gamma function
$\Gamma(a,z)=\int_{z}^{\infty}{{t^{a-1}}{e^{-t}}dt}$, Euler’s constant
$\gamma_{\text{E}}\approx 0.57712$, and $p_{m}=\sqrt{\omega(2m+\omega)}$ is
the extreme point of the exponent.

###  C.1 Alternative derivation of the local transverse DSF

Another way to obtain the result in Eq. (12) is to focus on the local DSF from
the start, defined as

$S^{\varepsilon}(\omega)=\int\mathrm{d}t\,e^{i\omega
t}\langle\varepsilon(0,t)\varepsilon(0,0)\rangle_{T}\,.$ (41)

This contains a disconnected piece proportional to $\delta(\omega).$ The first
nontrivial term in the expansion of the connected part is

$C_{11}(\omega)=\int\frac{\mathrm{d}\theta}{2\pi}\int\frac{\mathrm{d}\theta^{\prime}}{2\pi}|F_{2}^{\varepsilon}(\theta|\theta^{\prime})|^{2}e^{-\beta
m\cosh\theta}2\pi\delta[\omega+m(\cosh\theta-\cosh\theta^{\prime})]\,.$ (42)

Let us assume that $\omega>0$ so the energy conservation condition has two
real solutions $\theta_{+}$ and $\theta_{-}=-\theta_{+}$ for all $\theta,$
where we denote the positive solution by

$\theta_{+}=\mathrm{arccosh}(\cosh\theta+\omega/m)>0\,.$ (43)

Now
$\delta[\omega+m(\cosh\theta-\cosh\theta^{\prime})]=\frac{1}{m\sinh\theta_{+}}\left[\delta(\theta^{\prime}-\theta_{+})+\delta(\theta^{\prime}+\theta_{+})\right],$
so

$C_{11}(\omega)=m\int\frac{\mathrm{d}\theta}{2\pi}\frac{e^{-\beta
m\cosh\theta}}{\sqrt{(\cosh\theta+\omega/m)^{2}-1}}[1+\cosh\theta(\cosh\theta+\omega/m)]\,,$
(44)

where we used the identity
$\cosh^{2}[(\theta-\theta_{+})/2]+\cosh^{2}[(\theta+\theta_{+})/2]=1+\cosh\theta\cosh\theta_{+}.$
Introducing the shorthand notation $\tilde{\omega}=\omega/m$ and changing
integration variable to $u=\cosh\theta,$

$C_{11}(\omega)=\frac{m}{\pi}\int_{1}^{\infty}\mathrm{d}u\,e^{-m\beta
u}\frac{u^{2}+\tilde{\omega}u+1}{\sqrt{(u-1)(u+1)(u-1+\tilde{w})(u+1+\tilde{w})}}\,.$
(45)

The integral is singular for $\omega=0.$ To extract the small-$\omega$
behavior, we can approximate the integral by

$\begin{split}C_{11}(\omega)&\approx\frac{m}{\pi}\int_{1}^{\infty}\mathrm{d}u\,e^{-m\beta
u}\frac{u^{2}+1}{(u+1)\sqrt{(u-1)(u-1+\tilde{w})}}\approx\frac{m}{\pi}\int_{0}^{\infty}\mathrm{d}v\,e^{-m\beta(1+v)}\frac{1+v/2}{\sqrt{v(v+\tilde{w})}}\\\
&=\frac{m}{\pi}e^{-m\beta}\left[e^{\beta\omega/2}K_{0}\left(\frac{\beta\omega}{2}\right)+\frac{\sqrt{\pi}}{4m\beta}U(1/2,0,\beta\omega)\right]\,,\end{split}$
(46)

where $K_{0}(x)$ is the modified Bessel function of the second kind and
$U(a,b,z)$ is the confluent hypergeometric function. Expanding the result for
small $\omega$ we find

$C_{11}(\omega)\approx\frac{m}{\pi}e^{-m\beta}\left[-\ln\left(\frac{\beta
w}{4}\right)+\frac{1}{2m\beta}-\gamma_{\text{E}}\right]$ (47)

which agrees with the result in Eq. (40).

##  Appendix D  Calculations of local transverse DSF in the spin chain:
Truncated form factor series method

In this appendix we present the form factor calculation, analogous to that in
App. C, in the spin chain. The transverse magnetization on the lattice is
given in terms of the Jordan–Wigner fermions as

$\sigma_{j}^{x}=1-2c^{\dagger}_{j}c_{j}\,.$ (48)

The Hamiltonian is quadratic in these fermionic operators and it is
diagonalized by going to momentum space and performing the Bogoliubov
transformation

$\displaystyle\frac{1}{\sqrt{L}}\sum_{j=1}^{L}c_{j}e^{ikj}=c(k)$
$\displaystyle=\cos(\vartheta_{k}/2)\alpha(k)+i\sin(\vartheta_{k}/2)\alpha^{\dagger}(-k)\,,$
(49a)
$\displaystyle\frac{1}{\sqrt{L}}\sum_{j=1}^{L}c^{\dagger}_{j}e^{ikj}=c^{\dagger}(-k)$
$\displaystyle=i\sin(\vartheta_{k}/2)\alpha(k)+\cos(\vartheta_{k}/2)\alpha^{\dagger}(-k)$
(49b)

with

$e^{i\vartheta_{k}}=\frac{g-e^{ik}}{\sqrt{1+g^{2}-2g\cos
k}}=\frac{2J(g-e^{ik})}{\varepsilon(k)}\,.$ (50)

The ground state expectation value of $\sigma^{x}_{j}$ is then easily
calculated to be

$\left\langle\sigma^{x}\right\rangle\equiv\langle
0|\sigma^{x}_{j}|0\rangle=1-2\frac{1}{L}\sum_{n=1}^{L}\sin^{2}(\vartheta_{k}/2)=\frac{1}{L}\sum_{n=1}^{L}\left(1-2\sin^{2}(\vartheta_{k}/2)\right)=\frac{1}{L}\sum_{n=1}^{L}\cos(\vartheta_{k})\longrightarrow\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\cos(\vartheta_{k})\,.$
(51)

We can also compute the matrix elements

$\displaystyle\langle 0|\sigma^{x}_{j}|k,k^{\prime}\rangle$
$\displaystyle\equiv\langle
0|\sigma^{x}_{j}\,\alpha^{\dagger}_{k}\alpha^{\dagger}_{k^{\prime}}|0\rangle=\phantom{-}\frac{2i}{L}e^{-i(k+k^{\prime})j}\sin\left(\frac{\vartheta_{k}-\vartheta_{k^{\prime}}}{2}\right)\,,$
(52) $\displaystyle\langle k|\sigma^{x}_{j}|k^{\prime}\rangle$
$\displaystyle\equiv\langle
0|\alpha_{k}\,\sigma^{x}_{j}\,\alpha^{\dagger}_{k^{\prime}}|0\rangle=-\frac{2}{L}e^{i(k-k^{\prime})j}\cos\left(\frac{\vartheta_{k}+\vartheta_{k^{\prime}}}{2}\right)+\delta_{k,k^{\prime}}\left\langle\sigma^{x}\right\rangle\,.$
(53)

The first contributions to the local DSF in the Lehmann representation are

$\displaystyle C_{00}(\omega)$
$\displaystyle=2\pi\delta(\omega)\left\langle\sigma^{x}\right\rangle^{2}\,,$
(54) $\displaystyle C_{11}(\omega)$
$\displaystyle=\sum_{n,m}e^{-\beta\varepsilon_{n}}2\pi\delta(\omega+\varepsilon_{n}-\varepsilon_{m})\left[\delta(\omega)\left\langle\sigma^{x}\right\rangle^{2}\delta_{n,m}+\left(\frac{2}{L}\right)^{2}\cos^{2}\left(\frac{\vartheta_{k_{n}}-\vartheta_{k_{m}}}{2}\right)-\delta_{n,m}\frac{4}{L}\cos\left(\frac{\vartheta_{k_{n}}-\vartheta_{k_{m}}}{2}\right)\right]\,,$
(55)

where $\varepsilon_{n}=\varepsilon(k_{n})$ with $k_{n}=2\pi n/L,$
$n=-L/2+1,\dots L/2$ (Ramond sector). When calculating $D_{11}(\omega),$ the
first term in $C_{11}(\omega)$ cancels exactly
$\mathcal{Z}_{1}C_{00}(\omega)$. The remaining terms containing $\delta_{n,m}$
contribute to the $\omega=0$ response. Let us focus on the first nontrivial
term contributing at finite $\omega,$

$F(\omega)\equiv\left(\frac{2}{L}\right)^{2}\sum_{n,m}e^{-\beta\varepsilon_{n}}2\pi\delta(\omega+\varepsilon_{n}-\varepsilon_{m})\cos^{2}\left(\frac{\vartheta_{k_{n}}-\vartheta_{k_{m}}}{2}\right)\\\
\longrightarrow\;\;4\int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi}\int_{-\pi}^{\pi}\frac{\mathrm{d}k^{\prime}}{2\pi}e^{-\beta\varepsilon(k)}2\pi\delta(\omega+\varepsilon(k)-\varepsilon(k^{\prime}))\cos^{2}\left(\frac{\vartheta(k)-\vartheta(k^{\prime})}{2}\right)\,.$
(56)

The energy Dirac-delta can only be satisfied if

$\varepsilon_{\text{min}}=2J|1-g|\leq\varepsilon(k)+\omega\leq
2J(1+g)=\varepsilon_{\text{max}}\,,$ (57)

which implies

$\cos
k_{0}\equiv-1+\frac{1+g}{g}\tilde{\omega}-\frac{\tilde{\omega}^{2}}{2g}\leq\cos
k\leq 1+\frac{|1-g|}{g}\tilde{\omega}-\frac{\tilde{\omega}^{2}}{2g}\,,$ (58)

where $\tilde{\omega}\equiv\omega/(2J).$ Note that for $\omega<0$ the first
inequality is automatic, while for $0<\omega<2\Delta$ the second inequality is
automatically satisfied. We focus on the latter case from now on. Then $\cos
k^{\prime}$ is given by

$\cos k^{\prime}=\frac{1+g^{2}-(\omega+\varepsilon(k))^{2}/(4J^{2})}{2g}=\cos
k-\frac{\tilde{\omega}}{g}\sqrt{1+g^{2}-2g\cos
k}-\frac{\tilde{\omega}^{2}}{2g}$ (59)

which has two solutions which we denote by $k^{\prime}_{1}>0$ and
$-k^{\prime}_{1}.$ Changing variables and performing the integration over
$k^{\prime}$ we obtain

$F(\omega)=\frac{2}{\pi}\int_{-k_{0}}^{k_{0}}\mathrm{d}k\,e^{-\beta\varepsilon(k)}\frac{\varepsilon(k^{\prime}_{1}(k))}{4J^{2}g\sin(k^{\prime}_{1}(k))}\left[\cos^{2}\left(\frac{\vartheta(k)-\vartheta(k^{\prime}_{1}(k))}{2}\right)+\cos^{2}\left(\frac{\vartheta(k)+\vartheta(k^{\prime}_{1}(k))}{2}\right)\right]\,.$
(60)

Now we use the identity
$\cos^{2}[(x-y)/2]+\cos^{2}[(x+y)/2]=1+\cos(x)\cos(y),$ the explicit forms of
$\varepsilon(k)$ and

$\cos\vartheta(k)=\frac{g-\cos k}{\sqrt{1+g^{2}-2g\cos k}}$ (61)

together with $\varepsilon(k^{\prime}_{1})=\varepsilon(k)+\omega$ to arrive at

$F(\omega)=\frac{1}{Jh\pi}\int_{-k_{0}}^{k_{0}}\mathrm{d}k\,e^{-\beta\varepsilon(k)}\frac{2\tilde{\varepsilon}(k)(\tilde{\varepsilon}(k)+\tilde{\omega})-\sin^{2}k-\tilde{\omega}/g\cos
k(\tilde{\varepsilon}(k)+\tilde{\omega}/2)+\tilde{\omega}^{2}/2}{\tilde{\varepsilon}(k)\sqrt{1-(\cos
k-\tilde{\omega}/g\tilde{\varepsilon}(k)-\tilde{\omega}^{2}/(2g))^{2}}}\,,$
(62)

where $\tilde{\varepsilon}(k)=\varepsilon(k)/(2J).$ We would like to obtain an
approximate analytical result in the $\omega\to 0$ limit. Then we can keep
only the $O(\tilde{\omega}^{0})$ terms in the numerator. Changing variables to
$u=\cos k,$

$\begin{split}F(\omega)&\approx\frac{2}{Jg\pi}\int_{-1+(1+g)/g\,\tilde{\omega}-\tilde{\omega}^{2}/(2g)}^{1}\mathrm{d}u\,e^{-\beta
2Je(u)}\frac{2e(u)^{2}-(1-u^{2})}{\sqrt{1-u^{2}}e(u)\sqrt{1-(u-\tilde{\omega}e(u)/g-\tilde{w}^{2}/(2g))^{2}}}\\\
&\sim\frac{4}{\pi}\frac{|1-g|}{2Jg}e^{-\beta\Delta}\ln\tilde{\omega}=\frac{1}{\pi}\frac{\Delta}{J^{2}g}e^{-\beta\Delta}\ln\tilde{\omega}\end{split}$
(63)

with $e(u)=\sqrt{1+g^{2}-2gu}.$ It’s not surprising that the same logarithmic
divergent behavior appears as in Eq. (76) and Eq.(12). Comparing to the field
theory we have to keep in mind that between $\varepsilon$ and $\sigma^{x}$
there is a rescaling factor of $-2a=-\hbar c/J$ which leads to a perfect match
of the prefactors of $\ln\omega.$

##  Appendix E Transverse DSF $S^{xx}(q,\omega)$ in the spin chain

In this section we discuss the transverse DSF in the spin chain. In this
section, we report the calculation of the exact transverse DSF and specify its
low temperature and low frequency behavior at the end of the calculation. The
analogous derivation in the scaling limit can be found in Appendix G.

The transverse DSF follows by

$\begin{split}S^{xx}(\omega,q)&=\sum_{l=1}^{N}\int^{\infty}_{-\infty}\mathrm{d}te^{i\omega
t-iqla}\cdot\left[\left\langle\sigma^{x}_{l}(t)\sigma^{x}_{0}\right\rangle_{T}-\langle\sigma^{x}_{l}(t)\rangle_{T}\langle\sigma^{x}_{0}\rangle_{T}\right]\,,\end{split}$
(64)

where $\hbar\omega$ and $\hbar q$ are transferred energy and momentum,
respectively. In the following we set the lattice spacing as $a=1$ and
$\hbar=1$. Starting from $\sigma^{x}_{i}=1-2c^{\dagger}_{i}c_{i}$, going to
momentum space and performing the Bogoliubov rotation, the thermal expectation
values can be calculated in a straightforward way using the Hamiltonian Eq.
(18) and Wick’s theorem. The resulting expression for $S^{xx}(\omega,q)$ reads

$\begin{split}S^{xx}(\omega,q)&=\frac{4}{N}\sum_{k}\left[2\pi\delta(\omega-\epsilon_{k+q/2}+\epsilon_{k-q/2})\cdot(u_{k+q/2}^{2}u_{k-q/2}^{2}-u_{k+q/2}v_{k+q/2}u_{k-q/2}v_{k-q/2})\cdot(1-n_{k+q/2})n_{k-q/2}\right.\\\
&\left.\quad+2\pi\delta(\omega+\epsilon_{k+q/2}+\epsilon_{k-q/2})\cdot(u_{k-q/2}^{2}v_{k+q/2}^{2}+u_{k+q/2}v_{k+q/2}u_{k-q/2}v_{k-q/2})\cdot
n_{k+q/2}n_{k-q/2}\right.\\\
&\left.\quad+2\pi\delta(\omega-\epsilon_{k+q/2}-\epsilon_{k-q/2})\cdot(u_{k+q/2}^{2}v_{k-q/2}^{2}+u_{k+q/2}v_{k+q/2}u_{k-q/2}v_{k-q/2})\cdot(1-n_{k+q/2})(1-n_{k-q/2})\right.\\\
&\left.\quad+2\pi\delta(\omega+\epsilon_{k+q/2}-\epsilon_{k-q/2})\cdot(v_{k+q/2}^{2}v_{k-q/2}^{2}-u_{k+q/2}v_{k+q/2}u_{k-q/2}v_{k-q/2})\cdot
n_{k+q/2}(1-n_{k-q/2})\right]\\\
&=\int^{\pi}_{-\pi}\mathrm{d}k\left[2({1+f(k,q)})\cdot\delta(\omega-\epsilon_{k-q/2}+\epsilon_{k+q/2})\cdot(1-n_{k-q/2})n_{k+q/2}+({1-f(k,q)})\right.\\\
&\left.\cdot(\delta(\omega-\epsilon_{k-q/2}-\epsilon_{k+q/2})\cdot(1-n_{k-q/2})(1-n_{k+q/2})+\delta(\omega+\epsilon_{k-q/2}+\epsilon_{k+q/2})\cdot
n_{k-q/2}n_{k+q/2})\right]\,\end{split}$ (65)

with Fermi distribution function
$n_{k}=\left[e^{\epsilon_{k}/(k_{\text{B}}T)}+1\right]^{-1}$ and

$\begin{split}f(k,q)&=4J^{2}\left({[g-\cos(k-q/2)][g-\cos(k+q/2)]}-{\sin(k-q/2)\sin(k+q/2)}\right)/(\epsilon_{k-q/2}\epsilon_{k+q/2})\,.\end{split}$
(66)

Here we have taken advantage of particle and energy conservation during the
derivation, for example $\int^{\infty}_{-\infty}\mathrm{d}te^{i\omega
t}\langle
e^{iHt}\gamma_{k-q/2}\gamma_{k+q/2}^{\dagger}e^{-iHt}\gamma_{k+q/2}\gamma_{k-q/2}^{\dagger}\rangle_{T}=2\pi\delta(\omega+\epsilon_{k+q/2}-\epsilon_{k-q/2})(1-n_{{k-q/2}})n_{{k+q/2}}.$
Focusing on $\omega>0$ region, after integration, we obtain

$\begin{split}S^{xx}(\omega,q)&=\sum_{s=\pm}\left[2({1+f(k_{s},q)})\cdot\frac{(1-n_{k_{s}-q/2})n_{k_{s}+q/2}}{|D^{-+}_{k_{s}}(\omega,q)|}+({1-f(k_{s},q)})\cdot\frac{(1-n_{k_{s}-q/2})(1-n_{k_{s}+q/2})}{|D^{--}_{k_{s}}(\omega,q)|}\right]\,,\end{split}$
(67)

where $k_{\pm}$ are solutions of energy conservation constraint

$\cos(k_{\pm})=\frac{\omega^{2}\cos\left(\frac{q}{2}\right)}{16J^{2}g\sin^{2}\left(\frac{q}{2}\right)}\pm\left(\frac{\omega^{2}}{16J^{2}\sin^{2}\left(\frac{q}{2}\right)}-1\right)^{1/2}\left(\frac{\omega^{2}}{16J^{2}g^{2}\sin^{2}\left(\frac{q}{2}\right)}-1\right)^{1/2}\,,$
(68)

and the Jacobians are

$|D^{--}_{k_{s}}(\omega,q)|=\left|4J^{2}g\left(\frac{\sin\left(k_{s}+\frac{q}{2}\right)}{\epsilon_{k_{s}+q/2}}+\frac{\sin\left(k_{s}-\frac{q}{2}\right)}{\epsilon_{k_{s}-q/2}}\right)\right|,\qquad
s=\pm$ (69)
$|D^{-+}_{k_{s}}(\omega,q)|=\left|4J^{2}g\left(\frac{\sin\left(k_{s}+\frac{q}{2}\right)}{\epsilon_{k_{s}+q/2}}-\frac{\sin\left(k_{s}-\frac{q}{2}\right)}{\epsilon_{k_{s}-q/2}}\right)\right|,\qquad
s=\pm\,.$ (70)

Figure 3:  (a) The transverse DSF ${(J/\hbar)}S^{xx}(\omega,q)$ for $\omega>0$
at $\Delta=0.1J$, $k_{\text{B}}T=0.01J$. $q_{{}_{\text{IR}}}/a=\Delta/(\hbar
c)=(g-1)/a$ is the infrared wave number scale. The continuum above $2\Delta$
is contributed from the second term of Eq. (67). (b) The enlarged view of the
small rectangular region in (a) exhibits $S^{xx}(\omega,q)$ in the small
momentum and low energy region which comes from the first term in Eq. (67).

The transverse DSF, measured in units of $\hbar/J$, is plotted in Fig. 3 as a
function of dimensionless variables using the infrared frequency and wave
number scales $\Delta/\hbar$ and $\Delta a/(\hbar c)$, respectively. The upper
threshold in Fig. 3 (a) is given by

$\omega_{\text{up}}(q)=2\epsilon_{q/2\pm\pi}=4J\sqrt{g^{2}+2g\cos\left(q/2\right)+1}\,,$
(71)

while the lower thresholds are

$\displaystyle\omega_{\text{low1}}(q)$
$\displaystyle=2\epsilon_{q/2}=4J\sqrt{g^{2}-2g\cos\left(q/2\right)+1}\,,$
(72) $\displaystyle\omega_{\text{low2}}(q)$
$\displaystyle=\epsilon_{0}+\epsilon_{q}=2J\left(\sqrt{g^{2}-2g\cos
q+1}+|g-1|\right)\,.$ (73)

In Fig. 3 (b), the thresholds are $\omega_{\text{low'}}(q)=\epsilon_{q/2+k}$
and $\omega_{\text{up'}}(q)=\epsilon_{q/2-k}$ with
$k=\arccos\left[{\cos\left(q/2\right)}/{g}\right].$

###  E.1 Low temperature behavior of the local transverse DSF

We now focus on the local transverse DSF in the quantum disordered region with
gap much larger than the temperature, i.e. $k_{\text{B}}T\ll\Delta$. The
leading contribution in Eq. (67) is of the order
$e^{-\Delta/(k_{\text{B}}T)}$, and is given by the first term of Eq. (67).
Then the local transverse DSF follows immediately,

$\displaystyle S^{xx}(\omega)$ $\displaystyle=$
$\displaystyle\int^{\pi}_{-\pi}\frac{dq}{2\pi}S^{xx}(\omega,q)=\int^{\pi}_{q_{c}}\frac{dq}{\pi}S^{xx}(\omega,q)\approx\int_{q_{c}}^{\pi}\frac{dq}{\pi}\sum_{s=\pm}2[(1+f(k_{s},q)]\frac{(1-n_{k_{s}-q/2})n_{k_{s}+q/2}}{|D^{-+}_{k_{s}}(\omega,q)|}$
(74) $\displaystyle\approx$ $\displaystyle\int_{q_{c}}^{\pi}\frac{dq}{\pi}\
\left|\frac{\sqrt{1+g^{2}-2g\cos\frac{q}{2}}}{gJ\sin\frac{q}{2}}\right|\,e^{-\epsilon_{k_{-}+q/2}/(k_{\text{B}}T)}\,,$
(75)

where $q_{c}\approx{\omega}/{2J}$ is the lower bound obtained from
$\omega=\omega_{\text{up}}(q)$ at $\omega\to 0$ limit. The asymptotic behavior
of the integral in the $\omega\ll k_{\text{B}}T\ll\Delta$ regime is determined
in Appendix H with the result

$S^{xx}(\omega)\approx\frac{e^{-\frac{\Delta}{k_{\text{B}}T}}}{\pi}\left[-\frac{2\Delta}{J(\Delta+2J)}\left(\ln\left(\frac{\omega}{4k_{\text{B}}T}\right)+\gamma_{E}\right)+\frac{\Delta^{2}+12J^{2}+6\Delta
J}{6J^{2}(\Delta+2J)^{2}}k_{\text{B}}T\right]\,.$ (76)

The asymptotic result Eq. (76) shows that finite temperature local transverse
DSF diverges logarithmically as $\omega\to 0$. The energy conservation
constraints in Eq. (67) implies only the first term of Eq. (67) can contribute
to such low-energy behavior. Furthermore, the energy conservation leads to a
constraint for the phase space. After integration of Eq. (65), the constraint
gives rise to the $1/q$ dependency in the integrand of Eq. (75) for small $q$,
which is just the case for the lower bound dependent on the frequency. This
finally results in the logarithmic behavior. The temperature dependence shows
an exponential decay together with a logarithmic correction in the prefactor.
In the scaling limit we obtain

$S^{xx}(\omega)\approx\frac{\Delta}{\pi
J^{2}}e^{-\frac{\Delta}{k_{\text{B}}T}}\left[-\ln\left(\frac{\omega}{4k_{\text{B}}T}\right)+\frac{k_{\text{B}}T}{2\Delta}-\gamma_{E}\right]\,.$
(77)

Using $J=c/(2a)$ and recalling the rescaling factor $2a$ between the
$\sigma^{x}$ and the field theory operator $\varepsilon,$ we find perfect
agreement with the result (12).

In the $k_{\text{B}}T\ll\omega\ll\Delta$ region, we can simply approximate the
integral by the steepest descent method and obtain the asymptotic result

$S^{xx}(\omega)\approx\frac{1}{\pi}e^{-\Delta/(k_{\text{B}}T)}\sqrt{\frac{\pi
k_{\text{B}}T}{\omega}}\frac{2\Delta}{J(2J+\Delta)}\,.$ (78)

Figure 4:  The local transverse DSF as a function of temperature at fixed
$\Delta=0.1J$, and $\omega=10^{-4}J$. In the high $T$ region it clearly
deviates from the exponentially decaying behavior.

##  Appendix F Detailed field theory calculation of longitudinal DSF

In this section we turn to the DSF of the order parameter field. This operator
is highly nonlocal in terms of the Jordan–Wigner fermions prohibiting an exact
calculation based on free fermion techniques. However, one can still use the
truncated form factor series approach. This approach has been used in the
study of local spin DSF and NMR relaxation rate $1/T_{1}$ in Refs. Wu  _et
al._ (2014); Steinberg  _et al._ (2019). The calculation of the form factors
of $\sigma^{z}$ are far from being trivial, but are known exactly even on the
finite spin chain Iorgov  _et al._ (2011). Here we perform the calculation in
the paramagnetic phase in the scaling limit, focusing on the DSF of the
continuum spin operator $\sigma(x)$ in Eq. (6). In the disordered phase, the
$\sigma(x)$ operator creates and destroys particles, so its only non-zero
matrix elements are between states with particle numbers of different parity,
that is, the total number of particles in the two states must be odd. The
vacuum form factors are given by Berg  _et al._ (1979)

$F^{\sigma}_{n}(\theta_{1},\dots,\theta_{n})=\langle
0|\sigma|\theta_{1},\dots,\theta_{n}\rangle=\bar{\sigma}\sum_{i<j}\tanh\left(\frac{\theta_{i}-\theta_{j}}{2}\right)$
(79)

with $\bar{\sigma}=\bar{s}m^{1/8}$ where $\bar{s}$ is defined in Eq. (6) and
we work with $\hbar=c=1$ for the following field theory calculation. All other
matrix elements can be obtained by the crossing relation. For example,

$\langle\theta|\sigma|\theta_{1},\dots,\theta_{n}\rangle=\langle
0|\sigma|\theta+i\pi,\theta_{2},\dots,\theta_{n}\rangle\,,$ (80)

whenever $\theta\neq\theta_{j}$ which will be the case in our calculations.

Thus the first contributions to the DSF come from $D_{01}(t,x)=C_{01}(t,x)$
and $D_{10}(t,x)=C_{10}(t,x)$, yielding

$S^{zz}_{01}(\omega,q)=\bar{\sigma}^{2}\int\frac{\mathrm{d}\theta}{2\pi}(2\pi)^{2}\delta(\omega-m\cosh\theta)\delta(q-m\sinh\theta)\,,$
(81)

and $S_{10}(\omega,q)=e^{\beta\omega}S_{01}(-\omega,-q).$ Due to the energy
conserving Dirac-delta, both $S_{10}$ and $S_{01}$ are zero for $\omega<m$. It
is clear that all $D_{0,s}$ and $D_{r,0}$ will also give zero contribution,
which reflects the fact that the zero temperature result is identically zero.

Energy conservation at small frequencies also leads to a great simplification
in higher orders, similarly to the case of the transverse magnetization in the
previous sections. Because of the Dirac-delta and $\omega\approx 0$ the two
states in each matrix element must have almost equal energies, so
$S_{r,s}(\omega\approx 0,q)\sim e^{-\mathrm{max}(r,s)\cdot\beta m}$. This
implies that the classification in terms of orders of $e^{-\beta m}$ is
simplified, because in every order there is only a finite number of terms. For
instance, in the second order one has $S_{12}+S_{21}$, in the third order
$S_{23}+S_{32}$, in the fourth $S_{14}+S_{41}+S_{34}+S_{43}$, and so on.

Thus up to the second order one needs only two terms, $S_{12}$ and $S_{21}$.
We use the expression for $D_{12}$ given in Ref. Essler and Konik (2009) that
can be shown to be equivalent to the more general formula in Ref. Pozsgay and
Takács (2010),

$\begin{split}D_{12}&=\frac{1}{2}\int\frac{\mathrm{d}\theta}{2\pi}\int_{C_{+}}\int_{C_{-}}\frac{\mathrm{d}\theta_{1}}{2\pi}\frac{\mathrm{d}\theta_{2}}{2\pi}|F_{3}^{\sigma}(\theta+i\pi,\theta_{1},\theta_{2})|^{2}e^{-\beta
m\cosh\theta}\\\ &\times
2\pi\delta[q+m(\sinh\theta-\sinh\theta_{1}-\sinh\theta_{2})]2\pi\delta[\omega+m(\cosh\theta-\cosh\theta_{1}-\cosh\theta_{2})]\\\
&-\bar{\sigma}^{2}\int\frac{\mathrm{d}\theta}{2\pi}e^{-m\beta\cosh\theta}2\pi\delta(\omega-m\cosh\theta)2\pi\delta(q-m\sinh\theta)\,,\end{split}$
(82)

where the contours $C_{\pm}$ are running above and below the real axis,
respectively, to avoid the kinematical poles of the form factors. But
$\theta=\theta_{i}$ $(i=1,2)$ is impossible for $\omega<m$, so the integrals
avoid the poles even for real rapidities and there is no need to shift the
contours off the real axis. The last term is proportional to $D_{01}$ so it
does not contribute for $\omega<m$ and we are left with

$\begin{split}S^{zz}_{12}(\omega,q)&=\frac{1}{2}\int\frac{\mathrm{d}\theta}{2\pi}\int\frac{\mathrm{d}\theta_{1}}{2\pi}\int\frac{\mathrm{d}\theta_{2}}{2\pi}|F_{3}^{\sigma^{z}}(\theta+i\pi,\theta_{1},\theta_{2})|^{2}e^{-\beta
m\cosh\theta}\\\ &\times
2\pi\delta[q-m(\sinh\theta_{1}+\sinh\theta_{2}-\sinh\theta)]2\pi\delta[\omega-m(\cosh\theta_{1}+\cosh\theta_{2}-\cosh\theta)]\,.\end{split}$
(83)

Exploiting the Dirac-deltas we perform the integrals over rapidities
$\theta_{1,2}.$ The Jacobian of the transformation is
$m^{2}\sinh(\theta_{1}-\theta_{2}).$ The set of two constraint equations
coming from the Dirac-deltas has two solutions,
$\\{\theta_{1},\theta_{2}\\}=\\{\theta_{+},\theta_{-}\\}$ and
$\\{\theta_{1},\theta_{2}\\}=\\{\theta_{-},\theta_{+}\\}$, where

$\theta_{\pm}=\log\left[\frac{1}{2}\left(\omega_{\theta}+q_{\theta}\pm\frac{\sqrt{(\omega_{\theta}^{2}-p_{\theta}^{2})(\omega_{\theta}^{2}-q_{\theta}^{2}-4)}}{\omega_{\theta}-q_{\theta}}\right)\right]$
(84)

with $q_{\theta}\equiv
q/m+\sinh\theta,\omega_{\theta}\equiv\omega/m+\cosh\theta.$ The rapidities
$\theta_{1,2}$ must be real which gives restrictions on the remaining rapidity
$\theta.$ The reality condition of $\theta_{1,2}$ is equivalent to the
condition that $e^{\theta_{1}}+e^{\theta_{2}}$ and
$e^{\theta_{1}}\,e^{\theta_{2}}$ must be positive, which gives
$\omega_{\theta}+q_{\theta}>0,\,\omega_{\theta}-q_{\theta}>0.$ Moreover, the
combination under the square root must also be positive,
$\omega_{\theta}^{2}-q_{\theta}^{2}>4.$ One of the first two conditions, e.g.
$\omega_{\theta}-q_{\theta}>0$ can then be dropped which leaves us with two
conditions. The solution of $\omega_{\theta}^{2}-q_{\theta}^{2}>4$ for
$|\omega|<m$ is the following:

$\displaystyle\theta$ $\displaystyle<\theta^{(-)}$
$\displaystyle\qquad|\omega|$ $\displaystyle<q\,,$ (85a) $\displaystyle\theta$
$\displaystyle>\theta^{(+)}$ $\displaystyle\qquad-|\omega|$
$\displaystyle>q\,,$ (85b) $\displaystyle\theta$
$\displaystyle<\theta^{(-)}\;\text{ or }\;\theta>\theta^{(+)}$
$\displaystyle\qquad-\omega$ $\displaystyle<q<\omega\,,$ (85c)

and for $\omega<q<-\omega$ ($\omega<0$) there is no solution. Here

$\theta^{(\pm)}=\log\left[\frac{q^{2}-\omega^{2}+3m^{2}\pm\sqrt{(q^{2}-\omega^{2}+m^{2})(q^{2}-\omega^{2}+9m^{2})}}{2m(\omega-q)}\right].$
(86)

It turns out that the other condition, $\omega_{\theta}+q_{\theta}>0$, is
automatically satisfied, so Eqs. (85) give the integration domain of $\theta$
in the various cases depending on $\omega$ and $q.$ Thus we find

$S^{zz}_{12}(\omega,q)=2\int_{D}\frac{\mathrm{d}\theta}{2\pi}\frac{e^{-\beta
m\cosh\theta}|F_{3}^{\sigma}[\theta+i\pi,\theta_{+}(\theta),\theta_{-}(\theta)]|^{2}}{m^{2}\sqrt{(\omega_{\theta}^{2}-q_{\theta}^{2})(\omega_{\theta}^{2}-q_{\theta}^{2}-4)}}\,,$
(87)

where $D$ denotes the domain given in Eqs. (85) and we used that
$\sinh(\theta_{1}-\theta_{2})=\sqrt{(\omega_{\theta}^{2}-q_{\theta}^{2})(\omega_{\theta}^{2}-q_{\theta}^{2}-4)}/2.$
It is easy to see that $S_{21}(\omega,q)=e^{\beta\omega}S_{12}(-\omega,q),$ so
we have the total leading $\mathcal{O}(e^{-2m\beta})$ contribution to the DSF.
The result $S^{zz}_{12}(\omega,q)+S^{zz}_{21}(\omega,q)$ is plotted in Fig. 5.

Figure 5:  The leading $\mathcal{O}[e^{-2\Delta/(k_{\text{B}}T)}]$
contribution to the longitudinal dynamic structure factor in the scaling
limit, $\Delta^{2}/(\hbar^{2}c\,\bar{\sigma}^{2})S^{zz},$ at
$k_{\text{B}}T=\Delta/10.$

The corresponding local DSF reads

$\begin{split}S^{zz}_{12}(\omega)&=\frac{1}{2}\\!\int\\!\frac{\mathrm{d}\theta}{2\pi}\\!\int\\!\frac{\mathrm{d}\theta_{1}}{2\pi}\\!\int\\!\frac{\mathrm{d}\theta_{2}}{2\pi}e^{-\beta
m\cosh\theta}|F_{3}^{z}(\theta+i\pi,\theta_{1},\theta_{2})|^{2}2\pi\delta[\omega-m(\cosh\theta_{1}+\cosh\theta_{2}-\cosh\theta)]\\\
&=\frac{1}{2}\int\frac{\mathrm{d}\theta_{1}}{2\pi}\int\frac{\mathrm{d}\theta_{2}}{2\pi}\frac{e^{\beta\omega-\beta
m(\cosh\theta_{1}+\cosh\theta_{2})}}{\sqrt{(m\cosh\theta_{1}+m\cosh\theta_{2}-\omega)^{2}-m^{2}}}\left[|F_{3}^{\sigma}(\theta_{0}+i\pi,\theta_{1},\theta_{2})|^{2}+\\{\theta_{0}\to-\theta_{0}\\}\right]\,,\end{split}$
(88)

where
$\theta_{0}=\mathrm{arccosh}(\cosh\theta_{1}+\cosh\theta_{2}-\omega/m)>0.$ At
the second order we also need $S_{21}(\omega)=e^{\beta\omega}S_{12}(-\omega).$

We can give approximate expressions for $S_{12}$ and $S_{21}$. For $\beta m\gg
1$ only a small region around the origin in the $(\theta_{1},\theta_{2})$
plane contributes, so we can expand both the exponent and the rest of the
integrand to second order using the explicit form factors. Performing the
resulting Gaussian integrals, and expanding the result in $\omega/m$ (with
$T/m\ll 1$) we obtain the result in Eq. (14):

$\begin{split}S^{zz}(\omega)\approx
S^{zz}_{12}(\omega)+S^{zz}_{21}(\omega)&\approx\frac{\bar{\sigma}^{2}}{m}\frac{3\sqrt{3}}{4\pi}\left(\frac{1}{m\beta}\right)^{2}e^{-2m\beta}\left[e^{\beta\omega}\left(1+2\frac{\omega}{m}\right)+\left(1-2\frac{\omega}{m}\right)\right]\\\
&\approx\frac{\bar{\sigma}^{2}}{m}\frac{6\sqrt{3}}{4\pi}\left(\frac{T}{m}\right)^{2}e^{-2\frac{m}{T}}=\frac{\bar{\sigma}^{2}}{\Delta}\frac{3\sqrt{3}}{2\pi}\left(\frac{k_{\text{B}}T}{\Delta}\right)^{2}e^{-2\frac{\Delta}{k_{\text{B}}T}}\,.\end{split}$
(89)

The correction terms to this result are the third order $S_{23}+S_{32}$.
However, these terms contain singularities for which the regularization has
not yet been worked out explicitly. But as we discussed, unlike the case of
the broadening of the Dirac-delta in the zero temperature DSF, there is no
physical reason why unexpected singularities should show up in the higher
terms, thus we stop at the second order.

##  Appendix G Exact transverse DSF in the field theory

Using the plane wave expansion Eq. (19) and $\varepsilon=i\bar{\psi}\psi,$ the
connected correlation function

$C^{\varepsilon}(x,t)=\langle\varepsilon(x,t)\varepsilon(0,0)\rangle_{T}-\langle\varepsilon(x,t)\rangle\langle\varepsilon(0,0)\rangle_{T}$
(90)

can be written as a four-fold rapidity integral of a linear combination of
thermal expectation values of products of four creation/annihilation
operators. Using the thermal Wick’s theorem,

$\langle
a^{\pm}(\theta)b_{1}b_{2}b_{3}\rangle_{T}=f^{\pm}(\theta)\Big{(}\\{a,b_{1}\\}\langle
b_{2}b_{3}\rangle_{T}-\\{a,b_{2}\\}\langle
b_{1}b_{3}\rangle_{T}+\\{a,b_{3}\\}\langle b_{1}b_{2}\rangle_{T}\Big{)}\,,$
(91)

where $a^{+}=a^{\dagger},a^{-}=a,$ and
$f^{+}(\theta)=(1+e^{\beta\varepsilon_{\theta}})^{-1}=f(\theta),$
$f^{-}(\theta)=(1+e^{-\beta\varepsilon_{\theta}})^{-1}=1-f(\theta),$ one
arrives at

$\begin{split}C^{\varepsilon}(x,t)&=\frac{m^{2}}{4}\int\frac{\mathrm{d}\theta}{2\pi}\int\frac{\mathrm{d}\theta^{\prime}}{2\pi}\left[f(\theta)f(\theta^{\prime})(e^{\theta-\theta^{\prime}}-1)e^{i(p+p^{\prime})\cdot
x}+(1-f(\theta))(1-f(\theta^{\prime}))(e^{\theta-\theta^{\prime}}-1)e^{-i(p+p^{\prime})\cdot
x}\right.\\\
&\left.+f(\theta)(1-f(\theta^{\prime}))(e^{\theta-\theta^{\prime}}+1)e^{i(p-p^{\prime})\cdot
x}+(1-f(\theta))f(\theta^{\prime})(e^{\theta-\theta^{\prime}}+1)e^{-i(p-p^{\prime})\cdot
x}\right]\,,\end{split}$ (92)

where we used the Lorentz product notation, $p\cdot x=\varepsilon t-px.$ At
zero temperature $T=0,$ $f(\theta)=0$ and we obtain the closed form result

$\begin{split}C(x,t)&=\frac{m^{2}}{4}\int\frac{\mathrm{d}\theta}{2\pi}\int\frac{\mathrm{d}\theta^{\prime}}{2\pi}(e^{\theta-\theta^{\prime}}-1)e^{-i(p+p^{\prime})\cdot
x}=\frac{m^{2}}{4}\int\frac{\mathrm{d}\theta}{2\pi}e^{\theta}e^{-ip\cdot
x}\int\frac{\mathrm{d}\theta^{\prime}}{2\pi}e^{-\theta}e^{-ip^{\prime}\cdot
x}-\frac{m^{2}}{4}\left(\int\frac{\mathrm{d}\theta}{2\pi}e^{-ip\cdot
x}\right)^{2}\\\
&=m^{2}K_{0}^{2}\left(m\sqrt{x^{2}-t^{2}}\right)-m^{2}K_{1}^{2}\left(m\sqrt{x^{2}-t^{2}}\right)\,.\end{split}$
(93)

Note that since $(p+p^{\prime})\cdot
x=(\varepsilon_{\theta}+\varepsilon_{\theta^{\prime}})t-(p_{\theta}+p_{\theta^{\prime}})x$
and $\varepsilon_{\theta}+\varepsilon_{\theta^{\prime}}\geq 2\Delta,$ after
Fourier transformation $S(\omega,q)|_{T=0}=0$ for $0<\omega<2\Delta.$ At low
temperature, the leading order can be obtained by approximating
$f(\theta)\approx e^{-\beta\varepsilon(\theta)}$ and keeping only first powers
of $f(\theta):$

$\begin{split}C(x,t)&\approx\frac{m^{2}}{4}\int\frac{\mathrm{d}\theta}{2\pi}\int\frac{\mathrm{d}\theta^{\prime}}{2\pi}\left[(1-f(\theta)-f(\theta^{\prime}))(e^{\theta-\theta^{\prime}}-1)e^{-i(p+p^{\prime})\cdot
x}\right.\\\
&\left.+f(\theta)(e^{\theta-\theta^{\prime}}+1)e^{i(p-p^{\prime})\cdot
x}+f(\theta^{\prime})(e^{\theta-\theta^{\prime}}+1)e^{-i(p-p^{\prime})x}\right]\,.\end{split}$
(94)

Taking the Fourier transform, for frequencies $|\omega|<m$ only the second
line gives nonzero contribution and it recovers the expression Eq. (11).

##  Appendix H Asymptotic analysis of the integral Eq. (75)

In this appendix we report the details of the asymptotic analysis of Eq. (75)
for the local transverse DSF. We approximate the integral by dividing it into
two integrals at the extreme point
$q_{m}=\sqrt{\frac{\omega\Delta}{J(\Delta+2J)}}$ of the exponent:

$\begin{split}S^{xx}(\omega)&\approx\int_{q_{c}}^{\pi}\frac{dq}{\pi}\
\left|\frac{\sqrt{1+g^{2}-2g\cos\frac{q}{2}}}{gJ\sin\frac{q}{2}}\right|\,e^{-\epsilon_{k_{-}+q/2}/(k_{\text{B}}T)}\\\
&\approx\int_{q_{c}}^{\pi}\frac{dq}{\pi}\
\frac{\frac{\Delta}{Jq}+q\left(\frac{J}{2\Delta}+\frac{\Delta}{24J}+\frac{1}{2}\right)}{\frac{\Delta}{2}+J}\cdot\exp\left\\{-\frac{1}{k_{\text{B}}T}\left[\Delta+\frac{\Delta}{4J(\Delta+2J)}\frac{\omega^{2}}{q^{2}}+q^{2}J\left(\frac{J}{2\Delta}+\frac{1}{4}\right)\right]\right\\}\\\
&\approx\int_{q_{c}}^{q_{m}}\frac{dq}{\pi}\
\frac{\frac{\Delta}{Jq}+q\left(\frac{J}{2\Delta}+\frac{\Delta}{24J}+\frac{1}{2}\right)}{\frac{\Delta}{2}+J}\cdot\exp\left\\{-\frac{1}{k_{\text{B}}T}\left[\Delta+\frac{\Delta}{4J(\Delta+2J)}\frac{\omega^{2}}{q^{2}}\right]\right\\}\\\
&+\int_{q_{m}}^{\pi}\frac{dq}{\pi}\
\frac{\frac{\Delta}{Jq}+q\left(\frac{J}{2\Delta}+\frac{\Delta}{24J}+\frac{1}{2}\right)}{\frac{\Delta}{2}+J}\cdot\exp\left\\{-\frac{1}{k_{\text{B}}T}\left[\Delta+q^{2}J\left(\frac{J}{2\Delta}+\frac{1}{4}\right)\right]\right\\}\\\
&=\frac{1}{\pi}\frac{e^{-\frac{1}{k_{\text{B}}T}\left(\Delta+\frac{\Delta
J}{\Delta+2J}+\frac{\omega}{4}\right)}}{96\Delta
J^{3}(\Delta+2J)^{2}}\left[\omega\left(\Delta^{2}+12J^{2}+6\Delta
J\right)\left(4\Delta Je^{\frac{\Delta J}{k_{\text{B}}T(\Delta+2J)}}-\omega
e^{\frac{\omega}{4k_{\text{B}}T}}(\Delta+2J)\right)\right.\\\ &\left.-\Delta
Je^{\frac{\omega}{4k_{\text{B}}T}+\frac{\Delta
J}{k_{\text{B}}T(\Delta+2J)}}\left(96\Delta
J(\Delta+2J)-\frac{\omega^{2}}{k_{\text{B}}T}\left(\Delta^{2}+12J^{2}+6\Delta
J\right)\right)\left(\Gamma\left(0,\frac{J\Delta}{k_{\text{B}}T(2J+\Delta)}\right)-\Gamma\left(0,\frac{\omega}{4T}\right)\right)\right]\\\
&+\frac{1}{\pi}\frac{e^{-\frac{1}{k_{\text{B}}T}\left(\Delta+\frac{\pi^{2}J(\Delta+2J)}{4\Delta}+\frac{\omega}{4}\right)}}{6J^{2}(\Delta+2J)^{2}}\left[6\Delta
J(\Delta+2J)e^{\frac{1}{k_{\text{B}}T}\left(\frac{\pi^{2}J(\Delta+2J)}{4\Delta}+\frac{\omega}{4}\right)}\left(\text{Ei}\left(-\frac{J\pi^{2}(2J+\Delta)}{4\Delta
T}\right)-\text{Ei}\left(-\frac{\omega}{4k_{\text{B}}T}\right)\right)\right.\\\
&\left.+k_{\text{B}}T\left(\Delta^{2}+12J^{2}+6\Delta
J\right)\left(e^{\frac{\pi^{2}J(\Delta+2J)}{4\Delta
k_{\text{B}}T}}-e^{\frac{\omega}{4k_{\text{B}}T}}\right)\right]\\\
&\approx\frac{1}{\pi}e^{-\frac{\Delta}{k_{\text{B}}T}}\left(-\frac{2\Delta}{J(\Delta+2J)}\left(\ln\left(\frac{\omega}{4k_{\text{B}}T}\right)+\gamma_{E}\right)+\frac{\Delta^{2}+12J^{2}+6\Delta
J}{6J^{2}(\Delta+2J)^{2}}k_{\text{B}}T\right)+\cdots\\\
&\approx\frac{\Delta}{\pi
J^{2}}e^{-\frac{\Delta}{k_{\text{B}}T}}\left[-\ln\left(\frac{\omega}{4k_{\text{B}}T}\right)+\frac{k_{\text{B}}T}{2\Delta}-\gamma_{E}\right]+\cdots\,,\end{split}$
(95)

with the incomplete gamma function $\Gamma(a,z)=\int_{z}^{\infty}dt\
t^{a-1}e^{-t}$ and the exponential integral function
$\text{Ei}(z)=-\int_{-z}^{\infty}\frac{e^{-t}}{t}\,dt$. The last line is
obtained by taking scaling limit $\Delta/J\rightarrow 0$, namely, $g\to
g_{c}$, and the result agrees with field theory result Eq.(12).

##  Appendix I NMR relaxation rates for large nuclear spin

For nuclear spin $I>1/2$, the nuclear quadrupole interaction splits the
nuclear spin energy levels, and Eq.(16) needs to be evaluated based on Bloch-
Wangsness-Redfield theory using the density matrix for nuclear spin
$\rho_{\alpha\alpha^{\prime}}$ Slichter (1990); Pennington and Slichter
(1991); Pennington (1989); Barrett (1992),

$\frac{d\rho_{\alpha\alpha^{\prime}}}{dt}=\sum\limits_{\beta\beta^{\prime}}R_{\alpha\alpha^{\prime},\beta\beta^{\prime}}\rho_{\beta\beta^{\prime}},$
(96)

where $\alpha$, $\alpha^{\prime}$, $\beta$ and $\beta^{\prime}$ specify the
nuclear spin energy levels, and
$R_{\alpha\alpha^{\prime},\beta\beta^{\prime}}$ is the element of the
relaxation matrix $R$. In this approach,
$1/T_{2}=R_{\alpha\alpha-1,\alpha\alpha-1}$ for the $I_{z}=\alpha$ to
$\alpha-1$ transition of a given nuclear spin $I$ Slichter (1990), and

$\frac{1}{T_{2}}=A~{}\frac{1}{T_{1}}+\gamma_{n}^{2}~{}\overline{h^{2}}~{}\tau_{o},$
(97)

where the pre-factor $A$ is a constant that depends on $I$, $\gamma_{n}$ is
the nuclear gyromagnetic ratio of the observed nuclear spin,
$\overline{h^{2}}$ represents the averaged fluctuating hyperfine magnetic
field along the direction of the external magnetic field (i.e. $x$-axis in the
present case of TFIC), $\tau_{\text{o}}$ is the correlation time
($\omega_{n}\tau_{\text{o}}\ll 1$), and the second term represents
$1/T_{2}^{\prime}$ within the framework of Redfield’s theory. In the case of
nuclear spin $I=1/2$ with no nuclear quadrupole splitting, $A=1/2$ Slichter
(1990). For the $I_{z}=+1/2$ to $-1/2$ central transition of $I=3/2$, earlier
work showed that $A=7/2$ Pennington and Slichter (1991); Pennington (1989). In
the case of $I=9/2$ at 93Nb sites in the TFIC candidate material CoNb2O6, the
calculations of $A$ are straightforward but rather tedious, and we obtained
$A=49/2$.

## References

  * Special Issue: Quantum Criticality and Novel Phases (2013) Special Issue: Quantum Criticality and Novel Phases, Phys. Status Solidi B 250, 417 (2013).
  * Paschen and Si (2021) S. Paschen and Q. Si, Nat. Rev. Phys. 3, 9 (2021).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Ramshaw  _et al._ (2015) B. J. Ramshaw, S. E. Sebastian, R. D. McDonald, J. Day, B. S. Tan, Z. Zhu, J. B. Betts, R. Liang, D. A. Bonn, W. N. Hardy, and N. Harrison, Science 348, 317 (2015).
  * Kirchner  _et al._ (2020) S. Kirchner, S. Paschen, Q. Chen, S. Wirth, D. Feng, J. D. Thompson, and Q. Si, Rev. Mod. Phys. 92, 011002 (2020).
  * Coleman and Schofield (2005) P. Coleman and A. J. Schofield, Nature 433, 226 (2005).
  * Löhneysen  _et al._ (2007) H. v. Löhneysen, A. Rosch, M. Vojta, and P. Wölfle, Rev. Mod. Phys. 79, 1015 (2007).
  * Dai  _et al._ (2009) J. Dai, Q. Si, J.-X. Zhu, and E. Abrahams, Proc. Natl Acad. Sci. USA 106, 4118 (2009).
  * Licciardello  _et al._ (2019) S. Licciardello, J. Buhot, J. Lu, J. Ayres, S. Kasahara, Y. Matsuda, T. Shibauchi, and N. E. Hussey, Nature 567, 213 (2019).
  * Niemeijer (1967) T. Niemeijer, Physica 36, 377 (1967).
  * Pfeuty (1970) P. Pfeuty, Ann. Phys. 57, 79 (1970).
  * Sachdev (2011) S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Barouch and McCoy (1971) E. Barouch and B. M. McCoy, Phys. Rev. A 3, 786 (1971).
  * Suzuki (1971) M. Suzuki, Prog. Theor. Phys. 46, 1337 (1971).
  * Suzuki (1976) M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976).
  * Kopp and Chakravarty (2005) A. Kopp and S. Chakravarty, Nature Physics 1, 53 (2005).
  * Coldea  _et al._ (2010) R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Science 327, 177 (2010).
  * Wang  _et al._ (2018) Z. Wang, T. Lorenz, D. I. Gorbunov, P. T. Cong, Y. Kohama, S. Niesen, O. Breunig, J. Engelmayer, A. Herman, J. Wu, K. Kindo, J. Wosnitza, S. Zherlitsyn, and A. Loidl, Phys. Rev. Lett. 120, 207205 (2018).
  * Cui  _et al._ (2019) Y. Cui, H. Zou, N. Xi, Z. He, Y. X. Yang, L. Shu, G. H. Zhang, Z. Hu, T. Chen, R. Yu, J. Wu, and W. Yu, Phys. Rev. Lett. 123, 067203 (2019).
  * Zhang  _et al._ (2020) Z. Zhang, K. Amelin, X. Wang, H. Zou, J. Yang, U. Nagel, T. Rõõm, T. Dey, A. A. Nugroho, T. Lorenz, J. Wu, and Z. Wang, Phys. Rev. B 101, 220411 (2020).
  * Zou  _et al._ (2021) H. Zou, Y. Cui, X. Wang, Z. Zhang, J. Yang, G. Xu, A. Okutani, M. Hagiwara, M. Matsuda, G. Wang, G. Mussardo, K. Hódsági, M. Kormos, Z. He, S. Kimura, R. Yu, W. Yu, J. Ma, and J. Wu, Phys. Rev. Lett. 127, 077201 (2021).
  * Wu  _et al._ (2018) J. Wu, L. Zhu, and Q. Si, Phys. Rev. B 97, 245127 (2018).
  * Wu  _et al._ (2014) J. Wu, M. Kormos, and Q. Si, Phys. Rev. Lett. 113, 247201 (2014).
  * unp (2013) Márton Kormos, Jianda Wu, and Qimiao Si, unpublished notes (2013).
  * Steinberg  _et al._ (2019) J. Steinberg, N. P. Armitage, F. H. L. Essler, and S. Sachdev, Phys. Rev. B 99, 035156 (2019).
  * Kinross  _et al._ (2014) A. W. Kinross, M. Fu, T. J. Munsie, H. A. Dabkowska, G. M. Luke, S. Sachdev, and T. Imai, Phys. Rev. X 4, 031008 (2014).
  * Jaccarino (1965) V. Jaccarino,  _Nuclear Resonance in Antiferromagnets_ , edited by G. T. Rado and H. Suhl, Vol. Magnetism IIA (Academic Press, 1965).
  * (28) There is no field corresponding to $\sigma_{y}$ in the Ising field theory. However, the DSF of $S^{yy}$ can be exactly established via $S^{zz}$ in the corresponding lattice model Wu  _et al._ (2014).
  * Essler and Konik (2009) F. H. L. Essler and R. M. Konik, J. Stat. Mech. 2009, P09018 (2009).
  * Pozsgay and Takács (2010) B. Pozsgay and G. Takács, J. Stat. Mech. 2010, P11012 (2010).
  * Iorgov  _et al._ (2011) N. Iorgov, V. Shadura, and Y. Tykhyy, J. Stat. Mech. 2011, P02028 (2011).
  * Berg  _et al._ (1979) B. Berg, M. Karowski, and P. Weisz, Phys. Rev. D 19, 2477 (1979).
  * Delfino and Mussardo (1995) G. Delfino and G. Mussardo, Nucl. Phys. B 455, 724 (1995).
  * Moriya (1956) T. Moriya, Prog. Theor. Physics (Kyoto) 16 (1956).
  * Wang  _et al._ (2021) X. Wang, H. Zou, K. Hódsági, M. Kormos, G. Takács, and J. Wu, Phys. Rev. B 103, 235117 (2021).
  * Szécsényi and Takács (2012) I. M. Szécsényi and G. Takács, J. Stat. Mech. 2012, P12002 (2012).
  * Slichter (1990) C. P. Slichter,  _Principles of magnetic resonance_ , 3rd ed. (Springer-Verlag, Berlin, 1990).
  * Pennington and Slichter (1991) C. H. Pennington and C. P. Slichter, Phys. Rev. Lett. 66, 381 (1991).
  * Pennington (1989) C. H. Pennington, (Ph.D. thesis, University of Illinois at Urbana-Champaign, 1989).
  * Barrett (1992) S. E. Barrett, (Ph.D. thesis, University of Illinois at Urbana-Champaign, 1992).