Introduction

Geometrical frustration is a central theme of condensed-matter physics because it can generate exotic magnetic states. These states can typically be divided into spin liquids, in which frustration inhibits long-range magnetic order, and spin solids, in which perturbations to the dominant frustrated interactions drive magnetic order1. Defying this classification, some frustrated magnets exhibit partial magnetic order2,3,4,5,6,7,8—the coexistence of order and disorder in the same magnetic phase. Magnetic partial order can be driven by fluctuations in an “order-by-disorder” scenario9, by interactions between emergent degrees of freedom in spin-fragmented states10,11,12, or by proximity to a quantum critical point13, while structural partial order can drive the behavior of materials such as fast-ion conductors14,15, Pb-based photovoltaics16,17, and high-pressure elemental phases18. To benchmark theories of partially ordered states9,19,20, experimental determination of the nature of partial magnetic ordering in real materials is crucial.

Materials in which magnetic ions occupy a pyrochlore lattice of corner-sharing tetrahedra provide opportunities for realizing exotic frustrated states21. The frustrated pyrochlore antiferromagnet Gd2Ti2O7 is a canonical partially ordered system in which magnetic Gd3+ ions (S = 7/2, L = 0) undergo two phase transitions at T1 = 1.1 K and T2 = 0.75 K22,23,24,25,26. Both low-temperature (T ≪ T2) and intermediate (T2 < T < T1) phases have magnetic propagation vector \({{{\bf{k}}}}=(\frac{1}{2}\frac{1}{2}\frac{1}{2})\)27,28, which is uncommon among rare-earth pyrochlores and suggests exchange interactions extend beyond nearest neighbors20. Both magnetic structures are also partially ordered, as shown by the coexistence of magnetic Bragg and diffuse scattering in polarized-neutron scattering measurements28. The presence of partial ordering in a magnetic ground state at millikelvin temperatures is highly unusual21,28. However, the low-temperature magnetic structure of Gd2Ti2O7 has not yet been conclusively solved, for two reasons. First, the large neutron absorption cross-section of natural Gd makes neutron-scattering experiments on large crystals challenging. Second, most experimental probes are unable to distinguish a magnetic structure that orders with a single \({{{\bf{k}}}}=(\frac{1}{2}\frac{1}{2}\frac{1}{2})\) wavevector (1-k structure) from structures that superpose symmetry-equivalent \({{{\bf{k}}}}\in \langle \frac{1}{2}\frac{1}{2}\frac{1}{2}\rangle\) (2-k, 3-k, and 4-k structures). This phenomenon is known as the “multi-k problem”29. Developing methodologies to solve the multi-k problem is of broad relevance, because multi-k structures often correspond to noncollinear spin textures with nontrivial topological properties, such as skyrmions, hedgehogs, and vortex crystals (see, e.g., refs. 30,31,32,33,34).

Figure 1 shows the 1-k, 2-k, 3-k, and 4-k magnetic-structure candidates and their space groups. All are partially ordered, but each has a different modulation of the ordered magnetic moment μord: the 1-k and 4-k structures have 25% interstitial paramagnetic sites, whereas 2-k and 3-k structures have more complicated μord modulations. It was proposed in ref. 28 that magnetic diffuse-scattering measurements in the low-temperature phase support a 4-k structure with cubic magnetic symmetry. However, this result was called into question by the observation of transverse magnetization in small applied magnetic fields H∥〈112〉 and 〈100〉, which is inconsistent with cubic symmetry35,36. Theory proposes that the ordering is driven by a subtle interplay of energetic and entropic terms9. For interaction parameters relevant to Gd2Ti2O7, 1-k and 4-k structures are degenerate, and the 4-k structure is stabilized at T1 by thermal fluctuations in an “order-by-disorder” mechanism. The same model predicts a second phase transition at T2 into a 2-k ground state9. These striking predictions have awaited a conclusive experimental test.

Fig. 1: Candidate magnetic structures.
figure 1

Basic candidate structures for Gd2Ti2O7, showing a 1-k, b 2-k, c 3-k, and d 4-k structures. The symmetry of each structure is labeled. A single crystallographic unit cell is shown in each case; spin orientations are reversed in adjacent unit cells because of the \({{{\bf{k}}}}=(\frac{1}{2}\frac{1}{2}\frac{1}{2})\) propagation vector. Spin directions are shown as arrows with lengths proportional to the ordered moment magnitude μord, and the spin planes of triangular plaquettes for the 1-k and 4-k structures are shaded gray. Arrows of different colors indicate symmetry-inequivalent magnetic sites, and paramagnetic sites with zero μord in a, c, and d are shown as gray spheres.

Full size image

In this article, we experimentally determine the nature of partial magnetic order in Gd2Ti2O7 using neutron-scattering measurements of isotopically enriched powder and single-crystal samples, combined with symmetry analysis and spin-wave calculations. We show that the low-temperature state of Gd2Ti2O7 is actually 2-k, in agreement with theory9 but in contradiction with the interpretation of previous experiments28. Remarkably, inelastic neutron-scattering (INS) experiments on a powder sample unambiguously distinguish between n-k structures, demonstrating that neither single crystals nor symmetry-breaking perturbations such as applied field are prerequisite to solve the multi-k problem.

Our paper is structured as follows. We first present single-crystal neutron-diffraction measurements in applied magnetic field that suggest non-cubic magnetic symmetry. We then perform a comprehensive symmetry analysis of candidate magnetic structures. Finally, we show that only a partially ordered 2-k structure is consistent with low-temperature INS data. We conclude by discussing the general implications of our study for understanding partially ordered states and solving multi-k structures.

Results

Evidence for field-dependent magnetic domain population

An example of the multi-k problem29 is that measurements of magnetic Bragg intensities in zero applied magnetic field do not directly distinguish the structures shown in Fig. 1, due to spherical averaging in powder samples or averaging over degenerate magnetic domains in single crystals. To address this problem, we performed single-crystal neutron-diffraction measurements and applied a weak magnetic field \({{{\bf{H}}}}\parallel [1\bar{1}0]\) to break the domain degeneracy at T = 0.07 K after zero-field cooling. Domains of the cubic 4-k structure are related only by translational and time-reversal symmetries and hence appear identical to neutrons, whereas domains of other n-k structures are related by rotational symmetries and hence have different diffraction patterns. A field-induced domain imbalance is therefore expected to leave the diffraction pattern unchanged only if the low-temperature structure is 4-k.

The magnetic field dependence of selected magnetic Bragg intensities is shown in Fig. 2. Magnetic Bragg peaks in the (hhl) plane disappear in small applied field 0.2 ≤ μ0H ≤ 0.5 T, while magnetic Bragg peaks outside the (hhl) plane become more intense. These observations are incompatible with the cubic 4-k structure, unless the applied field actually causes a magnetic phase transition rather than a domain imbalance. This scenario occurs in Er2Ti2O737, but appears unlikely in Gd2Ti2O7, in which no experimental signature of such a phase transition is observed in either specific heat25 or torque magnetometry35 measurements at base temperature and for μ0H\(\parallel\) 〈110〉 of less than 2 T. The field-induced uniform magnetization is also too small to suppress antiferromagnetic Bragg peaks significantly (M = 0.2 μB for μ0H = 0.2 T26,36) and so need not be considered. Moreover, our conclusion is unaffected by the significant neutron absorption of our single-crystal sample, because it is indicated by the field dependence of the peak intensities, and not by comparing the intensities of different peaks at the same field. These observations provide evidence against the cubic 4-k structure; however, they do not distinguish between non-cubic structures such as 1-k and 2-k candidates.

Fig. 2: Field dependence of single-crystal magnetic Bragg scattering.
figure 2

Relative intensities of selected single-crystal magnetic Bragg peaks at different values of applied magnetic field \({{{\bf{H}}}}\parallel [1\bar{1}0]\) at T = 0.07 K. Magnetic Bragg peaks are labeled in each panel and applied fields μ0H are as labeled above the panels. To emphasize the field dependence, intensities are normalized such that the maximum intensity in each panel is unity. Error bars represent one standard deviation. For the \((\frac{3}{2}\frac{3}{2}\frac{7}{2})\) peak, the integrated intensity is 0.11(2) units for H = 0 and 0.22(2) units for μ0H = 0.5 T, consistent with the expected doubling; the estimated background is shown as a light blue line.

Full size image

Magnetic symmetry analysis

To constrain further the low-temperature magnetic structure, we reinterpret published powder neutron-diffraction data measured at 0.25 K (from ref. 28) using a comprehensive symmetry analysis. Figure 3 shows the magnetic diffraction data, with nuclear Bragg peaks subtracted, as a function of wavevector magnitude Q = ∣Q∣ = 2π/d. The possible magnetic irreducible representations for the \(Fd\bar{3}m\) space group (with origin choice 2) are denoted L1+, L2+, L3+, L1−, and L3− in Miller and Love’s notation38. Importantly, this analysis makes no assumptions about the nature of the phase transition at T2 (continuous or discontinuous), because we do not assume an intermediate-temperature parent structure, but instead consider the paramagnetic space group as the parent symmetry. Previous studies have shown that nearly all the features of the data can be modeled using a single irrep, L1+, which yields the fit shown in Fig. 3a27,28. The L1+ model is the best currently available and generates the four basic n-k structures shown in Fig. 1. Crucially, however, the \((\frac{1}{2}\frac{1}{2}\frac{1}{2})\) magnetic Bragg peak observed in the 0.25 K data is absent for the L1+ model (inset to Fig. 3a). This peak appears at T228, and is accompanied by a large specific-heat anomaly observed in several different samples24,25, suggesting that it is intrinsic and sample-independent. Hence, while the L1+ irrep is the main contributor to the low-temperature magnetic structure, at least one other irrep must also be present. We therefore investigated the four possible combinations of L1+ with one other irrep. Only the (L1+, L3−) and (L1+, L3+) irrep pairs allow nonzero intensity of the \((\frac{1}{2}\frac{1}{2}\frac{1}{2})\) peak and so are candidates. For each irrep pair, a given set of magnetic distortion-mode amplitudes39,40 yields several magnetic structures with identical Bragg profiles. Accordingly, we treated the magnetic distortion-mode amplitudes as free parameters that we optimized against our diffraction data. Figure 3 compares fits to diffraction data for the single-irrep L1+ model with the two-irrep (L1+, L3−) and (L1+, L3+) models. The two-irrep models yield nonzero intensity of the weak \((\frac{1}{2}\frac{1}{2}\frac{1}{2})\) peak—and slightly improved overall fits to the data—by increasing the number of free parameters as indicated in Fig. 3. Unfortunately, the inclusion of an additional irrep also increases the number of candidate structures from four to 32, all of which are consistent with the diffraction data. The magnetic space groups, mode amplitudes, and ordered-moment lengths for these structures are given in Supplementary Tables 1–3. To compress this search space, we apply a physical criterion that no site should have μord > 7.0 μB—the maximum value for S = 7/2 Gd3+ ions—which reduces the number of candidate structures to eight. Importantly, this criterion rules out all 3-k structures, so we do not consider these further. Four of the eight remaining structures are monoclinic (L1+, L3−) variants of the 1-k structure, in which paramagnetic spins order with a small μord28. However, these structures are disfavored by symmetry because the L3− irrep is not a symmetry-allowed secondary order parameter (SOP) of the 1-k structure39,40. The other four structures comprise three 2-k and one 4-k. Of these, the second irrep is a symmetry-allowed SOP in only one candidate—a (L1+, L3+) 2-k structure. This structure satisfies the moment-length criterion, unlike the single-irrep 2-k structure that was therefore discounted in previous work28,41. Hence, while the (L1+, L3+) 2-k structure has the same magnetic Bragg profile as other candidates, physical and symmetry arguments favor it.

Fig. 3: Powder neutron-diffraction data and magnetic Rietveld refinements.
figure 3

Experimental powder neutron-diffraction data (from ref. 28, λ = 2.42 Å) in the low-temperature phase (0.25 K) and Rietveld fits for a the L1+ irrep, b the (L1+, L3−) irrep pair, and c and the (L1+, L3+) irrep pair. Experimental data are shown as black points, fits as red lines, and data—fit as blue lines. For each model, the number of free parameters (magnetic distortion modes) n and the goodness-of-fit metric Rwp are shown. The insets show the \((\frac{1}{2}\frac{1}{2}\frac{1}{2})\) magnetic Bragg peak on an expanded scale.

Full size image

Inelastic neutron scattering and spin-wave analysis

While powder-averaged Bragg scattering cannot directly distinguish n-k structures, this limitation need not apply to the powder-averaged excitation spectrum. We therefore turn to powder INS experiments with high energy resolution (≈0.025 meV FWHM). Figures 4a, b show background-subtracted powder INS data in the intermediate phase (0.77 K) and the low-temperature phase (~0.05 K), respectively. The magnetic scattering at 0.77 K is broad in Q and E. By contrast, the low-temperature data show two relatively flat modes at energies of approximately 0.06 meV and 0.17 meV. Additional INS measurements on a thin piece of our single crystal show that the background-subtracted single-crystal scattering integrated over \((h,k,l)=(0\pm 1,0\pm 1,\frac{3}{2}\pm \frac{1}{2})\) qualitatively resembles the powder data (Fig. 4c).

Fig. 4: Inelastic neutron-scattering data and spin-wave models.
figure 4

a Inelastic neutron scattering data collected on a powder sample in the intermediate phase (0.77 K). b Powder INS data collected in the low-temperature phase (~0.05 K). c Single-crystal INS data integrated over \((h,k,l)=(0\pm 1,0\pm 1,\frac{3}{2}\pm \frac{1}{2})\) in the low-temperature phase (0.05 K, solid red diamonds) and intermediate phase (0.85 K, empty gray circles). Error bars represent one standard deviation. d Optimized 2-k magnetic structure. e Powder-averaged linear spin-wave theory (LSWT) calculation for the L1+ 4-k structure. f Powder LSWT calculation for the optimized (L1+, L3+) 2-k structure discussed in the text. g Powder LSWT calculation for the L1+ 1-k structure. h Powder INS data and spin-wave calculations at 0.05 K. Experimental data are given in absolute intensity units by normalization to the nuclear Bragg profile. Gray crosses show the measured elastic line (∣E∣ < 0.03 meV) with the incoherent scattering cross-section subtracted, empty blue diamonds show the measured low-E mode (0.03 ≤ E < 0.12 meV), and filled red circles show the measured high-E mode (0.12 ≤ E ≤ 0.21 meV). Spin-wave calculations for low-E and high-E modes are shown as blue and red lines, respectively, and have been vertically scaled by the same factor to match the experimental data.

Full size image

We use linear spin-wave theory to test candidate magnetic structures against the low-temperature excitation spectrum. The minimal spin Hamiltonian for Gd2Ti2O7 is given by

$$H={J}_{1}\mathop{\sum}\limits_{\langle i,j\rangle }{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}+{J}_{2}\mathop{\sum}\limits_{\langle \langle i,j\rangle \rangle }{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}+D\mathop{\sum}\limits_{i}{({S}_{i}^{z})}^{2}+{H}_{{{{\rm{dip}}}}},$$
(1)

where J1 and J2 denote Heisenberg exchange interactions between nearest neighbor and next-nearest neighbor spin pairs, which are denoted by angle brackets 〈 〉 and 〈〈 〉〉, respectively; D is a single-ion anisotropy term that arises from mixture of the excited 6P7/2 atomic state into the 8S7/2 ground state42; and

$${H}_{{{{\rm{dip}}}}}={D}_{{{{\rm{dip}}}}}{r}_{{{{\rm{nn}}}}}^{3}\mathop{\sum}\limits_{i,j}\frac{{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}-3({{{{\bf{S}}}}}_{i}\cdot {\hat{{{{\bf{r}}}}}}_{ij})({{{{\bf{S}}}}}_{j}\cdot {\hat{{{{\bf{r}}}}}}_{ij})}{{r}_{ij}^{3}}$$
(2)

is the long-ranged magnetic dipolar interaction including contributions from all spin pairs i, j, whose energy scale DdipS(S + 1) = 0.84 K is fixed by S and the nearest-neighbor distance rnn = 3.60 Å22. We assume the literature value of J1S(S + 1) = 4.8 K22 and include a small ferromagnetic J2 = − 0.04J1 to stabilize \({{{\bf{k}}}}=(\frac{1}{2}\frac{1}{2}\frac{1}{2})\) ordering9,20,22; however, the spin-wave spectrum is not strongly affected by choosing different small values of J2. Electron-spin resonance (ESR) experiments find an easy-plane anisotropy that favors spin alignment perpendicular to local 〈111〉 axes42; we take DS2 = 1.5 K to match our INS data optimally. As an independent check on these interaction parameters, we performed measurements of the paramagnetic (T > T1) diffuse scattering using the D4 diffractometer at the ILL43,44, and find these are consistent with Monte Carlo simulations45 of the paramagnetic diffuse scattering for our parameter values (Supplementary Fig. 1).

A prerequisite for spin-wave modeling is that the magnetic structure is a local energy minimum of Eq. (1). We therefore tested which of the eight candidate structures are proximate to energy minima by iteratively aligning each spin with its mean field and checking for stability via the absence of imaginary spin-wave modes. Two candidate structures—one 2-k and one 4-k—are locally stable; both derive from the (L1+, L3+) irrep pair that yields the best fit to diffraction data (Fig. 3c). By contrast, all candidate 1-k structures with nonzero interstitial μord are unstable. Figure 4d shows the optimized (L1+, L3+) 2-k structure, which resembles the refined structure (Supplementary Fig. 2 and Supplementary Table 4). Compared to the L1+ 2-k structure shown in Fig. 1b, it has canted magnetic moments with more uniform magnitudes, μord/μB ∈ {6.1(1), 4.6(1), 6.2(1)} in a 1:1:2 ratio. However, the suppression of μord compared to its theoretical value of 7.0 μB indicates partial ordering. Figure 4e–g shows the calculated spin-wave spectra for the L1+ 4-k, (L1+, L3+) 2-k, and L1+ 1-k structures, respectively. The (L1+, L3+) 4-k structure involves weak ordering of the 25% of spins that are paramagnetic in the L1+ 4-k structure; however, the spectra for both these structures are very similar and we show only the latter. The 1-k and 4-k calculations strongly disagree with the 0.05 K experimental data. By contrast, the 2-k calculation reproduces well the experimental data, most importantly the prominent low-energy (~0.06 meV) mode that is absent for the other candidates. This key qualitative observation is not affected by different choices of interaction parameters in the physically appropriate regime J > D ≫ J2. The Q dependences of the low-E (0.03 ≤ E < 0.12 meV) and high-E (0.12 ≤ E ≤ 0.21 meV) modes shown in Fig. 4h confirm agreement between the 2-k calculation and the experimental data. We therefore conclude that the (L1+, L3+) 2-k structure is the correct low-temperature model.

Investigation of spin-lattice coupling

The 2-k structure has orthorhombic symmetry (magnetic space group Camcm), which is therefore expected to drive a crystallographic distortion via spin-lattice coupling. To test this hypothesis, we performed high-resolution powder neutron-diffraction measurements with ΔQ/Q = 5 × 10−5 on the HRPD instrument at ISIS46,47. However, our data do not show visible peak splitting [Supplementary Fig. 2], and while a statistically significant rhombohedral distortion could be refined, orthorhombic refinements were inconclusive due to their increased number of parameters. This suggests that spin-lattice coupling in Gd2Ti2O7 is too weak to yield an observable orthorhombic distortion in our measurements. Further experiments, such as NMR or Mössbauer spectroscopy, may allow the expected distortion to be observed, but the requirement for <0.7 K temperatures presents experimental challenges.

Discussion

Our experimental result that the low-temperature structure of Gd2Ti2O7 is 2-k confirms state-of-the-art theoretical predictions9. However, it contradicts a previous experimental study28, which proposed a 4-k low-temperature structure based on analysis of low-temperature magnetic diffuse scattering. This study did not consider 2-k structures, because the single-irrep 2-k structure is unphysical and two-irrep 2-k structures were not identified. It also assumed that ordered sites contribute no diffuse scattering. This assumption is incorrect, however, because spin-wave scattering from ordered sites contributes to the energy-integrated diffuse intensity. In general, the magnetic Bragg intensity is proportional to \({({\mu }_{{{{\rm{ord}}}}}/{\mu }_{{{{\rm{B}}}}})}^{2}\) and the total magnetic intensity is proportional to g2S(S + 1) for a spin-only ion. For a fully ordered system, μord/μB = gS, so that spin-wave scattering would comprise 1/(S + 1) = 22% of the total intensity for S = 7/2. By contrast, partially ordered Gd2Ti2O7 has μord/μB < gS. Combining the refined values of μord with the total-moment sum rule, we estimate that low-temperature diffuse scattering comprises 46(2)% of the total intensity. The inelastic scattering (E > 0.03 meV) shows clear Q dependence that is fully reproduced by the spin-wave calculation, whereas the elastic scattering (∣E∣ < 0.03 meV) is weak and essentially flat away from Bragg peaks. Therefore, the low-temperature magnetic diffuse scattering arises primarily from enhanced spin-wave fluctuations, and not from static spin disorder or quasielastic (paramagnetic) spin fluctuations.

Our neutron-scattering experiments, magnetic symmetry analysis, and spin-wave modeling reveal that the low-temperature magnetic structure of Gd2Ti2O7 is 2-k, solving a longstanding problem in the field of frustrated pyrochlore oxides. The 2-k structure does not contain any entirely paramagnetic sites, and the increase in uniformity of the ordered-moment magnitudes likely stabilizes it over higher-symmetry 1-k and 4-k candidates. Nevertheless, the 2-k structure remains partially ordered, because magnetic-moment magnitudes are modulated and suppressed by up to 35% compared to their fully ordered value, while spin-wave fluctuations are enhanced. The magnetic ground state thus combines partial ordering with a multi-k structure—an exotic combination that, to the best of our knowledge, has not been observed in any other material with localized magnetic moments. Its existence in Gd2Ti2O7 is especially remarkable because the large value of the spin quantum number and absence of orbital angular momentum for Gd3+ ions suggests that quantum effects are weak. Instead, the enhancement of spin fluctuations may be related to the near-degeneracy of 2-k and 4-k states that results from the interplay of frustrated exchange interactions with magnetic anisotropy and long-ranged dipolar interactions.

While our study has focused on the low-temperature ground state, our results also shed light on the intermediate phase. Theory predicts that a 4-k structure is stabilized by thermal fluctuations at T1 for our interaction parameters, before transitioning to a 2-k ground state at T29. This prediction is consistent with the experimental observation that the \((\frac{1}{2}\frac{1}{2}\frac{1}{2})\) magnetic Bragg peak is absent in the intermediate phase28, because the irreps required to generate this peak are not symmetry-allowed SOPs for the 4-k structure. Moreover, our INS data in the intermediate phase show only broad inelastic features (Fig. 4a), consistent with the observation of quasielastic scattering in neutron-spin-echo experiments in the intermediate phase48. A reduction in spin-wave scattering intensity is expected due to elevated temperature (T = 0.75T1), which suppresses the refined value of the ordered magnetic moment to 2.65(3) μB per ordered Gd at 0.77 K, for a 4-k model [Supplementary Fig. 4]. However, the near-complete suppression of spin-wave intensity of the intermediate phase may hint at the presence of purely paramagnetic sites in a 4-k structure. Such sites would contribute a broad continuum of magnetic scattering, as observed in our INS data (Fig. 4a), and may also suppress propagating spin-wave excitations, similar to the effect of paramagnetic impurities49. Modeling these effects in the intermediate phase would provide a possible avenue for future work.

Perhaps the most general implication of our study is that single-k and multi-k structures can be unambiguously distinguished using inelastic neutron-scattering experiments on powder samples in zero applied field. In contrast to previous approaches50,51,52, this approach requires neither single crystal samples nor the application of external fields that explicitly break the symmetry. We anticipate that this approach will be widely applicable, and will help identify and understand multi-k structures where noncoplanar spin textures have nontrivial topologies30,31,32,33,34.

Methods

Single-crystal sample growth

A single crystal was prepared by the floating-zone image furnace method53,54 and was 99.4% enriched with 160Gd to minimize neutron absorption by 155Gd and 157Gd.

Neutron diffraction experiments

Single-crystal neutron-diffraction measurements were performed using the WISH time-of-flight diffractometer at ISIS55 on a = 10 mm3 piece cut from our single crystal. This sample was aligned with the (hhl) scattering plane horizontal, and loaded within a dilution refrigerator and magnet, with the magnetic field along the vertical \([1\bar{1}0]\) direction. Details of the previously published powder neutron-diffraction measurements are given in ref. 28.

Inelastic neutron-scattering experiments

Powder inelastic neutron-scattering experiment measurements were performed using the DCS spectrometer at NIST56 with an incident wavelength of 8.0 Å. An approximately 0.2 g portion of the same isotopically enriched powder sample studied in ref. 28 was loaded in a Cu foil and wrapped around the circumference of a cylindrical Cu container, to minimize neutron absorption. Single-crystal inelastic neutron-scattering experiments measurements were performed on a thin piece of our crystal using the CNCS instrument at ORNL with an incident wavelength of 9.0 Å. The samples were cooled using dilution refrigerators.

Magnetic structure refinement

Magnetic refinements to the powder neutron-diffraction data shown in Fig. 3 were performed using the Topas Academic software57. Lattice and profile parameters (peak shape, zero offset, and background) were obtained by Pawley refinement to the experimental magnetic Bragg profile, and the intensity scale factor was obtained from refinement to the nuclear Bragg profile. Refinement of an overall Debye–Waller parameter yielded a positive value, indicating that absorption is not a major problem for the powder sample. For each irrep or pair of irreps, the magnetic distortion-mode amplitudes were then refined to the experimental magnetic Bragg profile, keeping lattice, scale, and profile parameters fixed. To identify different sets of mode amplitudes with identical fit quality (degenerate solutions), multiple refinements were performed with different initial values of the mode amplitudes. For the (L1+, L3−) irrep pair, we identified two degenerate sets of mode amplitudes, A and B, given in Supplementary Table 1. For the (L1+, L3+) irrep pair, we identified the four degenerate sets of mode amplitudes C–F given in Supplementary Table 2. By combining these sets of mode amplitudes with the n-k structures using the Isodistort software39,40, we identified 32 candidate magnetic structures. The magnetic space groups, mode amplitudes, and ordered-moment magnitudes for these structures are given in Supplementary Table 3.

Spin-wave calculations

Spin-wave calculations shown in Fig. 4e–h were performed using the SpinW software package58. For each candidate magnetic structure, we generated a magnetic unit cell containing 2 × 2 × 2 crystallographic unit cells (128 ordered magnetic moments); this “supercell” approach is necessary to calculate spin-wave spectra of multi-k structures. The candidate structures were then optimized by iteratively aligning the orientation of each ordered magnetic moment with its mean field. The magnetic dipolar interaction was evaluated by direct summation to a maximum interatomic distance of 40 Å. We checked that the optimized magnetic structure and excitation spectrum remained stable to summations over longer distances (up to 100 Å), indicating that satisfactory convergence was achieved. The optimized and refined (L1+, L3+) 2-k structures are given in Supplementary Table 4, and their calculated magnetic neutron-diffraction patterns match closely as shown in Supplementary Fig. 3. The calculated spin-wave spectra were powder averaged numerically using Fibonnaci spherical integration with 987 points on the sphere, multiplied by the squared magnetic form factor of Gd3+, and convolved with a Gaussian with FWHM of 0.025 meV to match the energy resolution of the experimental data.