Order out of a Coulomb phase and Higgs transtion: frustrated transverse interactions of Nd2Zr2O7

TL;DR

This study reveals that Nd$_2$Zr$_2$O$_7$ transitions from an all-in-all-out magnetic order to a quantum spin ice state above its Néel temperature, driven by frustrated transverse interactions and characterized by Coulombic phase behavior.

Contribution

It demonstrates the existence of a Coulombic phase and Higgs transition in Nd$_2$Zr$_2$O$_7$ driven by frustrated transverse interactions, supported by neutron scattering and theoretical calculations.

Findings

Disappearance of gapped magnon excitations above $T_N$

Presence of strong spin ice correlations at zero energy

Support for Coulombic phase above $T_N$

Abstract

The pyrochlore material Nd$_2$Zr$_2$O$_7$ with an "all-in-all-out" (AIAO) magnetic order shows novel quantum moment fragmentation with gapped flat dynamical spin ice modes. The parameterized spin Hamiltonian with a dominant frustrated ferromagnetic transverse term reveals a proximity to a U(1) spin liquid. Here we study magnetic excitations of Nd$_2$Zr$_2$O$_7$ above the ordering temperature ($T_\text{N}$) using high-energy-resolution inelastic neutron scattering. We find strong spin ice correlations at zero energy with the disappearance of gapped magnon excitations of the AIAO order. It seems that the gap to the dynamical spin ice closes above $T_\text{N}$ and the system enters a quantum spin ice state competing with and suppressing the AIAO order. Classical Monte Carlo, molecular dynamics and quantum boson calculations support the existence of a Coulombic phase above $T_\text{N}$. Our findings relate the magnetic ordering of Nd$_2$Zr$_2$O$_7$ with the Higgs mechanism and provide explanations for several previously reported experimental features.

Figures

Equations

$${\cal H}_{XYZ}=\sum_{\langle ij\rangle}\big{[}\tilde{J}_{\tilde{x}}\tilde{\tau}^{\tilde{x}}_{i}\tilde{\tau}^{\tilde{x}}_{j}+\tilde{J}_{\tilde{y}}\tilde{\tau}^{\tilde{y}}_{i}\tilde{\tau}^{\tilde{y}}_{j}+\tilde{J}_{\tilde{z}}\tilde{\tau}^{\tilde{z}}_{i}\tilde{\tau}^{\tilde{z}}_{j}\big{]}.$$

$$\text{gap}=\sqrt{(3|{\tilde{J}}_{\tilde{z}}|-{\tilde{J}}_{\tilde{x}})(3|{\tilde{J}}_{\tilde{z}}|-{\tilde{J}}_{\tilde{y}})}.$$

$$\frac{\text{d}{\mathbf{\tau }}_i(t)}{\text{d}t}={\mathbf{\tau }}_i\times \sum _j{J}_{ij}{\mathbf{\tau }}_j(t)$$

$$\begin{gathered} S(\mathbf{q},\omega )=\frac{1}{2\pi }\sum _{ij}\sum _{\mu \nu } \\ \int dt{e}^{i\omega t}\left({\delta }_{\mu \nu }-\frac{q_{\mu }{q}_{\nu }}{q^2}\right)⟨m_i^{\mu }(\mathbf{q},t)m_i^{\mu }(-\mathbf{q},0)⟩ \end{gathered}$$

$${\bf m}_{i}(\mathbf{q})=\\ \frac{1}{\sqrt{N_{\sf uc}}}\sum_{\bf R}\exp(-i\mathbf{q}\cdot(\mathbf{R}+\mathbf{u}_{i}/2)){\bf m}_{i}(\mathbf{R}+\mathbf{u}_{i}/2)\$$

$${\mathbf{u}}_0=(1,1,1),{\mathbf{u}}_1=(1,-1,-1),$$

$${\mathbf{u}}_2=(-1,1,-1),{\mathbf{u}}_3=(-1,-1,1).$$

$$\displaystyle{\bf m}_{i}=g_{zz}\mu_{B}\left(\cos(\vartheta)\tau^{\tilde{z}}_{i}+\sin(\vartheta)\tau^{\tilde{x}}_{i}\right)\hat{{\bf z}}_{i}\$$

$$g_{zz}=5.0,\vartheta =0.98$$

$${\hat{\mathbf{z}}}_i=\frac{1}{\sqrt{3}}{\mathbf{u}}_i.$$

$$S(\mathbf{q},\omega)=g_{zz}^{2}\mu_{B}^{2}\frac{1}{2\pi}\sum_{ij}\left(\hat{\bf z}_{i}\cdot\hat{\bf z}_{j}-\frac{(\hat{\bf z}_{i}\cdot\mathbf{q})(\hat{\bf z}_{j}\cdot\mathbf{q})}{q^{2}}\right)\\ \int dt\ e^{i\omega t}\bigg{(}\cos^{2}(\vartheta)\langle\tau^{\tilde{z}}_{i}(\mathbf{q},t)\tau^{\tilde{z}}_{i}(-\mathbf{q},0)\rangle+\\ \cos(\vartheta)\sin(\vartheta)\langle\tau^{\tilde{z}}_{i}(\mathbf{q},t)\tau^{\tilde{x}}_{i}(-\mathbf{q},0)\rangle+\\ \cos(\vartheta)\sin(\vartheta)\langle\tau^{\tilde{x}}_{i}(\mathbf{q},t)\tau^{\tilde{z}}_{i}(-\mathbf{q},0)\rangle+\\ \sin^{2}(\vartheta)\langle\tau^{\tilde{x}}_{i}(\mathbf{q},t)\tau^{\tilde{x}}_{i}(-\mathbf{q},0)\rangle\bigg{)}.$$

$$S(\mathbf{q},\omega )=g_{zz}^2{\mu }_B^2\left({\cos }^2(\vartheta )S^{\tilde{z}\tilde{z}}(\mathbf{q},\omega )+{\sin }^2(\vartheta )S^{\tilde{x}\tilde{x}}(\mathbf{q},\omega )\right)$$

$$S^{\tilde{z}\tilde{z}}(\mathbf{q},\omega )=\frac{1}{2\pi }\sum _{ij}\left({\hat{\mathbf{z}}}_i\cdot {\hat{\mathbf{z}}}_j-\frac{({\hat{\mathbf{z}}}_i\cdot \mathbf{q})({\hat{\mathbf{z}}}_j\cdot \mathbf{q})}{q^2}\right)$$

$$\int dt{e}^{i\omega t}⟨{\tau }_i^{\tilde{z}}(\mathbf{q},t){\tau }_i^{\tilde{z}}(-\mathbf{q},0)⟩$$

$$S^{\tilde{x}\tilde{x}}(\mathbf{q},\omega )=\frac{1}{2\pi }\sum _{ij}\left({\hat{\mathbf{z}}}_i\cdot {\hat{\mathbf{z}}}_j-\frac{({\hat{\mathbf{z}}}_i\cdot \mathbf{q})({\hat{\mathbf{z}}}_j\cdot \mathbf{q})}{q^2}\right)$$

$$\int dt{e}^{i\omega t}⟨{\tau }_i^{\tilde{x}}(\mathbf{q},t){\tau }_i^{\tilde{x}}(-\mathbf{q},0)⟩$$

$$S^{\tilde{x}\tilde{x}}(\mathbf{q},\omega )=\frac{1}{2}\delta (\omega )\sum _{ij}\left({\hat{\mathbf{z}}}_i\cdot {\hat{\mathbf{z}}}_j-\frac{({\hat{\mathbf{z}}}_i\cdot \mathbf{q})({\hat{\mathbf{z}}}_j\cdot \mathbf{q})}{q^2}\right)P_{ij}(\mathbf{q})$$

$$P_{ij}(\mathbf{q})={\delta }_{ij}-\sum _{a=0,1,b=0,1}{E}_{ia}(\mathbf{q})M_{ab}(\mathbf{q})E_{bj}^{†}(\mathbf{q}).$$

$$E(\mathbf{q})=\left(\begin{array}{cc}\exp (-i\mathbf{q}\cdot {\mathbf{u}}_0/2)& \exp (i\mathbf{q}\cdot {\mathbf{u}}_0/2)\\ \exp (-i\mathbf{q}\cdot {\mathbf{u}}_1/2)& \exp (i\mathbf{q}\cdot {\mathbf{u}}_1/2)\\ \exp (-i\mathbf{q}\cdot {\mathbf{u}}_2/2)& \exp (i\mathbf{q}\cdot {\mathbf{u}}_2/2)\\ \exp (-i\mathbf{q}\cdot {\mathbf{u}}_3/2)& \exp (i\mathbf{q}\cdot {\mathbf{u}}_3/2)\end{array}\right)$$

$$M(\mathbf{q})=(E^{†}(\mathbf{q})\cdot E(\mathbf{q}){)}^{-1}.$$

$$\delta (\omega )\to G(\omega )=\frac{2}{f}\sqrt{\frac{\log (2)}{\pi }}{e}^{-\frac{4\log (2){\omega }^2}{f^2}}$$

$${\tau }_i^{\tilde{\pm }}={\tau }_i^{\tilde{y}}\pm i{\tau }_i^{\tilde{z}}$$

$$\displaystyle\mathcal{H}=\sum_{\langle ij\rangle}\bigg{[}\tilde{J}^{\tilde{x}}\tau^{\tilde{x}}_{i}\tau^{\tilde{x}}_{j}-\tilde{J}^{\tilde{\pm}}\left(\tau^{\tilde{+}}_{i}\tau^{\tilde{-}}_{j}+\tau^{\tilde{-}}_{i}\tau^{\tilde{+}}_{j}\right)+$$

$$\displaystyle\qquad\qquad\tilde{J}^{\tilde{\pm}\tilde{\pm}}\left(\tau^{\tilde{+}}_{i}\tau^{\tilde{+}}_{j}+\tau^{\tilde{-}}_{i}\tau^{\tilde{-}}_{j}\right)\bigg{]}$$

$$Q_{\mathbf{R}}={ϵ}_{\mathbf{R}}\sum _{i=0}^3{\tau }_i^{\tilde{x}}(\mathbf{R}+{\mathbf{u}}_i/2)$$

$${\tau }_i^{\widetilde{+}}(\mathbf{R}+{\mathbf{u}}_i/2)=\frac{1}{2}{\psi }_{\mathbf{R}}^{†}\exp (i{A}_{\mathbf{R},\mathbf{R}+{\mathbf{u}}_i}){\psi }_{\mathbf{R}+{\mathbf{u}}_i}$$

$${\tau }_i^{\widetilde{-}}(\mathbf{R}+{\mathbf{u}}_i/2)=\frac{1}{2}{\psi }_{\mathbf{R}}\exp (-i{A}_{\mathbf{R},\mathbf{R}+{\mathbf{u}}_i}){\psi }_{\mathbf{R}+{\mathbf{u}}_i}^{†}.$$

$$\frac{1}{2}⟨\exp (i{A}_{\mathbf{R},\mathbf{R}+{\mathbf{u}}_i})⟩\to \kappa .$$

$${\tau }_i^{\widetilde{+}}(\mathbf{R}+{\mathbf{u}}_i/2)=\kappa {\psi }_{\mathbf{R}}^{†}{\psi }_{\mathbf{R}+{\mathbf{u}}_i}$$

$${\tau }_i^{\widetilde{-}}(\mathbf{R}+{\mathbf{u}}_i/2)=\kappa {\psi }_{\mathbf{R}}{\psi }_{\mathbf{R}+{\mathbf{u}}_i}^{†}.$$

$${\psi }_{\mathbf{R}}^{†}{\psi }_{\mathbf{R}}=1.$$

$${\psi }_{\mathbf{R}}=\frac{1}{\sqrt{1+d_{\mathbf{R}}^{†}{d}_{\mathbf{R}}+b_{\mathbf{R}}^{†}{b}_{\mathbf{R}}}}(d_{\mathbf{R}}+b_{\mathbf{R}}^{†})$$

$$Q_{\mathbf{R}}=d_{\mathbf{R}}^{†}{d}_{\mathbf{R}}-b_{\mathbf{R}}^{†}{b}_{\mathbf{R}}.$$

$$\left(d_{\mathbf{R}}^{†}{d}_{\mathbf{R}}\right)\left(b_{\mathbf{R}}^{†}{b}_{\mathbf{R}}\right)=0.$$

Full text

Competing interactions and geometrical frustration support highly degenerate states which suppress conventional magnetic order and lead to novel emergent states Lacroix2011book . Classical spin ice (CSI) is a prominent example which is realized in (Dy/Ho)2Ti2O7 pyrochlores consisting of a network of corner-sharing tetrahedra Gardner2010rev ; Fennell2009 ; Morris2009 ; Harris1997 ; Henley2005 . In the CSI, the single-ion Ising anisotropy due to the crystal electric field (CEF) interactions frustrates the ferromagnetic (FM) interactions between the spins. This creates the “2-in-2-out” (2I2O) local constraint (ice rule) on the spin configuration leading to infinite degeneracy on the pyrochlore lattice Harris1997 . However, an antiferromagnetic (AFM) interaction is not frustrated resulting in a long-range “all-in-all-out” (AIAO) order Lacroix2011book .

With the introduction of transverse spin couplings, CSI transforms to quantum spin ice (QSI) allowing quantum tunnelling between the degenerate ice-rule states which realizes a type of U(1) quantum spin liquid Hermele2004 ; Onoda2010 ; Onoda2011 ; Shannon2012 ; Savary2012 ; Lee2012 ; Benton2012 ; Kato2015 ; Huang2018 . Recently, there has been a tremendous effort aimed at finding materials supporting QSI Gingras2014 . Several materials have been examined in the search for QSI including, for example, Yb2Ti2O7 Ross2011 , Pr2(Zr/Hf)2O7 Kimura2013 ; Sibille2018 , and Tb2Ti2O7 Gardner1999 , but the evidence so far is ambiguous, complicated by multi-phase competitions, structural defects and low-lying crystal field levels Jaubert2015 ; Yan2017 ; Martin2017 ; Wen2017 ; Hamid2007 ; Princep2015 .

As a QSI candidate, Nd2Zr2O7 is an ideal material for modelling, having well-isolated, Kramers, Ising anisotropic, dipolar-octupolar CEF ground doublets and a clean, well-ordered structure and it has been under intensive study Huang2014 ; Hatnean2015 ; Lhotel2015 ; Xu2015 ; Petit2016 ; Benton2016 ; Xu2016 ; Opherden2017 ; Lhotel2018 ; Xu2018 ; Xu2019 . Although it has an AIAO order as the ground state below $T_{\text{N}}\sim 0.4$ K Lhotel2015 ; Xu2015 , it shows remarkably persistent spin dynamics Xu2016 , quantum moment fragmentation, gapped dynamical spin ice Petit2016 ; Xu2019 , gapped kagome spin ice in (111) fields Lhotel2018 and quantum spin-1/2 chains in (110) fields Xu2018 . The parameterized anisotropic pseudospin-1/2 $XYZ$ Hamiltonian indicates that the un-frustrated AFM $\tilde{J}^{\tilde{z}}\approx-0.046$ meV induces the AIAO order, though the FM transverse $\tilde{J}^{\tilde{x}}\approx 0.09$ meV is approximately twice as strong as $|\tilde{J}^{\tilde{z}}|$ Xu2019 ; Lhotel2018 ; Benton2016 ; Petit2016 ; Huang2014 . This is a result of the frustration for the FM $\tilde{J}^{\tilde{x}}$ term. According to linear spin wave theory, the gap to the flat spin ice modes closes at $\tilde{J}^{\tilde{x}}\//|\tilde{J}^{\tilde{z}}|=3$ Benton2016 where spin ice configurations and AIAO order have the same energy. Classically, this signals the formation of an extensive ground-state manifold with icelike character for the $\tilde{x}$ component of the spin and the mixing of these states by quantum fluctuations caused by $\tilde{J}^{\tilde{y}}$ and $\tilde{J}^{\tilde{z}}$ stabilizes a U(1) spin liquid with dynamic emergent gauge fields Benton2016 . It was pointed out also that if there is a gapless Coulomb phase above $T_{\text{N}}$, the ordering transition would be a candidate for a Higgs transition Benton2016 .

In this paper, we study the magnetic excitations of Nd2Zr2O7 above $T_{\text{N}}$ using high-energy-resolution inelastic neutron scattering. We find that the gapped magnon excitations disappear and the pinch point pattern characteristic for spin ice correlations is still present but at zero energy. The gapless nature of the strong spin ice correlations points to a QSI state above $T_{\text{N}}$, which supports the theoretical speculations Benton2016 . Classical Monte Carlo simulations (MC) indicate that the system does not simply become paramagnetic with short-range AIAO correlations but enters a intriguing state with ice correlations for the $\tau^{\tilde{x}}$ spin components. Calculations for a finite temperature QSI with propagating monopole excitations, are qualitatively consistent with the scattering data. Our findings shed light on the anomalously slow spin dynamics probed by muon spin relaxation commonly seen in frustrated magnets Xu2016 and explains the puzzling temperature dependence of polarized neutron scattering data reported in Ref. Petit2016 .

The Nd2Zr2O7 single crystal ($\sim 2.5$ gram) was grown by using an optical floating zone furnace in the Core Laboratory for Quantum Materials (QM Core Lab) in Helmholtz-Zentrum Berlin (HZB) and characterized using X-ray powder diffraction and Laue diffraction Xu2019 . Inelastic neutron scattering experiments were conducted on the direct-geometry time-of-flight (tof) spectrometer CNCS at SNS in Oak Ridge National Lab and on the indirect-geometry tof spectrometer Osiris at ISIS Neutron Source in Rutherford Appleton Lab Ehlers2016 ; Telling2016 . For the CNCS measurement, the sample was mounted on a 3He insert which cooled the sample down to 240 mK. Neutrons of incident wavelength 4.98 Å (3.315 meV) were used in the high-flux mode of the instrument (energy resolution $\sim 0.1$ meV). Data were collected at 240 mK, 450 mK and 20 K with a 360-degree sample rotation at a step of one degree. The large reciprocal space coverage provides an overview of the spin dynamics. To resolve the spin dynamics better, we performed the experiment on Osiris with a higher energy resolution 25 $\mu$eV using the PG(200) analyser which analyses scattered neutrons of energy 1.84 meV. The crystal was mounted onto a dilution refrigerator and data were collected at 30 mK, 450 mK and 20 K. The data at 20 K were used as the background for the two experiments. The software packages Dave dave , Mantid mantid and Horace horace were used for data processing. Some data were averaged based on the symmetry of the scattering plane in order to improve the statistics. The classical Monte Carlo simulations and semiclassical molecular dynamics calculations were done using the MATLAB package SPINW spinw , implementing our own codes (see supplementary materials SM ).

Figure 1(a) shows the constant energy slice at 0.025 meV of the CNCS data at 450 mK (the 20 K data is subtracted as background). It is a well-formed strong pinch point pattern but at a much lower energy than in the ordered state ($\sim 0.075$ meV) (see supplementary materials SM ). Fig. 1(b) presents the background-subtracted data measured on Osiris at 30 mK and 450 mK with a much higher energy resolution 25 $\mu$eV. Integrating over [-0.02, 0.02] meV, we see clearly a strong scattering arm along (111) of the pinch point pattern at 450 mK, consistent with the CNCS data, whereas the data at 30 mK does not show any signals except for the (220) magnetic Bragg peak. Conversely, in the data with integration over [0.06, 0.08] meV, the 450 mK data does not show a clear pattern while the data at 30 mK shows the expected pinch-point spinwave modes from the AIAO order Petit2016 ; Xu2019 . In Fig. 1(c), the $E-Q$ slices along the (22L) direction show that the gapped magnon excitations vanish above $T_{\text{N}}$ and a strong scattering appears around the elastic line. In addition, a weak continuum at finite energy around (220) at 450 mK is also a new feature which has a similar intensity with the dispersing magnon modes at 30 mK.

In Fig. 2, we show the one-dimentional cuts through high-energy-resolution data at $Q=(0.9,0.9,0.9)$ and $(220)$ measured at 30 mK, 450 mK and 20 K. Comparing with the 20 K background data, we see pronounced gapped inelastic signals for 30 mK [besides the magnetic (220) Bragg peak] and strong elastic signals for 450 mK [Fig. 2(a) and (c)]. After background subtraction as shown in Fig. 2(b) and (d), we see that with the disappearance of the magnon excitations, most of the spectral weight changes to be elastic at 450 mK (there is a minor over-subtraction at negative energy transfers).

Therefore, we can conclude that above $T_{\text{N}}$, there appears a non-trivial paramagnetic phase with significant spin-ice correlations. It is quite surprising that the strong 2I2O correlations appears just above the AIAO ordering temperature though the broad feature around the AIAO Bragg peaks at [220] and [113] may indicate the existence of short-range AIAO order.

To investigate the nature and origin of the high-temperature phase, we did classical Monte Carlo simulations based on the anisotropic spin Hamiltonian determined in Ref. Xu2019 ,

[TABLE]

where $\tau^{\alpha}_{i}$ ($\alpha=x,y,z$) is the $\alpha$ component of the pseudospin-1/2 at site $i$ defined in the rotated local frames and ${\tilde{J}}_{\tilde{\alpha}}$ is the corresponding nearest-neighbor exchange constant Huang2014 ; Benton2016 . The calculations were done with a supercell of $6\times 6\times 6$ cubic unit cells, which has 3456 spins in total (for details see the supplementary materials SM ). The simulated specific heat indicates that the system enters the AIAO phase at $T_{\text{N}}^{\prime}\approx 0.18$ K. Above $T_{\text{N}}^{\prime}$ at 0.25 K, the calculated neutron scattering structure factor [Fig. 3(a)] shows a pinch point pattern with broad signals around the AIAO Bragg peak wavevectors [220] and [113], which is quite consistent with the experiment.

To identify the spin correlations responsible for this novel scattering pattern, we calculated the thermal average of the amplitudes of the pseudospin components, $\langle|\tau^{\tilde{\alpha}}|\rangle$ ($\alpha=x,y,z$) and the probability distribution functions (pdfs) of $\tau^{\tilde{\alpha}}$ and the average over tetrahedra $1/4\sum_{\boxtimes}\tau^{\tilde{\alpha}}$. Fig. 3(b) shows the temperature dependence of $\langle|\tau^{\tilde{\alpha}}|\rangle$. At base temperature, only $\langle|\tau^{\tilde{z}}|\rangle$ are significant consistent with the AIAO order. With increasing the temperature, $\langle|\tau^{\tilde{z}}|\rangle$ decreases while $\langle|\tau^{\tilde{x}}|\rangle$ and $\langle|\tau^{\tilde{y}}|\rangle$ increase due to enhanced thermal fluctuations and all approach $\tau/2=0.25$ ($\tau=1/2$ is the amplitude of pseudospin) expected for completely random spin configuration. Remarkably, above $T_{\text{N}}^{\prime}$, $\langle|\tau^{\tilde{x}}|\rangle$ gets higher than the thermal average $\tau/2$ and $\langle|\tau^{\tilde{z}}|\rangle$ and then decays slowly with increasing temperature exhibiting a maximum around $T_{\text{N}}^{\prime}$. By contrast, $\langle|\tau^{\tilde{z}}|\rangle$ decreases quickly and even gets slightly lower than $\tau/2$, revealing the breakdown of the AIAO correlations.

The correlations of $\tau^{\tilde{x}}$ is the origin of the pinch point pattern. We characterize the correlations of $\tau^{\tilde{x}}$ with the probability distribution functions mentioned above. As shown in Fig. 3(c), at temperatures below $T_{\text{N}}^{\prime}$, pdf$(\tau^{\tilde{x}})$ is a Gaussian function centered at zero and its width increases with raising temperature. Above $T_{\text{N}}^{\prime}$, it surprisingly turns to be two peaks at $\pm\tau$ decaying with further increasing temperature, which means that $\tau^{\tilde{x}}$ points into the tetrahedron for half of the spins and out of the tetrahedron for the other half. The pdf of the average on tetrahedra $1/4\sum_{\boxtimes}\tau^{\tilde{x}}$ shows how the half/half in/out $\tau^{\tilde{x}}$ is distributed on the tetrahedra which is always a Gaussian function centered at zero (shown in Ref. SM ). The above two statistical quantities indicate that ice-rule correlations appear for the $\tau^{\tilde{x}}$ component on the tetrahedra which is consistent with the FM $\tilde{J}^{\tilde{x}}$. On the other hand, pdf$(\tau^{\tilde{z}})$ and pdf$(1/4\sum_{\boxtimes}\tau^{\tilde{z}})$ change from peaks at either $\tau$ or $-\tau$ (depending on the AIAO domain type) to be a broad Gaussian peak centered at zero indicating the loss of the AIAO correlations of $\tau^{\tilde{z}}$ above $T_{\text{N}}^{\prime}$ SM .

Both experiments and classical MC simulations reveal strong spin ice correlations in the system above $T_{\text{N}}$. Assuming the existence of a Coulombic phase with respect to $\tau^{\tilde{x}}$ above $T_{\text{N}}$, we have calculated the correlations based on the bosonic many-body theory of quantum spin ice SM ; Hao2014 . As shown in Fig. 4, the calculated dynamical structure factor at 450 mK (integrated over [0, 0.05] meV) exhibits a pinch point pattern with additional broad scattering around [220] and [113] in good agreement with experimental data. According to the theory, besides the pinch point pattern at zero energy due to the scattering of the Coulombic phase, monopole creation and hopping cause broad scattering around [220] and [113]. Neutrons can be scattered by monopoles via two different processes: (i) the incoming neutron flips a spin belonging to an ice-rule tetrahedron creating a pair of monopoles, which gives a continuum at finite energy above a small gap; (ii) at finite temperature where there are a finite density of monopoles already in the system, the incoming neutron can flip a spin belonging to a monopole tetrahedron causing this monopole to hop which gives a continuum of scattering around zero energy. These scattering features are shown in Fig. 4(b). This agrees with our data which exhibits strong broad signals around [220] and [113] at zero energy and a continuum at finite energy. However, the gapped feature is not clear in the data which could be attributed to possible fast decay of the coherent monopole excitations due to, for example, strong thermal fluctuations. The signal also may be contaminated by the scattering from possible short-range AIAO correlations. More experiments are needed to clarify this.

The existences of strong gapless spin-ice correlations above $T_{\text{N}}$ and the gapped flat dynamical spin-ice modes in the ordered state Petit2016 ; Xu2019 make the ordering a candidate for a Higgs transition where the emergent gauge field of the Coulomb phase is gapped by the condensation of emergent gauge charges (monopoles) Chang2012 ; Powell2011 . Our semiclassical molecular dynamics calculations (Fig. 3(d) and Ref. SM ) show that the gap closes with raising temperature which supports this picture. A Higgs transition was also proposed for pyrochlore Yb2Ti2O7 by demonstrating the first order nature of the ordering transition and the sudden suppression of the intensity of the pinch point pattern below the transition Chang2012 . The gapped and gapless pinch point patterns below and above $T_{\text{N}}$ is another important evidence for a Higgs transition.

Our results provide an explanation for the seemingly contradicting temperature dependences of the intensities of the AIAO Bragg peak and the pinch point pattern in the energy-integrated polarized neutron scattering data reported for Nd2Zr2O7 Petit2016 and very recently for Nd2ScNbO7 Mauws2019 . It was shown that with increasing temperature, the magnetic Bragg peak weakens and disappears at $T_{\text{N}}$ while the pinch point intensity maximizes around $T_{\text{N}}$ and persists to much higher temperatures ($\sim 1$ K) for both Nd2Zr2O7 and Nd2ScNbO7. This cannot be rationalized if the pinch point scattering is only a feature of the magnons in the ordered state. We argue that at $T_{\text{N}}$ the gapped pinch point pattern is replaced by a pinch point pattern at zero energy due to the Coulombic phase built on ${\tilde{\tau}}^{\tilde{x}}$ which has the strongest ice correlations around $T_{\text{N}}$ and gets weaker slowly with increasing temperature, similar to the temperature evolution of $\langle|\tau^{\tilde{x}}|\rangle$ in the MC simulations [Fig. 3(b)]. The temperature range where the pinch point pattern presents is also comparable with $\tilde{J}^{\tilde{x}}\sim 1$ K. In addition, the Coulombic phase with strongly correlated spins could induce slow spin dynamics which supports the observed anomalously slow paramagnetic spin dynamics in the muon spin relaxation experiments Xu2016 .

What is more, it was reported recently that a spin ice model with frustrated transverse terms exhibits competing phases and nematicity Mathieu2017 ; Benton2018 . Our results provide a concrete experimental and theoretical example of an Ising antiferromagnet with frustrated transverse terms which also could show interesting physics. Our further MC simulations in Ref. SM show that the AIAO ordering temperature is suppressed largely with increasing $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|$ which should be constant according to the mean field theory and on the other hand, the ordered phase invades the spin ice phase at finite temperature for $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|>3$ similar to the phase diagram in Ref. Mathieu2017 which is surprising because the spin ice phase should be more stable due to the higher entropy. This suggests that further theoretical study is needed.

In summary, we used high-energy-resolution inelastic neutron scattering, classical MC simulations, semi-classical molecular dynamics simulations and a bosonic theory of quantum spin ice to show that Nd2Zr2O7 has a non-trivial paramagnetic state with 2I2O spin-ice correlations, which is possibly a magnetic Coulombic phase, despite a long-range AIAO ordered ground state. We attributed it to the dominant frustrated transverse $\tilde{J}^{\tilde{x}}$ term in the spin Hamiltonian and related the ordering transition to the Higgs mechanism. Our results indicate that the paramagnetic phase of an ordered system may host unconventional spin correlations different from the ground-state order in nature due to competition and frustration among different terms of anisotropic exchange interactions. This expands the field for searching for quantum spin ice to the ordered systems with frustrated terms in the spin Hamiltonian and makes it interesting to examine several similar QSI candidates with AIAO order, such as Nd2(Hf/Sn/Pb)2O7, Nd2ScNbO7 and Sm2(Ti/Sn)2O7 Anand2015 ; Bertin2015 ; Hallas2015 ; Mauws2018 ; Mauws2019 ; Viviane2019 .

I INS data measured on CNCS and Osiris

Figure S1 shows constant energy slice of the (HHL) reciprocal plane (integrated over [0, 0.03] meV) without symmetrization where a pinch point pattern is clearly shown.

Figure S2 shows the $E-Q$ slices along the (HHH) direction (an arm of the pinch point pattern) of the Osiris data measured at 450 mK and 30 mK which show that above $T_{\text{N}}$, the gapped pinch point pattern disappears and a pinch point pattern at zero energy appears.

Fig. S3 shows the constant-energy slices in the (HHL) plane of the Osiris data measured at 30 mK and 450 mK. Fig. S4 shows the constant-energy slices in the (HHL) plane of the CNCS data measured at 240 mK and 450 mK. Both show that the spinwave excitations disappeared or are extremely weak above $T_{\text{N}}$ and a pinch point pattern at zero energy shows up. Fig. S5 shows the disappearance of the magnetic Bragg intensity at (220) and (-2-20) at 450 mK in the CNCS data.

II Quantum fluctuations in linear spin wave theory

Figure S6 shows the calculated reduced ordered spin for $0<\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|<3$ due to zero-point quantum fluctuations based linear spin wave theory using SPINW MATLAB package [43]. The gap to flat modes are calculated from the equation [33]

[TABLE]

They are calculated by changing $\tilde{J}^{\tilde{x}}$ with fixing $\tilde{J}^{\tilde{y}}$ and $\tilde{J}^{\tilde{z}}$ of the spin Hamiltonian for Nd2Zr2O7 [35]. With increasing $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|$, the gap decreases and the quantum fluctuations are enhanced drastically for $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|>2$ due to the proximity to a U(1) spin liquid. This can be related to the persistent spin dynamics in the muon spin relaxation experiments of Nd2Zr2O7 [31].

III Monte Carlo simulations

The classical single-spin-flip Metropolis Monte Carlo simulations are performed using SPINW Matlab package by implementing our own codes [43]. We used a $6\times 6\times 6$ supercell with periodic boundary conditions. At each temperature, $5\times 10^{4}$ Monte Carlo steps per spin were used for thermal equilibrium and $5\times 10^{5}-1\times 10^{6}$ Monte Carlo steps per spin for data collection. To avoid freezing of the dynamics at low temperatures, the spins are moved in a cone around its present direction with the cone angle as a free parameter chosen to keep the acceptance rate above $30$%.

Figure S7 shows the simulated temperature dependence of specific heat $C_{\text{v}}$ using the spin Hamiltonian of Nd2Zr2O7 with setting $\tilde{J}_{\tilde{x}}=\{1,\sim 2,3,3.2,3.3,3.5\}|\tilde{J}_{\tilde{z}}|$ and Fig. S8 shows $\langle|\tau^{\tilde{\alpha}}|\rangle$ ($\alpha=x,y,z$) as a function of temperature for different $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|$. They show that spin-ice correlations become stronger above $T_{\text{N}}$ with increasing $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|$ and the AIAO ordering temperature is suppressed largely from $\sim 0.20$ K for $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|=1$ to $\sim 0.13$ K for $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|=3$ which should be constant within the mean-field theory. It is very interesting that the system still shows a long-range AIAO order at the critical point $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|=3$ where spin ice state is expected which has the same energy with the AIAO order but a lower free energy because of its higher entropy. For $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|=3.2$ and $3.3$, $C_{\text{v}}$($T$) show double peaks which are related to tendency of AIAO ordering and a re-entrant transition out of AIAO order to the spin ice regime, as shown in Fig. S8. The interesting competition between the AIAO order and spin ice will be studies further.

Figure. S9 shows the probability distribution functions of $\tau^{\tilde{\alpha}}$ and the average over tetrahedra $1/4\sum_{\boxtimes}\tau^{\tilde{\alpha}}$ for $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|\sim 2$ (for Nd2Zr2O7) and $3.2$. For Nd2Zr2O7, strong spin-ice correlations of $\tau^{\tilde{x}}$ are built up with decreasing temperature and then suddenly suppressed below $T_{\text{N}}$. For $\tilde{J}^{\tilde{x}}/|\tilde{J}^{\tilde{z}}|=3.2$, the system first shows spin ice correlations at high temperature, then tendency of AIAO ordering and finally spin ice phase at low temperatures.

IV Spin dynamics at finite temperatures

The spin dynamics at finite temperature were calculated using a semiclassical molecular dynamics. An ensemble of configurations is first obtained from Monte Carlo simulations with a supercell of $10\times 10\times 10$ cubic unit cells (16000 spins in total) with periodic boundary conditions. The ensemble is then evolved in time according to the equation of motion,

[TABLE]

which is integrated numerically using eight-order Runge-Kutta method with an adaptive step-size control. In the simulation, $5\times 10^{5}$ Monte Carlo steps were used for thermalization at a fixed temperature and measurements were taken every 1000 steps in $5\times 10^{5}\sim 2.5\times 10^{6}$ Monte Carlo steps.

Fig. S10 shows the simulated $E-Q$ slices along (HH0) direction at different temperatures which show the gap decreases with increasing temperature and most of the spectral weight changes to be elastic.

V Theoretical calculation of structure factor in finite temperature quantum spin ice

Here we give details of theoretical calculations of the structure factor $S(\mathbf{q},\omega)$ for $T>T_{N}$. These calculations are made by modelling the system above $T_{N}$ as a spin ice with bosonic, propagating, quantum monopole (spinon) excitations. We first give some necessary definitions before discussing the low energy spin ice correlations [Sec. V.1] and the non-interacting description of propagating monopoles [Sec. V.2].

The structure factor of interest is:

[TABLE]

where $i,j=0,1,2,3$ index the four sublattices of the pyrochlore lattice and $\mu,\nu$ are spatial directions in the crystal coordinate frame. ${\bf m}_{i}(\mathbf{q})$ is the Fourier transform of the magnetic moment configuration on sublattice $i$:

[TABLE]

where $N_{\sf uc}$ is the number of unit cells in the system, ${\bf R}$ are the centers of the unit cells, $\frac{1}{2}\mathbf{u}_{i}$ is the relative position of site $i$ within a unit cell and ${\bf m}_{i}(\mathbf{R}+\mathbf{u}_{i}/2)$ is the magnetic moment on that site. We take a coordinate frame where the vectors ${\bf u}_{i}$ are

[TABLE]

The relationship between the magnetic moments ${\bf m}_{i}$ and the pseudospin operators $\tau^{\tilde{\alpha}}_{i}$ which appear in the theoretical description of the magnetism [27, 33] is

[TABLE]

where $\mu_{B}$ is the Bohr magneton. The parameters $g_{zz}$ and $\vartheta$ were estimated for Nd2Zr2O7 in Ref. [35]:

[TABLE]

and the local axes $\hat{\bf z}_{i}$ are

[TABLE]

Inserting Eq. (S6) into Eq. (S3), and defining lattice Fourier transforms of $\tau^{\tilde{\alpha}}_{i}$ analogously to Eq. (S4) gives:

[TABLE]

The symmetries of the XYZ Hamiltonian for dipolar-octupolar pyrochlores [27] guarantee that the cross terms $\langle\tau^{\tilde{z}}\tau^{\tilde{x}}\rangle$ vanish in the high-temperature phase. Using this, we can express $S(\mathbf{q},\omega)$ as a sum of two distinct contributions:

[TABLE]

where

[TABLE]

Our model of the phase just above $T_{N}$ has spin ice correlations with respect to the $\tilde{x}$ axis of pseudospin space. In this picture the correlator $S^{\tilde{x}\tilde{x}}(\mathbf{q},\omega)$ measures processes within the ice manifold where there are two spins with $\tau^{\tilde{x}}=\frac{1}{2}$ and two with $\tau^{\tilde{x}}=-\frac{1}{2}$ on each pyrochlore tetrahedron. The correlator $S^{\tilde{z}\tilde{z}}(\mathbf{q},\omega)$ then measures excitations out of that manifold: i.e. monopoles.

The plot in Fig. 4 of the main text is made using Eq. (S10), calculating $S^{\tilde{x}\tilde{x}}(\mathbf{q},\omega)$ as described in Section V.1 and $S^{\tilde{z}\tilde{z}}(\mathbf{q},\omega)$ as described in Section V.2.

V.1 Low energy spin-ice correlations

The calculation of $S^{\tilde{x}\tilde{x}}(\mathbf{q},\omega)$ assumes that at $T_{N}<T<\tilde{J}^{\tilde{x}}$, the system fluctuates among states obeying an ice rule for $\tau^{\tilde{x}}$.

In a quantum spin ice there will be some energy scale, $g$, for quantum tunnelling between these states [7]. At temperatures, $T\lesssim g$, the correlations will be significantly modified from their classical form with the signature pinch point singularities being suppressed, and a photon dispersion being visible in the inelastic spectrum [12]. At higher temperatures $T>g$, the correlations of $\tau^{\tilde{x}}$ will essentially take on their classical form [12, 14, 15] and will be centred on $\omega=0$. In our case, since to study the Coulomb phase we are restrictred to $T>T_{N}>g$, we assume the system to be in the latter regime where the correlations of $\tau^{\tilde{x}}$ are essentially classical.

The structure factor of a classical spin ice is well known and has the form

[TABLE]

where $P_{ij}(\mathbf{q})$ is the projection matrix given by Henley [6]:

[TABLE]

$E_{ia}(\mathbf{q})$ is a $4\times 2$ matrix

[TABLE]

and $M_{ab}(\mathbf{q})$ is a $2\times 2$ matrix

[TABLE]

To account for finite energy resolution in the experiment, we replace the delta function in Eq. (S13) with a Gaussian:

[TABLE]

with full width at half maximum (FWHM) $f=0.1$meV for Fig. 4(a) and $f=0.02$meV for Fig. 4(b).

V.2 Non-interacting model for propagating monopoles

The calculation of $S^{\tilde{z}\tilde{z}}$ is based on the intuition that if the system obeys an ice rule with respect to $\tau^{\tilde{x}}$, the action of $\tau^{\tilde{z}}$ will create pairs of violations of the ice rule: namely, monopoles. Our treatment of the monopoles follows the non-interacting theory introduced in Ref. [45].

Introducing ladder operators with respect to the $\tilde{x}$ axis of spin space:

[TABLE]

we can write the Hamiltonian as

[TABLE]

In Nd2Zr2O7 we have [35] $\tilde{J}^{\tilde{x}}\approx 0.091$meV, $\tilde{J}^{\tilde{\pm}}=-\frac{1}{4}\left(\tilde{J}^{\tilde{y}}+\tilde{J}^{\tilde{z}}\right)\approx 0.008$meV, $\tilde{J}^{\tilde{\pm}\tilde{\pm}}=\frac{1}{4}\left(\tilde{J}^{\tilde{y}}-\tilde{J}^{\tilde{z}}\right)\approx 0.015$meV.

The $\tilde{J}^{\tilde{x}}$ term in Eq. (S19) exerts an energy cost $\tilde{J}^{\tilde{x}}/2$ per monopole. The $\tilde{J}^{\tilde{\pm}}$ term generates monopole hopping. $\tilde{J}^{\tilde{\pm}\tilde{\pm}}$ is an interaction term for the monopoles.

The monopoles are defined on the pyrochlore tetrahedra. The centers of these tetrahedra form a diamond lattice, which may be divided into two interpenetrating FCC sublattices A and B. The monopole charge at a tetrahedron ${\bf R}$ is given by

[TABLE]

where $\epsilon_{\bf R}=+1$ for ‘A’ sublattice tetrahedra and $\epsilon_{\bf R}=-1$ for ‘A’ sublattice tetrahedra.

Following Ref. [45], we can relate the spin ladder operators to raising and lowering operators of the charge $\psi^{\dagger}_{\bf R},\psi_{\bf R}$

[TABLE]

$A_{{\bf R}{\bf R}^{\prime}}$ is the gauge field to which the charges are coupled. In the spirit of the mean field approximation of Ref. [9] we decouple the gauge field from the monopoles and replace

[TABLE]

to write

[TABLE]

For our present purposes we will determine $\kappa$ later from the sum rule on the structure factor. The charge raising and lowering operators obey a constraint

[TABLE]

As in Ref. [45] we write $\psi_{\bf R}$, $Q_{\bf R}$ in terms of boson operators $d_{\bf R},b_{\bf R}$ which respectively carry postive and negative charge.

[TABLE]

The $d$ and $b$ bosons obey a constraint, that they cannot occupy the same site at the same time

[TABLE]

Assuming low-density of the bosons we can also expand the denominator of Eq. (S27) to write

[TABLE]

We can then write $\mathcal{H}_{\tau}$ in terms of spinon operators, keeping only terms up to bilinear order:

[TABLE]

The first term in Eq. (S32) is simply a chemical potential term for the spinons while the second incorporates both hopping and charge conserving pair creation/annihilation processes. The coefficients of these terms may be renormalized from their bare values by both fluctuations of the gauge field and by interactions. To take this into account we will take

[TABLE]

and consider $\mu$ and $\eta$ as phenomenological parameters.

We are then left with two copies (one on the $A$ sublattice, one on the $B$ sublattice) of a non-interacting boson Hamiltonian on the FCC lattice.

To diagonalize Eq. (S32) we Fourier transform to reciprocal space, and then do a Bogoliubov transformation

[TABLE]

where

[TABLE]

This brings the Hamiltonian into the form

[TABLE]

with the spinon dispersion $\omega_{\mathbf{q}}$ given by Eq. (S39).

We then seek to turn the non-interacting treatment of the spinons into a calculation of $S^{\tilde{z}\tilde{z}}(\mathbf{q},\omega)$. We first write

[TABLE]

and note that the spinon hopping Hamiltonian conserves the total spinon charge on each FCC sublattice independently. The correlators $\langle{\tau}^{+}_{i}(-\mathbf{q},t){\tau}^{+}_{j}(\mathbf{q},0)\rangle$ and $\langle{\tau}^{-}_{i}(-\mathbf{q},t){\tau}^{-}_{j}(\mathbf{q},0)\rangle$ appearing in Eq. (S41) do not conserve these charges independently and must therefore vanish within this approximation. We are therefore left with

[TABLE]

The spin operators relate to the charge raising and lowering operators through Eqs. (S24)- (S25). Using these relationships, leads to the correlation function of interest:

[TABLE]

The expectation values in Eq. (S43) can be calculated from the non-interacting theory, as products of bilinear expectations values:

[TABLE]

$n_{B}(\omega)$ is the Bose-Einstein distribution at temperature $T$.

Using these we obtain the final result for $S^{\tilde{z}\tilde{z}}(\mathbf{q},\omega)$:

[TABLE]

where

[TABLE]

The first term in Eq. (S45) comes from the creation of monopoles pairs by an incoming neutron, the final one from the annihilation of monopole pairs and the middle two from processes where a neutron flips a spin to hop an already-present monopole.

For the purposes of the plot in Fig. 4 of the manuscript Eq. (S45) is evaluated with Monte Carlo integration over the Brillouin zone for each value of ${\bf q}$. In the calculations we set $\mu=0.05$meV, $\eta=0.002$meV, $T=450$mK. $\kappa$ is fixed by requiring agreement with the sum rule:

[TABLE]

Finite energy resolution is accounted for by replacing the delta functions in Eq. (S45) with Gaussians of FWHM $f=0.1$meV for Fig. 4(a) and $f=0.02$meV for Fig. 4(b) [Eq. (S17)].