# Spin dynamics and unconventional Coulomb phase in Nd2Zr2O7

M. Léger  Institut Néel, CNRS and Université Grenoble Alpes, 38000 Grenoble,
France  Laboratoire Léon Brillouin, Université Paris-Saclay, CNRS, CEA, CE-
Saclay, F-91191 Gif-sur-Yvette, France  E. Lhotel<EMAIL_ADDRESS>Institut Néel, CNRS and Université Grenoble Alpes, 38000 Grenoble, France  M.
Ciomaga Hatnean  Department of Physics, University of Warwick, Coventry, CV4
7AL, United Kingdom  J. Ollivier  Institut Laue Langevin, F-38042 Grenoble,
France  A. R. Wildes  Institut Laue Langevin, F-38042 Grenoble, France  S.
Raymond  Université Grenoble Alpes, CEA, IRIG, MEM, MDN, 38000 Grenoble,
France  E. Ressouche  Université Grenoble Alpes, CEA, IRIG, MEM, MDN, 38000
Grenoble, France  G. Balakrishnan  Department of Physics, University of
Warwick, Coventry, CV4 7AL, United Kingdom  S. Petit<EMAIL_ADDRESS>Laboratoire Léon Brillouin, Université Paris-Saclay, CNRS, CEA, CE-Saclay,
F-91191 Gif-sur-Yvette, France

###### Abstract

We investigate the temperature dependence of the spin dynamics in the
pyrochlore magnet Nd2Zr2O7 by neutron scattering experiments. At low
temperature, this material undergoes a transition towards an “all in - all
out” antiferromagnetic phase and the spin dynamics encompass a dispersion-less
mode, characterized by a dynamical spin ice structure factor. Unexpectedly,
this mode is found to survive above $T_{\rm N}\approx 300$ mK. Concomitantly,
elastic correlations of the spin ice type develop. These are the signatures of
a peculiar correlated paramagnetic phase which can be considered as a new
example of Coulomb phase. Our observations near $T_{\rm N}$ do not reproduce
the signatures expected for a Higgs transition, but show reminiscent features
of the “all in - all out” order superimposed on a Coulomb phase.

Geometrical frustration is well known to be one of the key ingredients leading
to unconventional states of matter, especially in magnetism [1, 2]. Among
them, spin ice and more generally Coulomb phases [3] have attracted
significant interest. These can be considered as an original state of matter
formed by disordered degenerate configurations where local degrees of freedom
remain strongly constrained at the local scale by an organizing principle. In
the case of spin ice, these degrees of freedom are Ising spins, sitting on the
sites of a pyrochlore lattice formed of corner sharing tetrahedra and aligned
along the axes which connect the corners of the tetrahedra to their center.
The organizing principle, the “ice rule”, states that each tetrahedron should
have two spins pointing in and two out, in close analogy with the rule which
controls the hydrogen position in water ice [4]. Importantly, the idea that
this local constraint can be considered as the conservation law of an
“emergent” magnetic flux (${\bf\nabla}\cdot{\bf B}=0$) was quickly imposed [5,
6, 7]. Quantum fluctuations can cause this flux to change with time, giving
rise to an emergent electric field, and eventually to an emergent quantum
electromagnetism [8, 9, 10]. This quantum spin ice state hosts spinon
(monopole in the spin ice language [11]) and photon like excitations. Despite
much work, however, experimental evidence for this enigmatic physics remains
elusive, with the possible exception of Pr2Hf2O7 [12]. Indeed, the conditions
for the realisation of this so-called quantum spin ice state are drastic:
transverse terms have to be sizable in the Hamiltonian to enable fluctuations
out of the local Ising axes, but should remain small enough to prevent the
stabilization of classical phases, called Higgs phases, characterized by
ordered components perpendicular to these axes [13, 11, 14].

The pyrochlore material Nd2Zr2O7 offers the opportunity to approach this
issue. Recent studies suggest that below 1 K this compound hosts a correlated
state, which could be a remarkable novel example of Coulomb phase [15, 16].
This phase would be described by a “two in – two out” rule as in spin ice, but
built on a pseudospin component different from the conventional $\langle
111\rangle$ Ising one. The “all in – all out” (AIAO) ordering previously
observed below $T_{\rm N}\approx 300$ mK [17, 18] would then correspond to the
pseudospin ordering in directions perpendicular to the components responsible
for the “high temperature” Coulomb phase. It was proposed that a Higgs
mechanism may account for this transition [16]. Such a process is invoked in
$U(1)$ quantum spin liquids when the deconfined spinon excitations undergo a
Bose-Einstein condensation, resulting in a Higgs phase along with a gapped
photon excitation [19, 13, 20].

In this letter, we show that the paramagnetic phase of Nd2Zr2O7 does carry
elastic spin ice-like correlations, and thus confirm the proposed Coulomb
phase picture above $T_{\rm N}$. We present a detailed study of the spin
dynamics as a function of temperature and explore the nature of this Coulomb
phase above and close to the transition. The spin excitations of Nd2Zr2O7 deep
in the AIAO phase include a peculiar spectrum with a flat band at the energy
$E_{0}\approx$ 70 $\mu$eV characterized by a spin ice-like ${\bf
Q}$-dependence [15, 21, 22]. Using neutron scattering experiments, we report
the temperature dependence of the gap $E_{0}$, and reveal that this gap
persists above $T_{\rm N}$. This result is robust, and withstands a small
substitution at the Zr site. The spectra recorded above $T_{\rm N}$ do not
show the spinon continuum expected in the Higgs scenario. Instead, we observe
dispersive features reminiscent of the AIAO ordered phase superimposed on the
Coulomb phase signal. This coexistence suggests that a strong exchange
competition is at work in this temperature range, emphasizing the originality
of the Coulomb phase above the transition.

The single crystal samples used in this work are the same as in our previous
studies (labeled #1 [17, 15, 21] and #2 [21]). In addition, results on a
single crystal of Nd2(Zr1-xTix)2O7, with $x=2.4$ % (Sample #3) (See
supplementary material [23]) are presented, not in order to analyze the role
of disorder but to illustrate the robustness of the results. Magnetic
properties were measured in very low temperature SQUID magnetometers developed
at the Institut Néel [29]. The composition and magnetic structure at low
temperature were determined using the D23 (CEA-CRG@ILL) neutron diffractometer
[23]. Polarized neutron scattering experiments were carried out at D7 (ILL) on
Sample #1. Inelastic neutron scattering (INS) experiments were carried out on
the IN5 (ILL) time of flight spectrometer on all samples and on the triple
axis spectrometer IN12 (CEA-CRG@ILL) for Sample #1. The INS data have been
analyzed using the cefwave software developed at LLB.

The XYZ Hamiltonian proposed to describe the properties of Nd based
pyrochlores due to the peculiar dipolar-octupolar character of the Nd3+ion
[30], writes:

${\cal H}=\sum_{\langle i,j\rangle}\left[{\sf
J}_{x}\tau^{x}_{i}\tau^{x}_{j}+{\sf J}_{y}\tau^{y}_{i}\tau^{y}_{j}+{\sf
J}_{z}\tau^{z}_{i}\tau^{z}_{j}+{\sf
J}_{xz}(\tau^{x}_{i}\tau^{z}_{j}+\tau^{z}_{i}\tau^{x}_{j})\right]$ (1)

In this Hamiltonian, $\tau_{i}$ is not the actual spin, but a pseudospin which
resides on the rare-earth sites of the pyrochlore lattice. Its $z$ component
relates to the usual magnetic moment and is directed along the local $\langle
111\rangle$ directions of the tetrahedra of the pyrochlore lattice. This
Hamiltonian can be rewritten by rotating the ${\bf x}$ and ${\bf z}$ axes in
the $({\bf x},{\bf z})$ plane by an angle $\theta$. In this
$({\bf\tilde{x}},{\bf\tilde{z}})$ rotated frame, the relevant parameters of
the Hamiltonian ${\cal H}$ are labeled $\tilde{{\sf J}}_{x,y,z}$, leading to
[30, 31]:

$\displaystyle{\cal H}_{\rm XYZ}=$ $\displaystyle\sum_{\langle
i,j\rangle}\left[{\tilde{\sf
J}_{x}}\tilde{\tau}^{\tilde{x}}_{i}\tilde{\tau}^{\tilde{x}}_{j}+{\tilde{\sf
J}_{y}}\tilde{\tau}^{\tilde{y}}_{i}\tilde{\tau}^{\tilde{y}}_{j}+{\tilde{\sf
J}_{z}}\tilde{\tau}^{\tilde{z}}_{i}\tilde{\tau}^{\tilde{z}}_{j}\right]$ (2)
$\displaystyle{\rm with}\quad\tan(2\theta)=\frac{2{\sf J}_{xz}}{{\sf
J}_{x}-{\sf J}_{z}}$

With time and maturation of the subject, the estimated parameters for Nd2Zr2O7
have evolved. Determinations of the $\tilde{\sf J}_{i}$ parameters are based
on the spin wave spectra measured at very low temperature in zero field [15,
31, 22] or applied field [21], while the angle $\theta$ is deduced from the
Curie-Weiss temperature [31] and/or the ordered AIAO magnetic moment [21, 22].
The sets of reported parameters are summarized in Table 1, where we have added
the parameters refined here for the Nd2(Zr1-xTix)2O7 sample (Sample #3) [23]
and have revisited the ones of Samples #1 and #2. From these values, two
interesting features stand out, which remain unexplained to date and should be
further explored to ascertain their relevance: (i) the larger the Néel
temperature, the larger the ordered moment along z is. (ii) very similar
$\tilde{\sf J}_{i}$ parameters are obtained for the various samples, despite
differences with regard to the amount of impurities or to the ordering
parameters.

The ${\sf J}$ parameters lead to an ordered AIAO ground state, where the
pseudospins point along the (local) direction ${\bf\tilde{z}}$, turned around
the ${\bf z}$-axis towards the ${\bf x}$-axis by the angle $\theta$ [31]. As
shown by INS experiments, peculiar excitations are associated with this ground
state. They manifest as an inelastic spin ice like flat mode at an energy
$E_{0}\approx$ 70 $\mu$eV, above which spin wave branches disperse (See Figure
4a for Sample #1) [15]. This excitation spectrum is understood in the
framework of the dynamic fragmentation [31, 32] as the sum of a dynamic
divergence-free contribution, giving rise to the flat mode at $E_{0}$ and of a
dynamic curl-free contribution, which takes the form of the dispersing
branches. These spin waves correspond to the propagation of magnetically
charged excitations and have a spectral weight made of half-moons in
reciprocal space [15, 33].

Sample / Ref. |  $m_{\rm ord}~{}(\mu_{\rm B})$ |  $T_{\rm N}$ (mK) | Hamiltonian parameters (K) |  $\theta$ (rad)
---|---|---|---|---
|  |  |  $\tilde{{\sf J}}_{x}$ |  $\tilde{{\sf J}}_{y}$ |  $\tilde{{\sf J}}_{z}$ |  $\tilde{{\sf J}}_{x}/|\tilde{{\sf J}}_{z}|$ |
#1 |  $0.8\pm 0.05$ |  285 |  1.18 |  -0.03 |  -0.53 |  2.20 |  1.23
#2 |  $1.1\pm 0.1$ |  340 |  1.0 |  0.066 |  -0.5 |  2.0 |  1.09
#3 |  $1.19\pm 0.03$ |  375 |  0.97 |  0.21 |  -0.53 |  1.83 |  1.08
[22] |  1.26 |  400 |  1.05 |  0.16 |  -0.53 |  1.98 |  0.98
[31] |  1.4 |  - |  1.2 |  0.0 |  -0.55 |  2.18 |  0.83
Table 1: Ordered moment $m_{\rm ord}$ along ${\bf z}$, transition temperature
$T_{\rm N}$ and Hamiltonian parametrization reported in different studies.
$\tilde{{\sf J}}_{i}$ parameters for Sample #1 and from Ref. 31 were obtained
from fits of the INS data reported in Ref. 15 and, for Sample #2 in Ref. 21.
$m_{\rm ord}$ from Ref. 31 is a calculated value. The total Nd3+ magnetic
moment is estimated to $\approx 2.4~{}\mu_{\rm B}$ [17, 18]. Figure 1:  (a-b)
Magnetic instantaneous correlations in Sample #1 as a function of temperature.
The 10 K dataset has been subtracted as a background reference. Measurements
in (a) were symmetrized. (c) “Spin ice” moment $m_{1}$ and AIAO ordered moment
$m_{2}$ along ${\bf z}$ as a function of temperature [23]. Lines are guides to
the eye.

Instantaneous spin-spin correlations $S({\bf Q})$ were measured in Sample #1
as a function of temperature between 60 mK and 1 K through polarized neutron
scattering experiments and are displayed in Figure 1 [23]. These measurements
integrate over the neutron energy loss up to 3.5 meV, and thus contain both
elastic and inelastic signals. At 1 K, a spin ice pattern can barely be
observed, revealing the onset of a Coulomb phase. Upon cooling, the spin ice
pattern becomes clearly visible below 600 mK. At 450 mK, the magnetic moment
$m_{1}$ responsible for the spin ice-like diffuse scattering is estimated to
$2.05\pm 0.3~{}\mu_{\rm B}$ [23], to be compared to the 2.4 $\mu_{\rm B}$ full
Nd moment [17, 18]. In addition to this signal, below 800 mK, magnetic diffuse
scattering spots appear around $(220)$, $(113)$ and symmetry related
positions. Intensity on these positions increases with cooling until they
transform into Bragg peaks below $T_{\rm N}$ (285 mK in this sample)
characteristic of the AIAO phase. At low temperature, the corresponding
ordered magnetic moment is $m_{2}=m_{\rm ord}=0.8\pm 0.05~{}\mu_{\rm B}$ (from
diffraction measurements) and the magnetic contribution to the spin ice like
diffuse scattering amounts to $m_{1}=2\pm 0.3~{}\mu_{\rm B}$ [23] (see Figure
1c). The moment embedded in the spin ice correlations is thus at maximum
around $T_{\rm N}$ and slightly decreases at lower temperature. The diffuse
scattering observed in the vicinity of the Bragg peak positions above $T_{\rm
N}$ might arise from AIAO diffuse scattering just above the ordering
transition, but could also be a signature of deconfined excitations, as
proposed in Ref. 16.

Figure 2:  Spectral function $S(E)$ at different temperatures [23] measured at
a wavelength $\lambda=$ 8.5 Å, hence an energy resolution of 20 $\mu$eV: (a)
in Sample #1, integrated around ${\bf Q}=(0.8~{}0.8~{}0.8)$. The grey and red
lines correspond to the fitted incoherent elastic $I_{c}$ and inelastic
$S_{0}$ contributions respectively. (b) and (c): integrated over the measured
${\bf Q}$ range in Samples #2 (b) and #3 (c).

To determine the spectral profile contained in those magnetic correlations,
and especially the elastic or inelastic nature of the spin ice correlations
associated to $m_{1}$, INS measurements have been carried out on the three
aforementioned samples (see Table 1) as a function of temperature. To
highlight the possible presence of an inelastic flat mode, the ${\bf
Q}$-integrated spectral function $S(E)=\int d{\bf Q}S({\bf Q},E)$ was
computed. As this quantity is akin to a density of states, it enhances the
contribution of the flat modes contained in the spectrum. Figure 2 displays
$S(E)$ at different temperatures. As previously shown [15], the inelastic flat
band is clearly seen at low temperature. It is still visible at finite energy
close to $T_{\rm N}$ (320 mK for Sample #2 and 315 mK for Sample #3) and above
$T_{\rm N}$ (340 mK for Sample #1), yet broadens significantly upon warming.
At the highest temperatures, the signal looks almost quasielastic. To obtain a
quantitative insight into the temperature evolution of the mode, data were
fitted for the three samples (as shown in Figure 2a for Sample #1) to the
following model [23]:

$S(E)=b+I_{c}(E)+F(E,T)\times\left[S_{0}(E)+S_{1}(E)\right]$ (3)

$b$ is a flat background, $I_{c}(E)$ is a Gaussian function centered at zero
energy to account for the elastic incoherent scattering. $F(E,T)=1+n(E,T)$ is
the detailed balance factor ($n$ is the Bose-Einstein distribution).
$S_{0}(E)$ and $S_{1}(E)$ are two Lorentzian profiles, centered on the energy
$E_{0,1}$ and of intensity $I_{0,1}$, which represent respectively the flat
band and the dispersive mode typical of the spin wave spectrum in Nd2Zr2O7.

The determined positions $E_{0}$ and intensities $I_{0}$ are shown in Figure 3
as a function of the temperature normalized to $T_{\rm N}$ for the three
samples. As anticipated from Figure 2, with increasing temperature, the band
at $E_{0}$ softens and broadens while its intensity decreases. Nevertheless,
$E_{0}$ is non-zero at $T_{\rm N}$ and a persistent dynamical behaviour is
observed in all samples at and above $T_{\rm N}$, up to about $2T_{\rm N}$.
Finally, the width of the features above the flat mode makes it hard to
extract quantitative information from $S_{1}$. However, close examination of
$S({\bf Q},E)$ measured for Sample #1 above $T_{\rm N}$ at 340 mK (see Figure
4) shows that, in all investigated directions, besides a strong quasielastic
contribution (the inelastic mode being hardly discernible due to the energy
resolution and the color scale), weak features are present close to the
position of the low temperature dispersions. These spin wave fingerprints,
highlighted by arrows on Figure 4(b) and which manifest as a broad signal in
${\bf Q}$-cuts (Figure 4(c-d)), are not compatible with the excitation
spectrum expected in the presence of monopole creation and hopping [16].

Several striking features emerge from these measurements. INS experiments
reveal that the intensity $I_{0}$ of the inelastic spin ice mode decreases
when increasing temperature. Since D7 polarized experiments show that the full
spin ice correlations, elastic and inelastic, are strongest around $T_{\rm
N}$, the spin ice pattern observed above $T_{\rm N}$ must contain a new spin
ice contribution, likely elastic, and different from the inelastic mode at
$E_{0}$. This is confirmed by magnetization measurements, which point to
ferromagnetic-like correlations, as expected for spin ice [23]. This elastic
signal could not be directly identified in the elastic line of the IN5 data
[23] certainly due to background issues, but we should stress that the D7
polarization analysis is definitely the most appropriate way to remove
properly nuclear contributions and visualize small magnetic contributions.
These results thus point to the coexistence of two spin ice-like
contributions, an elastic and an inelastic one with different origins, and
different temperature dependences.

Figure 3:  (a) $E_{0}$ and (b) $I_{0}$ as a function of reduced temperature
$T/T_{\rm N}$ obtained from measurements on IN5 (dots - see Figure 2) and IN12
(triangles), together with results from Monte-Carlo (MC) calculations from
Ref. 16 (red dots) [23]. Lines are guides to the eye. The large $I_{0}$
experimental value when $E_{0}=0$ is the signature of the persistent
quasielastic contribution above $T_{\rm N}$.

These two contributions can be understood as the manifestation of the strong
competition at play between the pseudospin components of Nd. The negative
value of $\tilde{\sf J}_{z}$ (see Table 1) promotes an AIAO phase built on
${\tilde{\tau}}^{\tilde{z}}$ while the positive $\tilde{\sf J}_{x}$ favors a
Coulomb phase, similar to a spin ice phase, but built on
${\tilde{\tau}}^{\tilde{x}}$. For $\tilde{\sf J}_{x}/|\tilde{\sf
J}_{z}|\approx 2$, the value determined for Nd2Zr2O7, the former is stabilized
at low temperature and the latter at finite temperature, due to the large
entropy associated to the Coulomb phase. In these two regimes, spin ice
contributions are expected, an elastic one in the Coulomb phase at “high”
temperature, and an inelastic one in the AIAO ordered phase (accompanied by
dispersive excitations). Remarkably, the observable $\tau_{z}$, which
corresponds to the magnetic dipolar moment along the local $\langle
111\rangle$ axes, is a combination of the ${\tilde{\tau}}^{\tilde{x}}$ and
${\tilde{\tau}}^{\tilde{z}}$ components of the pseudospin. It thus holds the
two competing contributions (AIAO and Coulomb), which contrasts with the
conventional spin ice case where the $z$ component carries elastic spin ice
correlations only.

Figure 4:  INS spectra of Sample #1 along several high symmetry directions at
60 mK (a) and 340 mK (b), measured on IN5 with $\lambda$=6 Å. Red arrows
highlight the dispersive modes and their fingerprints above $T_{\rm N}$. (c-d)
Constant ${\bf Q}$-cuts at these two temperatures, integrated (c) along
$(hh0)$ and (d) along $(hh2)$.

The present results shed light on the manner in which the system evolves from
the “high” temperature Coulomb phase to the low temperature AIAO ordered
phase. At high temperature, around 1 K, the elastic spin ice signal
characteristic of the ${\tilde{\tau}}^{\tilde{x}}$ Coulomb phase appears
first. Upon cooling, the inelastic spin ice contribution along with dispersive
spin wave branches emerge above $T_{\rm N}$ and coexist with the elastic one.
They can naturally be considered as excitations stemming from the short-range
AIAO correlations of the ${\tilde{\tau}}^{\tilde{z}}$ component observed below
800 mK (see Figure 1).

The system enters the long-range AIAO ordered state at a temperature $T_{\rm
N}\approx 300$ mK. It corresponds to about $|\tilde{\sf J}_{x}|/4$, thus to a
temperature scale far above the one obtained theoretically for the
stabilization of the quantum regime of spin ice, which is estimated to a few
percents of the characteristic exchange interaction [34, 35]. This indicates
that the Coulomb phase remains in its thermal regime down to $T_{\rm N}$.
Surprisingly, the ordering temperature is larger than semi-classical Monte-
Carlo calculations predictions [16]. At $T_{\rm N}$, the excitation spectrum
is gapped, with the coexistence between the elastic spin ice component and the
inelastic spectrum typical of AIAO ordering. The lack of a spinon continuum
which would condense at $T_{\rm N}$ seems to preclude a transition driven by a
Higgs mechanism.

Deeper in the AIAO phase, the inelastic component - together with the Bragg
peaks - develops at the expense of the elastic component. The weak maximum of
the spin ice $m_{1}$ moment around $T_{\rm N}$ can thus be interpreted as due
to the rise of the inelastic spin ice mode along with the persistence of the
elastic contribution of the Coulomb phase. The coexistence of the elastic and
inelastic signals is consistent with MC calculations [16], even if, close to
$T_{\rm N}$, the two modes are less distinguishable in the experiments than in
the calculations due to the strong broadening of the inelastic mode. Although
some distribution is observed between the samples, the measured temperature
dependence of the inelastic spin ice mode, described by the energy $E_{0}(T)$
and intensity $I_{0}(T)$, is also consistent with calculations [23], despite a
slightly stronger inelastic component in experiments above $T_{\rm N}$ (see
Figure 3).

In summary, we find that with increasing temperature, the now well-established
flat spin ice band characteristic of the AIAO ground state in Nd2Zr2O7,
softens while its intensity decreases. The energy of this mode remains however
finite at and above $T_{\rm N}$ and becomes overdamped with increasing the
temperature further. At the same time, a new elastic spin ice component
appears. The nature of the correlated phase above $T_{\rm N}$ is thus highly
unconventional with the coexistence of an (elastic) Coulomb phase and
fragmented excitations, resulting from the competition between the different
terms of the Hamiltonian. Our observations support a picture where the AIAO
ordering arises from a thermal spin ice phase, a scenario which is well
accounted for by semi-classical MC calculations from Ref. 16, and is different
from the proposed Higgs transition. When increasing the ratio $\tilde{\sf
J}_{x}/|\tilde{\sf J}_{z}|$, reentrant behaviors are predicted [16] while the
system approaches a quantum spin liquid ground state [31]. Tuning the
parameters of the Hamiltonian (2) with novel materials would thus be of high
interest to understand the unusual behavior of Nd2Zr2O7 and explore the
frontiers between thermal and quantum regimes.

###### Acknowledgements.

The work at the University of Warwick was supported by EPSRC, UK through Grant
EP/T005963/1. M. L. and S.P. acknowledge financial support from the French
Federation of Neutron Scattering (2FDN). M. L. acknowledges financial support
from Université Grenoble-Alpes (UGA). M.L., E.L. and S.P. acknowledge
financial support from ANR, France, Grant No. ANR-19-CE30-0040-02. S.P. and
E.L. acknowledge F. Damay for helpful remarks and J. Xu for providing the data
of his calculations. E.L. acknowledges C. Paulsen for the use of his
magnetometers.

## References

  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila, eds.,  _Introduction to Frustrated Magnetism_ (Springer-Verlag, Berlin, 2011).
  * Gardner  _et al._ [2010] J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Magnetic pyrochlore oxides, Rev. Mod. Phys. 82, 53 (2010).
  * Henley [2010] C. L. Henley, The “Coulomb phase” in frustrated systems, Annu. Rev. Condens. Matter Phys. 1, 179 (2010).
  * Harris  _et al._ [1997] M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Geometrical frustration in the ferromagnetic pyrochlore Ho2Ti2O7, Phys. Rev. Lett. 79, 2554 (1997).
  * Isakov  _et al._ [2004] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Dipolar spin correlations in classical pyrochlore magnets, Phys. Rev. Lett. 93, 167204 (2004).
  * Henley [2005] C. L. Henley, Power-law spin correlations in pyrochlore antiferromagnets, Phys. Rev. B 71, 014424 (2005).
  * Castelnovo  _et al._ [2008] C. Castelnovo, R. Moessner, and S. L. Sondhi, Magnetic monopoles in spin ice, Nature 451, 42 (2008).
  * Hermele  _et al._ [2004] M. Hermele, M. P. A. Fisher, and L. Balents, Pyrochlore photons: The $U(1)$ spin liquid in a $S=1/2$ three-dimensional frustrated magnet, Phys. Rev. B 69, 064404 (2004).
  * Shannon  _et al._ [2012] N. Shannon, O. Sikora, F. Pollmann, K. Penc, and P. Fulde, Quantum ice: A quantum Monte Carlo study, Phys. Rev. Lett. 108, 067204 (2012).
  * Benton  _et al._ [2012] O. Benton, O. Sikora, and N. Shannon, Seeing the light: Experimental signatures of emergent electromagnetism in a quantum spin ice, Phys. Rev. B 86, 075154 (2012).
  * Gingras and McClarty [2014] M. J. P. Gingras and P. A. McClarty, Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets, Rep. Prog. Phys. 77, 056501 (2014).
  * Sibille  _et al._ [2018] R. Sibille, N. Gauthier, H. Yan, M. Ciomaga Hatnean, J. Ollivier, B. Winn, G. Balakrishnan, M. Kenzelmann, N. Shannon, and T. Fennell, Experimental signatures of emergent quantum electrodynamics in a quantum spin ice, Nature Phys. 14, 711 (2018).
  * Savary and Balents [2012] L. Savary and L. Balents, Coulombic quantum liquids in spin-1/2 pyrochlores, Phys. Rev. Lett. 108, 037202 (2012).
  * Hao  _et al._ [2014] Z. Hao, A. G. R. Day, and M. J. P. Gingras, Bosonic many-body theory of quantum spin ice, Phys. Rev. B 90, 214430 (2014).
  * Petit  _et al._ [2016] S. Petit, E. Lhotel, B. Canals, M. Ciomaga Hatnean, J. Ollivier, H. Mutka, E. Ressouche, A. R. Wildes, M. R. Lees, and G. Balakrishnan, Observation of magnetic fragmentation in spin ice, Nature Phys. 12, 746 (2016).
  * Xu  _et al._ [2020] J. Xu, O. Benton, A. T. M. N. Islam, T. Guidi, G. Ehlers, and B. Lake, Order out of a Coulomb phase and Higgs transition: Frustrated transverse interactions in Nd2Zr2O7, Phys. Rev. Lett. 124, 097203 (2020).
  * Lhotel  _et al._ [2015] E. Lhotel, S. Petit, S. Guitteny, O. Florea, M. Ciomaga Hatnean, C. Colin, E. Ressouche, M. R. Lees, and G. Balakrishnan, Fluctuations and all-in–all-out ordering in dipole-octupole Nd2Zr2O7, Phys. Rev. Lett. 115, 197202 (2015).
  * Xu  _et al._ [2015] J. Xu, V. K. Anand, A. K. Bera, M. Frontzek, D. L. Abernathy, N. Casati, K. Siemensmeyer, and B. Lake, Magnetic structure and crystal-field states of the pyrochlore antiferromagnet Nd2Zr2O7, Phys. Rev. B 92, 224430 (2015).
  * Pekker and Varma [2015] D. Pekker and C. M. Varma, Amplitude / Higgs modes in condensed matter physics, Annu. Rev. Condens. Matter Phys. 6, 269 (2015).
  * Chang  _et al._ [2012] L.-J. Chang, S. Onoda, Y. Su, Y.-J. Kao, K.-D. Tsuei, Y. Yasui, K. Kakurai, and M. R. Lees, Higgs transition from a magnetic Coulomb liquid to a ferromagnet in Yb2Ti2O7, Nature Commun. 3, 992 (2012).
  * Lhotel  _et al._ [2018] E. Lhotel, S. Petit, M. Ciomaga Hatnean, J. Ollivier, H. Mutka, E. Ressouche, M. R. Lees, and G. Balakrishnan, Evidence for dynamic kagome ice, Nature Commun. 9, 3786 (2018).
  * Xu  _et al._ [2019] J. Xu, O. Benton, V. K. Anand, A. T. M. N. Islam, T. Guidi, G. Ehlers, E. Feng, Y. Su, A. Sakai, P. Gegenwart, and B. Lake, Anisotropic exchange hamiltonian, magnetic phase diagram, and domain inversion of Nd2Zr2O7, Phys. Rev. B 99, 144420 (2019).
  * [23] See Supplemental Material for details on single crystal growth, neutron diffraction, polarized neutron experiments, inelastic neutron experiments, spin dynamics in Ti doped sample, analysis of classical dynamics results, magnetization, which includes Refs. [24, 25, 26, 27, 28].
  * Ciomaga Hatnean  _et al._ [2015] M. Ciomaga Hatnean, M. R. Lees, and G. Balakrishnan, Growth of single-crystals of rare-earth zirconate pyrochlores, ${Ln}_{2}$Zr2O7 (with ${Ln}=$La, Nd, Sm, and Gd) by the floating zone technique, J. Cryst. Growth 418, 1 (2015).
  * Ciomaga Hatnean  _et al._ [2016] M. Ciomaga Hatnean, C. Decorse, M. R. Lees, O. A. Petrenko, and G. Balakrishnan, Zirconate pyrochlore frustrated magnets: crystal growth by the floating zone technique, Crystals 6, 79 (2016).
  * Rodríguez-Carvajal [1993] J. Rodríguez-Carvajal, Recent advances in magnetic structure determination by neutron powder diffraction, Physica B 192, 55 (1993).
  * Arnold and et al. [2014] O. Arnold and et al., Mantid - data analysis and visualization package for neutron scattering and $\mu$SR experiments, Nucl. Instrum. Methods Phys. Res. Sect. A 764, 156 (2014).
  * Ewings  _et al._ [2016] R. A. Ewings, A. Buts, M. D. Lee, J. van Duijn, I. Bustinduy, and T. G. Perring, Horace: Software for the analysis of data from single crystal spectroscopy experiments at time-of-flight neutron instruments, Nucl. Instrum. Methods Phys. Res. Sect. A 834, 132 (2016).
  * Paulsen [2001] C. Paulsen, Dc magnetic measurements, in  _Introduction to Physical Techniques in Molecular Magnetism: Structural and Macroscopic Techniques - Yesa 1999_ , edited by F. Palacio, E. Ressouche, and J. Schweizer (Servicio de Publicaciones de la Universidad de Zaragoza, 2001) p. 1.
  * Huang  _et al._ [2014] Y.-P. Huang, G. Chen, and M. Hermele, Quantum spin ices and topological phases from dipolar-octupolar doublets on the pyrochlore lattice, Phys. Rev. Lett. 112, 167203 (2014).
  * Benton [2016] O. Benton, Quantum origins of moment fragmentation in Nd2Zr2O7, Phys. Rev. B 94, 104430 (2016).
  * Brooks-Bartlett  _et al._ [2014] M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, Magnetic-moment fragmentation and monopole crystallization, Phys. Rev. X 4, 011007 (2014).
  * Yan  _et al._ [2018] H. Yan, R. Pohle, and N. Shannon, Half moons are pinch points with dispersion, Phys. Rev. B 98, 140402(R) (2018).
  * Savary and Balents [2013] L. Savary and L. Balents, Spin liquid regimes at nonzero temperature in quantum spin ice, Phys. Rev. B 87, 205130 (2013).
  * Huang  _et al._ [2018] C.-J. Huang, Y. Deng, Y. Wan, and Z. Y. Meng, Dynamics of topological excitations in a model quantum spin ice, Phys. Rev. Lett. 120, 167202 (2018).

Spin dynamics and unconventional Coulomb phase in Nd2Zr2O7
Supplementary Material

##  I Single crystal growth

Single crystals of Nd2(Zr1-xTix)2O7 ($x=0$ and 0.025) were grown by the
floating zone method, using a four-mirror xenon arc lamp optical image furnace
[24, 25]. A summary of the conditions used for each crystal growth is given in
Table S1.

Crystal | Sample label |  Lattice parameter | Growth rate |  Growth atmosphere, | Feed / seed
---|---|---|---|---|---
|  Å | (mm/h) |  pressure | rotation rate (rpm)
Nd2Zr2O7 | Sample #1 |  $10.66\pm 0.02$ | 12.5 |  Air, ambient | 15 / 30
Nd2Zr2O7 | Sample #2 |  $10.66\pm 0.04$ | 15 |  Air, ambient | 20 / 25
Nd2(Zr1-xTix)2O7 | Sample #3 |  $10.65\pm 0.02$ | 10 |  Air, ambient | 15 / 5
Table S1:  Summary of the samples with their crystal growth conditions. The
lattice parameters were obtained at 6 K on neutron diffractometers (Sample #1
and #3) and triple axis spectrometers (all samples).

Two different pure Nd2Zr2O7 samples had to be used in inelastic neutron
scattering experiments, because the first one broke when warming up the
dilution fridge after an experiment.

##  II Characterization of the Ti substituted sample (Sample #3)

Figure S1:  (a) Refinement of the crystal neutron structure factors at 6 K,
giving a refined Ti content equal to 2.4 %. (b) Measured intensity on the
magnetic peaks (220), (113), (351) and (260), and symmetry related peaks at 60
mK, obtained from the difference with the 6 K data, and compared to the
refined intensity.

We have studied a substituted sample, in which a small content of Zr is
replaced by Ti, slightly shrinking the structure. The nominal composition of
the studied sample is 2.5 % of Ti atoms. As shown below, this substitution
only slightly affects the magnetic properties of the magnetic Nd3+ sublattice
and the low temperature properties are qualitatively the same.

The value of the Ti content was refined by neutron diffraction, thanks to the
significant contrast between Zr and Ti. A series of Bragg peak intensities was
collected at 6 K on the single crystal neutron diffractometer D23 (CEA CRG-
ILL). The data are in agreement with the pyrochlore structure ($Fd{\bar{3}}m$
space group), with a lattice parameter of 10.65 Å and the 48f oxygen atoms at
the position $x_{\rm 48f}=0.336$. The Ti content is found to be 2.4 %. The
Fullprof refinement [26] of the structure factor is shown on Fig. S1(a).

The Néel temperature was determined from very low temperature magnetization
measurements, and found to be $T_{\rm N}=375$ mK.

The magnetic contribution raises below $T_{\rm N}$ in neutron diffraction
measurements on top of the crystalline peaks. The Fullprof refinement below
$T_{\rm N}$ confirms the same “all in - all out” (AIAO) magnetic structure as
in the pure sample (Fig. S1(b)). At 60 mK, the refined ordered Nd3+ magnetic
moment is $1.19\pm 0.03~{}\mu_{\rm B}$.

##  III Measurements and analysis of polarized neutron scattering experiments

Figure S2:  D7 data and analysis. Panels (a,b,c) display respectively the raw
data at 60 mK, the symmetrized and noise filtered data (10 K data have been
subtracted). Panel (d) shows the spin ice magnetic scattering function
calculated for a lattice containing 5488 spins. Panel (e) shows the same
calculation assuming random $\pm 1$ Ising spins. It serves as a background
reference, and is subtracted from (d), just as the 10 K data are subtracted
from the low temperature data.

In polarized neutron experiments carried out at D7 (ILL, France), we used the
$P_{z}$ polarization mode, ${\bf z_{Q}}$ being the axis normal to the
scattering plane and parallel to $[1\bar{1}0]$. We measured $N+I^{(z)}$ in the
NSF channel, and $I^{(y)}$ in the SF channel. Here we use conventional
notations: $N$ is the crystalline structure factor, while $I^{(y)}$ and
$I^{(z)}$ denote the spin-spin correlation functions between spin components
parallel to ${\bf y_{Q}}$ and ${\bf z_{Q}}$ respectively. The ${\bf y_{Q}}$
axis lies within the scattering plane, perpendicular both to ${\bf Q}$ and
${\bf z_{Q}}$. We used a wavelength $\lambda=4.85$ Å. The sample was rotated
by steps of 1 degree, and 2 positions of the detector bank have been combined.
Standard corrections (vanadium and quartz) have been processed. Finally, in
order to eliminate any background contribution, the data recorded at 10 K have
been subtracted from the data taken at lower temperatures.

In Figure 1 of the main article, the data have been symmetrized (for the three
lowest temperatures) while the noise was reduced by a mean filtering. This
image processing based treatment tends to reduce the variation between one
pixel and the next. The idea of mean filtering is to replace each pixel value
with the average value of its neighbors, including itself. This has the effect
of eliminating pixel values which are unrepresentative of their surroundings.
For the sake of illustration, Figure S2(a-c) shows the different steps of this
processing for the 60 mK data.

Unfortunately, it was not possible to determine intensities in absolute units
from the D7 measurements. To determine the magnetic moment responsible for the
spin ice-like diffuse scattering, we had to proceed in an alternative manner.
To this end, we carried out a series of calculations, assuming a “theoretical
sample crystal” consisting of Ising spins (of length unity) located at the
rare earth sites of a pyrochlore lattice of size $L$. We have considered $n$
spin ice configurations generated on this pyrochlore lattice (via a Monte-
Carlo algorithm) and computed the average structure factor from the obtained
magnetic moment ${\bf m}_{i,a=x,y,z}$ at each site $i$. The total magnetic
neutron intensity, proportional to the spin-spin correlation, is calculated
as:

$I^{(y)}+I^{(z)}=\sum_{i,j}\sum_{a,b=x,y,z}m_{i,a}\left(\delta_{ab}-\frac{{\bf
Q}_{a}{\bf Q}_{b}}{{\bf Q}^{2}}\right)m_{j,b}~{}e^{i{\bf Q}.({\bf R}_{i}-{\bf
R}_{j})}$

and the intensity in the SF $P_{z}$ mode, corresponding to the data, is given
by $\displaystyle I^{(y)}=\sum_{i,j}{\bf m}_{i}.{\bf y}_{\bf Q}\ {\bf
m}_{j}.{\bf y}_{\bf Q}\ e^{i{\bf Q}.({\bf R}_{i}-{\bf R}_{j})}$.

Noteworthy, the Monte-Carlo sampling was checked using the analytical method
proposed by C.L. Henley [6]. The spin-spin correlation function (per unit
cell) is written as:

$I^{(y)}=4~{}t_{0}^{2}~{}\sum_{a=x,y,z}M^{T}_{a}\left[I-E(E^{+}E)^{-1}E^{+}\right]M_{a}$

$E$ is a 2-column matrix and $M_{a=x,y,z}$ is a collection of 1-column vectors
containing the coordinates $a=x,y,z$ of the four magnetic moments ${\bf
m}_{i}$ belonging to a given tetrahedron:

$\displaystyle E$ $\displaystyle=\begin{pmatrix}e^{-i\pi{\bf
Q.u_{1}}}&e^{i\pi{\bf Q.u_{1}}}\\\ e^{-i\pi{\bf Q.u_{2}}}&e^{i\pi{\bf
Q.u_{2}}}\\\ e^{-i\pi{\bf Q.u_{3}}}&e^{i\pi{\bf Q.u_{3}}}\\\ e^{-i\pi{\bf
Q.u_{4}}}&e^{i\pi{\bf Q.u_{4}}}\\\ \end{pmatrix}$ $\displaystyle M_{a}$
$\displaystyle=\begin{pmatrix}M_{1,a}\\\ M_{2,a}\\\ M_{3,a}\\\
M_{4,a}\end{pmatrix}\quad{\rm where}\quad M_{i,a}=m_{i,a}-\frac{{\bf
m}_{i}.{\bf Q}}{Q^{2}}Q_{a}-\left(\sum_{b}\left(m_{i,b}-\frac{{\bf m}_{i}.{\bf
Q}}{Q^{2}}Q_{b}\right).z_{{\bf Q},b}\right)z_{{\bf Q},a}\quad{\rm
for}\quad{i=1,2,3,4}$

The four moments are defined as ${\bf m}_{1}=(c,c,c)$, ${\bf
m}_{2}=(-c,-c,c)$, ${\bf m}_{3}=(-c,c,-c)$, ${\bf m}_{4}=(c,-c,-c)$ with
$c=1/\sqrt{3}$ and attached to a tetrahedron. The ${\bf u_{i}}$ are vectors
pointing towards the four corners of a tetrahedron: ${\bf u_{1}}=(d,d,d)$,
${\bf u_{2}}=(-d,-d,d)$, ${\bf u_{3}}=(-d,d,-d)$, ${\bf u_{4}}=(d,-d,-d)$ and
$d=1/4$. Proper normalization condition imposes $t_{0}^{2}=2$.

In the same way, we computed the structure factor $I^{(y)}_{\rm rd}$ assuming
that the Ising spins have purely random $\pm 1$ values. This quantity was used
as a background reference, just as the 10 K measurement was used as described
above. We eventually considered the case where the spins are arranged in an
“all in – all out” (AIAO) ordering, leading to magnetic Bragg peaks with a
structure factor denoted hereafter $I^{(y)}_{\rm AIAO}$.

First we calculated the integrated intensity around $(11\bar{3})$ from the
experimental data at different temperatures, yielding $I^{\rm exp}_{\rm
AIAO}(T)$. On the other hand, the same quantity was determined from
$I^{(y)}_{\rm AIAO}$, yielding $I^{\rm c}_{\rm AIAO}$. Since the actual value
of the ordered AIAO moment at low temperature $m_{\rm AIAO}$ (called $m_{2}$
in the main text) is precisely known from diffraction measurements (D23), we
introduced the normalization factor $c$:

$c=\frac{I^{\rm exp}_{\rm AIAO}}{m_{\rm AIAO}^{2}~{}I^{\rm c}_{\rm AIAO}}$

In a second step, we have computed the experimental integrated intensity
within a box delineating the arm along $(hh{\bar{h}})$, yielding $\Delta
I^{\rm exp}_{\rm arm}$. To obtain an estimate of the moment $m_{\rm SI}$
(called $m_{1}$ in the main text) involved in the spin ice component, which
reflects the evolution seen on the maps presented in Figure 1 of the main
text, we proposed to compare $\Delta I^{\rm exp}_{\rm arm}$ to $\Delta I^{\rm
c}_{\rm arm}=m_{\rm SI}^{2}(I^{(y)}-I^{(y)}_{\rm rd})$. The estimation is then
made quantitative by looking for $m_{\rm SI}$ such that :

$\Delta I^{\rm exp}_{\rm arm}(T)=c\times\Delta I^{\rm c}_{\rm arm}$

The obtained values are listed in Table S2.

Other calculation methods have been tested and lead to similar results in
terms of absolute values and evolution with temperature. Furthermore, using
$m_{\rm AIAO}=0.8\pm 0.05~{}\mu_{\rm B}$, the accuracy on $m_{\rm SI}$ is
estimated to $\pm 0.3~{}\mu_{\rm B}$. This analysis confirms that the spin ice
pattern has a weak maximum close to $T_{\rm N}$ and persists up to 600 mK,
i.e. far above $T_{\rm N}$.

$T$ (mK) |  $I^{\rm exp}_{\rm arm}(T)$ |  $m_{\rm SI}$ ($\mu_{\rm B}$)
---|---|---
60 |  0.027 | 1.97
235 |  0.035 | 2.25
450 |  0.030 | 2.05
600 |  0.020 | 1.70
800 |  0.015 | 1.46
1000 |  0.0095 | 1.16
Table S2:  Values of the spin moment $m_{\rm SI}$ vs temperature, determined
from the procedure described in the text. Note that the temperature of 235 mK
was estimated from the amplitude of the magnetic Bragg peaks (the thermometer
indicated 300 mK). For 450 mK, no precise determination of the sample
temperature could be done, but, in the absence of magnetic Bragg peaks, the
temperature was definitely above $T_{\rm N}$.

##  IV Time of flight inelastic scattering measurements

Inelastic neutron scattering experiments were carried out on the IN5 disk
chopper time of flight spectrometer (ILL, France). A good compromise between
flux, energy resolution and accessible ${\bf Q}$ space was obtained with a
wave length $\lambda=6$ or 6.5 Å. However, to ensure a better energy
resolution $\Delta E=20~{}\mu$eV, necessary to fully resolve the dynamic spin
ice mode at $E_{0}$, experiments were also conducted with $\lambda=8.5$ Å. The
data were processed with the Mantid [27] and horace [28] softwares,
transforming the recorded time of flight, sample rotation and scattering angle
into energy transfer and ${\bf Q}$-wave vectors. The offset of the sample
rotation was determined based on the Bragg peak positions. In all the
experiments, the sample was rotated in steps of 1 degree and the counting time
was about 10 minutes per sample position.

It should be noticed that a very long thermalization time was systematically
necessary to cool down the sample to the lowest temperature. In addition, we
realized that when warming up from the lowest temperature, the sample
temperature was not necessarily the same as the temperature indicated by the
thermometer. For this reason, when possible (depending on the ratio between
the resolution and the temperature), we have refined the “true” temperature by
fitting the negative energy part of the spectra. It leads to the temperatures
indicated on Figures 2, 3 and 4 of the main text, which are quite different
from the thermometer temperatures. These temperatures are summarized in Table
S3.

Sample |  Thermometer |  Estimated
---|---|---
|  temperature |  temperature
Sample #1 |  450 mK |  $341\pm 100$ mK
Sample #2 |  60 mK |  $323\pm 78$ mK
|  300 mK |  $313\pm 68$ mK
|  450 mK |  $444\pm 111$ mK
Sample #3 |  275 mK |  $242\pm 35$ mK
|  350 mK |  $317\pm 38$ mK
Table S3: Estimated effective temperatures in the different experiments
performed on IN5. Figure S3:  Constant ${\bf Q}$-cuts at two temperatures, the
base temperature of 60 mK (blue) and above $T_{\rm N}$ (red) and which clearly
show the vestiges of spin waves.

Constant ${\bf Q}$-cuts from the data have been performed at the base
temperature (typically 60 mK) and above $T_{\rm N}$ (340 mK), to clearly show
the persistence of the spin wave signal above $T_{\rm N}$. These cuts,
displayed in Figure S3 are along $(1,1,\ell)$, $(2,2,\ell)$, $(h,h,0)$ and
$(h,h,2)$. Two of them are reproduced in the main text.

This residual spin wave signal above $T_{\rm N}$ is not expected in
conventional three dimensional paramagnets, in the absence of magnetic
frustration. It thus would not be observed in a standard “all in – all out”
antiferromagnet, which is predicted to behave classically close to the
antiferromagnetic transition. The persistence of the spin wave signal (and of
“all in – all out” diffuse scattering) quite far above $T_{\rm N}$ in Nd2Zr2O7
thus points out the unconventional nature of the magnetism in this compound
and is likely related to the strong competition at play with the Coulomb phase
observed above $T_{\rm N}$.

##  V Inelastic scattering measurements on a triple axis spectrometer

The temperature dependence of the spin dynamics in Sample #1 was also
investigated on the cold TAS spectrometer IN12 (ILL, France). Scans at
specific ${\bf Q}$ positions (0.5 0.5 2), (0 0 2.5) and (1.8 1.8 0) have been
performed at different temperatures ranging from 50 up to 800 mK. Those
positions were chosen since they probe different regions with respect to the
dispersion. (0.5 0.5 2) essentially probes the flat spin ice band, (1.8 1.8 0)
is sensitive to the zone boundary dispersive spin wave mode and (0.5 0.5 2) is
somehow intermediate. A final wave vector $k_{f}=1.05$ Å-1 was used (in
combination with nitrogen cooled Be filter) to ensure the best energy
resolution, $\Delta E=50~{}\mu$eV. A magnetic field was also applied along
$[1\bar{1}0]$. After correction from the detailed balance factor, we computed
the difference between data taken a given temperature $T$ and the 800 mK data.
Where applicable, we subtracted the data obtained at the same temperature but
under a 1 T magnetic field. We could then extract the energy and intensity of
the inelastic mode in the same way as for TOF measurements. The temperatures
below $T_{\rm N}$ were estimated from the intensity of the (220) magnetic
Bragg peak.

##  VI Spin dynamics in the Ti substituted sample (Sample #3)

###  VI.1 Determination of the parameters

Inelastic neutron scattering data carried out at IN5 (ILL) on a single crystal
sample show little evolution compared to the pure sample. The inelastic flat
spin ice mode is observed at $E_{0}\approx 70~{}\mu$eV, while the dispersing
mode stemming from the pinch point positions unfolds towards the zone centers,
for instance $(220)$ or $(113)$. This is illustrated in Figure S4, which shows
the dispersion along several reciprocal directions at 45 mK.

Figure S4:  Top: INS data taken at IN5 at 45 mK on the Sample #3 along high
symmetry directions. Black and white dots are the energies $E_{0}$ and $E_{1}$
respectively, fitted according to the procedure described in the main text
(see also equation (S1). Bottom: Spin wave calculations performed with the
parameters given in Table 1 (main text).

To determine the parameters of the XYZ Hamiltonian

${\cal H}_{\rm XYZ}=\sum_{\langle i,j\rangle}\left[{\tilde{\sf
J}_{x}}\tilde{\tau}^{\tilde{x}}_{i}\tilde{\tau}^{\tilde{x}}_{j}+{\tilde{\sf
J}_{y}}\tilde{\tau}^{\tilde{y}}_{i}\tilde{\tau}^{\tilde{y}}_{j}+{\tilde{\sf
J}_{z}}\tilde{\tau}^{\tilde{z}}_{i}\tilde{\tau}^{\tilde{z}}_{j}\right]$

(see also equation (2) of the main text), we use analytic calculations giving
the energy of the spin ice band [31]:

$E_{0}=\sqrt{(3|\tilde{{\sf J}}_{z}|-\tilde{{\sf J}}_{x})(3|\tilde{{\sf
J}}_{z}|-\tilde{{\sf J}}_{y})}$

as well as the energy of the dispersive modes at some high symmetry ${\bf Q}$
vectors [22]:

$\displaystyle{\bf Q}=(110),(112)$ $\displaystyle,~{}$
$\displaystyle\Delta_{2}=\sqrt{(3|\tilde{{\sf J}}_{z}|+\tilde{{\sf
J}}_{x})(3|\tilde{{\sf J}}_{z}|+\tilde{{\sf J}}_{y})}$ $\displaystyle{\bf
Q}=(220),(113)$ $\displaystyle,~{}$
$\displaystyle\Delta_{3}=3\sqrt{(|\tilde{{\sf J}}_{z}|+\tilde{{\sf
J}}_{x})(|\tilde{{\sf J}}_{z}|+\tilde{{\sf J}}_{y})}$

Simulations have then been performed to reproduce the data with the cefwave
software developed at LLB using the Hamiltonian (1):

${\cal H}=\sum_{\langle i,j\rangle}\left[{\sf
J}_{x}\tau^{x}_{i}\tau^{x}_{j}+{\sf J}_{y}\tau^{y}_{i}\tau^{y}_{j}+{\sf
J}_{z}\tau^{z}_{i}\tau^{z}_{j}+{\sf
J}_{xz}(\tau^{x}_{i}\tau^{z}_{j}+\tau^{z}_{i}\tau^{x}_{j})\right]$

The ground state configuration is first determined by solving this Hamiltonian
at the mean field level, where the expectation values
$\langle\tau^{x,y,z}_{j}\rangle$ are determined in a self-consistent manner.
Spin wave calculations are performed using a generalized susceptibility
approach out of the obtained configurations. Finally, the neutron cross
section is calculated from $\tau^{z}_{i}\tau^{z}_{j}$ correlations. Notably,
the simulations performed with the parameters of Table 1 (main text) reproduce
quite well the data, as shown in Figure S4.

###  VI.2 Temperature dependence of the spin dynamics

Inelastic data in Sample #3 were fitted in the whole measured ${\bf Q}$ space
using the model described in the main text:

$S({\bf Q},E)=b+I_{c}(E)+F(E,T)\times\left[S_{0}(E)+S_{1}(E)\right]$ (S1)

$b$ is a flat background (wavelength dependent), $I_{c}$ is a Gaussian
function centered at zero energy to represent the elastic incoherent
scattering. $F(E,T)=(1+n(E))$ is the detailed balance factor, and $S_{0}$ and
$S_{1}$ are two Lorentzian profiles which represent respectively the flat band
and the dispersive mode typical of the spin wave spectrum in Nd2Zr2O7.

Figure S5 shows the energy $E_{0}$ of the flat band at different temperatures
in the form of a map over the sector probed by TOF measurements. Figure S6
displays the intensities $I_{0}$ (panel a, upper row) and $I_{1}$ (panel b,
lower row). The map on the right of the same figure shows the energy $E_{1}$
of the dispersive spin wave mode. To check the overall consistency of the
fitting procedure, dashed lines visualize the directions of the scans reported
in Figure S4.

Note that for Sample #1, the fit was carried out at selected ${\bf Q}$ values
(${\bf Q}=(0.8\ 0.8\ 0.8)$, $(1.1\ 1.1\ 1.1)$, $(1/2\ 1/2\ 1/2)$, $(1/2\ 1/2\
3/2)$ and $(3/4\ 3/4\ 3/2)$).

Figure S5:  Temperature dependence of the flat spin ice band at $E_{0}$
deduced from the fit in Sample #3, as described in the main text. The portion
of $({\bf Q},E)$ space corresponds to the sector probed by TOF measurements
with $\lambda$=8.5 Å. Figure S6:  (a) Temperature dependence of the intensity
$I_{0}$ of the flat spin ice band from the fit in Sample #3, as described in
the main text. (b) shows the intensity $I_{1}$ of the dispersive mode, and (c)
shows its energy $E_{1}$. Dashed lines correspond to the directions of the
cuts shown in Figure S4. The portion of $({\bf Q},E)$ space corresponds to the
sector probed by TOF measurements with $\lambda=8.5$ Å.

Figure S7:  Temperature dependence of the intensity of the incoherent
scattering $I_{c}(E=0)$ in Sample #3. The portion of $({\bf Q},E)$ space
corresponds to the sector probed by TOF measurements with $\lambda=8.5$ Å.

Finally, aiming at identifying a possible elastic contribution with the spin
ice structure factor, Figure S7 shows the temperature evolution of the
intensity of the incoherent elastic contribution, i.e. the dominant
contribution $I_{c}(E=0)$ in the spectrum, to which the 45 mK map was
subtracted. Within experimental uncertainties, these maps are featureless and
no spin ice pattern can be clearly distinguished.

##  VII Analysis of classical dynamics results

Figure S8:  (a) From Ref. 16 (courtesy of J. Xu): Evolution of the gapped flat
mode for several temperatures (0.05, 0.1, 0.125, 0.15, 0.165, 0.175, 0.18,
0.185, 0.2 K) simulated using semi-classical molecular dynamics averaging over
Q from (0.1 0.1 0) to (0.9 0.9 0). (b-c) Temperature dependences fitted from
(a) (see equation S2): (b) Temperature dependence of the intensities $I_{0}$
and $I^{\prime}_{0}$ of the flat spin ice band and of the elastic
contributions respectively. (c) Temperature dependence of the energy $E_{0}$
of the flat spin ice band.

The main text of the present work compares the measured temperature dependence
of $E_{0}$ and $I_{0}$ in our three samples with Monte Carlo calculations
reported in Ref. 16. These calculations use effective exchange parameters from
Ref. 22, which are detailed in Table 1 of the main text. They give a Néel
temperature of 0.18 K.

From these calculations, as illustrated in the Figure 3 of Ref. 16, two
contributions are obtained: a spin ice elastic contribution which projects
onto the ${\bf z}$ axis with a factor $\sin^{2}\theta$, as well as spin waves
features characteristic of the AIAO phase, with especially a flat spin ice
band. To compare quantitatively these results with our data, those theoretical
curves have been fitted to two modes, following:

$I(E)=I_{0}~{}e^{-4\log
2(\frac{E-E_{0})}{\delta_{0}})^{2}}+\frac{I^{\prime}_{0}}{1+(\frac{E}{\delta^{\prime}_{0}})^{2}}$
(S2)

The result of this fit is illustrated in Figures S8(b) and S8(c), which
display respectively the intensity of the modes ($I_{0}$ and $I^{\prime}_{0}$)
and the position $E_{0}$ of the flat spin ice band. Interestingly, this energy
remains finite even above the calculated critical temperature $T_{\rm N}=0.18$
K.

##  VIII Correlations in magnetization measurements

Figure S9:  Magnetization $M$ vs $H/T$ measured for Sample #1 with the field
applied along: (a) [100], (b) [110] and (c) [111] at 1, 1.8 and 4.2 K.

In a paramagnet, isothermal magnetization curves scale as a function of the
variable $H/T$. The deviations to this scaling give insight into the nature of
the correlations that develop in the system. Upon cooling, if the
magnetization curve increases faster (slower) than the higher temperature
curve, it is the signature of the development of ferromagnetic-like
(antiferromagnetic-like) correlations.

We have plotted the magnetization as a function of $H/T$ for Nd2Zr2O7,
measured in Sample #1. As shown in Figure S9, the $M(H/T)$ curves rise above
the 4.2 K curves upon cooling down to 1 K, which indicates the development of
ferromagnetic correlations, consistent with the elastic spin ice picture
inferred from our neutron scattering measurements.

At 500 mK, the curves lie between the 1 and 4.2 K curves (not shown on the
figure for clarity), showing the development of antiferromagnetic correlations
compared to 1 K, but the persistence of global ferromagnetic correlations.
These antiferromagnetic correlations will end in the “all in – all out”
ordering at about 300 mK.