Emergent order in the kagome Ising magnet Dy3Mg2Sb3O14

Abstract

The Ising model--in which degrees of freedom (spins) are binary valued (up/down)--is a cornerstone of statistical physics that shows rich behaviour when spins occupy a highly frustrated lattice such as kagome. Here we show that the layered Ising magnet Dy₃ Mg₂ Sb₃ O₁₄ hosts an emergent order predicted theoretically for individual kagome layers of in-plane Ising spins. Neutron-scattering and bulk thermomagnetic measurements reveal a phase transition at ∼0.3 K from a disordered spin-ice-like regime to an emergent charge ordered state, in which emergent magnetic charge degrees of freedom exhibit three-dimensional order while spins remain partially disordered. Monte Carlo simulations show that an interplay of inter-layer interactions, spin canting and chemical disorder stabilizes this state. Our results establish Dy₃ Mg₂ Sb₃ O₁₄ as a tuneable system to study interacting emergent charges arising from kagome Ising frustration.

ARTICLE

Received 14 Jul 2016 | Accepted 4 Nov 2016 | Published 20 Dec 2016

DOI: 10.1038/ncomms13842 OPEN

Emergent order in the kagome Ising magnet Dy3Mg2Sb3O14

Joseph A.M. Paddison1,2, Harapan S. Ong1, James O. Hamp1, Paromita Mukherjee1, Xiaojian Bai2, Matthew G. Tucker3,4, Nicholas P. Butch5, Claudio Castelnovo1, Martin Mourigal2 & S.E. Dutton1

The Ising modelin which degrees of freedom (spins) are binary valued (up/down)is a cornerstone of statistical physics that shows rich behaviour when spins occupy a highly frustrated lattice such as kagome. Here we show that the layered Ising magnet Dy3Mg2Sb3O14 hosts an emergent order predicted theoretically for individual kagome layers of in-plane Ising spins. Neutron-scattering and bulk thermomagnetic measurements reveal a phase transition at B0.3 K from a disordered spin-ice-like regime to an emergent charge ordered state, in which emergent magnetic charge degrees of freedom exhibit three-dimensional order while spins remain partially disordered. Monte Carlo simulations show that an interplay of inter-layer interactions, spin canting and chemical disorder stabilizes this state. Our results establish Dy3Mg2Sb3O14 as a tuneable system to study interacting emergent charges arising from kagome Ising frustration.

1 Department of Physics, Cavendish Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge CB3 0HE, UK. 2 School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332, USA. 3 ISIS Neutron and Muon Source, Rutherford Appleton Laboratory, Harwell Campus, Didcot OX11 0QX, UK.

4 Spallation Neutron Source, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA. 5 NIST Center for Neutron Research, National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA. Correspondence and requests for materials should be addressed to J.A.M.P. (email: mailto:[email protected]

Web End [email protected] ) or to S.E.D. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842

The kagome latticea two-dimensional (2D) arrangement of corner-sharing trianglesis at the forefront of the search for exotic states generated by magnetic

frustration. Such states have been observed experimentally for Heisenberg14 and planar57 spins. If Ising spins lie within kagome planes and point either towards or away from the centre of each triangle, the potential for emergent behaviour is shown by considering a spin (magnetic dipole) as two separated and magnetic charges: the emergent charge T of a triangle is

dened as the algebraic sum over the three charges it contains (Fig. 1a)8. Ferromagnetic nearest-neighbour interactions favour T 1 states, yielding six degenerate states on each

triangle. This macroscopic ground-state degeneracy leads to a zero-point entropy S0E13 ln 92R per mole of Dy (where R is the molar gas constant), and suppresses spin order9, in analogy to three-dimensional (3D) spin-ice materials10,11. The long-range magnetic dipolar interaction generates an effective Coulomb interaction between emergent charges, driving a transition to an emergent charge ordered (ECO) state that is absent for nearest-neighbour interactions alone8,12. In this state, and charges

alternate, but the remaining threefold degeneracy of spin states for each charge means that spin order is only partial (Fig. 1b). The ECO state has two bulk experimental signatures: non-zero entropy S0E0.11R per mole of Dy12, and the presence of both

Bragg and diffuse magnetic scattering in neutron-scattering measurements13,14. Experimentally, kagome ECO states have been observed in spin-ice materials under applied magnetic eld15,16 and nano-fabricated systems in the 2D limit14,1719. However, a crucial experimental observation has remained elusivenamely, observation of the spatial arrangement of emergent charges in a bulk kagome material.

In this article, we show that an ECO state exists at low temperature in the recently-reported bulk kagome magnet Dy3Mg2Sb3O14 (ref. 20). Our experimental evidence derives from neutron-scattering and thermodynamic measurements, while Monte Carlo (MC) simulations reveal that this ECO state is stabilized by a combination of interactions between kagome layers, spin canting out of kagome layers and chemical disorder.

ResultsStructural and magnetic characterization. Structural and magnetic characterization suggests that Dy3Mg2Sb3O14 (ref. 20) is an ideal candidate for an ECO state. The material crystallizes in a variant of the pyrochlore structure (space group R 3m20) in which kagome planes of magnetic Dy3 alternate with triangular layers of non-magnetic Mg2 (Fig. 1c). X-ray and neutron powder diffraction measurements conrm the absence of a structural phase transition to t0.2 K (Supplementary Figs 1 and 2 and

Supplementary Tables 1 and 2) and reveal a small amount of site disorder in our sample, with 6(2)% of Dy kagome sites occupied by Mg (and 18(6)% of Mg sites occupied by Dy). Curie-Weiss ts to the magnetic susceptibility (Fig. 1d) yield a Curie-Weiss constant yCW 0.1(2) K for tting range 5rTr50 K,

consistent with ref. 20 (however, the value depends strongly on tting range). Demagnetization effects may also be signicant increasing yCW by 1.4 K in spin-ice materials21but cannot be quantitatively determined for a powder sample. The local Dy environment in Dy3Mg2Sb3O14 is similar to the cubic spin ice

Dy2Ti2O7 (ref. 22) (Supplementary Fig. 3), suggesting that Dy3 spins have an Ising anisotropy axis directed into or out of the kagome triangles with an additional component perpendicular to

+

a b

=

=

+

+

+3

+

=

=

3

c e

d

8

3

2

1

0 0 20 40

T (K)

4 6

5

4

3

2

1

0 0 2 4 6 8 10 12 14

6

(cm3 mol1 )

Dy

4

1 (mol Dy cm 3)

M( B per Dy)

2

60 80 1000

c

0H (T)

b

a

Figure 1 | Ising spins on the kagome lattice. (a) Relationship between spin vectors (arrows), magnetic dipoles (connected red and blue spheres) and emergent charge T of a triangle (labelled or 3). (b) Example of a microstate showing emergent charge order (ECO). (c) Partial crystal structure of

Dy3Mg2Sb3O14, showing kagome Dy1 xMgx site (blue spheres) and triangular Mg1 3xDy3x site (orange spheres), where x 0.06(2) for the sample of

Dy3Mg2Sb3O14 studied here. (d) Magnetic susceptibility data w(T) measured in an applied eld m0H 0.01 T after zero-eld cooling (left axis; black

squares), inverse magnetic susceptibility data w 1 (right axis; orange circles) and Curie-Weiss t over the range 5rTr50 K (blue line). (e) Dependence of magnetization M on applied magnetic eld m0H at different temperatures (labelled above each curve) and ts to the paramagnetic Ising model. Data are shown as solid coloured lines and ts as white dashed lines (note the nearly perfect agreement: as plotted the t lines are indistinguishable from the data).

In d,e, standard errors are derived from ts to the magnetization and are smaller than the symbol size or line width in the plots.

2 NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842 ARTICLE

the kagome planes. Experimentally, we conrm Ising anisotropy at low temperatures using isothermal magnetization measurements, which are ideally described by paramagnetic Ising spins with magnetic moment m 10.17(8) mB per Dy

(Fig. 1e). Moreover, our inelastic neutron-scattering measurements show that the ground-state Kramers doublet is separated from the rst excited crystal-eld state by at least 270 K (Supplementary Fig. 4), indicating that crystal-eld excitations are negligible at the low temperatures (r50 K) we consider.

Low-temperature spin correlations. The magnetic specic heat Cm(T) shows that spin correlations start to develop below 5 K and culminate in a large anomaly at T* 0.31(1) K that we attribute

to a magnetic phase transition (Fig. 2a and Supplementary Fig. 5). Below 0.20 K, the spins fall out of equilibrium, as is also reported in spin-ice materials23. In zero applied eld, the entropy change DSm(T) from 0.2 K to T 10 K is slightly less than the expected

Rln2 for random Ising spins; however, the full Rln2 entropy is recovered in a small applied eld of 0.5 T. The 0.05(3)R difference between DSm(10 K) in zero eld and in a 0.5 T eld could be explained either by ECO (with entropy 0.11R in the 2D case12), or by the B6% randomly-oriented orphan Dy spins on the Mg site (with entropy 0.06 Rln2). Neutron-scattering experiments on a powder sample of 162Dy3Mg2Sb3O14 distinguish these two scenarios by revealing the microscopic processes at play across T*.

Figure 2b shows magnetic neutron-scattering data at 0.5 K (above T*) and at the nominal base temperature of 0.03 K (below T*). At 0.5 K, our data show magnetic diffuse scattering only, with a broad peak centred at E0.65 1 that is characteristic of ice-rule correlations in structurally related pyrochlore magnets24. In contrast, at 0.03 K, strong magnetic diffuse scattering is observed in addition to magnetic Bragg peaks. These peaks develop at Tr0.35 K; that is, as T* is crossed. No additional peaks are observed on further cooling and the magnetic scattering does not change between 0.1 and 0.03 K. Between 0.03 and 50 K, the scattering is purely elastic within our maximum experimental resolution of E17 meV (Supplementary

Fig. 6), indicating that the spins uctuate on a timescale longer than B0.2 ns. Our 0.03 K data suggest two immediate conclusions. First, the magnetic Bragg peaks are described by the propagation vector k 0; that is, order preserves the

crystallographic unit cell below T*. Second, a large fraction of the magnetic scattering is diffuse; hence, correlated spin disorder persists below T* and involves the majority of spins. These results cannot be explained by only a small fraction of orphan spins, but are consistent with an ECO state13,14.

Average magnetic structure. We use reverse Monte Carlo (RMC) renement25,26 to t spin microstates to data collected between0.03 and 4 K. A single RMC microstate can capture both the

a b c

14

12

10

8

6

4

14

12

10

8

6

4

2

7

6

5

4

3

2

1

0

14

12

0 T

0.5 T

avg = 2.82(4) B

1 Dy1 )

2 mol Dy )

1

1 mol Dy )

1

1 Dy1 )

10

8

6

4

C m/T (J K

S m(J K

I T 50 K (barn sr

I 0.03 0.5 K (barn sr

4

4

2

0

2

0

2

2

0 0.1

1

10 0.5

0

1 1

Q (1)

1.5 1.5

2 2

2.5

T (K)

Q (1)

d

1 3

+ + + +

Figure 2 | Low-temperature magnetism of Dy3Mg2Sb3O14. (a) Magnetic heat capacity divided by temperature Cm/T (left axis; black points) and magnetic entropy change DSm(T) (right axis; orange curves). Zero-eld data and data measured in applied eld m0H 0.5 Tare shown (elds labelled on each curve).

Error bars represent the addition of statistical and systematic uncertainties, where statistical uncertainty is calculated from a least-squares t of the measured data to a two-timescale relaxation model, and systematic uncertainty is calculated assuming a 5% error on the sample mass. (b) Magnetic neutron-scattering data (black circles) at T 0.03 K and 0.5 K obtained by subtracting a high-temperature (50 K) measurement as background, ts from

reverse Monte Carlo (RMC) renements (red lines) and data t (blue lines). The 0.5 K curves are vertically shifted by 10 barn sr 1 Dy 1 for clarity. Error

bars on neutron-scattering data indicate one standard error propagated from neutron counts. (c) Magnetic Bragg scattering obtained as the difference between 0.03 and 0.5 K data (black circles), t from Rietveld renement (red line) and difference (blue line). The inset shows the model of the average magnetic structure obtained from Rietveld renement. (d) The vector average of the three microstates that are equally occupied in a ECO state yields an average all-in/all-out structure with ordered moment mavg m/3, consistent with experimental observations.

NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842

average spin structure responsible for Bragg scattering and the local deviations from the average responsible for diffuse scattering (Fig. 2b and Supplementary Fig. 7). We determine the average spin structure by two methods: rst, by averaging rened RMC microstates onto a single unit cell; second, by using a combination of symmetry analysis and Rietveld renement to model the magnetic Bragg prole (obtained as the difference between 0.03 and 0.5 K data) (Fig. 2c). Details of the Rietveld renements are given in Supplementary Note 1. Both approaches yield the same all-in/all-out average spin structure (inset to Fig. 2c and Supplementary Fig. 8). The ordered magnetic moment at 0.03 K, mavg

2.82(4) mB per Dy, is much less than the total moment of mE10 mB. These results are consistent with ECO: Fig. 2d shows that averaging over the three possible ECO microstates for a given triangle generates an all-in/all-out average structure, as observed experimentally; moreover, the expected ordered moment for ECO, m/3E3.3 mB per Dy13, is in general agreement with the measured value of 2.82(4) mB per Dy.

Evidence for emergent charge order. To look for signatures of ECO in real space, we compare the temperature evolution of mavg

with the percentage of T 3 charges (Fig. 3a). The latter

quantity, f3, takes a value of 25% for random spins, 100% for an all-in/all-out microstate, and 0% for a microstate that fully obeys the T 1 ice rule. The value of f3 extracted from RMC

renements decreases with lowering temperature to a minimum value of o5% below 1 K; these values represent upper bounds because RMC renements were initialized from random micro-states. Crucially, below T*, the T 1 rule is obeyed while mavg

is non-zero (Fig. 3a); this coexistence of ice-rule correlations with an all-in/all-out average structure is a dening feature of the ECO state13,14. We conrm ECO by calculating the charge-correlation function T 0

T rab

h i, the average product of charges separated

by radial distance rab on the honeycomb lattice formed by the triangle midpoints. At 0.5 K, this function decays with increasing rab, indicating that T 1 charges are disordered (Fig. 3b). At

0.03 K, T 0

T rab

h i shows two key features that indicate an

ECO state: a diverging correlation length, and an alternation in sign with a negative peak at the nearest-neighbour distance (Fig. 3c). The magnitude of T 0

T rab

h i found experimentally

(E0.6 (0.94 3mavg/m)2) is smaller than the value of unity

corresponding to an ideal ECO state, which indicates that the alternation of charges contains some errors; we show below this is probably due to the presence of site disorder.

Explanation of emergent charge order. Why does Dy3Mg2 Sb3O14 show fundamentally the same ECO as predicted for a 2D kagome system of in-plane Ising spins? This is far from obvious, because the real material differs from the existing model8 in three respects: (i) the spins are canted at an angle of 26(2) to the kagome planes, (ii) the planes are layered in 3D and (iii) there is Dy/Mg site disorder (Fig. 2c). This puzzle is elucidated by Monte Carlo simulations for a minimal model containing the nearest-neighbour exchange interaction J 3.72 K determined for

structurally-related Dy2Ti2O7 (refs 22,27), and the long-range magnetic dipolar interaction D 1.28 K calculated from experi

mentally determined DyDy distances. In 2D, spin canting interpolates between two limitsan ECO transition followed by lower-temperature spin ordering for in-plane spins8, and a single spin-ordering transition for spins perpendicular to kagome planes28and hence destabilizes ECO compared with the 2D in-plane limit. In contrast, the stacking of kagome planes stabilizes 3D ECOuniquely minimizing the effective Coulomb interaction between emergent chargesbut leaves the spin-ordering transition temperature essentially unchanged. The

effect of random site disorder is shown in Fig. 3d. Disorder broadens the specic-heat anomalies and suppresses the ECO transition temperature. In spite of this, we nd that a distinct ECO phase persists for 6% Mg on the Dy site; that is, the estimated level of disorder present in our sample of Dy3Mg2Sb3O14. Moreover, simulated magnetic specic-heat (Fig. 3d) and powder neutron-scattering (Supplementary Fig. 9) curves with B4 to 6% Mg on the Dy site show remarkably good agreement with experimental data, especially given that J is not optimized for Dy3Mg2Sb3O14 .

Implications of emergent charge order. An ECO microstate can be coarse-grained into a magnetization eld with two components: the all-in/all-out average spin structure with non-zero divergence, and the local uctuations from the average that are captured by (divergence-free) dimer congurations on the dual honeycomb lattice13. These two components are independent, which leads to descriptions of the ECO state in terms of spin fragmentation13,14. Without site disorder, the uctuating component yields pinch-point features in single-crystal diffuse-scattering patterns, the signature of a Coulomb phase13,29. Figure 3e shows that the introduction of site disorder blurs the pinch points and reduces the magnitude of the ordered moment in the ECO phase. We nd good overall agreement between patterns from model simulations with B4 to 6% Mg on the Dy site and from RMC microstates rened to powder data (Fig. 3e). These results suggest that pinch-point scattering could be observed in single-crystal samples of Dy3Mg2Sb3O14 with low levels of disorder. Our simulations also suggest why a transition from ECO to spin ordering is not observed experimentally: single-spin-ip dynamics (arguably more appropriate to real materials) become frozen in the ECO state and non-local (loop) dynamics are required to observe the spin-ordering transition in Monte Carlo simulations.

DiscussionThe ECO state in Dy3Mg2Sb3O14 is the rst realization of ordering of emergent degrees of freedom in a solid-state kagome material. Phase transitions driven by emergent excitations are rarerelated examples being the critical end-point in spin ice11,30,31 and the recent report of spin fragmentation in pyrochlore Nd2Zr2O7 (ref. 32). Moreover, the unusually slow spin dynamics offer the exciting possibility of measuring nite-time (Kibble-Zurek) scaling at the ECO critical point31. The ECO state in Dy3Mg2Sb3O14 presents an intriguing comparison with other partially ordered magnets. In Gd2Ti2O7, symmetry breaking yields two inequivalent Gd sites, only one of which orders33,34; in contrast, in the ECO state, all spins possess both ordered and disordered components. In Ho3Ga5O12, local antiferromagnetic correlations coexist with average antiferromagnetic order35, whereas in the ECO state, the average order is antiferromagnetic (all-in/all-out) while the local correlations are ferromagnetic (two-in/one-out or vice versa). Whether the predicted spin-ordering8 eventually occurs in Dy3Mg2Sb3O14 remains to be seen: spin freezing36,37 or site disorder may prevent its onset. We expect physical and/or chemical perturbations to control the properties of Dy3Mg2Sb3O14 ; for example, application of magnetic eld slightly tilted from the c-axis should drive a Kastelyn transition towards spin-ordering15,16; modied synthesis conditions may allow the degree of site mixing to be controlled20; and application of chemical pressure may alter the spin-canting angle and/or the distance between kagome layers, potentially generating a novel spin-ordering phase instead of ECO for sufciently large canting28. Substitution of Dy3 by other lanthanide ions20,3840 may increase the ratio of exchange to

4 NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842 ARTICLE

a

1

1

Ordered moment, avg ( Bper Dy)

C m/T (J K

4 25

20

15

10

5

0

3

2

b

0.5 K

0.03 K

% Triangles with 3 charge, 3

(0)(r)

0

c

1

0

40

30

20

10

0 0 0.1 0.2 0.3 0.4 0.5

RMC

0

0

0.1

1

10 10 20 30

T (K)

T *

r ()

d e

BraggMax

Max Diffuse

0

2 mol Dy )

1

0

0

T (K)

2

Figure 3 | Emergent charge order in Dy3Mg2Sb3O14. (a) Temperature evolution of the ordered magnetic moment per Dy, mavg (left axis) and the number of triangles for which T 3, f3 (right axis). Values of mavg from Rietveld renements to experimental data are shown as lled black circles, and values of

mavg from Monte Carlo (MC) simulations (with 4% Mg on the Dy site) are shown as a black dotted line. Upper bounds on f3 from reverse Monte Carlo (RMC) renements to experimental data are shown as hollow orange squares, and values from MC simulations as an orange dashed line. The location of

T* is shown by a vertical grey line, and the background is shaded blue below T* and green above T*. Throughout, error bars for results from Rietveld renements indicate one standard error from least-squares tting, and error bars from RMC renements are derived by assuming 10% uncertainty on the absolute intensity normalization of the magnetic scattering data. (b) Charge-correlation function T 0

T rab

h i obtained from RMC renements at 0.5 K,

and (c) T 0

T rab

h i from RMC at 0.03 K. Solid bars show correlation magnitudes, with positive correlations shown in red and negative correlations in

blue. (d) Magnetic heat capacity from MC simulations (system size N 7,776 spins) for different amounts of random site disorder (the % Mg on the Dy

site is labelled above each curve). The uncertainty in the MC results was assessed by computing the standard deviation of statistically-independent simulations; standard errors are smaller than the symbols in the gures. (e) Single-crystal neutron scattering calculations in the (hk0) plane from MC simulations at T 0.2 K for different amounts of random site disorder (the % Mg on the Dy site is labelled on each segment of the plot). The single-crystal

calculation from RMC renement to 0.03 K powder data (for 6% Mg on the Dy site) is shown for comparison. Separate colour scales are used for the intensity of the diffuse scattering and the {110} Bragg peaks, and the location of a pinch point is indicated by small white arrows.

dipolar interactions, offering promising routes towards exotic spin-liquid behaviour: dimensionality reduction by effective layer decoupling (when exchange dominates over dipolar interactions), and realization of quantum kagome systems with local spin anisotropies.

Methods

Sample preparation. Powder samples of Dy3Mg2Sb3O14 were prepared from a stoichiometric mixture of dysprosium (III) oxide (99.99%, Alfa Aesar*), magnesium oxide (99.998%, Alfa Aesar*) and antimony (V) oxide (99.998%, Alfa Aesar*). For neutron-scattering experiments a B5 g sample isotopically enriched with

162Dy (94.4(2)% 162Dy2O3, CK Isotopes*) was prepared. For all samples, starting materials were intimately mixed and pressed into pellets before heating at 1,350 C for 24 h in air. This heating step was repeated until the amount of impurity phases as determined by X-ray diffraction was no longer reduced on heating. The enriched sample contained impurity phases of MgSb2O6 (6.4(5) wt%) and Dy3SbO7(0.97(8) wt%), the latter of which orders antiferromagnetically at TE3 K (ref. 41). *The name of a commercial product or trade name does not imply endorsement

or recommendation by the National Institute of Standards and Technology (NIST).

X-ray diffraction measurements. Powder X-ray diffraction was carried out using a Panalytical Empyrean* diffractometer with Cu Ka radiation (l 1.5418 ).

Measurements were taken between 5r2yr120 with D2y 0.02.

*The name of a commercial product or trade name does not imply endorsement or recommendation by NIST.

Neutron-scattering measurements. Powder neutron diffraction measurements were carried out on the General Materials (GEM) diffractometer at the ISIS Neutron and Muon Source, Harwell, UK42, at T 0.50, 0.60, 0.90, 2.0, 4.0, 25 and

300 K. For T 25 and 300 K measurements, around 4.2 g of isotopically enriched

powder was loaded into a f 6 mm vanadium can and cooled in a ow cryostat.

For measurements at Tr25 K, the same sample was loaded into a f 6 mm

vanadium can, which was attached directly to a dilution refrigerator probe and loaded within a ow cryostat. Inelastic neutron-scattering experiments were carried out on the Disk Chopper Spectrometer (DCS) at the NIST Center for Neutron Research, Gaithersburg MD, USA43, at T 0.03, 0.10, 0.20, 0.30, 0.35, and 0.50 K.

Around 1.1 g of isotopically enriched powder was loaded into a f 4.7 mm copper

can and mounted at the base of a dilution refrigerator. The temperature was measured at the mixing chamber and does not necessarily reect the sample temperature for 0.1 and 0.03 K, as the spins progressively fall out of equilibrium. On DCS, data were measured with incident wavelengths of 1.8, 5 and 10 . The

NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842

1.8 data were used to look for crystal-eld excitations (Supplementary Fig. 4). The 10 data were used to look for low-energy quasi-elastic scattering (Supplementary Fig. 6). The 5 data were integrated over the energy range

0.15rEr0.15 meV to obtain the total scattering (Supplementary Fig. 10). Data reduction was performed using the MANTID and DAVE44 programs. All data were corrected for detector efciency using a vanadium standard, normalized to beam current (GEM) or incident beam monitor (DCS), and corrected for absorption by the sample.

Crystal-structure renements. Combined Rietveld analysis of the 300 K X-ray and neutron (GEM) diffraction data was carried out using the FULLPROF suite of programs45. The individual patterns were weighted so that the total contribution from X-ray and neutron diffraction was equal; that is, data from each of the ve detector banks on GEM was assigned 20% of the weighting of the single X-ray pattern. The neutron scattering cross-section for Dy was xed to bDy 0.6040 fm, to reect the isotopic composition as determined by inductively

coupled plasma mass spectrometry. Peak shapes were modelled using a pseudo-Voigt function, convoluted with an Ikeda-Carpenter function or an axial divergence asymmetry function for neutron and X-ray data, respectively. Backgrounds were tted using a Chebyshev polynomial function. At 25 K, Rietveld analysis of only the neutron diffraction data was carried out. In addition to the impurity phases observed in X-ray diffraction, a small amount (o1 wt%) of vanadium (IV) oxide from corrosion of the vanadium sample can was also observed in the neutron-diffraction data. The t to 300 K data is shown in Supplementary Fig. 1, rened values of structural parameters are given in Supplementary Table 1, and selected bond lengths are given in Supplementary Table 2.

Magnetic measurements. Magnetic susceptibility measurements, w(T) M(T)/H,

were made using a Quantum Design* Magnetic Properties Measurement System (MPMS) with a superconducting interference device (SQUID) magnetometer. Measurements were made after cooling in zero eld (ZFC) and in the measuring eld (FC) of m0H 0.1 T over the temperature range 2rTr300 K. Isothermal

magnetization M(H) measurements were made using a Quantum Design* Physical Properties Measurement System (PPMS) at selected temperatures 1.6rTr80 K between 14rm0Hr14 T. A global t to the M(H) data for TZ5 K (Fig. 1e) was

performed using the powder-averaged form for free Ising spins,

Mpowder

Ising

1. A random site-disorder model with 6% non-magnetic Mg on the Dy site was assumed, and Si 0 for atomic positions occupied by Mg. Ising variables were

initially assigned at random, and then rened against experimental data in order to minimize the sum of squared residuals,

w2W

XQIcalc Q Iexpt Q s Q

2

; 2

where I(Q) is the magnetic total-scattering intensity at Q, subscripts calc and expt denote calculated and experimental intensities, respectively, s(Q) is an experimental uncertainty, and W is an empirical weighting factor. For data collected on GEM, a rened at-in-Q background term was included inthe calculated I(Q). For data collected at Tr0.35 K, we obtain

Icalc(Q) I

Bragg(Q) I

diffuse(Q) I

random(Q), where subscripts Bragg, diffuse and random indicate magnetic Bragg, magnetic diffuse and high-temperature contributions, respectively. Here, Irandom(Q) 23C[mf(Q)/mB]2, where the constant

C (gnre/2)2 0.07265 barn and f(Q) is the Dy3 magnetic form factor54. The

Bragg and diffuse contributions were separated by applying the identity Si hSii DSi to each atomic position55, where the average spin direction hSii

is obtained by vector averaging the supercell onto a single unit cell, and the local spin uctuation DSi Si hSii. The Bragg contribution is given by

IBragg Q C

f Q

mB

m

2

Z

; 1

where H is applied magnetic eld, and magnetic moment m is the only tting parameter21. The tted value m 10.17(8) mB per Dy is in close agreement with the

expected value of 10.0 mB for a Kramers doublet ground state with g 4/3 and

mJ 15/2; in particular, the reduced value of the saturated magnetization,

MsatEm/2, is as expected for powder-averaged Ising spins21.

*The name of a commercial product or trade name does not imply endorsement or recommendation by NIST.

Heat-capacity measurements. Heat-capacity measurements were carried out on a Quantum Design* Physical Properties Measurement System instrument using dilution fridge (0.07rTr4 K) and standard (1.6rTr250 K) probes in a range of measuring elds, 0rm0Hr0.5 T. To ensure sample thermalization at low temperatures, measurements were made on pellets of Dy3Mg2Sb3O14 mixed with an equal mass of silver powder, the contribution of which was measured separately and subtracted to obtain Cp. The magnetic specic heat Cm was obtained by subtracting modelled lattice Cl and nuclear Cn contributions from Cp. We obtained Cl by tting an empirical Debye model to the 10oTo200 K data, with yD

1 tanh

mH cos y kBT

d cos y

22p2N

c

2

G2 R Q G

; 3

in which G is a reciprocal lattice vector with length G, V is the volume of the unit cell, Nc is number of unit cells in the supercell, R(Q G) is the resolution

function determined from Rietveld renement56. The magnetic structure factor F>(G)

Pi

NV

XGF? G

hSii> exp(iG ri), where supercript > indicates projection

perpendicular to G, and the sum runs over all atomic positions in the unit cell. The diffuse contribution is given by

Idiffuse Q C

f Q

mB

2

1 N

8 <

:

2

3

XiDSij j2

Xj 6 iAij sinQrijQrij Bijsin Qrij Qrij

3

!

" #

cos Qrij

Qrij

2

;

4

9 =

;

where sums run over all atomic positions in the supercell, rij is the radial

distance between positions i and j, and the correlation coefcientsAij DSi DSj (DSi rij)(DSj rij)/r2ij and Bij 3(DSi rij)(DSj rij)/r2ij DSi DSj

(refs 53,57). For data collected at TZ0.5 K, which show no magnetic Bragg scattering, we obtain Icalc(Q) I

diffuse(Q) I

272(13) K. To obtain a lower bound on the contact hyperne and electronic quadrupolar contributions to Cp23,46, we used previous experimental results on dysprosium gallium garnet47, a related material for which these contributions are known down to T 0.037 K. Correcting for the larger static electronic moment

E4.2 mB of dysprosium gallium garnet compared with hmiZ2.5 mB below 0.2 K for

Dy3Mg2Sb3O14, we obtained the high-temperature tail of the nuclear hyperne contributions as Cp A/T2 with A 0.0032 J K mol 1Dy (Supplementary Fig. 5).

*The name of a commercial product or trade name does not imply endorsement or recommendation by NIST.

Average magnetic structure analysis. The magnetic Bragg prole was obtained by subtracting data collected at To0.5 K from the 0.5 K data. Renements were carried out using the Rietveld method within the FULLPROF suite of programs45, as described above. For the magnetic-structure renement shown in Fig. 2c, candidate magnetic structures were determined using symmetry analysis48 via the SARAH49 and ISODISTORT50 programs, as described in Supplementary Note 1. The average magnetic structure is described by the irreducible representation G3, in

Kovalevs notation51. The basis vectors of the magnetic structure are given in Supplementary Table 3 and rened values of structural parameters are given in Supplementary Table 4.

Magnetic total scattering. To isolate the total magnetic contribution to the neutron-scattering data, data collected at a high temperature Thigh44yCW

were subtracted from the low-temperature data of interest, where Thigh 25 K

(GEM data) or 50 K (DCS data). For the data obtained below the magnetic ordering temperature of the Dy3SbO7 impurity phase (E3 K (ref. 41)), a rened model of the magnetic Bragg scattering of Dy3SbO7 was subtracted, as described in Supplementary Note 2 (we note that the orthorhombic crystal structure of Dy3SbO7 (ref. 52) allowed the impurity Bragg peaks to be readily distinguished from sample peaks). The t to neutron data of the Dy3SbO7 magnetic-structure model is shown in Supplementary Fig. 11, the magnetic basis vectors are givenin Supplementary Table 5, and rened values of structural parameters aregiven in Supplementary Table 6. The data were placed on an absolute intensity scale (barn sr 1 Dy 1) by normalization to the calculated nuclear Bragg prole at

Thigh.

Reverse Monte Carlo renements. Renements to the total (Bragg diffuse)

magnetic scattering were performed using a modied version of the SPINVERT program53 available from J.A.M.P. In these renements, a microstate was generated as a periodic supercell containing N 7776 Dy3 spin vectors Simsi^ei, where

m 10.0 mB is the xed magnetic moment length, the unit vector ^ei species the

local Ising axis determined from Rietveld renement, and the Ising variable si

random(Q), where Si replaces DSi everywhere. All renements employed the Metropolis algorithm with single-spin ip dynamics, and were performed for 200 proposed ips per spin, after which no signicant reduction in w2 was observed. Fits-to-data at T 0.03, 0.20, 0.50,

0.60, 0.90, 2.0 and 4.0 K are shown in Supplementary Fig. 7.

Monte Carlo simulations. Simulations were performed for the dipolar spin ice model27,58, extended to the geometry of interest in this work. The model is dened for Ising spins Simsi^ei, which are constrained to point along the

local easy-axis directions ^ei and can thus be described by the Ising pseudospin variables, si 1. The Hamiltonian comprises an exchange term of strength

J between nearest-neighbour spins hi, ji, and long-range dipolar interactions

of characteristic strength D (m0/4p)m2/r3nn between all pairs of spins, where

mE10 mB is the magnitude of the Dy3 spin and rnn is the nearest-neighbour distance of the lattice. The Hamiltonian is thus given by

H J

Xi;j h i

sisj ^ei ^ej

Dr3nn

Xi4j

sisj ^ei ^ej

r3ij

!

; 5

3 ^ei rij

r5ij

6 NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842 ARTICLE

where rij is the vector of length rij connecting spins i and j. We use D 1.28 K as

calculated from experimentally determined Dy-Dy distances, and J 3.72 K

from Dy2Ti2O7 (ref. 27), which has a similar Dy environment to Dy3Mg2Sb3O14 (ref. 22) (Supplementary Fig. 3). We treat the long-range dipolar interactions using Ewald summation58,59 with tinfoil boundary conditions at innity. In simulations including site disorder, non-magnetic ions are simulated by setting the corresponding si to zero. Our unit cell comprises three stacked kagome layers, each layer made from four kagome triangles. The whole system comprisesN 7776 spins in total, commensurate with the possible3

p

p spin-ordered state found in 2D (ref. 8). We use both single-spin ip and loop dynamics58,60, with Metropolis weights. Loop dynamics are necessary to ensure ergodicity at low temperatures and explore possible long-range spin-ordered states. We use the short loop algorithm58,60. One Monte Carlo sweep is dened as N single spin-ip attempts, followed by the proposal of loop moves until the cumulative number of proposed spin-ips (in the loops) is at least N. We use an annealing protocol, initializing the system at high temperature with B104N single spin-ip attempts, then decrease the temperature incrementally. After each temperature decrement, the system is updated with B103 Monte Carlo sweeps to ensure equilibration before collecting data every B10 Monte Carlo sweeps. Powder-averaged magnetic neutron-scattering patterns calculated from Monte Carlo are shown in Supplementary Fig. 9.

Data availability. The underlying research materials can be accessed at the following location: http://dx.doi.org/10.17863/CAM.4902

Web End =http://dx.doi.org/10.17863/CAM.4902 .

3

References

1. de Vries, M. A. et al. Scale-free antiferromagnetic uctuations in the s 1/2

kagome antiferromagnet herbertsmithite. Phys. Rev. Lett. 103, 237201 (2009).2. Han, T.-H. et al. Fractionalized excitations in the spin-liquid state of a kagome-lattice antiferromagnet. Nature 492, 406410 (2012).

3. Fk, B. et al. Kapellasite: A kagome quantum spin liquid with competing interactions. Phys. Rev. Lett. 109, 037208 (2012).

4. Mendels, P. & Bert, F. Quantum kagome frustrated antiferromagnets: One route to quantum spin liquids. C. R. Phys. 17, 455470 (2016).

5. Zhou, H. D. et al. Partial eld-induced magnetic order in the spin-liquid kagom Nd3Ga5SiO14. Phys. Rev. Lett. 99, 236401 (2007).

6. Zorko, A., Bert, F., Mendels, P., Marty, K. & Bordet, P. Ground state of the easy-axis rare-earth kagome langasite Pr3Ga5SiO14. Phys. Rev. Lett. 104, 057202 (2010).

7. Cairns, A. B. et al. Encoding complexity within supramolecular analogues of frustrated magnets. Nat. Chem. 8, 442447 (2016).

8. Chern, G.-W., Mellado, P. & Tchernyshyov, O. Two-stage ordering of spins in dipolar spin ice on the kagome lattice. Phys. Rev. Lett. 106, 207202 (2011).

9. Wills, A. S., Ballou, R. & Lacroix, C. Model of localized highly frustrated ferromagnetism: The kagom spin ice. Phys. Rev. B 66, 144407 (2002).

10. Bramwell, S. T. & Gingras, M. J. Spin ice state in frustrated magnetic pyrochlore materials. Science 294, 14951501 (2001).

11. Castelnovo, C., Moessner, R. & Sondhi, S. L. Magnetic monopoles in spin ice. Nature 451, 4245 (2008).

12. Mller, G. & Moessner, R. Magnetic multipole analysis of kagome and articial spin-ice dipolar arrays. Phys. Rev. B 80, 140409 (2009).

13. Brooks-Bartlett, M. E., Banks, S. T., Jaubert, L. D. C., Harman-Clarke, A. & Holdsworth, P. C. W. Magnetic-moment fragmentation and monopole crystallization. Phys. Rev. X 4, 011007 (2014).

14. Canals, B. et al. Fragmentation of magnetism in articial kagome dipolar spin ice. Nat. Commun. 7, 11446 (2016).

15. Fennell, T., Bramwell, S. T., McMorrow, D. F., Manuel, P. & Wildes, A. R. Pinch points and Kasteleyn transitions in kagome ice. Nat. Phys. 3, 566572 (2007).

16. Matsuhira, K., Hiroi, Z., Tayama, T., Takagi, S. & Sakakibara, T. A new macroscopically degenerate ground state in the spin ice compound Dy2Ti2O7 under a magnetic eld. J. Phys.: Condens. Matter 14, L559 (2002).

17. Mengotti, E. et al. Real-space observation of emergent magnetic monopoles and associated Dirac strings in articial kagome spin ice. Nat. Phys. 7, 6874 (2011).

18. Zhang, S. et al. Crystallites of magnetic charges in articial spin ice. Nature 500, 553557 (2013).

19. Anghinol, L. et al. Thermodynamic phase transitions in a frustrated magnetic metamaterial. Nat. Commun. 6, 8278 (2015).

20. Dun, Z. L. et al. Magnetic ground states of the rare-earth tripod kagome lattice Mg2RE3Sb3O14 (RE Gd,Dy,Er). Phys. Rev. Lett. 116, 157201 (2016).

21. Bramwell, S. T., Field, M. N., Harris, M. J. & Parkin, I. P. Bulk magnetization of the heavy rare earth titanate pyrochlores - a series of model frustrated magnets.J. Phys.: Condens. Matter 12, 483 (2000).22. Farmer, J. M. et al. Structural and crystal chemical properties of rare-earth titanate pyrochlores. J. Alloys Comp. 605, 6370 (2014).

23. Pomaranski, D. et al. Absence of Paulings residual entropy in thermally equilibrated Dy2Ti2O7. Nat. Phys. 9, 353356 (2013).

24. Hallas, A. M. et al. Statics and dynamics of the highly correlated spin ice Ho2Ge2O7. Phys. Rev. B 86, 134431 (2012).

25. McGreevy, R. L. & Pusztai, L. Reverse Monte Carlo simulation: A new technique for the determination of disordered structures. Mol. Simul. 1, 359367 (1988).

26. Paddison, J. A. M. & Goodwin, A. L. Empirical magnetic structure solution of frustrated spin systems. Phys. Rev. Lett. 108, 017204 (2012).

27. den Hertog, B. C. & Gingras, M. J. P. Dipolar interactions and origin of spin ice in Ising pyrochlore magnets. Phys. Rev. Lett. 84, 34303433 (2000).

28. Chioar, I. A., Rougemaille, N. & Canals, B. Ground-state candidate for the classical dipolar kagome Ising antiferromagnet. Phys. Rev. B 93, 214410

2016:

29. Fennell, T. et al. Magnetic coulomb phase in the spin ice Ho2Ti2O7. Science

326, 415417 (2009).

30. Sakakibara, T., Tayama, T., Hiroi, Z., Matsuhira, K. & Takagi, S. Observation of a liquid-gas-type transition in the pyrochlore spin ice compound Dy2Ti2O7 in a magnetic eld. Phys. Rev. Lett. 90, 207205 (2003).

31. Hamp, J., Chandran, A., Moessner, R. & Castelnovo, C. Emergent Coulombic criticality and Kibble-Zurek scaling in a topological magnet. Phys. Rev. B 92, 075142 (2015).

32. Petit, S. et al. Observation of magnetic fragmentation in spin ice. Nat. Phys. 12, 746750 (2016).

33. Champion, J. D. M. et al. Order in the heisenberg pyrochlore: the magnetic structure of Gd2Ti2O7. Phys. Rev. B 64, 140407 (2001).

34. Stewart, J. R., Ehlers, G., Wills, A. S., Bramwell, S. T. & Gardner, J. S. Phase transitions, partial disorder and multi-k structures in Gd2Ti2O7. J. Phys.:

Condens. Matter 16, L321 (2004).35. Zhou, H. D. et al. Intrinsic spin-disordered ground state of the Ising garnet Ho3Ga5O12. Phys. Rev. B 78, 140406 (2008).

36. Snyder, J. et al. Low-temperature spin freezing in the Dy2Ti2O7 spin ice. Phys. Rev. B 69, 064414 (2004).

37. Matsuhira, K. et al. Spin dynamics at very low temperature in spin ice Dy2Ti2O7. J. Phys. Soc. Jpn. 80, 123711 (2011).

38. Sanders, M. B., Krizan, J. W. & Cava, R. J. RE3Sb3Zn2O14 (RE La, Pr, Nd, Sm,

Eu, Gd): a new family of pyrochlore derivatives with rare earth ions on a 2D kagome lattice. J. Mater. Chem. C 4, 541550 (2016).39. Sanders, M. B., Baroudi, K. M., Krizan, J. W., Mukadam, O. A. & Cava, R. J. Synthesis, crystal structure, and magnetic properties of novel 2D kagome materials RE3Sb3Mg2O14 (RE La, Pr, Sm, Eu, Tb, Ho): comparison to

RE3Sb3Zn2O14 family. Phys. Stat. Solidi B 253, 20562065 (2016).40. Scheie, A. et al. Effective spin-1/2 scalar chiral order on kagome lattices in Nd3Sb3Mg2O14. Phys. Rev. B 93, 180407 (2016).

41. Fennell, T., Bramwell, S. T. & Green, M. A. Structural and magnetic characterization of Ho3SbO7 and Dy3SbO7. Can. J. Phys. 79, 14151419 (2001).

42. Hannon, A. C. Results on disordered materials from the GEneral Materials diffractometer, GEM, at ISIS. Nucl. Instr. Meth. Phys. Res. A 551, 88107 (2005).

43. Copley, J. & Cook, J. The Disk Chopper Spectrometer at NIST: a new instrument for quasielastic neutron scattering studies. Chem. Phys. 292, 477485 (2003).

44. Azuah, R. T. et al. DAVE: A comprehensive software suite for the reduction, visualization, and analysis of low energy neutron spectroscopic data. J. Res. Natl. Inst. Stan. Technol. 114, 341 (2009).

45. Rodrguez-Carvajal, J. Recent advances in magnetic structure determination by neutron powder diffraction. Physica B 192, 5569 (1993).

46. Henelius, P. et al. Refrustration and competing orders in the prototypical Dy2Ti2O7 spin ice material. Phys. Rev. B 93, 024402 (2016).

47. Filippi, J., Lasjaunias, J., Ravex, A., Tchou, F. & Rossat-Mignod, J. Specic heat of dysprosium gallium garnet between 37 mK and 2 K. Solid State Commun. 23, 613616 (1977).

48. Wills, A. S. Magnetic structures and their determination using group theory.J. Phys. IV France 11, 133158 (2001).49. Wills, A. S. A new protocol for the determination of magnetic structures using simulated annealing and representational analysis (SARAh). Physica B 276-278, 680681 (2000).

50. Campbell, B. J., Stokes, H. T., Tanner, D. E. & Hatch, D. M. ISODISPLACE: a web-based tool for exploring structural distortions. J. Appl. Crystallogr. 39, 607614 (2006).

51. Kovalev, O. V. Representations of the Crystallographic Space Groups (Gordon and Breach Science Publishers, 1993).

52. Siqueira, K. P. F. et al. Crystal structure of uorite-related Ln3SbO7(Ln La-Dy) ceramics studied by synchrotron X-ray diffraction and Raman

scattering. J. Solid State Chem. 203, 326332 (2013).53. Paddison, J. A. M., Stewart, J. R. & Goodwin, A. L. Spinvert: a program for renement of paramagnetic diffuse scattering data. J. Phys.: Condens. Matter 25, 454220 (2013).

54. Brown, P. J. International Tables for Crystallography. vol. C, 454460 (Kluwer Academic Publishers, 2004).

NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13842

55. Frey, F. Diffuse scattering from disordered crystals. Acta Crystallogr. B 51, 592603 (1995).

56. Mellergrd, A. & McGreevy, R. L. Reverse Monte Carlo modelling of neutron powder diffraction data. Acta Crystallogr. A 55, 783789 (1999).

57. Blech, I. A. & Averbach, B. L. Spin correlations in MnO. Physics 1, 3144 (1964).

58. Melko, R. G. & Gingras, M. J. P. Monte Carlo studies of the dipolar spin ice model. J. Phys.: Condens. Matter 16, R1277 (2004).

59. de Leeuw, S. W., Perram, J. W. & Smith, E. R. Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants. Proc. R. Soc. London A 373, 2756 (1980).

60. Melko, R. G., den Hertog, B. C. & Gingras, M. J. P. Long-range order at low temperatures in dipolar spin ice. Phys. Rev. Lett. 87, 067203 (2001).

Author contributions

H.S.O., P.M. and S.E.D prepared the samples. H.S.O., P.M., X.B., M.M. and S.E.D. performed and analysed the thermo-magnetic measurements. J.A.M.P., P.M., X.B., M.G.T., N.P.B. and S.E.D. performed the neutron-scattering measurements and J.A.M.P., M.M. and S.E.D. analysed the data. J.A.M.P. carried out the RMC renements. J.O.H. and C.C. carried out the Monte Carlo simulations. C.C. and S.E.D. conceived the project, which was supervised by C.C., M.M. and S.E.D. J.A.M.P. wrote the paper with input from all authors.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Paddison, J. A. M. et al. Emergent order in the kagome Ising magnet Dy3Mg2Sb3O14. Nat. Commun. 7, 13842 doi: 10.1038/ncomms13842 (2016).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

Acknowledgements

Work at Cambridge was supported through the Winton Programme for the Physics of Sustainability. The work of J.A.M.P., X.B. and M.M. and facilities at Georgia Tech were supported by the College of Sciences through M.M. start-up funds. J.A.M.P. gratefully acknowledges Churchill College, Cambridge for the provision of a Junior Research Fellowship. H.S.O. acknowledges a Teaching Scholarship (Overseas) from the Ministry of Education, Singapore. J.O.H. is grateful to the Engineering and Physical Sciences Research Council (EPSRC) for funding. C.C. was supported by EPSRC Grant No. EP/G049394/1, and the EPSRC NetworkPlus on Emergence and Physics far from Equilibrium. Experiments at the ISIS Pulsed Neutron and Muon Source were supported by a beamtime allocation from the Science and Technology Facilities Council. This work utilized facilities at the NIST Center for Neutron Research. Monte Carlo simulations were performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service (http://www.hpc.cam.ac.uk/

Web End =http://www.hpc.cam.ac.uk/) and the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk/

Web End =http://www.archer.ac.uk/, for which access was provided by an ARCHER Instant Access scheme). We thank G.-W. Chern, J. Goff,A. L. Goodwin, G. Lonzarich, G. Moller, D. Prabhakaran, J. R. Stewart and A. Zangwill for valuable discussions, and M. Kwasigroch for preliminary theoretical work.

8 NATURE COMMUNICATIONS | 7:13842 | DOI: 10.1038/ncomms13842 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications