Evidence for dynamic kagome ice

The search for two-dimensional quantum spin liquids, exotic magnetic states remaining disordered down to zero temperature, has been a great challenge in frustrated magnetism over the last few decades. Recently, evidence for fractionalized excitations, called spinons, emerging from these states has been observed in kagome and triangular antiferromagnets. In contrast, quantum ferromagnetic spin liquids in two dimensions, namely quantum kagome ices, have been less investigated, yet their classical counterparts exhibit amazing properties, magnetic monopole crystals as well as magnetic fragmentation. Here, we show that applying a magnetic field to the pyrochlore oxide Nd2Zr2O7, which has been shown to develop three-dimensional quantum magnetic fragmentation in zero field, results in a dimensional reduction, creating a dynamic kagome ice state: the spin excitation spectrum determined by neutron scattering encompasses a flat mode with a six arm shape akin to the kagome ice structure factor, from which dispersive branches emerge.

They are written in the cubic F d3m structure of the pyrochlore lattice, for the site dependent local a i and b i vectors spanning the local (x, y) anisotropy planes. The CEF axes z i of the rare earth ions are perpendicular to these planes.
Supplementary Figure 1. Sketch of a tetrahedron in the pyrochlore lattice. The black arrows show the local zi anisotropy (crystal field -CEF) axes and the green disks represent the local (ai, bi) planes. Further details are given in Supplementary Table 1.
Supplementary Figure 2. Field dependence of the magnetic moments in the unit cell for H [111]. y (in blue) gives the moment of the apical spin, with its CEF axis along the field. x (in green) gives the moment of the three remaining spins (also called kagome). The magenta curve displays the magnetization obtained from macroscopic measurements. Red points correspond to the magnetization calculated using the refinements results and projected onto the field direction, leading to M (H) = (x − y)/4. The error bars are provided by the Rietveld refinement made using Fullprof.

Supplementary Note 1. Neutron diffraction
Experimental details: The neutron diffraction data were collected using the D23 single crystal diffractometer (CEA-CRG, ILL France) operated with a copper monochromator and using λ = 1.28Å. The single crystal of Nd 2 Zr 2 O 7 was glued on the Cu finger of a dilution insert and placed in a cryomagnet. The experiments have been conducted with the (vertical) field H parallel to the [111] direction. We first checked that the Bragg intensities remain zero on the forbidden q vectors of the F centred lattice. This implies that the field induced structures are described by a k = 0 propagation vector. The additional magnetic intensity is then observed on top of the nuclear Bragg peaks. A series of integrated intensities at 10 K and H = 0 was measured (about 60 reflections) as a reference. Refinements were performed using the Fullprof sofware suite 1 , working on the differences with the T = 10 K data set. Finally, for a set of chosen Bragg peaks, we ramped the field back and forth, up to 1 and back to -1 T, yielding a precise evolution of the magnetic intensities. The 1 T field was chosen to saturate the sample.
Model of the magnetic structure in a [111] field: Since the propagation vector is k = 0, the four tetrahedra of the crystalline unit cell host identical arrangements of the spins. Furthermore, because of the strong Ising anisotropy of the Nd 3+ ions, it is reasonable to consider that the magnetic moments are forced to align along the local CEF axis (labelled z i , see Supplementary Figure 1 and Table 1). In our refinement, we thus assume the following model (see Table 1): where x and y are the fitted parameters. Note that m 2 is parallel to [111], hence to the applied field. From a physical point of view, situations where x and y have the same sign correspond to a generalized all in -all out structure. When x and y have different signs, however, the model describes a structure which resembles the 1 out -3 in, but m 2 and m 1,3,4 may have different amplitudes. The magnetization per Nd 3+ ion is given by: Projected along the [111] field, it simply becomes M = (x − y)/4. The results obtained increasing the field from 0 to 3 T along with the calculated magnetization are shown in Supplementary Figure 2. The excellent agreement with the macroscopic magnetization is also shown for comparison. It is worth noting that the magnetic structure may also be described in terms of generalized charges Q living on the dual (diamond) lattice defined by the centres of the tetrahedra. Q is a vector with three components defined by: In zero field, the all in -all out structure can be considered as a Q z charged staggered pattern.
Analysis of the ramps: We now turn to the analysis of the intensities measured while ramping the magnetic field (see Supplementary Figure 3). To this end, we first write the magnetic structure factor for a wavevector q : where R i denotes the position of the i th spin in the unit cell and f i is the form factor. The magnetic Bragg intensity is then given by: For q = (2, 0, 0), (1,1, 1) and (1, 3,1), analytic calculations show that: where V is the volume, I N q denote the corresponding nuclear contributions and the coefficients are given by: A = 174.55, B = −11.64 and C = 34.91. In these expressions, the form factor has been neglected since it is close to 1 for those q vectors. It is then convenient to consider the new intermediate variables e and t defined as: hence: Note that e = 0 at 0 T since the structure is all in -all out (x = y). We now consider reduced quantities where the o and 1 superscripts denote respectively the values at 0 and 1 T. We then have = η e 2 t 2 (A + 2B + C) + 2(B + C)t + Ct 2 − ζ .
As a result, the variable e is readily obtained from the field evolution of the (2, 0, 0) or (1,1, 1) intensity: Note that the data are consistent, i.e. i 200 and i 111 give the same information. Finally, we need to solve the following equation, deduced from the expression for i 131 : For each field value, or in other words, for each e value implicitly determined from i 200 , there exist two solutions t 1 and t 2 , hence two solutions for the y/x = (1 + t) ratio. These solutions take the form of the black and grey branches shown in Supplementary Figure 4a. The final choice between the two possibilities relies on physical grounds and is shown by the red points in Supplementary Figure 4a: when the field is strongly negative or positive, a 3 in -1 out structure is expected, i.e. y/x ≤ 0. This corresponds to the black branch in Supplementary Figure 4a. In zero field, the all in -all out structure is stabilized, which corresponds to y/x = 1, hence to the point where the black and grey solutions coincide in Supplementary Figure 4b. Since x and y should be regular functions of the field, except close to the transition at 0.08 T where a discontinuity can be expected, the physical solution stays on the black branch. This holds up to the transition at 0.08 T. Above, there is a jump to join the grey solution. This gives the field evolution for the x and y parameters shown in Supplementary Figure 4c. The blue arrow marks the discontinuity at the transition. As a very good energy resolution is needed (about 20 µeV), we used a wavelength λ = 8.5Å. The data were then processed with the horace software 2 , transforming the recorded time of flight, sample rotation and scattering angle into energy transfer and 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. We finally prepared a set of constant q scans by integrating over a small q range with ∆h = ∆k = 0.05. We also prepared a series of constant energy maps (ranging from E = 0.05 up to 0.25 meV) that were symmetrized on the basis of the expected 6-fold symmetry of the intensity in this scattering plane. These maps are displayed in Supplementary Figure 7 and 8, giving an overview of the spin dynamics for a series of applied magnetic fields. Note that those maps were obtained by integrating over a range ∆ω = 10 µeV, which is about half the energy resolution. Supplementary Figure 6 shows the corresponding q-integrated spectra. Finally, Supplementary Figure 9  Above µ 0 H = 0.25 T, a different picture sets in. First of all, the structure factor of the flat mode evolves: the intensity close to q = 0 strongly decreases while 6 arms take shape. With increasing field, the intensity of the arms strongly decreases (see the first row of Supplementary Figure 7; note also the change in the color-scale). In addition, two narrow and dispersive branches, labelled A and B, can be observed (see Supplementary Figures 9 and  10). The A branch is descended from the zero field dispersive branch. For larger fields, this branch is still there but is quite difficult to observe since it becomes extremely weak. However, by taking the q-average of the scattering (see Supplementary Figure 6), which basically gives the density of states (enhancing flat modes), it becomes clearly visible. From this data, we can conclude that the A branch bandwidth becomes narrower, and progressively shifts towards higher energies.
The second dispersing B branch emerges from the kagome ice pinch point positions ((−4/3, 2/3, 2/3)-like points) above the 6 arms mode and closes at (−2, 2, 0) (and symmetry related points). With increasing the field, this B branch develops up to about 0.10 meV at µ 0 H = 0.25 T. The top of the band slightly shifts to higher energies with further increasing the energy. At µ 0 H = 0.75 T, for instance, it shifts up to 0.13 meV.
In this expression, τ i is a pseudo-spin that resides on the sites of the pyrochlore lattice. The (x, y, z) coordinates refer to the site dependent frames, where z is the local 111 axis (see also Supplementary Table 1). Note that the z component is related to the physical magnetic moment m = g τ z z, while the x and y components are non observable quantities. They respectively transform as dipolar and octupolar moments under symmetries. A rotation by an angle θ in the (x, z) plane allows one to define a new (x,z) frame, new pseudo-spin components (τx,τỹ,τz) and new coupling constants (J x ,J y ,J z ) so that the Hamiltonian becomes diagonal: The rotation matrix is defined as: With this definition, the angle between the z andz axes is θ. The following relations are especially useful: Importantly, the energies of the spin wave modes solely depend on the values of the three parametersJ x ,J y andJ z . This does not mean that J xz or the angle θ do not play any role; the latter especially appears when projecting the pseudo-spin onto the local z magnetic direction. As a result, it enters the spin correlation functions (i.e. the inelastic neutron scattering cross section) as well as the magnetization along with its field dependence.
In Supplementary Reference 4, we assumed J x = J xz = 0, following arguments based on the analysis of the crystal field. Furthermore, we found that: allow us to reproduce quantitatively the INS spectra. These parameters also led to an all in -all out ordering along the y octupolar axis. Note that at the classical level, this particular ordering is stabilized provided that: Spin wave calculations showed then that the spectrum encompasses two contributions, a flat mode E o at: along with dispersive branches, in excellent agreement with INS results. It is worth emphasizing that this analysis, however, does not explain the origin of the magnetic all in -all out ordering (along the z axis). It was actually designed in the framework of the fragmentation scenario, following the idea that the magnetization can be decomposed into two independent fragments, the divergence full and the divergence free emergent fields of a Helmholtz-Hodge decomposition. The above description, and especially the condensation of the octupolar ordering, was thus thought to apply to the divergence free part only.
In Supplementary Reference 5, O. Benton followed a different route. Taking advantage of a global symmetry of the Hamiltonian, he proposed to inter-change the axes compared to Supplementary Reference 4 so that: and introduced a non zero angle θ = 0.83(48 • ) chosen to reproduce the Curie-Weiss temperature. This set of parameters corresponds to: For this choice, the classical ground state is an all in -all out configuration with respect to the rotatedz axis. The spins τ thus point along thez axes. This rotated all in -all out order has projections onto the original magnetic z axis: and along the x axis: (note again that the latter cannot be observed directly with neutrons). INS results are equally well reproduced with these parameters since the spin wave energies do not depend on θ.

Supplementary Note 4. Fitting the INS data
In the following, we adopt the scheme discussed immediately above and in Supplementary Reference 5, but revisit the values of the parameters to fit the data taken in zero and in applied field. Our initial choice of θ is different, and is based on the following relation: where m ord = 0.8µ B is the experimental all in -all out ordered moment and m sat = 2.28 is the calculated saturation moment determined from crystal field coefficients 6 . This corresponds to an angle: To ensure that an all in -all out state is stabilized along thez axis, we assume that the coupling constants belong to the domain D defined byJ For couplings within D, the spin wave calculations 4,5 show that the spectrum encompasses a flat mode E o at: along with dispersive branches. To determine the values compatible with our experiments, we let the three parameters vary within D, for the fixed value of θ given above. Our goodness criterion is based on a χ 2 defined by the squared difference between the calculated and measured spin wave energies. In our definition, χ 2 takes into account q points along (−2, h, 2 − h) at H = 0 and 0.75 T. The inverse χ 2 is shown in Supplementary Figure 11. We find that the spin wave spectrum is well reproduced for a sub-domain D of D with (units are in K): but one is left with a quite large uncertainty regardingJ x andJ y . We note that the best χ 2 is obtained for J x ∼ 0.16,J y ∼ 0.97 and θ = 69 • , yet this choice does not predict the correct neutron intensities.
To solve this problem and resolve this uncertainty, we investigate more closely the sets of parameters found in D . To this end, we consider, for each set of coupling constants in D , a trial θ such that 0 ≤ θ ≤ π/2, and compare the calculated and measured magnetic moments for the different ions in the unit cell (see Supplementary Note 1 which discusses the neutron diffraction). A new χ 2 is defined as the sum of squared differences between those quantities for a series of fields between -1 and 1 T. For:  Figure 14).
Furthermore, the calculated magnetic structure experiences a field evolution which resembles the diffraction data. Supplementary Figures 12(b) and (c) compare the calculated magnetic moment of the apical Nd 3+ ion (y) and of the three kagome ions (x) with diffraction data, at 50 mK and 0.5 K respectively. While the overall agreement is good, there is a significant discrepancy between the experimentally observed field at which the abrupt transition occurs and the field value produced by the calculation at 50 mK. In the experiment, it is due to the transition from the all in -all out configuration to the 3 in -1 out structure and occurs at about 0.08 T. In the calculation, a similar abrupt transition also occurs, but at about 1 T at 50 mK and 0.35 T at 0.5 K. The origin of this transition is better understood when looking at the τ x components of the pseudo-spin. Supplementary Figure 12(d) shows that the transition occurs when the τ x component of the apical spins changes sign, becoming positive. In other words, it corresponds to the field where the apical spins are fully aligned along the z axis. Because of the primarily antiferromagneticJ z interactions, an all in -all out configuration is always favoured. Here the field already imposes a 3 in -1 out configuration for the z components, which is not favourable. A gain in exchange energy can be however achieved if the x components keep the same sign, hence precipitating the transition. Supplementary Note 5. Decomposition of the spectrum into divergence free and divergence full dynamical fields The decomposition of the spectrum in terms of (lattice) divergence free and divergence full dynamical fields has first been shown in Supplementary Reference 5. Here, we follow a similar route, expanding the calculations by taking into account the effect of the magnetic field. We first write the equations of motion in zero field taking advantage of the fact that the XYZ model is written using site dependent local bases: where the sum over j sites is restricted to nearest neighbours andJ is the exchange tensor. Here, m i is written in the local basis attached to site i, and not in the global Cartesian frame.
Our derivation assumes that the ground state is uniform, i.e. that the spins locally point along the same local direction, which is nothing but an all in -all out like configuration. In other terms, we assume here that the system is essentially an antiferromagnet. Calling m o the local equilibrium magnetization, we linearize those equations to obtain: Note that m o is indeed identical, whatever the site, since it is written in local bases.
Since m o is an equilibrium solution, we also note that: hence, m o andJ m o must be collinear. Writing the exchange tensor in the diagonal form (see section III): we observe that if m o is along one of these axes, for instancez, then In other words, m o andJ m o are indeed collinear. The energy is finally minimized when m o is along the axis corresponding to the smallest eigenvalue among (J x ,J y ,J z ).
Thanks to the pyrochlore structure, the sum over the neighbours sites j can be written in terms of the components of the generalized charge: defined in each tetrahedron (see Supplementary Figure 15). Note that by definition, those charges reside on the dual lattice of the pyrochlore, which is the diamond lattice. Since the latter is bipartite, we shall consider two dual sub-lattices, A and B, with: for A diamond sites and for B diamond sites. Interestingly, the uniform static solution m o can be understood as a staggered pattern, from the point of view of the charge, with Q = ±4m o . Finally, using these definitions of Q A,B , we shall write the equation of motion as: where the indexes i A and i B denote the A and B tetrahedra the site i belongs to. This form of the linearized equations of motion is the starting point to derive a peculiar partition of the solutions.
Charged solutions: We start with the subset of charged solutions. Forming the equations of motion for Q A,B using the previous equations, one obtains: Those solutions correspond to the propagation of the charge throughout the diamond lattice. We shall also make several points: • the A ↔ B symmetry is clear from the above equation.
• The sum over B refers to the four neighbouring B tetrahedra of a given A tetrahedron.
• Since for a lattice containing N site, there are (N/4 + N/4) A and B tetrahedra, the number of these modes is N/2. In the Helmholtz-Hodge decomposition of the m i field, these modes are obviously associated with the dynamical divergence full field.
• Finally, we write the exchange tensor in its diagonal formJ =J xxx +J yỹỹ +J zzz (assuming that m o is alongz, while (x,ỹ) span the plane perpendicular to m o ) and introduce the coordinates of the charge in the corresponding (x,ỹ,z) frame: The Qz A components remain constant (consistent with the all in -all out structure), while the Qx A and Qỹ A components evolve in time following: Alternate loops weathervane modes: The second set of solutions is constructed by introducing loops L consisting of N sites in the pyrochlore lattice. Note that the shortest of these loops are the hexagons that form the kagome planes. We consider: Forming the equation of motion for A, we obtain: We note that two successive sites belong, for example, to the same A tetrahedron but to different B ones: 1 A = Furthermore, since L is a loop, we have N B = 1 B, hence: Assuming alternate coefficients u k+1 = −u k , and u N = −u 1 , we obtain: This shows that the alternate loop solutions correspond to a resonating dispersionless field, decoupled from the charged solutions. To determine the corresponding energy, we write the exchange tensor in its diagonal form:J = J xxx +J yỹỹ +J zzz , and assume thatz is along m o . We thus have: The eigenvalues of above matrix are of the form ±i∆ with: and ∆ corresponds to the frequency of those dispersionless modes. Since we consider only the cases whereJ z < 0, we writeJ z = −|J z | and finally: which is the energy of the flat mode described in Supplementary notes 3 and 4.
Influence of a [111] magnetic field: Applying a magnetic field along one of the 111 high symmetry directions is especially interesting as the sites of the pyrochlore lattice separate in kagome (K) and triangular (T ) planes forming the subset of apical spins. The latter are rapidly polarized by the field since the Ising direction is precisely along H. The equations of motion can now be written as: and we still consider a linearized version around a uniform solution. However, as emphasized above, the field has different influences depending on the site: for the apical sites, g i H = −g H z while for the kagome sites, g i H = g H/3 z.
Our aim is to show that once the apical spins are fully polarized, the kagome ice displays, on one hand charged solutions, and on the other, flat alternate loop modes. Owing to the dimensional reduction, the charges are, however, kagome charges, and the loops are confined within the kagome layers. To this end, we introduce the kagome charge Q u=x,y,z = i∈ m u i of a given kagome triangle . It is worth noting that the same bipartite property holds in this 2D case. As a result, A and B triangles can be defined, in close analogy with the pyrochlore case (see Supplementary  Figure 16). Let us consider, on one hand, a site i where the i index runs in the kagome plane. This site belongs to an A and to a B kagome triangle, and has two neighbours sites in the T triangular planes, one in the top (j = t i ) plane and one in the bottom (j = b i ) plane. The equation of motion can be written as: (50) On the other hand, let's consider a site j where j runs in the triangular plane. This site is also the apical site of A and B kagome triangles, denoted j A and i B. With this notation, the equation of motion becomes: If the field is strong enough to fully polarize the apical sites, we may neglect the term j=ti,biJ δm T j , and hence we obtain: In this limit, we obtain a decomposition of the excitations similar to what we found in zero field, in terms of divergence full and divergence free fields, but within the kagome layers. The charged solutions are given by: and the dispersionless solutions by hexagonal loops of alternate spins: Supplementary Figure 15. Partial view of the pyrochlore lattice in the all in -all out magnetic structure depicted by the black arrows. The staggered charged pattern is shown using blue and red tetrahedra. The nodes of the diamond lattice, dual of the pyrochlore, are the centres of the tetrahedra. This dual lattice hosts the charge. It is bipartite, which allows one to distinguish A and B charges. Clearly, each pyrochlore site belongs to an A and to a B tetrahedron.