Anisotropic exchange Hamiltonian, magnetic phase diagram and domain inversion of Nd$_2$Zr$_2$O$_7$

TL;DR

This study refines the anisotropic exchange Hamiltonian for Nd$_2$Zr$_2$O$_7$, explores its magnetic phase diagram, and investigates domain inversion mechanisms, providing insights into its complex magnetic behavior and dynamical kagome ice phenomena.

Contribution

The paper introduces a refined anisotropic exchange Hamiltonian for Nd$_2$Zr$_2$O$_7$ and combines experimental and theoretical methods to analyze its magnetic phases and domain inversion.

Findings

Qualitative agreement between calculated and experimental phase diagrams.

Identification of metastable states and domain inversion mechanisms.

Explanation of dynamical kagome ice in [111] fields.

Abstract

We present thermodynamic and neutron scattering measurements on the quantum spin ice candidate Nd$_2$Zr$_2$O$_7$. The parameterization of the anisotropic exchange Hamiltonian is refined based on high-energy-resolution inelastic neutron scattering data together with thermodynamic data using linear spin wave theory and numerical linked cluster expansion. Magnetic phase diagrams are calculated using classical Monte Carlo simulations with fields along \mbox{[100]}, \mbox{[110]} and \mbox{[111]} crystallographic directions which agree qualitatively with the experiment. Large hysteresis and irreversibility for \mbox{[111]} is reproduced and the microscopic mechanism is revealed by mean field calculations to be the existence of metastable states and domain inversion. Our results shed light on the explanations of the recently observed dynamical kagome ice in Nd$_2$Zr$_2$O$_7$ in \mbox{[111]} fields.

Tables

Table 1. Table 1: Comparison of the Hamiltonian parameterisations and the calculated T CW subscript 𝑇 CW T_{\text{CW}} of Nd 2 Zr 2 O 7 .

${\tilde{J}}_{\tilde{x}}$ (meV)	${\tilde{J}}_{\tilde{y}}$ (meV)	${\tilde{J}}_{\tilde{z}}$ (meV)	$\vartheta$ (rad.)	$g_{zz}$	$T_{\text{CW}}$(K)	Refs.
-0.047	0	0.103	0	4.5	-	Petit2016
0.103	0	-0.047	0.83	-	-	Benton2016
0.086(4)	0.006(20)	-0.043(4)	1.26(17)	4.55	0.43(8)	Lhotel2018
0.091(9)	0.014(6)	-0.046(2)	0.98(3)	5.0(1)	0.28(4)	This work

Equations

$${\mathcal{H}}_{\text{ex}}=\sum _{⟨ij⟩}{J}_x{\tau }_i^x{\tau }_j^x+J_y{\tau }_i^y{\tau }_j^y+J_z{\tau }_i^z{\tau }_j^z+J_{xz}({\tau }_i^x{\tau }_j^z+{\tau }_i^z{\tau }_j^x),$$

$${\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{]}.$$

$$J_x={\tilde{J}}_{\tilde{x}}{\cos }^2\vartheta +{\tilde{J}}_{\tilde{z}}{\sin }^2\vartheta ,$$

$$J_z={\tilde{J}}_{\tilde{z}}{\cos }^2\vartheta +{\tilde{J}}_{\tilde{x}}{\sin }^2\vartheta ,$$

$$\tan 2\vartheta =\frac{2J_{xz}}{J_x-J_z},$$

$$J_{xz}=({\tilde{J}}_{\tilde{x}}-{\tilde{J}}_{\tilde{z}})\sin \vartheta \cos \vartheta ,$$

$${\tilde{\tau }}_i^{\tilde{x}}$$

$${\tilde{\tau }}_i^{\tilde{y}}$$

$${\tilde{\tau }}_i^{\tilde{z}}$$

$$m_z$$

$$={\mu }_{\text{B}}{g}_{zz}({\tilde{\tau }}^{\tilde{z}}\cos \vartheta +{\tilde{\tau }}^{\tilde{x}}\sin \vartheta ),$$

$${\Delta }_1=\sqrt{(3|{\tilde{J}}_{\tilde{z}}|-{\tilde{J}}_{\tilde{x}})(3|{\tilde{J}}_{\tilde{z}}|-{\tilde{J}}_{\tilde{y}})},$$

$${\Delta }_2$$

$${\Delta }_3$$

$${\Delta }_1=0.075(4),{\Delta }_2=0.186(4),{\Delta }_3=0.272(5)\text{meV}.$$

$${\tilde{J}}_{\tilde{x}}=0.091(9),{\tilde{J}}_{\tilde{y}}=0.014(6),{\tilde{J}}_{\tilde{z}}=-0.046(2)\text{meV}.$$

$$C_{\text{p}}\approx \frac{\alpha }{T^2}$$

$$\alpha =\frac{3R}{16k_{\text{B}}^2}({\tilde{J}}_{\tilde{x}}^2+{\tilde{J}}_{\tilde{y}}^2+{\tilde{J}}_{\tilde{z}}^2),$$

$$m_{\text{ord}}=\left(\frac{1}{2}-\delta {\tilde{\tau }}^{\tilde{z}}\right)\cos \vartheta g_{zz}{\mu }_{\text{B}},$$

$$g_{zz}=5.0(1),\vartheta =0.98(3)\text{radian}.$$

$$T_{\text{CW}}=\frac{{\tilde{J}}_{\tilde{z}}{\cos }^2\vartheta +{\tilde{J}}_{\tilde{x}}{\sin }^2\vartheta }{2k_B}.$$

$$J_x$$

$$J_z$$

$$m_{\text{AIAO}}=\frac{1}{N}\sum {\tilde{\tau }}_i^{\tilde{z}},$$

$$\varepsilon=\big{[}3g_{zz}H\tau\cos(\vartheta+\phi_{1})\cos\alpha+g_{zz}H\tau\cos(\vartheta+\phi_{2})\cos\beta\big{]}\\ +6\big{[}\tilde{J}_{\tilde{x}}\tau^{2}\sin\phi_{1}\sin\phi_{2}+\tilde{J}_{\tilde{z}}\tau^{2}\cos\phi_{1}\cos\phi_{2}\\ +\tilde{J}_{\tilde{x}}\tau^{2}\sin^{2}\phi_{1}+\tilde{J}_{\tilde{z}}\tau^{2}\cos^{2}\phi_{1}\big{]}$$

$$\frac{1}{N}⟨\mathcal{O}⟩=\sum _{c}M(c)W(c).$$

$$W(c)=⟨\mathcal{O}{⟩}_c-\sum _{s\subset c}W(s)$$

$$⟨\mathcal{O}{⟩}_{Euler3}=\frac{1}{2}⟨\mathcal{O}{⟩}_{NLC2}+\frac{1}{2}⟨\mathcal{O}{⟩}_{NLC3}$$

$$⟨\mathcal{O}{⟩}_{Euler4}=\frac{1}{4}⟨\mathcal{O}{⟩}_{NLC2}+\frac{1}{2}⟨\mathcal{O}{⟩}_{NLC3}+\frac{1}{4}⟨\mathcal{O}{⟩}_{NLC4}$$

Full text

I Introduction

Recently magnetic pyrochlore oxides have attracted a lot of attention because of geometrical frustration on the network of corner-sharing tetrahedra in the crystal structure Gardner2010rev . A rich spectrum of magnetic behaviors have been observed in rare-earth pyrochlores. For example, (Ho/Dy)2Ti2O7 show a classical spin ice state and monopole-like excitations originating from the Ising-anisotropy and effective ferromagnetic spin interactions Harris1997 ; Fennell2009 ; Morris2009 . Classical spin ice can be melted by quantum fluctuations coming from transverse terms in the spin Hamiltonian. This leads to quantum spin ice which is a type of long-sought quantum spin liquid with massive many-body entanglement and fractionalized excitations Hermele2004 ; Shannon2012 ; Savary2012 ; Gingras2014 . The quantum effect in the classical spin ice is generally ignored but could be significant in pyrochlores with Ising-anisotropic light rare earth elements, e.g. Ce, Pr, Sm, Nd Onoda2010 ; Onoda2011 ; Rau2015 .

Among the candidates for quantum spin ice, Nd-containing pyrochlores are quite special due to the peculiar symmetry of the dipolar-octupolar crystal field ground state doublet of the Nd3+ ion Huang2014 ; Li2016 ; Li2017 . It was shown theoretically that the pseudospin-1/2 Hamiltonian has the form of the $XYZ$ model which supports different types of symmetry-enriched quantum spin ice phases, e.g., dipolar and octupolar quantum spin ice states Lee2012 ; Li2016 ; Li2017 ; Benton2018 .

On the experimental side, the study of Nd2Zr2O7 is rather intensive and fruitful compared to the other Nd-containing pyrochlores Blote1969 ; Lutique2003 ; Hatnean2015 ; Lhotel2015 ; Xu2015 ; Xu2016 ; Petit2016 ; Benton2016 ; Opherden2017 ; Xu2018 ; Lhotel2018 ; Anand2015 ; Anand2017 ; Opherden2018 ; Bertin2015 . The single-ion crystal field ground state was confirmed to be an Ising-anisotropic dipolar-octuplar doublet and the collective ground state was found to be an “all-in-all-out” (AIAO) antiferromagnetic order (Néel temperature $T_{\text{N}}\approx 0.4$ K) with persistent spin dynamics Lhotel2015 ; Xu2015 ; Xu2016 . Single crystal inelastic neutron scattering (INS) in zero field indicated dynamical spin ice correlations and quantum moment fragmentation, and a parameterized spin Hamiltonian was proposed Petit2016 ; Benton2016 . Very recently, quantum spin-1/2 chains and dynamical kagome ice were detected in Nd2Zr2O7 with external fields along [110] and [111] directions and the observed neutron scattering spectra permit a complete determination of the exchange parameters Xu2018 ; Lhotel2018 . In addition, the susceptibility and magnetization measurements reveal a field-induced transition (critical field $H_{\text{c}}\sim 0.1$ T) with a large hysteresis in [111] fields Opherden2017 ; Lhotel2015 .

Although the measurements and analyses on Nd2Zr2O7 are quite comprehensive, a qualitative understanding of the observations based on a spin Hamiltonian is still lacking Petit2016 ; Benton2016 ; Lhotel2018 . The recently determined exchange parameters still cannot describe the INS data in several aspects (e.g. the transition field and spin excitation energies in fields). In this paper, we first refined the spin Hamiltonian by a comprehensive analysis of the high-energy-resolution INS data and thermodynamic properties of Nd2Zr2O7 using linear spin wave theory, high temperature series expansion and numerical linked cluster expansion (NLCE). Second, we calculated the phase diagrams for external fields along [100], [110] and [111] directions using classical Monte Carlo simulations finding qualitative agreement with experiment. Third, we found that the large hysteresis for the [111] field is caused by an irreversibility of the magnetisation or domain inversion due to the existence of a metastable state based on mean field calculations which is critical for explaining the observed dynamical kagome modes in [111] fields Lhotel2018 .

II Methods

II.1 Single crystal growth and structural and magnetic characterizations

The Nd2Zr2O7 single crystals were grown by using the optical floating zone furnace in the Core Laboratory for Quantum Materials (QM Core Lab) at Helmholtz-Zentrum Berlin (HZB) and characterized by X-ray powder diffraction and Laue diffraction. The susceptibility, magnetization and specific heat measurements above 2 K were performed on MPMS (SQUID) and PPMS (VSM) (Quantum Design) also at the QM Core Lab at HZB. Specific heat measurements below 1 K with fields along the [111] crystallographic direction were performed using the thermal relaxation method at Augsburg University (see Appendixes A, B, and C for more details).

II.2 Neutron scattering

The single crystal inelastic neutron scattering experiments were conducted on the direct-geometry time-of-flight (tof) spectrometer Cold Neutron Chopper Spectrometer (CNCS) at Spallation Neutron Source (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, a single crystal of $\sim 2.5$ g 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 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 which cooled the sample down to 30 mK and data were collected at 30 mK and 20 K. The data at 20 K were used as the background for both experiments. The software packages Dave dave , Mantid mantid and Horace horace were used for data processing.

Polarized neutron diffuse scattering with $Z$-polarization analysis was performed on the DNS diffractometer at FRM2, Munich. The single crystal was mounted on a dilution refrigerator. The incident neutron wavelength used was 3.3 Å. Data were collected at $\sim 300$ mK and 23 K (background reference), rotating the sample through 160 degrees in steps of 1 degree.

II.3 Spin wave and thermodynamic property simulations

The spin waves are calculated based on the linear spin wave theory with Holstein-Primakoff and Bogoliubov transformations Benton2016 . The Matlab package SpinW is also used for the spin wave calculation spinw .

Susceptibility and specific heat were analysed using the numerical linked cluster expansion (NLCE), which is a method of generating a series expansion for quantities in the thermodynamic limit from exact diagonalisation of small clusters. The series is guaranteed to converge at high temperature and for high magnetic fields but can be useful outside these limits (Appendix E) Tang2013 ; Singh2012 ; Applegate2012 .

The temperature-field phase diagram was simulated using the classical Monte Carlo method. The simulation is done on a cubic cluster of spins with periodic boundary conditions and the spins are treated as classical vectors of fixed length $|\bm{\tau}|=1/2$ (see Appendix F for details).

III Results of inelastic neutron scattering

Figure 1 shows colour-coded intensity maps of constant energy slices through the INS data measured on CNCS in the (HHL) reciprocal plane at 240 mK in the ordered phase (after background subtraction). Consistent with the data in Ref. Petit2016 , we see a highly structured pattern with the symmetry of the scattering plane which evolves with increasing energy transfer: from a gapped pinch point pattern at low energy to a pattern with intensity at the wavevectors [220] and [113] at high energy. The intensity is strongest at around $0.075$ meV where a pinch point pattern is observed. The energy-integrated polarized neutron diffraction data (Fig. 2) shows a pinch point pattern in agreement with the INS data and the data in Ref. Petit2016 . This pinch point pattern is gapped which is clearly shown in Ref. Petit2016 and in our high-resolution data measured on Osiris at 30 mK (Fig. 3 and 4). The gapped pinch point pattern is related to divergence-free spin dynamics Benton2016 ; Yan2018 ; Tomonari2018 .

Above the energy of the pinch point mode, there are modes with dispersions starting from the arms of the pinch point pattern and ending at the wavevectors of the AIAO Bragg peaks at $\sim 0.27\,$meV. These dispersing modes are recently recognized as dispersing pinch points due to curl-free spin dynamics Yan2018 ; Tomonari2018 . An integration of the INS data over energy range [0.09, 0.30] meV above the flat pinch point mode, also shows a pinch point pattern with a geometry orthogonal to the flat pinch point pattern at $\sim 0.075$ meV (Fig. 3). This confirms the theoretical prediction Yan2018 ; Tomonari2018 .

IV Spin Hamiltonian and spin wave simulations

The symmetry-allowed nearest-neighbor spin Hamiltonian for the dipolar-octupolar doublet on the pyrochlore lattice is given as Huang2014

[TABLE]

where $\tau^{\alpha}_{i}$ is the $\alpha$ component of the pseudospin-1/2 at site $i$ and $J_{\alpha}$ is the corresponding nearest-neighbor exchange constant ($\alpha=x,y,z$). The pseudospins are defined in the local coordinate frames with the local [111] crystallographic directions as the $z$ axes Huang2014 ; Benton2016 . Due to the peculiar symmetry of the dipolar-octupolar doublet, the interactions are uniform for every bond and the cross-coupling $J_{xz}$ can be removed by a pseudospin rotation around the local $y$ axes by a angle $\vartheta$ which leads to the $XYZ$ model Huang2014 ; Benton2016

[TABLE]

The relations between the exchange parameters and the pseudospin $\bm{\tau}_{i}$ in the original and in the rotated local frames are Benton2016 :

[TABLE]

[TABLE]

As pointed out in Ref. Benton2016 , the phase diagram and spin dynamics in zero field are determined only by the exchange parameters $\{\tilde{J}_{\tilde{x}},\tilde{J}_{\tilde{y}},\tilde{J}_{\tilde{z}}\}$ and the rotation angle $\vartheta$ controls how the pseudospin couples to an external probe because the magnetic moment (along the local [111] directions) is given by

[TABLE]

where $g_{zz}$ is the only non-zero component of the Ising anisotropic $g$ tensor. Although the system has a strong Ising anisotropy (only $g_{zz}$ is non-zero), the exchange interaction is non-Ising. As a result, transverse spin fluctuations (magnons) with respect to the AIAO ordering direction $\tilde{z}$ may exist, and are visible to neutrons as their projections to the local $z$ axes are non-zero Benton2016 . As such, despite the presence of strong Ising anisotropy in the system, the neutron scattering data could be analyzed within the framework of spin-wave theory. The neutron scattering data in zero magnetic field in arbitrary units can only determine the three exchange parameters. $\vartheta$ can be determined from the overall neutron scattering intensities in absolute units, the Curie-Weiss temperature, the ordered magnetic moment or the field responses of the system (see Sec. V).

Single crystal inelastic neutron scattering data of Nd2Zr2O7 show gapped magnon excitations in the AIAO ordered phase. However, due to the low energy scale of the system, it is difficult to determine the magnon dispersions accurately and thus extract the exchange parameters. According to Ref. Benton2016 , the gap to the flat mode is given by

[TABLE]

which provides a constraint for the data analysis. We further extracted another two equations for two characteristic energies which have a simple relation with the exchange parameters, namely the energy $\Delta_{2}$ at Brillouin zone boundary (e.g. $\bm{Q}=[110],[112]$ or $[221]$) and the highest energy of the dispersion $\Delta_{3}$ at the Brillouin zone center,

[TABLE]

With Eqs. 6 to 8, it is convenient to extract $\{\tilde{J}_{\tilde{x}},\tilde{J}_{\tilde{y}},\tilde{J}_{\tilde{z}}\}$ unambiguously, provided that the three energies can be determined accurately. Fitting our high-energy-resolution data measured on Osiris (Fig. 4), we have

[TABLE]

Accordingly, the extracted exchange parameters are

[TABLE]

Another solution is obtained by swapping the values of $\tilde{J}_{\tilde{x}}$ and $\tilde{J}_{\tilde{y}}$ but it produces the wrong neutron scattering intensity. As shown in Fig. 1, Fig. 4(c) and in Appendix D, the calculated spin wave scattering patterns agree well with the data in zero field.

Furthermore, the high-temperature series expansion shows that the high-temperature specific heat can be described by Baker1967

[TABLE]

with

[TABLE]

where $R$ is the molar gas constant and $k_{\text{B}}$ is Boltzmann’s constant. As shown in Fig. 5, the calculated heat capacity (without any fitting parameters) agrees well with the experimental data above 4 K (the nuclear lattice contribution has been subtracted) which validates the refined spin Hamiltonian.

V Determining $g_{zz}$ and $\vartheta$

The determination of $g_{zz}$ and $\vartheta$ is carried out with $\tilde{J}_{\tilde{x},\tilde{y},\tilde{z}}$ fixed to the values in Eq. 10. We establish the values of these two parameters by analysing two quantities:

• •

the inverse magnetic susceptibility $1/\chi(T)$,

• •

the ground state ordered moment $m_{\text{ord}}$.

The inverse susceptibility is calculated using numerical linked cluster expansion. For the calculation we have used NLCE up to third order. Fig. 6 shows the total squared error between the experimental and calculated susceptibility as a function of $g_{zz}$ and $\vartheta$ (a constant van Vleck term was included in the fit and optimized independently for each value of $g_{zz}$ and $\vartheta$). There is an extended minimum in the total squared error. To fix the parameters we use a further piece of information: the ground state ordered moment.

The ordered moment for the dipolar AIAO order (with pseudospins oriented along the local $\tilde{z}$ axes) is given by

[TABLE]

where $\delta\tilde{\tau}^{\tilde{z}}$ is the moment reduction from zero point quantum fluctuations. We calculate $\delta\tilde{\tau}^{\tilde{z}}\approx 0.05$ from linear spin wave theory using the exchange parameters in Eq. 10.

Experimentally, we have the ordered moment at 0.1 K $m_{\text{ord}}=1.26(2)\mu_{\text{B}}$ according the neutron diffraction experiments in Ref. Xu2015 . The region of parameter space which is consistent with this is marked by the dashed white lines in Fig. 6. The overlap between this region and the extended minimum in the squared error between experimental and theoretical inverse susceptibility gives our estimates:

[TABLE]

The resulting fit to the inverse susceptibility is shown in Fig. 7. The obtain $g_{zz}$ accounts for the measured saturated magnetization for the single crystal (Appendix B). Table 1 presents a comparison of the determined spin Hamiltonian with the ones in previous works. Our parameters are consistent with Ref. Lhotel2018 published recently. The major difference comes from $g_{zz}$ and $\vartheta$ which is because a fixed $g_{zz}$ was used in Ref. Lhotel2018 rather than being fitted.

In addition, the Curie-Weiss (CW) temperature is given in Ref. Benton2016

[TABLE]

It yields $T_{\text{CW}}=0.28(4)$ K which is consistent with the experimental data 0.31(1) K (Appendix B), solving the puzzle that Nd2Zr2O7 has an antiferromagnetic order but positive Curie-Weiss temperature Hatnean2015 ; Lhotel2015 ; Xu2015 ; Benton2016 .

Finally, the exchange parameters in the original local coordinate frames are

[TABLE]

VI Heat capacity in [111] fields

We have used the model parameterisation to calculate the heat capacity in a magnetic field oriented along the [111] direction using fourth order NLCE. The terms in the series expansion oscillate strongly as the temperature is decreased, but convergence is improved by using Euler transformations Applegate2012 . Fig. 8 shows the comparison between calculations and the experimental data at external fields of $H=0,0.2,0.5,1,2$ T.

The calculations at $H=0,0.2$ T show a sharp increase of the heat capacity as the temperature is lowered, consistent with the approach to a phase transition. As the phase transition at $T\approx 0.4$ K is approached, the agreement between NLCE calculation and experiment becomes worse, which is an expected consequence of the correlation length exceeding the size of the largest cluster used in NLCE. For higher fields $H=0.5,1,2$ T, there is smooth crossover from high to low temperature in both calculation and experiment. Thus the NLCE calculation reproduces the qualitative behaviour seen in experiment for all fields.

VII Monte Carlo phase diagram

In this section we present results for the magnetic phase diagram in an applied magnetic field determined using classical Monte Carlo simulations of the Hamiltonian determined above. We consider three directions of magnetic field: along the [100], [110] and [111] crystal axes. All-in-all-out order is detected using the order parameter

[TABLE]

where $N=65536$ which is the number of the spins in the simulation.

The phase diagrams as a function of field strength $H$ and temperature $T$ are shown in Fig. 9 where the value of the order parameter in the simulation is indicated by the color scale. The phase transition temperature for each value of $H$ is estimated by the position of a sharp peak in the susceptibility of $m_{\text{AIAO}}$. For fields along the [111] direction, simulations show a sharp peak in the order parameter susceptibility for fields up to $H\sim 0.3$ T, with the peak temperature being essentially field independent. For $H>0.3$ T, there is a smooth crossover into the low temperature phase. This is consistent with the experimental heat capacity data which shows that the sharp behavior at low field turns into a smooth crossover somewhere between 0.2 T and 0.5 T (Fig. 8).

For fields along the [100] direction, simulations show a sharp phase transition, with the transition temperature going to $T=0$ at $\sim 0.25$ T. For fields larger than this there is no phase transition. The measured magnetization and susceptibility in Refs. Lhotel2015 ; Opherden2017 indicate a qualitatively similar behavior but with the phase transition vanishing at the lower $H_{\text{c}}\sim 0.1$ T.

For fields along the [110] direction, simulations again show a sharp phase transition, with the transition temperature going to zero at 0.6 T. There is actually a small region of re-entrance, where there are two phase transitions for fields slightly above 0.6 T. It would appear, based on the neutron scattering data in a [110] field in Ref. Xu2018 that experimentally the AIAO order is removed by a field 0.25 T, therefore the classical simulations overestimate the stability of the low temperature ordered state against a [110] field.

Finally, we consider the issue of hysteresis for the [111] field. We have studied the magnetisation curves at $T=0.05$ K, by initializing the system from a polarised state and then equilibriating at $H=1$ T. The final configuration in the simulation is then used to initialise the simulation at the next (lower) value of field. In this way the simulation is able to mimic the history dependence of a hysteresis loop and we can compare the magnetisation curves obtained when sweeping the field upwards and downwards. For fields along the [110] and [100] directions [Figs. 10(a) and (b)] there is no difference between the up and down sweeps, and thus no hysteresis. For fields along the [111] direction [Figs. 10(c)] there is a large hysteresis loop which we will explain in the next section. This result compares reasonably well with experiments which also show pronounced hysteresis for fields along the [111] direction, but only very weak hysteresis for fields along the [100] direction Opherden2017 .

VIII Hysteresis and domain dependence in [111] fields

The magnetization and AC susceptibility in magnetic fields along the [111] direction show abrupt changes as a function of the field, which were associated with a field-induced transition and changes in the domain structure Lhotel2015 ; Opherden2017 . Here we present the microscopic mechanism of the field-induced transition which reveals the existence of metastable states and domain dependence behaviors of the system in [111] fields.

First, the pyrochlore lattice can be viewed as stacking of kagome and triangle planes along the [111] direction as shown in Fig. 11. Second, there are two degenerate all-in-all-out orders, namely all-in-all-out (AIAO) and all-out-all-in (AOAI), leading to the possibility for two types of magnetic domain in a sample cooled in zero field. The two types of domains are not equivalent with respect to a field along the [111] direction. As shown in Fig. 12, a [111] field flips the moments on the kagome lattice for the AIAO domain but flips the moments on the triangular lattice for the AOAI domain.

In mean field theory, assuming site-independent ansatz within the two types of layers in [111] fields at zero temperature, the energy on the pyrochlore lattice can be simplified to be based on a single tetrahedron given by

[TABLE]

where the spins are taken as classical vectors $\tau[\cos(\phi_{i}),0,\sin(\phi_{i})]$ with $\tau=1/2$ and $\phi_{1}$ and $\phi_{2}$ are the spin canting angles with respect to the ${\tilde{z}}$ axes in the ${\tilde{x}}$-${\tilde{z}}$ planes for the kagome and triangular lattices and $\cos\alpha=1/3$ and $\cos\beta=1$ with $\alpha$ and $\beta$ being the angles between the local [111] axes and the [111] field for the kagome and triangular lattices. Here we have only two angular variables for the magnetic lattice with four sublattices because there only two types of sites (in the kagome and triangle planes) in the [111] field. The non-magnetic $\tilde{\tau}^{\tilde{y}}$ component is excluded for the current situation because $\tau_{\tilde{y}}$ is not coupled to the field directly ($g_{yy}=0$) and $\tilde{J}_{\tilde{y}}$ is small compared with the other terms.

Fig. 12(a) and (b) shows the calculated magnetisation in increasing and decreasing [111] fields obtained by optimizing the energy (Eq. 18) locally for the two types of domains (the optimized state in the previous field is used as the start for the current field). The calculation is consistent qualitatively with the susceptibility data in Ref. Opherden2017 and the Monte Carlo simulations in Sec. VII. We found that the hysteresis only happens to the AIAO domain. Fig. 12(c) shows the moment configurations in the progress of changing the field for the AIAO and AOAI domains. It reveals that the flip of the moments on the kagome layers for the AIAO domain causes the sudden jump in the magnetization and in decreasing field the moments on the triangular layers reverse their direction smoothly recovering to an AOAI state (not the initial AIAO state) leading to the irreversibility. On the contrary, the AOAI domain always has the moments on the triangular layers flipped and reversed smoothly for increasing and decreasing fields. This result is consistent with the recent single crystal neutron diffraction experiments in magnetic fields Lhotel2018 and supports the proposed domain inversion explanation for the field dependence of the magneto-resistivity in the pyrochlores containing 5$d$ elements Tardif2015 ; Fujita2016 ; Tian2016 ; Fujita2018 ; Ma2015 .

To investigate the irreversibility further, we calculated the energy landscape for field along the [111] direction. Fig. 13 shows the energy landscape as a function of the spin orientations on the kagome and triangular layers in the [111] field of 0.8 T (lower than the mean-field $H_{\text{c}}\sim 1.2$ T) with the exchange Hamiltonian determined above. As we can see, the AOAI spin configuration is smoothly connected to the global energy minimum corresponding to a field polarized 3-in-1-out state and there is an energy barrier in-between for the AIAO domain. Therefore the AIAO state is a metastable state in the [111] field. In an increasing field, when the energy barrier vanishes at the critical field, the AIAO state transforms to the ground state leading to a sudden change of the spin configuration. When decreasing the field, the polarized state returns to AOAI state which is smoothly connected to the polarized state in energy.

The recovery to the AOAI single-domain state from the field-polarized 3-in-1-out state also can be understood by considering the exchange fields on the kagome and triangle layers. In the field-polarized state, the spins on the triangle layers are subjected to a stronger molecular field against the applied magnetic field because they have six broken bonds with their neighbors from the kagome layers while the spins on the kagome layers only have two.

The domain-dependent response to the [111] field has significance for the spin dynamics, i.e. dynamical kagome spin ice observed in Nd2Zr2O7 in [111] fields Lhotel2018 . First, the two domains have different spin configurations in fields below the critical field and are expected to show different spin dynamics. Our spin wave calculation shows that the gap to the dynamical kagome ice modes exhibits opposite behaviors for the two types of domains: decreasing for the AIAO domain and increasing for the AOAI domain. This may contribute to the broadening of the line width in the INS data Lhotel2018 if the sample is cooled in zero field and contains two types of domains. Second, above $H_{\text{c}}$ only the AOAI domain exists in the sample based on which the spin wave simulation should be done specifically though the AIAO domain is still present in mean-field theory. As shown in the Fig. 14, the spin wave of the AOAI domain shows a better fit to the data in Ref. Lhotel2018 (the experimental $H_{\text{c}}\sim 0.1$ T). However, there are still discrepancies especially in scattering intensity between the data and theory which is attributed to the limitation of mean field theory and the possible quantum fluctuations in the system. Third, it was shown that the two-dimensional dynamical kagome ice is quite robust appearing both below and above the critical field Lhotel2018 . This is because the spin configuration on the kagome lattice is always 3-in-3-out in [111] fields.

IX Conclusions

We refined the spin Hamiltonian by the combined analyses of the high-energy-resolution inelastic neutron scattering data, susceptibility, magnetisation and specific heat of Nd2Zr2O7 using linear spin wave theory and numerical linked cluster expansion. Using classical Monte Carlo simulation, we have calculated the magnetic phase diagram in the presence of applied magnetic fields in the [100], [110] and [111] directions which qualitatively agrees with the experimental phase diagram and reproduce the large hysteresis seen for fields along [111]. However, simulations overestimate the stability of the low temperature all-in-all-out order against external magnetic field, compared to experimental results. This may be related to quantum effects as shown in Yb2Ti2O7 Hitesh2017 .

The microscopic mechanism for the hysteresis is revealed in our mean-field calculation to be the existence of a metastable state and the AIAO-AOAI domain inversion in [111] field. Including the domain-dependent response to the [111] field, the observed dynamical kagome ice can be well described though there are still disagreements Lhotel2018 . In addition, the microscopic mechanism for the domain inversion could also be applied to the Nd-containing 5$d$ pyrochlores which shows hysteresis in the magneto-resistivity as a function of the field Tardif2015 ; Fujita2016 ; Tian2016 ; Fujita2018 ; Ma2015 and may also be related to the anomaly in the AC susceptibility of Nd2Hf2O7 Opherden2018 .

Acknowledgements.

We thank Y.-P. Huang, M. Hermele, S. T. Bramwell and A. T. Boothroyd for helpful discussions on the related theory. We acknowledge Helmholtz Gemeinschaft for funding via the Helmholtz Virtual Institute (Project No. VH-VI-521). This research used resources at the Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory. Experiments at the ISIS Neutron and Muon Source were supported by a beamtime allocation RB1810504 from the Science and Technology Facilities Council (DOI: 10.5286/ISIS.E.92924095).

Appendix A X-ray diffraction and crystallography

The structure of the single crystal was characterized using powder X-ray diffraction (XRD) (Bruker-D8, Cu-$K_{\alpha}$) and Rietveld refinements using the software Fullprof Suite fullprof . The XRD pattern and refinement are shown in Fig. 15. The crystallographic parameters are $a=10.651(1)$ Å  and $x_{\text{O1}}=0.3347(1)$ which are consistent with the previous reports for the powder and single crystal samples Xu2015 ; Lhotel2015 ; Hatnean2015

Appendix B Susceptibility and magnetization

Figure. 16 shows the susceptibility of the single crystal sample of Nd2Zr2O7 above 2 K. The susceptibility increases with decreasing temperature and shows no anomaly above 2 K. The $\chi(T)$ data were fitted by the Curie-Weiss law with the Van Vleck term $\chi(T)=\chi_{0}+C/(T-T_{\text{cw}})$ for the temperature range 10 $\leq T\leq 30$ K. This yields $T_{\text{CW}}=0.31(1)$ K, \mu_{\text{eff}}=2.47(1)\,$$\mu_{\text{B}}/Nd and $\chi_{0}=3.5(1)\times 10^{-3}\,$emu/mol Nd, which are consistent with those reported for powder and single crystal samples Xu2015 ; Lhotel2015 ; Hatnean2015 . The positive $T_{\text{CW}}$ indicates an effective ferromagnetic interaction between the Nd3+ moments though the sample shows an “all-in-all-out” antiferromagnetic order. This is a consequence of the dipolar-octupolar nature of the effective spin of the Nd3+ ion in pyrochlores Benton2016 .

Figure. 17 shows the magnetization $M(H)$ data for the single crystal with field along the three main cubic directions at 2 K. The magnetization is highly reduced (comparable to the free-ion value) and anisotropic, consistent with the previous reports Lhotel2015 ; Hatnean2015 . The magnetizations with fields applied along the [100] and the [110] directions are the highest and the lowest, respectively, reminiscent of the spin ice materials, consistent with the Ising anisotropy of the Nd3+ moment found in the crystal field analyses Xu2015 ; Lhotel2015 .

After subtracting the Van Vleck susceptibility, the data shows saturated magnetizations $M_{\text{s}}$ 1.39, 1.02 and 1.21 $\mu_{\text{B}}$/Nd for the three directions. For a pyrochlore with local [111] Ising anisotropy, $M_{\text{s}}$ for the three directions should be $\mu(1/\sqrt{3})$ , $\mu(\sqrt{2/3}\times 2)/4$ and $\mu(1+1/3\times 3)/4$, respectively where $\mu$ is the magnetic moment of Nd3+ Fukazawa2002 . With the refined \mu=g_{zz}\mu_{\text{B}}/2=2.5(1)\,$$\mu_{\text{B}}, we have 1.44(5) 1.02(4) and 1.25(5) $\mu_{\text{B}}$ which is consistent with the data within 5% deviation.

Appendix C Specific heat

Figure 18 shows the specific heat $C_{\text{p}}$ data of single crystal Nd2Zr2O7 and powder La2Zr2O7 samples and the data for powder Nd2Zr2O7 in Refs. Blote1969 ; Lutique2003 . Above 10 K, the specific heat of the two compounds are nearly the same because of their similar structures and thus the La2Zr2O7 data was used as a non-magnetic background (at higher temperatures they are different due to the CEF effect) Xu2015 . At low temperatures, Nd2Zr2O7 shows a $\lambda$-shape peak at T_{\text{N}}$$\,\approx 0.4\,K similar to the powder sample Blote1969 ; Lutique2003 . The upturn below 0.1 K is caused by the nuclear hyperfine interactions whose contribution is $\sim T^{-2}$ in the high temperature region Blote1969 . The magnetic 4$f$ electrons contribute to the $C_{\text{p}}$ peak, which originates from the magnetic correlations and excitations.

After subtracting the phonon and nuclear contributions, the calculated magnetic entropy is $\sim 0.93\,R\ln(2)$/mol Nd below 10 K. The entropy released due to the phase transition is close to $R\ln(2)$/mol Nd, indicating the establishment of long-range order in a spin-1/2 system, consistent with the neutron diffraction and crystal field analysis in Refs. Lhotel2015 ; Xu2015 .

For a normal antiferromagnet with a linear magnon dispersion at low energy, the C_{\text{p}}$$(T) should show a $T^{3}$ temperature dependence at sufficient low temperature below $T_{\text{N}}$. It was reported that the $T^{3}$ law applies to the low temperature $C_{\text{p}}$ of Nd2Sn2O7 which has the same magnetic order as for Nd2Zr2O7 Bertin2015 ; Xu2015 . However, the inelastic neutron scattering data of Nd2Zr2O7 clearly shows gapped magnon excitations. Therefore, the $T^{3}$ law is invalid for Nd2Zr2O7. The model $\sim T^{-2}\exp(-\Delta/T)$ ($\Delta$ is the magnon gap size) is normally used to describe the temperature dependence of the specific heat of gapped magnon excitations Quilliam2007 . It would be a straight line in the $\log(C_{\text{p}}T^{2})$-$1/T$ plot. The data below 0.18 K is fitted to $A\,T^{-2}+B\,T^{-2}\exp(-\Delta/T)$ (Fig. 18) which yields $\Delta=0.75(3)\,$K. The obtained gap is quite close to the measured one ($\approx 0.075\,\text{meV}\,=\,0.87\,\text{K}$) as shown in Sec. IV (see Ref. Xu2017 for details).

Appendix D Spin wave fitting

Figure 19 shows the INS data measured on CNCS comparing with the calculated scatering pattern using linear spin wave theory. Fig. 20 shows the magnon dispersion plotted over the INS data of Ref. Petit2016 .

Appendix E Numerical Linked Cluster Expansion (NLCE)

Figure 21 shows the clusters used in the NLCE calculation. We have used Numerical Linked Cluster Expansion (NLCE) to calculate the magnetic susceptibility and heat capacity for comparison with experiment. An introduction to the NLCE method is given in [Tang2013, ].

NLCE calculates extensive quantities per site $\frac{\langle\mathcal{O}\rangle}{N}$ as sums over contributions from clusters $c$ which can be embedded in the lattice

[TABLE]

$M(c)$ is the number of times $c$ can be embedded in the lattice, divided by the number of sites $N$. $W(c)$ is the weight of the cluster:

[TABLE]

where $\langle\mathcal{O}\rangle_{c}$ is the expectation value of $\mathcal{O}$ calculated from exact diagonalization of cluster $c$ with open boundary conditions. The second term in Eq. (20) is a sum over the weights of all subclusters of $c$.

We use a series of clusters beginning with a single site (zeroth order NLC, NLC0) and then all subsequent clusters are constructed from full tetrahedra. The $n^{\text{th}}$ order of NLC (NLCn), calculates the sum in Eq. (19) up to clusters of $n$ tetrahedra. Our calculations of the magnetic susceptibility (Fig. 7) were made up to third order (NLC3) and calculations of the heat capacity (Fig. 7) were made up to fourth order (NLC4).

In the absence of a magnetic field, there is only one kind of cluster at each order up to NLC3, and two distinct clusters at NLC4. These are illustrated in Fig. 21. In the presence of an applied magnetic field along the $111$ direction there are two distinct types of cluster which must be taken into account in each of NLC2 and NLC3 and six distinct types of cluster in NLC4.

To improve convergence, we have used Euler transformation on the heat capacity calculations Applegate2012 . The Euler transformed result at third order is:

[TABLE]

where $\langle\mathcal{O}\rangle_{\sf NLCn}$ is the estimate of $\langle\mathcal{O}\rangle$ up to $n^{\text{th}}$ order in NLC. The Euler transformed result at fourth order is

[TABLE]

Appendix F Monte Carlo Simulations

Classical Monte Carlo simulations have been used to calculate the magnetic phase diagrams in applied field (Fig. 9) and hysteresis curves (Fig. 10), for the [111], [110] and [100] field directions. In these simulations the spins are treated as classical vectors with fixed length $|{\bf S}|=1/2$. The simulations use the standard Metropolis algorithm, with single spin updates using the Marsaglia method marsaglia72 . Simulations are performed on cubic clusters of $N=16L^{3}$ sites, with $L$ being the number of 16-site cubic unit cells along each side of the cluster. We define one Monte Carlo step (MCS) as one sweep of the whole lattice of $N$ sites, attempting an update at each site.

Details of the simulations for the phase diagrams are given in Appendix F.1 and for the hysteresis loops in Appendix F.2.

F.1 Determination of phase diagrams

The phase diagrams in Fig. 9 were determined using simulations on an $L=16$ ($N=65536$) cluster. The simulations were performed separately for each value of applied field and followed the following protocol:

The system was initalized to a random spin configuration. 2. 2.

The system was then equilibriated at $T=1$ K using 5000 MCS without measurements being taken. 3. 3.

Observables are then averaged over the course of 50000 MCS, with each observable measured every 50 MCS. 4. 4.

The temperature is then reduced by $\delta T=0.0202$ K, and the system is equilibriated at the new temperature using 500 MCS. 5. 5.

The cycle is repeated down to $T=0.01$ K.

F.2 Hysteresis

The magnetic curves in Fig. 10 were determined using simulations on an $L=4$ ($N=1024$) cluster. The simulations for the down sweep followed the following protocol:

The system was initalized to a spin configuration fully polarized along the direction of the field. 2. 2.

The system was then equilibriated with an external field of $H=1.0$ T and temperature $T=0.05$ K using 500 MCS. 3. 3.

Observables are then averaged over the course of 100000 MCS, measuring every 100 MCS. 4. 4.

The magnetic field is then decreased by $\delta H=0.0202$ T, and the system is equilibriated at the new value of external field using 500 MCS. 5. 5.

The cycle is repeated down to $H=-1.0$ T.

The simulations for the up sweep are exactly the same, but initialising from a negatively polarised state and increasing the field from $H=-1.0$ T.