# Dynamical structure factor of the three-dimensional quantum spin liquid
candidate NaCaNi2F7

Shu Zhang  Department of Physics and Astronomy, Johns Hopkins University,
Baltimore, MD 21218  Institute for Quantum Matter, Johns Hopkins University,
Baltimore, MD 21218  Hitesh J. Changlani  Department of Physics, Florida State
University, Tallahassee, Florida 32306  Department of Physics and Astronomy,
Johns Hopkins University, Baltimore, MD 21218  Institute for Quantum Matter,
Johns Hopkins University, Baltimore, MD 21218  Kemp W. Plumb  Department of
Physics, Brown University, Providence, RI 02912  Oleg Tchernyshyov  Department
of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218
Institute for Quantum Matter, Johns Hopkins University, Baltimore, MD 21218
Roderich Moessner  Max-Planck Institute for the Physics of Complex Systems,
01187 Dresden, Germany

(August 27, 2024)

We study the spin-1 pyrochlore material NaCaNi2F7 with a combination of
molecular dynamics simulations, stochastic dynamical theory and linear spin
wave theory. The dynamical structure factor from inelastic neutron scattering
is well described with a near-ideal Heisenberg Hamiltonian incorporating small
anisotropic terms and weak second-neighbor interactions. We find that all
three approaches reproduce remarkably well the momentum dependence of the
scattering intensity as well as its energy dependence with the exception of
the lowest energies. These results are notable in that (i) the data show a
complete lack of sharp quasiparticle excitations in momentum space over much,
if not all, of the energy range; (ii) linear spin-wave theory appears to apply
in a regime where it would be expected to fail for a number of reasons. We
elucidate what underpins these surprises, and note that basic questions about
the nature of quantum spin liquidity in such systems pose themselves as a
result.

Introduction: Quantum spin liquids Anderson (1973) are enigmatic phases of
matter characterized by the absence of symmetry breaking and conventional
quasiparticles (magnons). The search for their realisation in actual magnetic
materials has targeted, but is not limited to, materials involving the
geometrically frustrated Lacroix  _et al._ (2011) triangular, kagome Zhou  _et
al._ (2017) and pyrochlore Gardner  _et al._ (2010) geometries with low spins.

Indeed, while there have been significant efforts to synthesize quantum spin
liquid materials in spin-1/2 systems in two dimensions, fewer efforts have
been devoted to three dimensions, see Ref. Zhou  _et al._ (2017) for a review.
This strategic choice is not without reason: high lattice coordination number
and high spin typically suppress quantum fluctuations, and thus in particular
favor conventional forms of magnetic order over quantum spin liquids in three
dimensions. However, it is now clear that this perspective is too pessimistic:
as a matter of principle, we know that certain types of spin liquid – in
particular, Coulombic U(1) spin liquids – can exist in $d=3$ but not in $d=2$
Moessner and Sondhi (2003); Hermele  _et al._ (2004); and as a matter of
practice, it looks as if quantum spin liquid phases need by no means be
restricted to $S=1/2$ exclusively Iqbal  _et al._ .

Despite several recent advances in the field, however, our understanding of
the actual properties of low-spin Heisenberg spin liquids in three dimensions
is very limited, as they are beyond the scope of practically all exact or
controlled approximate theoretical schemes. We are at a loss to describe
either their ground states or excitation spectra, unlike Ising models like
spin ice, where the simplest quantum versions Moessner and Sondhi (2003);
Hermele  _et al._ (2004) are amenable to quantum Monte Carlo (QMC) simulations
Banerjee  _et al._ (2008); Benton  _et al._ (2012). Experimental data is
therefore a particularly indispensable guide for our understanding of these
magnets, for an early review see Harris and Zinkin (1996).

In the quest to identify quantum spin liquids in real materials Knolle and
Moessner , one relies heavily on characteristic signatures in the magnetic
excitation spectra, as their ground states are often largely featureless. By
contrast, the excitations of spin liquids can be downright spectactular,
including in particular fractionalised Rajaraman (2002) and other unusual
emergent quasiparticles such as spinons in the spin-1/2 Heisenberg
antiferromagnet chain Bethe (1931); Lake  _et al._ (2013); Mourigal  _et al._
(2013), Majorana fermions in the Kitaev honeycomb model Kitaev (2006);
Banerjee  _et al._ (2017); Kasahara  _et al._ (2018), and magnetic monopoles
Castelnovo  _et al._ (2008) and photons in the U$(1)$ spin liquid Moessner and
Sondhi (2003); Hermele  _et al._ (2004); Benton  _et al._ (2012).

The dual challenge is thus to identify novel behaviour in experimental data on
candidates quantum spin liquid materials in $d=3$, and to devise a theoretical
framework for understanding the underlying behaviour. Here, we report progress
for the fluoride pyrochlore NaCaNi2F7 Krizan and Cava (2015), a prime $S=1$
quantum spin liquid candidate.

In NaCaNi2F7, the magnetic Ni2+ ions reside on the three-dimensional
pyrochlore lattice (Fig. 1), where what little is known theoretically about
quantum Heisenberg models for $S=1/2$ and $S=1$ points towards quantum spin
liquid behavior Canals and Lacroix (1998); Iqbal  _et al._ , while the
classical case is well-established to be a Coulomb spin liquid Villain (1979);
Moessner and Chalker (1998a, b), whose many-body dynamics is by now fairly
well understood Moessner and Chalker (1998a, b); Conlon and Chalker (2009).

We analyze the magnetic excitation spectrum obtained by inelastic neutron
scattering on NaCaNi2F7 in Ref. Plumb  _et al._ , which we supplement with new
data from a different experiment. We present three tractable complementary
theoretical approaches that reproduce the dynamical structure factor
$\mathcal{S}(\mathbf{q},\omega)$ for all momenta $\mathbf{q}$ and for a broad
range of energies $\omega$. At the highest energies, the quality of the
agreement differs between models (and becomes harder to assess on account of a
considerable phonon background). At low energies, we find the well-known
pinch-point motifs, while at intermediate energies, characteristic structures
complementary to the pinch points appear Yan  _et al._ ; Mizouguchi  _et al._
. Overall, the main disagreement between experiment and theory appears at the
lowest energies, as discussed below.

In the light of the abovementioned challenges posed by three-dimensional
quantum spin liquids, the capacity of our relatively simple approaches to
yield a wide-ranging account of the observed dynamics is as striking as it is
encouraging for the study of other yet unexplored systems and models in this
class. We therefore include a discussion of the broader implications of our
results about the nature of the quantum dynamics in such a setting, which we
believe may be of importance well beyond the material studied here.

From a standpoint of basic phenomenology, we particularly emphasize that none
of the employed theoretical approaches relies on the existence of branches of
spin waves or other quasiparticles with well-defined wavevector $\mathbf{q}$
and frequency $\omega$, nor do they require the presence of delicate quantum
coherence. At the same time, linear spin wave theory is among the methods
which is successful in this context!

We next introduce our model, and the main results on the dynamical structure
factor of NaCaNi2F7are presented in the form of a comparison of scattering
intensities as a function of momentum (Fig. 2) and energy (Figs. 3, 4).
Besides their interpretation and discussion, for the methodologically
interested reader, we collate all necessary technical information in a set of
self-contained technical appendices.

Figure 1: (a). Pyrochlore lattice in one cubic unit cell. (b). Nearest-
neighbor and next-nearest-neighbor interactions.

Figure 2: (Color) Momentum dependent dynamical structure factor in [HHL] plane
(a) and [H0L] plane (b) at constant energies: inelastic neutron experiment
compared to linear spin wave theory, molecular dynamics, and the stochastic
model. The data in panel (a) were collected on CNCS and the data in panel (b)
were collected on MACS. Raw neutron scattering intensity has been corrected by
the magnetic form-factor for Ni2+ Brown (2006). To focus on the wavevector
dependence, the data are rescaled for each value of energy, for the
experimental data, by the maximum magnetic scattering intensity; and by the
maximum intensity in the MD simulations for the six theory panels, with an
additional factor $\beta\omega$ between MD/SLN and LSWT (Eq. (23), see text),
${\mathcal{N}}(\omega)=\beta\omega{\mathcal{S}}^{\textrm{MD}}_{\textrm{max}}(\omega)$,
where $\beta=1/k_{B}T$.

Model and methods: We use the Hamiltonian (Fig. 1)

$H\\!=\\!\frac{1}{2}\sum_{ij}\sum_{\mu\nu}J^{\mu\nu}_{ij}s^{\mu}_{i}s^{\nu}_{j},$
(1)

where subscripts $i$ and $j$ refer to lattice sites and superscripts $\mu$ and
$\nu$ refer to Cartesian components of spins in the global frame. The
interaction matrix is parameterized by four exchange parameters between
nearest neighbors
$J_{01}=(J_{2},J_{4},J_{4};-J_{4},J_{1},J_{3};-J_{4},J_{3},J_{1})$ with
$J_{1}\\!=\\!J_{2}\\!=\\!3.2(1)$ meV, $J_{3}\\!=\\!0.019(3)$ meV,
$J_{4}\\!=\\!-0.070(4)$ meV, and isotropic between next-nearest-neighbors
$J_{NNN}\\!=\\!-0.025(5)$ meV, obtained by fitting to the equal-time
correlations by some of us in a previous study Plumb  _et al._ . The
interaction matrices for other pairs follow from appropriate symmetry
transformations.

The methods utilized are, firstly, molecular dynamics (MD) simulations of the
pyrochlore magnet Moessner and Chalker (1998a) where the semiclassical Landau-
Lifshitz equations of motion for the spins are integrated numerically,
averaged over initial conditions obtained from Monte Carlo simulations of our
Hamiltonian at the temperature $T\\!=\\!1.8$ K. Secondly, we employ a self-
consistent Gaussian approximation adapted to frustrated magnets Garanin and
Canals (1999) and extended into a stochastic model for their dynamics by
Conlon and Chalker Conlon and Chalker (2009), which we refer to as stochastic
large-$N$ (SLN). Thirdly, we use the linear spin-wave theory (LSWT) to
describe spin dynamics near a low-energy state (again averaged over an
ensemble obtained form Monte Carlo simulations). Details of these are provided
in the appendix.

The central object of investigation are the dynamical spin correlations as
captured by the structure factor in energy-momentum space

$\displaystyle{\mathcal{S}}({\mathbf{q}},\omega)=\sum_{\mu\nu}(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^{2}})$
(2) $\displaystyle\qquad\times\frac{1}{2\pi
N}\sum_{i,j=1}^{N}\int_{-\infty}^{\infty}dt\,e^{-i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})+i\omega
t}\langle s_{i}^{\mu}(t)s_{j}^{\nu}(0)\rangle.$

The classical expression is given above and the quantum expression is
sensitive to time order of the spin operators (see Appendix A.2). Throughout
this work, we rescale the two "classical" approaches MD and SLN by
$\beta\omega$ to make comparison with LSWT. The factor essentially arises from
the classical equipartition-based suppression of the low-$T$ intensity, as
opposed to the non-vanishing matrix elements in the quantum case, see Eq.
(23).

Results: The dynamical structure factor obtained experimentally and by the
three theories is depicted as a function of wavevector in a set of cuts at
various energies, Fig. 2, and as a function of energy along a set of paths
through reciprocal space, Fig. 3.

Fig. 2 displays normalized momentum cuts in the [HHL] and [H0L] planes at
energies 0.5, 2, 4 and 8 meV. At low energies, the pinch points characteristic
of Heisenberg pyrochlore magnets are clearly visible. This is in itself
interesting, as the presence of well-defined pinch points implies that each
tetrahedron has vanishing total magnetization Isakov  _et al._ (2004); Henley
(2010). In general, however, adjacent tetrahedra cannot both be in spin
singlet states, as their total spin operators do not commute. Therefore, while
for the classical theories, the pinch points sharpen as $\sqrt{T}$ as $T$ is
lowered Zinkin (1996), for $S=1/2$ they were found to be quite smeared out
Canals and Lacroix (1998), becoming sharper as $S$ increases. For $S=1$, a
theoretical prediction for the full-width at half maximum of the pinch point
in the static correlations at [002] (located at $(0,0,4\pi)$ in reciprocal
space) of $\delta q^{\mathrm{FWHM}}_{\mathrm{PP}}=4\pi/3$ Iqbal  _et al._ is
comparable to the value $\approx\pi$ extracted form the low-$T$ experimental
data.

As the energy increases, the overall intensity distribution changes little
initially, but whatever sharp features were present wash out; e.g., the
intensity minimum in the scattering rhombus around [202] is slowly filled in
and the pinch points broaden. At higher energies, the experimental signal is
increasingly polluted by the phonon background at large $q$, but it is still
possible to identify a qualitative rearrangement of the weight, especially in
the [HHL] data, with the area around the pinch-points growing into prominent
pairs of "half-moons" features at 8 meV. These have recently been identified
as a dispersing complement to the pinch points Yan  _et al._ ; Mizouguchi  _et
al._ . This feature is present in MD and LSWT, but not in SLN, which is
relaxational and does not capture the spin precession at high frequencies.

We next turn to the energy dependence of the data, depicted in Fig. 3, with
additional cuts from a different neutron instrument (see Appendix C) presented
in Fig. 4. The general shapes of experiment and MD/LSWT are very similar—a
broad signal with a vertical appearance reminiscent of a fountain. SLN fails
to capture the high-energy structure, which can therefore be ascribed to the
precessional spin dynamics not captured by this method; otherwise, the theory
plots essentially agree with one another.

The largest disagreement between theory and experiments occurs at low
frequencies, especially around [220], where a large increase of the
experimental signal below $\omega=1$ meV is not reflected in theory, Fig. 4.
We return to the issue of the low-frequency regime in the discussion below.

Figure 3: (Color) Energy dependence of dynamical structure factor along
momentum cuts [22L] and [HH2]. Neutron scattering intensity is in absolute
units (see Appendix C). Rescaled MD and LSWT in particular reproduce well the
shape of the broad dispersive curve, disagreeing mainly at the lowest
energies, while SLN fails to capture high-energy structure.

Figure 4: (Color) Comparison of dynamical structure factor between experiment
and MD/LSWT. (a) Energy-dependence at $\mathbf{q}=$ [220], [22$\frac{1}{2}$]
and [221]. Log scale is used for the y-axis to include the quasi-elastic
signals. (b) Momentum dependence along [22L] at $\omega=$ 2,4,6 and 8meV. The
neutron scattering intensity $I^{\mathrm{exp}}(\mathbf{q},\omega)$ is
background subtracted and normalized by the total spectral weight $\int
d\omega d^{3}{\mathbf{q}}\,I^{\mathrm{exp}}(\mathbf{q},\omega)$. MD (solid
lines) and LSWT (dashed lines) data
$\mathcal{S}^{\mathrm{calc}}(\mathbf{q},\omega)$ are normalized by
$(2/3)S(S+1)$, under an isotropic approximation to the sum rule.

Discussion: We next address several more general questions arising from our
central observation of what we believe is remarkable agreement between theory
and experiment at all but the lowest energies on one hand, and between MD and
LSWT on the other. We start with the latter, which is quite unexpected: the
$S=1$ Heisenberg pyrochlore antiferromagnet unifies several reasons why LSWT
should break down. Instead, it works (un)reasonably well, as evidenced in the
comparison by eye with experiment as well as in the detailed quantitative
agreement with MD, Fig. 4. The inauspicious ingredients are, firstly, the
absence of a state with long-range order around which to perturb, the
existence of which would have guaranteed a Goldstone mode as long-lived magnon
excitation. Other settings which lack long-range order, such as the $S=1/2$ or
$S=1$ Heisenberg chain, instead show a breakdown of LSWT, as their respective
low-energy descriptions involve not the gapless magnons but rather
fractionalized $S=1/2$ spinons and Haldane’s famous gap. Secondly, the spin
length, $S=1$, really is not particularly large in our setting, so that one
would generically expect at least considerable quantum renormalization
effects, all the more so since the classical local exchange field, a central
feature in suppressing fluctuations, is reduced as a result of geometric
frustration from $6S$ in a ferromagnet to $2S$. Finally, a finite fraction of
the spin-wave modes live at or near zero energy in LSWT, which implies the
onset of the many-particle continuum already at the bottom of the single-
particle spectrum. Above this onset, spin waves are generally expected to
cease to be a useful description of the excitation spectrum Chernyshev and
Zhitomirsky (2009).

So why does linear spin wave theory nonetheless work so well? LSWT actually
finds another route to work—it is not a theory of universal low-energy
hydrodynamic excitations, but more a description of the statistically typical
short-to-intermediate time behavior, which in fact does not do a good job
precisely at the lowest energies; thus in the end conforming to at least a
subset of the above expectations.

To see this, think of the (near-)zero frequency modes responsible for motion
between (near-)degenerate ground states, and oscillatory excitations around
these as driving this motion Moessner and Chalker (1998a, b); Conlon and
Chalker (2009). The latter have finite frequency and finite scattering rates.
The central ingredient is that, statistically, there appears to be no
difference between the fast spectra of states visited as the slow modes
evolve, so that such motion is not reflected in the broad spectra we consider
here. One may expect sample-to-sample fluctuations due to the disorder
inherent in the randomly sampled starting configuration to be small, not least
because the almost uniform exchange field implies that disorder is mainly in
off-diagonal, terms of the dynamical matrix. In keeping with this, we find
self-averaging in practice as only a few configurations are needed to obtain
smooth spectra for large system sizes (see Appendix B.3); and as expected for
weak disorder in three dimensions, the spin wave modes away from the band
edges are delocalized, as diagnosed by the scaling of their inverse
participation ratio with system size (Fig. 5).

,

Figure 5: (Color) Inverse participation ratio (IPR) of the normalized real-
space spin wave modes $\psi(\mathbf{r},\omega)$ as a function of frequency for
different system sizes. See Appendix B.3 for the calculation details. The
modes are delocalized everywhere except for possibly at the edges of the
spectrum.

This also resolves the conundrum why scattering of the spin waves—unavoidable
as the slow modes evolve Moessner and Chalker (1998b)—does not invalidate the
spin wave picture. Given LSWT finds a broad continuum in frequency space to
begin with, any further broadening of an individual mode due to its limited
lifetime will be small in temperature $T$, and therefore parametrically
smaller than the total (largely $T$-independent) bandwidth. Therefore, unlike
in the case of an initially sharp mode, lifetime broadening is insignificant.

Regarding the low-$T$ limit, the zero modes mentioned above have no dynamics
in LSWT (that is just the statement that their frequency is zero). The motion
along the ground-state manifold is thus essentially frozen out, and LSWT in
fact fails completely to capture their motion arising from scattering high-
energy excitations, which is present in (not rescaled) MD and SLN theories.

From the preceding paragraph, it is clear that our comparison is not
particularly sensitive to the detailed nature of the low-frequency behavior.
Indeed, it has been a common theme of several recent studies of exotic
magnetic dynamics that scattering away from low energies are most instructive.
While this part of the spectrum is not universal, it may permit simple models,
e.g. in terms of deconfined spinons in the case of the Heisenberg chain Lake
_et al._ (2013) which in its detailed agreement with experiment may be more
convincing than the relatively featureless, and fragile, low-energy universal
features.This is all the more so since, like here, this portion is often
experimentally harder to access. Furthermore, the most characteristic aspects
of the spin liquid ground states are topological in nature and as such
invisible to experimental probes that couple to local correlations anyway
Knolle and Moessner .

Further complicating the low-energy analysis is the presence of disorder and
(partial) freezing Krizan and Cava (2015); Cai  _et al._ , which will need to
be included in a separate nontrivial modeling effort Saunders and Chalker
(2007). Also, further small terms in the Hamiltonian to which the previous
fitting procedure may be insensitive, can additionally lead to shifts of
weight on a scale which is small in absolute units but nonetheless notable at
low energies. Further, to accurately model a low-energy window comparable in
size to the temperature, a more detailed correspondence between classical and
quantum calculations than our simple rescaling ansatz would be needed.

While bearing all of this in mind, we emphasize the complete absence of sharp
quasiparticle peaks characteristic of magnons with well-defined momenta and
energies both in theory and experiment. This reflects the spatially disordered
nature of the spin configurations in our classical theory, while posing the
question about the appropriate description of the corresponding low-
temperature quantum state. In particular, it will be interesting to know if
the small-spin pyrochlore Heisenberg antiferromagnet exhibits no well-defined
quasiparticle excitations at all.

The final basic issue raised by our study is the role of the "quantumness" in
this compound. The relative success of fully classical modeling across a broad
range of energies, at temperatures far below the Curie-Weiss scale, is rather
unexpected. The low-energy discrepancies discussed above seem like a small
price to pay for the huge simplicity of the theoretical approaches we have
employed. This calls for an experiment on analogous compounds with larger
spin, to investigate whether the low-energy regime will be better modeled
while retaining the other features already successfully accounted for.

Employing semi-classical modelling for what "ought to be" a quantum spin
liquid is not without precedent. This was done for the Kitaev honeycomb model
Kitaev (2006), which has the benefit that the availability of an exact
solution of the dynamical structure factor Knolle  _et al._ (2014) of the spin
liquid allows for a reliable comparison in detail. There Samarakoon  _et al._
(2017), the high-frequency portion of the response was accounted for modulo a
reasonable amount of data post-processing, while the physics related to the
emergent fluxes at low energies—the most direct manifestation of
fractionalization—remained inaccessible.

This of course suggests a similar scenario here, namely that qualitative
signatures of a quantum spin liquid are visible only at the lowest energies,
perhaps even only below the scope of the experimental data. In this case, the
challenge is to identify a framework which can account for such a rapid
crossover into a classical regime, where quantum mechanics mainly enters in
the mode occupation numbers. An alternative would be the absence of a
qualitatively distinct low-frequency quantum spin liquid regime altogether.
This could either happen intrinsically, if the emergent low-energy description
is amenable to a semi-classical description; or extrinsically, in that the
quantum spin liquid behavior is so fragile in practice that disorder or
coupling to phononic degrees of freedom destroys it entirely. There are many
tantalizing open questions. The minimal next step for which experimental input
would be most valuable would be to consider materials with other values of
spin—ideally both $S=1/2$ and higher spin values—as well as extending the
experimental window further down towards the asymptotic low-frequency
behavior, if possible in a sample including minimal disorder.

Acknowledgements: We thank John Chalker, Chris Laumann and Johannes Reuther
for helpful discussions; and Collin Broholm for collaboration on the
experimental side of the project. This work was supported through the
Institute for Quantum Matter at Johns Hopkins University, by the U.S.
Department of Energy, Division of Basic Energy Sciences, Grant DE-
FG02-08ER46544, and by the Deutsche Forschungsgemeinschaft via grant SFB 1143.
HJC also thanks Florida State University for start-up funds. We gratefully
acknowledge the Johns Hopkins Homewood High Performance Cluster (HHPC) and the
Maryland Advanced Research Computing Center (MARCC), funded by the State of
Maryland, for computing resources. A portion of this research used resources
at the Spallation Neutron Source, a DOE Office of Science User Facility
operated by the Oak Ridge National Laboratory.

##  Appendix A Correspondence between classical and quantum spin wave theories

We discuss linear spin wave theory from the classical and quantum
perspectives. We begin with a classical ground state (or generally a local
minimum of energy), so that every spin $\mathbf{s}_{i}$ is in a state of
equilibrium. It is convenient to define a local frame with three mutually
orthogonal unit vectors $\mathbf{u}_{i}$, $\mathbf{v}_{i}$, and
$\mathbf{w}_{i}$, where $\mathbf{u}_{i}\\!=\\!\mathbf{s}_{i}/S$ points along
the equilibrium direction. Small deviations of $\mathbf{s}_{i}$ from its
equilibrium position can be parametrized in terms of two coordinates $x_{i}$
and $y_{i}$ as follows:

$\displaystyle\mathbf{s}_{i}$ $\displaystyle=$
$\displaystyle\sqrt{S^{2}-S(x_{i}^{2}+y_{i}^{2})}\,\mathbf{u}_{i}+\sqrt{S}(x_{i}\mathbf{v}_{i}+y_{i}\mathbf{w}_{i})$
(3) $\displaystyle\approx$
$\displaystyle\left(S-\frac{x_{i}^{2}+y_{i}^{2}}{2}\right)\mathbf{u}_{i}+\sqrt{S}(x_{i}\mathbf{v}_{i}+y_{i}\mathbf{w}_{i}).$

The dynamics of variables $x_{i}$ and $y_{i}$ is governed by the Lagrangian

$L=\sum_{i}\frac{1}{2}\left(y_{i}\frac{dx_{i}}{dt}-x_{i}\frac{dy_{i}}{dt}\right)-U,$
(4)

where $U$ is the potential energy encoding the spin interactions. Upon
expanding it to the second order in the deviations from equilibrium, we obtain

$L=\frac{1}{2}z^{T}\Gamma\frac{d}{dt}z-\frac{1}{2}z^{T}\mathcal{H}z,$ (5)

where $\mathcal{H}$ is a symmetric matrix, $\Gamma$ is a skew-symmetric matrix
and $z$ is a column vector:

$z\equiv\left(\begin{array}[]{c}x_{1}\\\ y_{1}\\\ \vdots\\\ x_{N}\\\
y_{N}\end{array}\right),\quad\Gamma=\left(\begin{array}[]{ccccc}0&-1&\ldots&0&0\\\
1&0&\ldots&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\ldots&0&-1\\\
0&0&\ldots&1&0\end{array}\right).$ (6)

###  A.1 Classical approach

In the classical approach, the momentum- and energy-dependent spin correlation
is defined as

$\displaystyle\mathcal{S}_{\text{classical}}^{\mu\nu}(\mathbf{q},\omega)=$
$\displaystyle\frac{1}{2\pi
N}\sum_{i,j=1}^{N}\int_{-\infty}^{\infty}dt\,e^{-i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})+i\omega
t}$ (7) $\displaystyle\qquad\qquad\qquad\times\langle
s_{i}^{\mu}(t)s_{j}^{\nu}(0)\rangle.$

Here $i,\,j\\!=\\!1,\dots,N$ are indices labeling lattice sites and
$\mu,\,\nu\\!=\\!x,y,z$ are indices for Cartesian spin components.

Tthe classical equations of motion for the deviations $z$ are given by the
Lagrangian (5), which correctly describe the precession of spins (24),

$\Gamma\frac{d}{dt}z=\mathcal{H}z.$ (8)

The solution to this equation (8)

$z(t)=\sum_{\alpha}c_{\alpha}\psi_{\alpha}e^{-i\omega_{\alpha}t}$ (9)

is a superposition with amplitude $c_{\alpha}$ of the normal modes
$\psi_{\alpha}$ of the eigenproblem

$(i\omega_{\alpha}\Gamma+\mathcal{H})\psi_{\alpha}=0.$ (10)

This eigenproblem has the following properties: the eigenvalues come in pairs
of real numbers $\omega_{-\alpha}\\!=\\!-\omega_{\alpha}$ and the eigenvectors
$\psi_{-\alpha}^{\dagger}\Gamma\psi_{-\alpha}\\!=\\!-\psi_{\alpha}^{\dagger}\Gamma\psi_{\alpha}$;
we thus choose the orthonormalization

$\psi_{\beta}^{\dagger}(-i\Gamma)\psi_{\alpha}=\text{sgn}(\omega_{\alpha})\delta_{\alpha\beta}.$
(11)

As a result,

$\psi_{\beta}^{\dagger}\mathcal{H}\psi_{\alpha}=\text{sgn}(\omega_{\alpha})\omega_{\alpha}\delta_{\alpha\beta}=|\omega_{\alpha}|\delta_{\alpha\beta}.$
(12)

and the potential energy is diagonalized in $c_{\alpha}$

$U=\frac{1}{2}\sum_{\alpha}\omega_{\alpha}c_{\alpha}^{\star}c_{\alpha}.$ (13)

The Boltzmann distribution gives the thermal average of the amplitudes
$\langle
c^{\star}_{\beta}c_{\alpha}\rangle\\!=\\!\delta_{\alpha\beta}/\beta|\omega_{\alpha}|$,
thus

$\langle
z_{k}(t)z_{l}(0)\rangle=\sum_{\alpha}\frac{[\psi_{\alpha}]_{k}\,[\psi_{\alpha}^{\dagger}]_{l}}{\beta|\omega_{\alpha}|}e^{-i\omega_{\alpha}t},$
(14)

where $k,\,l\\!=\\!1,\dots,2N$. Denote
$\psi_{\alpha}^{i}\equiv([\psi_{\alpha}]_{2i},[\psi_{\alpha}]_{2i+1})^{T}$ and
$\eta^{\mu}_{i}\equiv(v_{i}^{\mu},w_{i}^{\mu})$, The spin correlation function
(7) can be expressed as

$\displaystyle\mathcal{S}_{\text{classical}}^{\mu\nu}(\mathbf{q},\omega)=\frac{S}{N}\sum_{i,j=1}^{N}\,e^{-i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}$
(15)
$\displaystyle\qquad\qquad\qquad\times\sum_{\alpha}\frac{(\eta^{\mu}_{i}\ldotp\psi_{\alpha}^{i})\,(\eta^{\mu}_{j}\ldotp\psi_{\alpha}^{j})^{\star}}{\beta|\omega_{\alpha}|}\delta(\omega-\omega_{\alpha}).$

###  A.2 Quantum statistics

In the quantum approach, the spin correlation can be computed as the imaginary
part of the retarded response function

$\displaystyle\mathcal{G}_{+}^{\mu\nu}=\frac{1}{N}\sum_{i,j=1}^{N}\int_{0}^{\infty}dt\,e^{-i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})+i\omega
t}(-i)\langle[s_{i}^{\mu}(t),s_{j}^{\nu}(0)]\rangle.$ (16)

As the temperature goes to zero,

$\mathcal{S}_{\text{quantum}}^{\mu\nu}(\mathbf{q},\omega)=\frac{1}{\pi}\lim_{T\rightarrow
0}(1-e^{-\beta\omega})^{-1}\text{Im}\mathcal{G}_{+}^{\mu\nu}=\frac{1}{\pi}\text{Im}\mathcal{G}_{+}^{\mu\nu}.$
(17)

To evaluate the functional average over the ensemble, we take
$t\rightarrow-i\tau$ in (5) to get the imaginary-time partition function

$\mathcal{Z}=\int\mathcal{D}z\;\exp\left[{-\frac{1}{2}\int_{0}^{\beta}d\tau\;z^{T}(-i\Gamma\frac{d}{d\tau}+\mathcal{H})z}\right].$
(18)

In the representation of Matsubara frequencies $\omega_{n}\\!=\\!2\pi n/\beta$
($n\in\mathbb{Z}$),

$z(\tau)=\frac{1}{\sqrt{\beta}}\sum_{\omega_{n}}\zeta_{n}e^{-i\omega_{n}\tau}.$
(19)

We diagonalize the action in the partition function (18) by decomposing
$\zeta_{n}$ into orthonormal vectors $\psi_{\alpha}$ (11) and calculate the
matrix elements of the propagator:

$\langle\zeta_{n}\zeta_{n}^{\dagger}\rangle_{kl}=\sum_{\alpha}\frac{[\psi_{\alpha}]_{k}\,[\psi_{\alpha}^{\dagger}]_{l}}{-i\omega_{n}+\omega_{\alpha}}\text{sgn}(\omega_{\alpha}).$
(20)

Replacing $i\omega_{n}$ by $\omega+i0^{+}$ gives the retarded response
function

$\langle
z_{k}(\omega)z_{l}(-\omega)\rangle_{+}=-\sum_{\alpha}\frac{[\psi_{\alpha}]_{k}\,[\psi_{\alpha}^{\dagger}]_{l}}{\omega-\omega_{\alpha}+i0^{+}}\text{sgn}(\omega_{\alpha}).$
(21)

The spin correlation function (17) is thus

$\displaystyle\mathcal{S}_{\text{quantum}}^{\mu\nu}(\mathbf{q},\omega)=\frac{S}{N}\sum_{i,j=1}^{N}\,e^{-i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}$
(22)
$\displaystyle\qquad\qquad\times\sum_{\alpha}(\eta^{\mu}_{i}\ldotp\psi_{\alpha}^{i})\,(\eta^{\mu}_{j}\ldotp\psi_{\alpha}^{j})^{\star}\,\delta(\omega-\omega_{\alpha})\text{sgn}(\omega_{\alpha}),$

with the same notation as in (15).

###  A.3 Correspondence

Comparing Eqs. (15) and (22), we arrive at a relation between the finite
temperature classical calculation and the zero temperature quantum calculation
for spin correlations under the linear spin wave framework,

$\beta\omega\mathcal{S}_{\text{classical}}^{\mu\nu}(\mathbf{q},\omega)=\mathcal{S}_{\text{quantum}}^{\mu\nu}(\mathbf{q},\omega),$
(23)

which is applicable at low temperature and positive energy transfer
($\beta\omega\gg 1$).

##  Appendix B Numerical methods

In this work, we have studied the dynamical structure factor of the spin-1
pyrochlore material NaCaNi2F7 with a combination of classical Monte Carlo and
molecular dynamics simulations (MD), stochastic "large-N" (SLN) and real-space
linear spin wave theory (LSWT). Here we provide additional details needed to
reproduce our data in the main text.

###  B.1 Monte Carlo-molecular dynamics

We study the dynamics of interacting classical spins by integrating the
Landau-Lifshitz equation Landau and Lifshitz (1975)

$\frac{d}{dt}\mathbf{s}_{i}=-\mathbf{s}_{i}\times\frac{\partial
H}{\partial\mathbf{s}_{i}},$ (24)

which describes the precession of the spin in the local exchange field. Here
time $t$ is in the unit of meV-1.

Following previous work Moessner and Chalker (1998a, b); Conlon and Chalker
(2009), we perform Monte Carlo plus molecular dynamics simulations with the
Hamiltonian $H$ derived previously by some of us for NaCaNi2F7 Plumb  _et al._
(see main text). The initial configuration (IC) of spins is drawn by a Monte
Carlo (MC) run from the Boltzmann distribution $\exp(-\beta H)$ at temperature
$T\\!=\\!1.8\leavevmode\nobreak\ K$. Then, for each starting configuration the
spins are deterministically evolved according to Eq. (24). This is done for
many independent initial configurations and the result is averaged,

$\displaystyle\langle s^{\mu}_{i}(t)s^{\nu}_{j}(0)\rangle=\sum_{\text{IC from
MC}}s^{\mu}_{i}(t)s^{\nu}_{j}(0)\big{|}_{\text{IC}}.$ (25)

We simulated spins on the pyrochlore lattice with $N\\!=\\!16L^{3}$ sites,
where $L^{3}$ is the number of cubic unit cells and we show the results for
$L\\!=\\!8$. $6000$ initial configurations were used from independent Monte
Carlo runs.

To compute the dynamical spin structure factor in practice, we evolve the
Landau-Lifshitz equation (24) for a long but finite time
$T_{s}\\!=\\!60\leavevmode\nobreak\ \mathrm{meV}^{-1}$ , using the fourth-
order Runge Kutta method with discretized time steps. The time step $\delta
t\\!=\\!0.02$ meV-1 is chosen to be large enough to allow us to reach large
$T_{s}$ efficiently, yet small enough to ensure that the energy is constant
during the entire time evolution (to an accuracy of roughly six digits in the
energy per site).

In order to compute energy-momentum spin correlations computationally
efficiently, we perform a Fourier transform of the spin configurations during
the time evolution

$\tilde{s}^{\mu}({\bf
q},\omega)=\frac{1}{\sqrt{N}}\sum_{i}\frac{1}{T_{s}}\sum_{t}e^{-i{\bf
q}\cdot{\bf r}_{j}+i\omega t}s^{\mu}({\bf r}_{i},t).$ (26)

Then (7) is equivalent to

$\displaystyle\mathcal{S}^{\mu\nu}({\bf{q}},\omega)$ $\displaystyle=$
$\displaystyle\frac{T_{s}}{2\pi}\langle\tilde{s}^{\mu}({\bf
q},\omega){\tilde{s}^{\nu}}(-{\bf q},-\omega)\rangle.$ (27)

As an overall check, we took the integration $\int
d\omega\,\mathcal{S}^{\mu\nu}({\bf{q}},\omega)$ to compare with the static
spin correlations $\mathcal{S}^{\mu\nu}({\bf{q}})$ of the spin configurations
sampled from Monte Carlo runs. They are consistent with each other.

###  B.2 Stochastic model

In the stochastic model, we study the hydrodynamic motion of the spins under
the large-$N$ approximation ($N$ as of spin components). The spin
configuration drifts under a generalized force, while a noise with Gaussian
distribution provides thermal fluctuations Conlon and Chalker (2009):

$\frac{d}{dt}s^{\mu}_{i}=\gamma\sum_{j}\Delta_{ij}\frac{\partial E}{\partial
s_{j}^{\mu}}+\xi^{\mu}_{i}(t),$ (28)

where $\gamma$ is a dynamical parameter to fit to the MD data (Fig. 6) and
$\Delta_{ij}\\!=\\!A^{(1)}_{ij}-z\delta_{ij}$ is the lattice Laplacian
($A^{(1)}_{ij}$ is the first-order adjacency matrix and $z$ is the
coordination number).

The Boltzmann factor is

$\beta E\equiv\sum_{ij}\sum_{\mu\nu}\frac{1}{2}(\beta
J^{\mu\nu}_{ij}+\lambda\delta_{ij}\delta^{\mu\nu})s^{\mu}_{i}s^{\nu}_{j},$
(29)

where $\lambda$ is the Lagrange multiplier fixing the average length of the
"soft" spins. As is the case in the self-consistent Gaussian approximation
Garanin and Canals (1999); Plumb  _et al._ , $\lambda$ is solved self-
consistently from

${}NS^{2}=\sum_{\mathbf{q}\in\text{BZ}}\sum_{\rho}\frac{1}{\beta\epsilon_{\rho}(\mathbf{q})+\lambda},$
(30)

where $\epsilon_{\rho}(\mathbf{q})$ is the eigenvalues of the interaction
matrix in the reciprocal space.

The noise variables follow the independent Gaussian distribution
$\langle\xi^{\mu}_{i}(t)\rangle\\!=\\!0$ and
$\langle\xi^{\mu}_{i}(t)\xi^{\nu}_{j}(t^{\prime})\rangle\\!=\\!-(2\gamma/\beta)\Delta_{ij}\delta^{\mu\nu}\delta(t-t^{\prime})$,
whose amplitude is determined by the fluctuation-dissipation theorem.

The lattice Fourier transform is performed on each of the four fcc sublattices
($a\\!=\\!1,\dots,4$),

${}\tilde{s}_{a}^{\mu}(\mathbf{q})=\frac{1}{\sqrt{N_{c}}}\sum_{i\in
a}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}}s_{i}^{\mu},$ (31)

where $N_{c}$ is the number of up- or down- tetrahedra and $4N_{c}\\!=\\!N$ is
the number of sites. The lattice Laplacian and the interaction matrix
transform accordingly into momentum space, denoted by $\Delta(\mathbf{q})$ and
$J(\mathbf{q})$ respectively. Define $P\\!\equiv\\!-\Delta(\mathbf{q})\otimes
I_{3}$ and $Q\\!\equiv\\!J(\mathbf{q})+(\lambda/\beta)I_{12}$, where $I_{n}$
is an $n\times n$ identity matrix; the equation of motion in the momentum-
energy space of the 12-component vector
$\tilde{S}\\!\equiv\\!(\tilde{s}_{1}^{x},\tilde{s}_{1}^{y},\tilde{s}_{1}^{z},\dots,\tilde{s}_{4}^{x},\tilde{s}_{4}^{y},\tilde{s}_{4}^{z})^{T}$
can be expressed by matrix multiplication:

$\tilde{S}(\mathbf{q},\omega)=G(\mathbf{q},\omega)\,\tilde{\xi}(\mathbf{q},\omega),$
(32)

where the Green’s function is

$G^{-1}(\mathbf{q},\omega)=-i\omega I_{12}+\gamma PQ.$ (33)

Because $P$ and $Q$ are real and symmetric, $Q$ is positive definite and $P$
is semi-positive definite, there exists a similarity transformation under
matrix $V$ to obtain a diagonal matrix $\Lambda$ with real and non-negative
entries (generally $PQ\neq QP$),

$QP=V\Lambda V^{-1}.$ (34)

It can then be derived

$\displaystyle\langle\tilde{S}_{\alpha}(\mathbf{q},\omega)\tilde{S}_{\beta}(-\mathbf{q},-\omega)\rangle=\frac{2\gamma}{\beta}(GPG^{\dagger})_{\alpha\beta}$
(35)
$\displaystyle\qquad\qquad=\frac{2\gamma}{\beta}\left[PV(\omega^{2}+\gamma^{2}\Lambda^{2})^{-1}V^{-1}\right]_{\alpha\beta}.$

The classical spin correlations (7) thus evaluate to be

$\mathcal{S}^{\mu\nu}(\mathbf{q},\omega)=\frac{1}{8\pi}\sum_{\alpha\beta}\kappa^{\mu}_{\alpha}\kappa^{\nu}_{\beta}\langle\tilde{S}_{\alpha}(\mathbf{q},\omega)\tilde{S}_{\beta}(-\mathbf{q},-\omega)\rangle,$
(36)

where $\kappa^{\mu}\\!=\\!(1,1,1,1)\otimes\hat{\mathbf{e}}^{\mu}$ and the
factor from $N\\!=\\!4N_{c}$ is accounted.

We further confirm the fluctuation-dissipation theorem is obeyed despite the
subtlety from the non-commutation of the matrices. The static spin
correlations upon integration of (36) over energies is analytically the same
as in our previous work Plumb  _et al._ .

Figure 6: (Color) The time-dependent dynamical structure factor
$\mathcal{S}(\mathbf{q},\omega)\\!=\\!\sum_{\mu\nu}(\delta_{\mu\nu}-q_{\mu}q_{\nu}/q^{2})\mathcal{S}^{\mu\nu}(\mathbf{q},\omega)$
from the MD (dots) and SLN (lines) along the path [000] $\rightarrow$ [222]
$\rightarrow$ [220] $\rightarrow$ [000] in momentum space at 4 different times
including the static structure factor ($t\\!=\\!0$). This was used to fit the
dynamical parameter $\gamma$ in the stochastic model, $\gamma\\!=\\!0.165$.

###  B.3 Real space linear spin wave theory

A typical application of linear spin wave theory is usually based on a
classical ground state candidate that is derived or postulated. Often, the
simplest ground states have a small unit cell of $n$ spins (say a pattern on a
single tetrahedron, $n=4$, for the pyrochlore lattice Ross  _et al._ (2011)),
which repeats in real space, so the wavevector is a good quantum number. The
Hamiltonian in momentum space is block diagonalized and gives $n$ bands with
sharp dispersions. However, the situation is markedly different for
inhomogeneous ground states or states with very large unit cells. It becomes
hard to track bands for large unit cells and for inhomogeneous states the loss
of translational invariance means momentum is no longer a good quantum number.

This latter situation is typical for the classical Heisenberg model on the
pyrochlore lattice. It has many ground states which satisfy the condition
$\sum_{i\in\boxtimes}\mathbf{s}_{i}\\!=\\!0$, and most of them are
inhomogeneous. (Small anisotropic and further-neighbor interactions lift this
degeneracy, but can still lead to the formation of many low-energy minima.)
Thus here we consider finite lattice clusters and perform linear spin wave
theory in real space directly using the formalism given by Eq. (22). A former
general consideration can be found in Ref. Maestro and Gingras (2004).

To closely mimic the situation in MD, we first assemble an ensemble of
classical ground (or metastable low-energy) states for the NaCaNi2F7
Hamiltonian by performing replica Monte Carlo runs ranging from very low
temperature ($T\\!=\\!0.1$ K) to high temperature ($T\\!=\\!10$ K). This
allows for good equilibration of spin configurations and largely prevents from
the simulation getting "stuck". (Formally, if all independent runs are run for
infinitely long they must all find the true ground state, this is not the case
in practical finite runs.) The last spin configuration encountered in each run
at the lowest temperature is used as the starting configuration for an
iterative algorithm. This algorithm works by aligning one spin with its local
exchange field; keeping all the other spins in the configuration fixed. One
sweep consists of $N$ such moves (one for each spin). Multiple sweeps are
performed until the spin directions stop changing completely, which guarantees
that a stable low-energy minimum has been achieved.

Figure 7: (Color) Real space linear spin wave theory results for different
lattice sizes $N\\!=\\!16L^{3}$ (a) and different numbers of samples $N_{s}$
with $L\\!=\\!8$ (b).

Then given a low-energy stable classical spin configuration, we construct and
then directly diagonalize the potential $U$ in real-space. To evaluate the
spin correlations (22) with discrete eigenvalues $\omega_{\alpha}$, we use a
very narrow Gaussian function as an approximation to the Dirac delta function,

$\delta(\omega-\omega_{\alpha})=\lim_{\epsilon\rightarrow
0}\frac{1}{\sqrt{2\pi}\epsilon}\exp\left[{-\frac{(\omega-\omega_{\alpha})^{2}}{2\epsilon^{2}}}\right],$
(37)

with $\epsilon=0.01$.

Finally, we average the dynamical structure factor over $N_{s}$ low-energy
spin configurations. The data in the main text is collected for
$N_{s}\\!=\\!250$ and $N\\!=\\!16L^{3}$ with $L\\!=\\!8$. Fig. 7(a) shows our
results for the dynamical structure factor along two different momentum cuts
as a function of energy for $L\\!=\\!4,6,8$. The finite size effects are
small. Fig. 7(b) shows the same cuts for different numbers of spin
configurations used in the average, all with $L\\!=\\!8$.

It can be shown that for a spin configuration with translational symmetry, our
real-space LSWT is equivalent to the result in momentum-space, which presents
sharp dispersions. Yet even for a single spin configuration (see Fig. 7(b)),
the branches are broad. This confirms that the classical local minima found in
the Monte Carlo simulation are indeed inhomogeneous.

We have also analyzed the inverse participation ratio (IPR) of the spin wave
modes to understand the localization and delocalization effects. In analogy to
the density distribution
$\rho_{e}(\mathbf{r})\\!=\\!|\psi_{e}(\mathbf{r})|^{2}$ of a given electron
wavefunction, which is normalized $\int
d^{3}\mathbf{r}\,\rho_{e}(\mathbf{r})\\!=\\!1$, we consider the normalization
(11) and define the density distribution of the spin wave modes at energy
$\omega_{\alpha}$ to be

$\rho_{\alpha}(\mathbf{r}_{i})=(\psi_{\alpha}^{i})^{\dagger}(-i\Gamma_{2})\psi_{\alpha}^{i},$
(38)

where $\Gamma_{2}\\!=\\!\left(\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right)$ and
$\psi_{\alpha}^{i}\equiv([\psi_{\alpha}]_{2i},[\psi_{\alpha}]_{2i+1})^{T}$,
corresponding to the two transverse spin deviations from the equilibrium
direction. The IPR is given by

$\text{IPR}(\omega_{\alpha})\equiv\sum_{i=1}^{N}|\rho_{\alpha}(\mathbf{r}_{i})|^{2}.$
(39)

For a delocalized mode, we expect $\text{IPR}(\omega_{\alpha})\sim 1/N$ while
for a localized mode $\text{IPR}(\omega_{\alpha})\sim O(1)$. In the main text
Fig. 5, we show the IPR multiplied by the number of sites, i.e.
$N\times\text{IPR}(\omega)$. The collapse of the values for different lattice
sizes (except perhaps at the edges of the spectrum, especially the upper edge)
indicates the delocalization of spin wave modes for a wide energy energy
range.

##  Appendix C Experimental data

The nonpolarized inelastic neutron scattering experiment probes the sum of the
components of the spin correlation function that is perpendicular to the
momentum transfer $\mathbf{q}$:

$\mathcal{S}({\bf{q}},\omega)=\sum_{\mu,\nu}\left(\delta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q^{2}}\right)\mathcal{S}^{\mu\nu}({\bf{q}},\omega).$
(40)

The measured intensity in the unit of scattering cross-section is given by

$I({\bf q},\omega)=\left(\frac{\gamma
r_{o}}{2}\right)^{2}|gf(q)|^{2}\mathcal{S}({\bf q},\omega),$ (41)

where $f(|{\bf q}|)$ is the magnetic form factor, $g$ is the g-factor, and
$\gamma r_{0}=0.539\times 10^{-12}$ cm is the neutron magnetic scattering
length. Throughout this work we use the dipole approximation to the Ni2+ form-
factor Brown (2006) and an estimation of $g=2.28$ Yamaguchi  _et al._ (1999).

Data presented in Fig. 2(b) and Fig. 3 is identical to the published data set
in Ref. Plumb  _et al._ . In addition to the previously published data set, we
have also collected a new complementary data set with finer momentum and
energy resolution. The inelastic neutron scattering data presented in Fig.
2(a) and Fig. 4 was collected on the cold neutron chopper spectrometer (CNCS)
at Oak Ridge National Lab using fixed incident neutron energies of 2.5 meV,
6.59 meV, and 12 meV providing energy resolution (FWHM) at the elastic lines
of 0.07 meV, 0.4 meV, and 0.95 meV respectively. The same single crystal
sample used for previous studies was mounted with the [HH0] and [00L]
directions in the horizontal scattering plane of the instrument. Data was
collected with [00L] initially directed along the incident neutron beam and
the sample rotated over 180∘. The sample was cooled to 350 mK for all
measurements. (No phase transition happens between 350 mK and 1.8 K and the
dynamical structure factor at finite energy ($\gtrsim 0.1$ meV) is mostly
temperature independent under the freezing temperature $T_{f}=3.6$ K.) Data in
Fig. 2(a) has been symmetrized by folding about the [HH0] and [00L] axes.

The data contains energy-dependent non-magnetic background intensity arising
from incoherent nuclear scattering and scattering from the sample environment.
This background contribution is estimated by taking a cut around the point
${\mathbf{q}}=[0\,0\,3.65]$ and subtracted from the data in Fig. 4.

## References

  * Anderson (1973) P. W. Anderson, Mater. Res. Bull. 8, 153 (1973).
  * Lacroix  _et al._ (2011) C. Lacroix, P. Mendels, and F. Mila,  _Introduction to Frustrated Magnetism_ (Springer, Berlin, Heidelberg, 2011).
  * Zhou  _et al._ (2017) Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Gardner  _et al._ (2010) J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Rev. Mod. Phys. 82, 53 (2010).
  * Moessner and Sondhi (2003) R. Moessner and S. L. Sondhi, Phys. Rev. B 68, 184512 (2003).
  * Hermele  _et al._ (2004) M. Hermele, M. P. A. Fisher, and L. Balents, Phys. Rev. B 69, 064404 (2004).
  * (7) Y. Iqbal, T. Müller, P. Ghosh, M. J. P. Gingras, H. O. Jeschke, S. Rachel, J. Reuther, and R. Thomale,  arXiv:1802.09546 .
  * Banerjee  _et al._ (2008) A. Banerjee, S. V. Isakov, K. Damle, and Y. B. Kim, Phys. Rev. Lett. 100, 047208 (2008).
  * Benton  _et al._ (2012) O. Benton, O. Sikora, and N. Shannon, Phys. Rev. B 86, 075154 (2012).
  * Harris and Zinkin (1996) M. J. Harris and M. P. Zinkin, Mod. Phys. Lett. B 10, 417 (1996).
  * (11) J. Knolle and R. Moessner,  arXiv:1804.02037 .
  * Rajaraman (2002) R. Rajaraman, “Fractional charge,” in  _Quantum [Un]speakables: From Bell to Quantum Information_ (Springer, Berlin, Heidelberg, 2002) pp. 383–399.
  * Bethe (1931) H. Bethe, Z. Physik 71, 205 (1931).
  * Lake  _et al._ (2013) B. Lake, D. A. Tennant, J.-S. Caux, T. Barthel, U. Schollwöck, S. E. Nagler, and C. D. Frost, Phys. Rev. Lett. 111, 137205 (2013).
  * Mourigal  _et al._ (2013) M. Mourigal, M. Enderle, A. Klöpperpieper, J.-S. Caux, A. Stunault, and H. M. Ronnow, Nat. Phys. 9, 435 (2013).
  * Kitaev (2006) A. Kitaev, Ann. Phys. 321, 2 (2006).
  * Banerjee  _et al._ (2017) A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Science 356, 1055 (2017).
  * Kasahara  _et al._ (2018) Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, Nature 559, 227 (2018).
  * Castelnovo  _et al._ (2008) C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).
  * Krizan and Cava (2015) J. W. Krizan and R. J. Cava, Phys. Rev. B 92, 014406 (2015).
  * Canals and Lacroix (1998) B. Canals and C. Lacroix, Phys. Rev. Lett. 80, 2933 (1998).
  * Villain (1979) J. Villain, Z. Phys. B. 33, 31 (1979).
  * Moessner and Chalker (1998a) R. Moessner and J. T. Chalker, Phys. Rev. Lett. 80, 2929 (1998a).
  * Moessner and Chalker (1998b) R. Moessner and J. T. Chalker, Phys. Rev. B 58, 12049 (1998b).
  * Conlon and Chalker (2009) P. H. Conlon and J. T. Chalker, Phys. Rev. Lett. 102, 237206 (2009).
  * (26) K. W. Plumb, H. J. Changlani, A. Scheie, S. Zhang, J. W. Kriza, J. A. Rodriguez-Rivera, Y. Qiu, B. Winn, R. J. Cava, and C. L. Broholm,  arXiv:1711.07509 .
  * (27) H. Yan, R. Pohle, and N. Shannon,  arXiv:1806.08520 .
  * (28) T. Mizouguchi, L. D. C. Jaubert, and M. Udagawa,  arxiv:1806.08534 .
  * Brown (2006) P. J. Brown,  _International Tables for Crystallography_ , Vol. C (Springer, Berlin, Heidelberg, 2006) Chap. 4.4.5, pp. 454 –461.
  * Garanin and Canals (1999) D. A. Garanin and B. Canals, Phys. Rev. B 59, 443 (1999).
  * Isakov  _et al._ (2004) S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 93, 167204 (2004).
  * Henley (2010) C. L. Henley, Annu. Rev. Condens. Matter Phys. 1, 179 (2010).
  * Zinkin (1996) M. P. Zinkin, Ph.D. thesis, University of Oxford (1996).
  * Chernyshev and Zhitomirsky (2009) A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. B 79, 144416 (2009).
  * (35) Y. Cai, M. N. Wilson, A. M. Hallas, L. Liu, B. A. Frandsen, S. R. Dunsiger, J. W. Krizan, R. J. Cava, O. Rubel, Y. J. Uemura, and G. M. Luke,  arXiv:1806.10970 .
  * Saunders and Chalker (2007) T. E. Saunders and J. T. Chalker, Phys. Rev. Lett. 98, 157201 (2007).
  * Knolle  _et al._ (2014) J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moessner, Phys. Rev. Lett. 112, 207203 (2014).
  * Samarakoon  _et al._ (2017) A. M. Samarakoon, A. Banerjee, S.-S. Zhang, Y. Kamiya, S. E. Nagler, D. A. Tennant, S.-H. Lee, and C. D. Batista, Phys. Rev. B 96, 134408 (2017).
  * Landau and Lifshitz (1975) L. D. Landau and E. M. Lifshitz,  _The Classical Theory of Fields_ , 4th ed. (Butterworth-Heinemann, Oxford, 1975).
  * Ross  _et al._ (2011) K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Phys. Rev. X 1, 021002 (2011).
  * Maestro and Gingras (2004) A. G. D. Maestro and M. J. P. Gingras, J. Phys. Condens. Matter 16, 3339 (2004).
  * Yamaguchi  _et al._ (1999) H. Yamaguchi, K. Katsumata, M. Hagiwara, M. Tokunaga, H. L. Liu, A. Zibold, D. B. Tanner, and Y. J. Wang, Phys. Rev. B 59, 6021 (1999).