Emergent order in the kagome Ising magnet Dy3Mg2Sb3O14

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 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 ∼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 Dy3Mg2Sb3O14 as a tuneable system to study interacting emergent charges arising from kagome Ising frustration.

Momentum and energy dependence of neutron-scattering intensity at nominal T = 0.03 K, measured using the DCS instrument at the National Institute of Standards and Technology (NIST) [2] with incident neutron wavelength λ = 1.8Å. An empty can measurement has been subtracted from the data in order to remove background scattering. These measurements reveal no crystal electric-field excitations to a maximum energy transfer 23 meV (270 K), demonstrating that the ground-state Kramers doublet is well separated from excited crystal-field states. showing the following contributions: total measured (C p , full blue circles), modelled lattice (C l , green dashed line), modelled nuclear (C n , solid orange line), and extracted magnetic (C m = C p − C n − C l , open blue circles). The main panel shows C/T and the inset shows C. The temperature region where sample coupling falls below 90% is shaded grey. The magnetic specific heat C m displays a very small peak at T i = 3.35(5) K, which is probably associated with the magnetic ordering of the ∼ 2 wt% Dy 3 SbO 7 impurity phase [3]. There is also a small peak in C m around 0.1 K, which is of uncertain origin given the poor sample coupling below approximately 0.2 K. Error bars represent the addition of statistical and systematic uncertainties, where statistical uncertainty is calculated from a least-squares fit of the measured data to a two-timescale relaxation model, and systematic uncertainty is calculated assuming a 5% error on the sample mass.
Supplementary Figure 6: Low-energy inelastic scattering for Dy 3 Mg 2 Sb 3 O 14 . Energy dependence of neutron-scattering intensity at 0.03 K (black circles) and 50 K (red squares), measured with incident neutron wavelength λ = 10Å using the DCS spectrometer at NIST [2]. The data are integrated over the range 15 ≤ 2θ ≤ 105 • , which contains no nuclear Bragg peaks but intense magnetic diffuse scattering. The black dotted vertical lines show the calculated energy resolution (FWHM) of the instrument, 16.9 µeV. The data are well described by Gaussian fits (solid lines), which yield fitted values of the FWHM of 13.1(2) µeV at 0.03 K and 13.5(4) µeV at 50 K; hence the energy line-width due to dynamical spin fluctuations is limited by the instrumental resolution at both temperatures. Error bars show one standard error propagated from neutron counts. Calculations are performed using the same approach as the RMC refinements (described in the Methods section), and a total magnetic moment µ = 10.0 µ B per Dy is assumed. No fitting parameters are included in order to match the experimental data. Good overall agreement is achieved between 0.2 K simulations and experimental data for between ∼4 and 6% Mg on the Dy site; values close to 4% yield best agreement with Bragg intensities, whereas values close to 6% yield best agreement with the diffuse-scattering profile. We anticipate that the agreement may be further improved by fitting the value of the nearest-neighbour exchange interaction and/or by considering possible short-range correlations of Mg occupancy of the Dy site.

Supplementary
Supplementary Figure 10: Validity of quasi-static approximation. Energy dependence of neutron-scattering intensity at 0.03 K (black circles) and 50 K (red squares), measured with incident neutron wavelength λ = 5Å using the DCS spectrometer at NIST [2]. The data are integrated over the range 15 ≤ 2θ ≤ 45 • , which contains no nuclear Bragg peaks but intense magnetic diffuse scattering. The black dotted vertical lines show the calculated energy resolution (FWHM) of the instrument, and the grey box shows the range of energy integration for the neutron-scattering data shown in Fig. 2b. The energy line-width due to dynamical spin fluctuations is limited by the instrumental resolution at both 0.03 and 50.0 K, which demonstrates that the quasi-static approximation is ideally satisfied. The decrease in peak intensity at 50 K occurs because diffuse intensity is redistributed to higher 2θ. Standard errors (propagated from neutron counts) are smaller than the symbol size in the plots.  determined from combined analysis of 300 K X-ray and neutron powder diffraction data, and 25 K neutron powder diffraction data. Fixed parameters are denoted by an asterisk ( * ).
Supplementary   Mg2 3b, (0, 0, 1 2 )    Fig. 2b]. To isolate this magnetic Bragg scattering, we subtracted the 0.5 K data from the 0.03, 0.10, 0.20, 0.30 and 0.35 K data. All magnetic Bragg peaks are indexed by the propagation vector k = (0, 0, 0) with respect to the hexagonal unit cell of Dy 3 Mg 2 Sb 3 O 14 . Symmetry-allowed magnetic structures were determined using the program SARAh [4] and verified using ISODISTORT [5,6]. The analysis first determines the group of symmetry elements that leave k invariant, and then decomposes the magnetic representations of Dy1 and Dy2 sites into irreducible representations. The decomposition of the magnetic representation for the Dy1 site is and the decomposition for the Mg2 site is where different irreps are labelled by subscript numbers (using the notation of Kovalev [7]), the dimensionality of each irrep is given by the superscript number, and the product of the superscript and the coefficient yields the number of basis vectors associated with the irrep. Rietveld refinements were performed to the 0.03 − 0.50 K data to determine the average magnetic structure. In these refinements, the occupancy of the Dy site was fixed at 94%, as estimated from the crystal-structure refinements [Supplementary Figure 1 and Supplementary Tables  1 and 2]. The magnetic peak-shape parameters were modelled as Gaussian and fixed to equal the nuclear peakshape parameters refined to 50 K data. The intensity scale factor was fixed at the value obtained from Rietveld refinement to the 50 K data, and the background was fitted by Chebychev polynomials. By testing each of the irreps Γ 1 , Γ 3 , and Γ 5 against the data, we found that only the Γ 3 irrep allowed a good fit. The orientation of the ordered moment of atom i is given by µ i,avg = ν C ν m ν,i ; here, the C v denotes the refined coefficient of the basis vector m ν,i = m ν,i aâ + m ν,i bb + m ν,i cĉ , whereâ,b,ĉ are unit vectors parallel to the unit-cell axes. The basis vectors for this magnetic structure are shown in Supplementary Table 3. The component of the ordered magnetic moment parallel to the c-axis on the Dy1 site, µ c,avg /µ avg = 0.44(3), shows that Dy spins are canted slightly more towards the c-axis than for cubic spin ices such as Dy 2 Ti 2 O 7 , where the equivalent projection equals 1/3.
Because Γ 3 occurs in the magnetic representation for both Dy1 and Mg2 sites, ordered moments may form on both sites in a single phase transition. The basis vectors of the Dy1 site describes an "all-in/all-out" structure, while basis vectors of the Mg2 site describe a uniform ferromagnetic component along the c-axis. Mindful of the fact that weak ferromagnetic components can be poorly determined from powder-diffraction data, we performed two separate refinements to the 0.03 − 0.50 K data. First, we refine the basis-vector coefficients on both sites ("2-site refinement"); second, we constrain the magnetic moment on the Mg2 site to equal zero ("1-site refinement"). Supplementary Table 4 shows the results from each refinement. Our data are consistent with a small ordered moment (∼ 1µ B ) on the Mg2 site; however, refining this moment yields an insignificant improvement in the fit, and does not change parameter values associated with the Dy1 site. In subsequent refinements, we therefore fix µ avg ≡ 0 on the Mg2 site for the sake of simplicity. The fits obtained to 0.10, 0.20, 0.30 and 0.35 K data are shown in Supplementary Figure 8. For these refinements, we fixed the magnetic structure to the result from the refinement to 0.03 − 0.50 K data, and refined only the background parameters and the overall scale factor in order to determine the temperature dependence of µ avg [Fig. 3a].
Supplementary Note 2: Average magnetic structure of Dy 3 SbO 7 Magnetic Bragg peaks appear in our neutron-diffraction data at T ≤ 2 K, but are absent in T ≥ 4 K data. These peaks cannot be indexed by high-symmetry propagation vectors of Dy 3 Mg 2 Sb 3 O 14 , but are indexed by the propagation vector k = (0, 0, 0) for the Dy 3 SbO 7 impurity phase [8]. We therefore identify these peaks with magnetic ordering of Dy 3 SbO 7 , which is reported to occur at 3.0(3) K [3].
A magnetic-structure model for Dy 3 SbO 7 has not been reported previously, so we turn to symmetry analysis to identify possible magnetic structures. The crystal structure of Dy 3 SbO 7 (space group Cmcm [8]) contains two inequivalent Dy sites, 4a (Dy1) and 8g (Dy2). The symmetry-allowed magnetic structures were determined by representational analysis using the program SARAh [4]. (irreps, Γ) of this group. The decomposition of the magnetic representation for the Dy1 site is and the decomposition for the Dy2 site is Because a single magnetic transition is reported in Dy 3 SbO 7 [3], we consider only the irreps that appear in the decomposition for both Dy1 and Dy2 sites-namely, Γ 1 , Γ 3 , Γ 5 , and Γ 7 . We used Rietveld refinement to test each of these irreps against the 0.60 − 4.0 K data collected on the GEM diffractometer at ISIS [9]. In these refinements, the magnetic peak-shape parameters were modelled as convoluted pseudo-Voigt and Ikeda-Carpenter functions, and the background was fitted by a linear interpolation between between fixed points (for which the background level was subsequently refined at lower Q). The Γ 3 irrep provided the best fit-todata, shown in Supplementary Figure 11a. This magnetic structure is shown in Supplementary Figure 11b and describes a non-collinear ferrimagnet. Because of the uncertainty associated with the wt% Dy 3 SbO 7 in our sample, it was not possible to determine accurately the values of the ordered magnetic moment on Dy1 and Dy2 sites; however, the data are consistent with moment lengths on Dy1 and Dy2 sites that are equal to within ∼20%. The basis vectors are given in Supplementary Table 5, and refined values of magnetic-structure parameters are given in Supplementary Table 6. Having determined a model of the magnetic structure of the Dy 3 SbO 7 impurity phase, we calculated the magnetic scattering from Dy 3 Mg 2 Sb 3 O 14 as I = I meas −I T high −I imp,Bragg +I imp,pm , where I T high denotes a hightemperature (25 K or 50 K) measurement which is subtracted to remove non-magnetic scattering, I imp,Bragg denotes the calculated magnetic Bragg profile for the Dy 3 SbO 7 impurity, and I imp,pm denotes the calculated paramagnetic intensity for the Dy 3 SbO 7 impurity.