Introduction

The magnetic skyrmion lattice (SkL) is a periodic array of spin vortices, which may be considered as an assembly of individual, tubular skyrmion particles protected against decay by their topological winding number1,2,3. Skyrmions hold significant technological promise, for example as tiny bits for data storage devices4, which are highly controllable by small applied electrical currents5. From a fundamental perspective, we may classify nearly all previously reported real-world realizations of the SkL state into two families: (1) thin ferromagnetic slabs, where dipolar interactions enable formation of topological bubbles with characteristic size of λ = 100 nm–10 μm6. (2) Magnets with broken inversion symmetry (i.e. chiral or polar structures)1,2,3,7, where competing Heisenberg and Dzyaloshinskii–Moriya interactions favor twisted spin structures with λ~10–200 nm. Magnetic interfaces, inversion-breaking by definition, are included in this second category8,9.

In the search for even smaller skyrmions (λ = 1–10 nm), which are expected to show giant responses to optical, electrical, and magnetic stimuli4, attention has turned to systems with competing exchange interactions, or generalized Ruderman–Kittel–Kasuya–Yosida (RKKY) interactions10,11,12,13,14,15. Noteworthy in this context is the paradigm shift concerning the role of spin–orbit coupling: centrosymmetric materials, where Dzyaloshinskii–Moriya interactions are absent or cancel out globally, have now moved into the cross-hairs of the search for skyrmion host compounds. They offer a path not only towards the miniaturization of spin textures, but also towards new physical properties related to the energetic near-degeneracy of left- and right-handed screws as well as Néel and Bloch-type spin helicities10,11,16. Amongst these predicted properties are time-dependent helicity changes in response to an applied current11, a modified excitation spectrum11, mixed phases and near-degeneracy of skyrmions and antiskyrmions10, background-free generation of second harmonic light due to a large toroidal moment17, and qualitatively new decay mechanisms18.

The recently discovered SkL phase with giant topological Hall effect (THE) in centrosymmetric Gd2PdSi3, with a triangular lattice of classical Gd3+ spins, represents a seminal step into the direction of realizing the theoretical predictions listed above19. However, Gd2PdSi3 suffers from unavoidable disorder in the (Pd,Si) sublattice and a related crystallographic superstructure, rendering it challenging to perform advanced real-space imaging techniques19.

In our present work, we reveal the SkL and competing magnetic orders in centrosymmetric Gd3Ru4Al12, a good metal (ρ300K/ρ2K~6) in which the number of magnetic moments per layer and skyrmion is even smaller than in Gd2PdSi3. Using a combination of real-space imaging, x-ray and neutron scattering, as well as measurements of the THE, the presence of the SkL is established unambiguously. On the basis of our experimental observations, we discuss how a combination of frustrated RKKY interactions, local ion anisotropy, and thermal fluctuations in Gd3Ru4Al12 provides the first opportunity to observe metastable skyrmions in a centrosymmetric material.

Results

Structural properties

Gd3Ru4Al12 crystallizes in hexagonal space group P63/mmc with weak, yet finite, anisotropy in the magnetization and transport properties20,21,22,23. We label unit vectors aligned with the principal crystallographic axes as a, b, and c (Fig. 1a, b). For illustrative purposes, the structure may be decomposed into buckled Ru4Al8 layers (containing the inversion center) and perfectly planar Gd3Al4 sheets stacked along the c-axis (Fig. 1a). The magnetic moments at the rare earth (Gd) site materialize a breathing kagomé network (Fig. 1b), equivalent to a triangular lattice of trimer plaquettes formed by sets of three Gd moments22,23. A quantitative measure of the distortion away from the ideal Kagome structure is the ratio between alternating (breathing) nearest (r) and next nearest (r’) neighbor Gd–Gd distances, r/r′ = 0.73 (Fig. 1b). Information about crystal growth and characterization, as well as a discussion of standard experimental techniques used in this work, is provided in the Methods section and in Supplementary Note 1.

Fig. 1: Crystal structure and zero-field magnetic order of a Gd-based breathing kagomé lattice.
figure 1

a Hexagonal unit cell of Gd3Ru4Al12, where a, b, and c are crystallographic lattice directions. b Within the Gd3Al4 layer, rare earth (Gd) atoms form a distorted kagomé net with alternating distances r, r’ between nearest neighbors. Al and Ru atoms are not shown. The black rhombus indicates the size of the primitive unit cell. c Magnetic susceptibility (blue, left axis) increases continuously in the paramagnetic state as temperature is lowered. The inverse susceptibility H/M (red, right axis) is fitted by the Curie–Weiss expression (dashed line) at high temperature. d, e Specific heat cP(T) and M/H show two phase transitions in zero magnetic field. f, g At the (7, 0, 0) + q3 = (7 + q, −q, 0) incommensurate reflection, resonant x-ray scattering with polarization analysis provides modulated moments within (m⊥ q,c, blue solid triangles) and perpendicular to (m//c, red open triangles) the hexagonal plane, as well as the magnitude of the ordering vector q. Inset of (c) six directions of qi are allowed by symmetry. The black hexagon indicates a conventional unit cell in real space. The transition temperatures TN2 > TN1 bound the red shaded area in (d–g).

Full size image

Magnetic phase diagram from magnetization and transport experiments

In good agreement with previous reports22,23, we observe Curie–Weiss behavior in the magnetic susceptibility M/H at temperature T > 200 K (Fig. 1c). The dominant ferromagnetic interaction, associated with magnetic coupling within the trimer plaquette, is evident in the high Curie–Weiss temperature TCW ≈ 70 K. Long-range magnetic order sets in far below TCW at TN2 = 18.6 K, underlining the importance of competing magnetic interactions in this material, likely RKKY couplings which oscillate as a function of Gd–Gd distance22. Two sharp anomalies in cP(T) at TN1 and TN2 (Fig. 1d), a double kink profile of the magnetic susceptibility M/H (Fig. 1e), as well as two changes of slope in the zero-field longitudinal resistivity ρxx(T) (Supplementary Fig. 5) all suggest successive evolution of order parameters, as is frequently the case in magnets with competing interactions24,25,26. Microscopic evidence for the phase transitions in zero field was obtained using elastic x-ray scattering (REXS), in resonance with the gadolinium L2 absorption edge and with polarization analysis at the detector. Sinusoidal magnetic order in the hexagonal a–b plane at TN2 = 18.6 K can be separated from the onset of three-dimensional helical order at TN1 = 17.2 K (Fig. 1f). As the incommensurate magnetic modulation vector q was found to be aligned within the hexagonal plane in the scattering experiments (Fig. 1g), six directions of q are equivalent by symmetry. We label three of these directions by q1, q2, and q3 (inset of Fig. 1c). The qi vectors are locked to the a* (and equivalent) directions, where a* and b* are unit vectors in reciprocal space. The helical pitch length corresponding to q is λ = 2.8 nm at T = 2.4 K, much smaller than in typical non-centrosymmetric B20 compounds such as MnSi (λ ≈ 19 nm)2.

In the following, we establish the magnetic phase diagram using comprehensive measurements of the magnetic susceptibility χDC = ∂M/∂H and Hall conductivity σxy. We have corrected the significant demagnetization effect for these datasets, and more generally for all data recorded under isothermal conditions (as well as all phase diagrams). In this spirit, isothermal data are plotted as a function of the internal magnetic field Hint = H-NM, where H is the externally applied magnetic field, N is the demagnetization factor, and M is the bulk magnetization (Methods and Supplementary Table 2).

In the configuration H//c, the degeneracy of the qi is maintained. Several magnetic phase boundaries are marked by open black symbols in the contour plot of χDC (Fig. 2a). In anticipation of the REXS results of Fig. 3, we label this cornucopia of competing magnetic states as helical (H), transverse conical (TC), fan-like (F), SkL, field-aligned ferromagnet (FA-FM), and the as-yet unidentified phase V. Raw data of χDC are presented in Supplementary Fig. 7. Some signatures of the bulk phase transitions and large magnetoresistance, reported for single crystals as part of our work, can be observed even in polycrystals23.

Fig. 2: Field-induced phase transitions in Gd3Ru4Al12 and topological Hall effect in the SkL phase.
figure 2

a Contour plot of magnetic susceptibility χDC = ∂M/∂H in which we distinguish helical (H), fan-like (F), transverse conical (TC), and skyrmion lattice (SkL) states, as well as the as-yet unidentified phase (V). FA-FM labels the field-aligned ferromagnetic state. Black open circles are phase boundaries obtained from χDC for ∂H/∂t, ∂T/∂t < 0 (small symbols bound phase V. c.f. Supplementary Note 7). The transition between SkL and TC is marked by black open squares, corresponding to the onset of hysteresis between FC and zero-field cooled (ZFC) cases in panel (d). b Isotherms of the Hall conductivity σxy at various temperatures. Dashed and solid lines mark increasing and decreasing magnetic fields, respectively. A polynomial background was subtracted to identify the topological Hall conductivity σxyTHE (gray shaded area). Curves shifted by offsets of Δσxy = 680 Ω−1 cm−1. c Large positive σxyTHE obtained from (b) is observed only in the SkL phase. Black open dots and squares indicate phase boundaries as in (a). d σxy(T) at fixed external magnetic field, for FC (bold, black line) and ZFC (thin, red line) conditions. Outside the TC phase, ZFC and FC curves overlap. In the TC phase, the FC signal is larger due to the additional σxyTHE of metastable skyrmions (m-SkL, gray cross-hatched area). In the SkL state itself, we have estimated the background of σxy with a polynomial (dark gray dashed line for μ0Hint = 1.22 and 1.42 T) and marked the σxyTHE by gray coloring (no hatching). Curves shifted by offsets of Δσxy = 250 Ω−1 cm−1 for clarity.

Full size image
Fig. 3: Resonant elastic x-ray scattering (REXS) with polarization analysis clarifies the magnetic order of bulk Gd3Ru4Al12.
figure 3

The scattering plane is (H, K, 0). Blue open symbols correspond to Iπσ scattering intensity, related to magnetization components in the hexagonal basal plane (m2//q, m2⊥q,c), while red solid symbols are for Iππ intensities, proportional to m2//c. Data in (d–i) were obtained in the zero-field cooled (ZFC) state, whereas (k–l) report field-cooled (FC) measurements. a–i The comparison of scattering intensities at (4-q, 4, 0) and (4, 4-q, 0) reflections allows us to separate the two in-plane components m//q and m⊥q,c for the three magnetic orders sketched in (a–c), observed in the sequence (a-b-c) as the magnetic field is increased starting from the ZFC state at T = 2.4 K. Only the helical state in (d, g) has m//c ≠ 0. k–m At T = 7 K and μ0Hint = 1.5 T, the fan-like state coexists with helical modulations under FC conditions. As the fan-like state has m//c = 0, Iππ scattering directly probes the helical component. Comparable magnitude of Iππ ~ m2//c at the three incommensurate reflections around (7, 0, 0) is consistent with multi-q helimagnetic order. The nano-scale skyrmion lattice (SkL) proposed here is sketched in (j), with each cone representing the magnetic moment on an individual Gd-site. Red (blue) correspond to positive (negative) mz component of the local moment. Due to the lack of full refinement of the Gd-site moment, (j) visualizes but one possible manifestation of the SkL order. Inset of (m), red circles of radius ri ~ √Iππ ~ m//c, illustrating approximately equal weight for the three directions q1, q2, and q3 marked by light pink arrows. Error bars correspond to one standard deviation.

Full size image

Out of this large number of magnetic phases, the SkL is distinguished by a large THE due to the non-zero integer winding number of the magnetic texture and the resulting Berry curvature of conduction electrons27,28. This transport signature provides direct evidence for the chiral nature of the magnetic order in the SkL phase. In Fig. 2b, a box-shaped and strongly field- and temperature-hysteretic Hall conductivity signal (shaded in gray) emerges on the back of a smooth background in an intermediate range of magnetic fields. We approximate the background by a low-order (odd) polynomial and extract the topological Hall conductivity σxyTHE. The topological signal as obtained from the isothermal field scans is confined within the boundaries of the SkL phase (Fig. 2c). Meanwhile, measurements of σxy(T) at fixed magnetic field and for increasing temperature (dT/dt > 0) show a large split between curves recorded under zero-field cooled (ZFC) and field-cooled (FC) sample conditions, exclusively at intermediate field values (Fig. 2d, μ0Hint = 1.22 and 1.42 T). The natural conclusion is that a metastable SkL state, with its largely enhanced σxy(T), can be sustained at the lowest temperatures in the FC experiment, where the SkL is absent under ZFC conditions. This behavior suggests the stabilization of the SkL by thermal fluctuations (c.f. Discussion section). The point of divergence between the ZFC and FC curves at T = 5–8 K in Fig. 2d marks the first order phase transition between the TC and SkL states in our phase diagram (labeled in Fig. 2a,c by black open squares). Detailed susceptibiltity measurements evidence that the boundaries of phases TC and SkL with all surrounding phases are also strongly of first order (Supplementary Fig. 7).

Resonant elastic x-ray scattering (REXS) and microscopic magnetic structure

We now proceed to study the field-induced magnetic phases using REXS (Fig. 3) and real-space imaging (Fig. 4, next section), before finally returning to a semiquantitative analysis of the Hall signal. For polarization analysis in REXS, three mutually orthogonal components of the q-modulated magnetic moment m(q) are separated viz.34

$${\mathbf{m}}\left( {\mathbf{q}} \right) = {\mathbf{c}}m_{//{\mathrm{c}}}\left( {\mathbf{q}} \right) + {\hat{\mathbf{q}}}m_{//{\mathrm{q}}}\left( {\mathbf{q}} \right) + \left( {{\mathbf{c}} \times {\hat{\mathbf{q}}}} \right)m_{ \bot {\mathrm{q}},{\mathrm{c}}}({\mathbf{q}}),$$
(1)

where \({\hat{\mathbf{q}}}\) is a vector of unit length along q. In our experiment, the incoming beam of x-rays is linearly polarized with electric field component Eω within the π-plane spanned by ki and kf, the wave-vectors of the incoming and outgoing beams (π-polarization). Two components of the scattered x-ray intensity are separated at the detector: Iππ, with Eω remaining within the π-plane, and Iπσ, with Eω now perpendicular to the π-plane. In the scattering geometry where ki, kf⊥c, we have Iππ ~ m2//c always (see Methods). We chose the incommensurate satellite reflections at (4 + q, 4, 0) and (4, 4-q, 0) so that Iπσ ~ m2⊥q,c and Iπσ ~ m2//q, respectively. Starting from the ZFC state at T = 2.4 K and increasing the magnetic field, this convenient experimental configuration allows us to identify the helical ground state (H) with m⊥q,c, m//c ≠ 0 and m//q = 0 (H = 0, Fig. 3a, d, g), the transverse conical (TC) state with m//c = 0 and finite values for both in-plane components of m(q) (μ0Hint = 1.5 T, Fig. 3b, e, h), as well as the fan-like (F) state, which has only m⊥q,c ≠ 0 (μ0Hint = 2.9 T, Fig. 3c, f, i). It was confirmed that the incommensurate reflections vanish in the field-aligned state (not shown). The TC ground state in finite field is likely stabilized by weak in-plane anisotropy of the local magnetic moment. Weak in-plane anisotropy was also observed in magnetization measurements (Supplementary Fig. 7).

Fig. 4: Real-space observations of helical stripes and SkL in a (001) thin-plate sample of Gd3Ru4Al12 using Lorentz-transmission electron microscopy (L-TEM).
figure 4

We report data obtained from two different areas A, B of the same thin plate. The magnetic field was applied perpendicular to the sample plane. All data were recorded under field cooling (FC) conditions. a* and b* are crystalline axes in reciprocal space. Fast Fourier transform patterns in b, d correspond to real-space images (underfocused L-TEM images) shown in a and c, respectively. The defocusing length of the electron beam was ld = −3 μm (a, b) and ld = −193 nm (c–h). In (b), the scale of the data was magnified as compared to (d), cutting the featureless high-q regime for clarity. Red (yellow) circles in (b, d) mark intensity due to the crystal lattice (the magnetic order). Focusing on representative parts of the real space image (yellow dashed box), filtered images are shown as insets in (a, c) (same length scale as main panel). Helical stripes (T = 8 K, μ0Hint = 0.59 T, μ0H = 0.7 T) and SkL (T = 8 K, μ0Hint = 1.53 T, μ0H = 1.95 T) are visible in the data. e–h T-dependence of the fast Fourier transform patterns of the magnetic and lattice images at 1.53 T. The magnetic contrast of the SkL is suppressed above T = 17 K. See Methods for filtering procedure and other experimental details. Scale bars in panels (a, c) correspond to 10 nm.

Full size image

In Methods (Fig. 5) and Supplementary Notes 2, 3 we present bulk neutron scattering data obtained on a 160Gd isotope-enriched single crystal. Firstly, we find excellent quantitative agreement of small-angle neutron scattering (SANS) and REXS, indicating that the REXS experiment is not seriously affected by surface strain and can be used to characterize the bulk properties of Gd3Ru4Al12. Secondly, our neutron experiment with H//a* confirms the multi-domain nature of the zero-field helical ground state. Thirdly, neutron scattering also provides proof that the magnetic modulations on the breathing kagomé layers are ferromagnetically stacked along the c-axis, by ruling out magnetic reflections at (q, 0, (2n−1)/2) and (q, 0, 2n−1) for n = 1 and 2. As compared to the triangular lattice, kagomé structures introduce additional complexity due to the larger number of atoms per crystallographic unit cell. While the scope of this work does not include a full refinement of the magnetic structure, the interesting question of the local spin alignment on the trimer plaquette remains to be resolved in future studies35.

Fig. 5: Small-angle neutron scattering (SANS) experiments as obtained at TAIKAN for 160Gd enriched Gd3Ru4Al12 in H = 0 and at low T.
figure 5

a Qz = 0 cut of scattering intensity in the SANS experiment. Three magnetic reflections corresponding to the multi-domain helical ground state are observed. A fourth reflection (marked by a white circle) originates from a minority crystallographic domain with different alignment of the c-axis. Scattered neutron intensity is shown on a logarithmic scale. b Integrated intensity of elastic neutron scattering onsets around TN2 = 18.6 K. c Weak temperature dependence of the absolute value of the magnetic modulation vector q = 2π/λ on the order of 2% was observed. The values of q obtained independently from REXS and SANS on different crystal batches are in good agreement. d Half width at half maximum (w//) of the magnetic reflection in elastic neutron scattering increases as TN2 is approached from below. Red shaded area in (b–d) is bounded by TN1, TN2.

Full size image

We have also performed REXS experiments in the SkL phase at T = 7 K, μ0Hint = 1.5 T under field cooling (Fig. 3j–l). In this experiment, the three reflections (7 + q, 0, 0), (7, q, 0), and (7 + q, −q, 0)—corresponding to q1, q2, and q3 in the inset sketch of Fig. 3l—were chosen. We find very strong Iπσ(q3) ~ m2⊥q,c but weaker Iπσ(q1) and Iπσ(q2), a telltale sign of the fan-like state. The large fan-like signal in this experiment likely arises due to the proximity to a first order phase transition and associated phase separation. Crucially, there is also significant Iππ ~ m2//c with comparable intensities for all the three qi. Our data, taken with an x-ray beam spot size of ~1 mm2, suggest about 20–50% volume fraction (fV) of the helical component m//c. Roughly equal population of the scattering signal related to helical order for the three qi is consistent with a topological multi-q ordered state, such as a lattice of Bloch skyrmions, in the SkL phase.

Spin-vortices in real space imaging

In the scattering study presented here, all information about the relative phase of the three helical modulations making up the SkL is lost. Thus, these experiments cannot confirm the topological nature of the SkL state. Meanwhile, imaging of the real-space spin structure using Lorentz-transmission electron microscopy (L-TEM) on a thin-plate sample29 provides unambiguous evidence for spin vortices in the SkL state, as shown in the following (Fig. 4). For skyrmions with Bloch-type character, clockwise or counter-clockwise helicity (in-plane rotation direction of spins) should appear as bright and dark dots in the underfocused L-TEM images, respectively3. On the other hand, the helical structure with in-plane q-vector can be visualized as alternating bright and dark stripes in the underfocused L-TEM image.

All our L-TEM data were recorded under field cooled conditions. At T = 8 K and μ0Hint = 0.59 T, we found stripe-like magnetic contrast of a single-q helical order with pitch λ = 2.8 nm (Fig. 4a, b). At the same temperature, but at higher field μ0Hint = 1.53 T, the vortex-like pattern of our real space image translates into six-fold symmetric incommensurate reflections in the Fourier transform (Fig. 4c, d). In T-dependent measurements the magnetic contrast vanishes at T = 17–20 K (Fig. 4e–h). The region of the B–T phase diagram occupied by the SkL phase was found to be slightly expanded in L-TEM measurements on thin plates of Gd3Ru4Al12, as compared to bulk samples. This behavior is consistent with previous work on other skyrmion hosting materials3,4. Although the values of λ from real space imaging and scattering experiments are in quantitative agreement, local lattice strain appears to rotate q towards the a-axis in the thin plate sample (see Methods for a discussion). The strain effect is also manifested in the fast Fourier transform data of our highest quality real space image, which is slightly distorted (area B, Fig. 4d).

With the aim of amplifying the weak contrasts of magnetic skyrmions and helical stripes in the L-TEM data, we cut background noise by preserving selected fast Fourier transform components as exemplified in Fig. 4b, d: (i) yellow circles mark in-plane q-vectors related to the helical structure and the SkL, while (ii) red circles mark the Fourier components related to the atomic crystal lattice of Gd3Ru4Al12. The filtered fast Fourier transforms are then converted back to filtered real-space images as shown in the insets of Fig. 4a, c. Note that (ii) are visible only in the case of the high-field data measured at very small defocusing length of the electron beam (c.f. Methods, where more experimental details are provided). In combination with the scattering and transport experiments, our L-TEM study firmly establishes the presence of the SkL phase in this compound.

Estimate of spin polarization in the conduction band from THE

Armed with microscopic knowledge of the magnetic order in the SkL phase, we proceed towards a semiquantitative analysis of the THE. In the continuum approximation, the emergent magnetic field from hexagonal lattice skyrmion textures is Bem = −(h/e)√3/(2λ2) ≈ −460 T27,30. This enormous effective field is related to the topological Hall resistivity through the normal Hall coefficient R0, the volume fraction of the skyrmion phase fV, and the effective spin polarization P of conduction electrons30

$$\rho _{yx} ^{\;\;\;\;\;{\rm{THE}}} = f_V \cdot P \cdot R_0 \cdot B_{em}$$
(2)

After extraction of the topological Hall conductivity σxyTHE (Fig. 2b), we estimate the topological Hall resistivity as ρyxTHE = σxyTHE ∙ ρxx2. Extrapolating the value of R0 from higher temperatures (see Supplementary Note 5) and using fV = 20–100% in the SkL phase, we arrive at P = 0.01–0.05 (higher fV corresponds to lower P). These values of P appear reasonable in comparison with related materials such as Gd2PdSi3 (P = 0.07)19. Note that the continuum approximation underlying Eq. (2) may be rendered inaccurate when λ becomes comparable to the crystallographic lattice spacing. Hence, the observed magnitude of the THE should be taken as a merely semiquantitative measure of Bem.

Discussion

Our combined experimental effort shows that a topological SkL phase is stabilized in the centrosymmetric breathing kagomé lattice Gd3Ru4Al12. This system charms with conceptual simplicity: Large local spins, whose moment size is affected little by thermal fluctuations, are coupled weakly to a sea of conduction electrons. In the SkL phase, the rare earth lattice imparts its scalar spin chirality onto the conduction electrons and drives a giant THE in a limited window of temperature and magnetic field.

Superficially, the magnetic phase diagram (Fig. 2) suggests similarities with Bloch skyrmions in chiral magnets; the SkL is thermodynamically stable only at elevated temperatures, and a metastable skyrmion state survives at low T under field-cooled conditions. Unlike in chiral magnets, however, thermal fluctuations are by no means necessary to stabilize equilibrium topological spin textures in centrosymmetric lattices. On the one hand, an equilibrium SkL was observed experimentally in triangular lattice Gd2PdSi3 down to very low temperature (at least T/TN = 0.1)19. On the other hand, numerical simulations have consistently shown extended parameter regimes with a ground state equilibrium SkL, both in the case of frustrated exchange interactions31 and RKKY systems with four-spin-two-sites (biquadratic) couplings32,33. In the present system, we have mapped the RKKY interaction to an effective spin model which shows frustrated antiferromagnetic couplings at the second nearest neighbor level (Supplementary Note 4). As this is a metallic system however, treatment as a Kondo lattice with RKKY couplings and possibly sizeable four-spin (biquadratic) couplings may be more appropriate for a quantitative description. A prominent role for thermal fluctuations is re-introduced in the case of Gd3Ru4Al12 due to a particularly delicate balance of terms in the free energy when the magnetic field is finite: rather weak in-plane anisotropy (Supplementary Note 7) of the magnetic moments favors transverse conical (TC) order at the lowest T, while thermal fluctuations promote the SkL state. In consequence, metastable skyrmions could be observed in transport experiments for the first time in a centrosymmetric magnet. We stress that whereas Gaussian fluctuations of the moment amplitude provide the route to the SkL in non-centrosymmetric magnets1, thermal fluctuations of the moment direction away from the easy plane of local anisotropy are expected to be essential in Gd3Ru4Al12.

Our real space imaging study provides compelling evidence for incommensurate, topological multi-q order, for the first time in a centrosymmetric crystal with tiny (a few nanometer sized) magnetic vortices. More broadly, our results show that the breathing kagomé lattice is a rich host for new topological magnetic phases, which may exert a colossal emergent magnetic field to produce unprecedented electrical, thermal, and thermoelectric phenomena, like the large THE reported here.

Methods

Synthesis and characterization

Large, cm-sized crystals of Gd3Ru4Al12 were grown under Ar gas flow in a floating zone (FZ) furnace equipped for high vacuum operation. The samples were characterized by powder x-ray diffraction (XRD), energy-dispersive x-ray spectroscopy (EDX), and examination under an optical microscope equipped for polarization analysis. Laue XRD was used to obtain samples with oriented surfaces. We also grew a crystal rod using the Czochralsky pulling technique, but found that the FZ approach is more stable and reproducibly yields high quality crystals. As Gd in natural abundance is a strong neutron absorber, an additional 160Gd isotope-enriched single crystal batch was grown for elastic neutron scattering experiments. Due to oxide impurities in the 160Gd raw material, this growth was of lower quality as characterized by transport (RRR = ρ(300 K)/ρ(2 K) ~ 2.6) and the sharpness of bulk phase transitions in the magnetization data. The final mass of the sample used for elastic neutron scattering was about 15 mg.

Bulk measurement techniques

Magnetization and heat capacity were measured using cube- or rectangular cuboid-shaped polished samples in commercial Quantum Design MPMS and PPMS cryostats, in an attempt to minimize adverse effects due to the demagnetization field. The magnetic field H was applied parallel to the crystallographic c-axis for all data shown in the main text. Due to experimental constraints, high-field vibrating sample magnetization (VSM) measurements were complemented by low-field (μ0H < 7 T) extraction magnetization data (DC-M) on the same sample in the same configuration, to arrive at a more reliable estimate for the absolute value of M. The VSM data were then scaled to the DC-M results (scaling factor ~1.06). These measurements were performed in a Quantum Design PPMS-14T (VSM) and a Quantum Design MPMS3 (DC-M) system, respectively.

Transport experiments

Transport measurements were carried out using polished and oriented plates of approximate dimensions 2.5 × 1 × 0.15 mm3 with electrical contacts attached by silver paste (Dupont) or H20E silver epoxy (Epo-Tek). For the data presented in Fig. 2, the face of the sample plate was perpendicular to the c-axis (H // c) and the charge current density was J // a*. Similar data were obtained for a sample with J // a (Supplementary Fig. 6). For the transport experiments, we used a Quantum Design PPMS cryostat for temperature and magnetic field control, but a custom arrangement of lockin- and pre-amplifiers replaced the PPMS measurement circuits. The applied current density was J ~ 3 × 104 A/m2 and the excitation frequency was 9–15 Hz. We calibrated the absolute value of ρxx using a long, bar-shaped crystal of dimensions 3 × 0.5 × 0.5 mm3 in order to reduce systematic errors arising from the measurement of the sample geometry when using a standard optical microscope. For J // a or a*, the resistivity was 110(5) μΩ cm at room temperature. The field dependence of ρyx and ρxx was calculated from anti-symmetrized and symmetrized voltage traces, respectively. Due to significant hysteresis at low temperature, we recorded full field ramps for both increasing and decreasing magnetic field and paired, for example, the data with dH/dt < 0, H < 0 with dH/dt > 0, H > 0 for the (anti-)symmetrization routine.

Resonant elastic x-ray scattering (REXS)

Magnetic x-ray scattering experiments were conducted in resonance with the Gd L2-edge on beamline BL-3A at Photon Factory in KEK, Japan with the scattering plane being (H, K, 0), i.e. the incoming and outgoing beams of wave vector ki and kf were both in the plane perpendicular to the crystallographic c-axis. The 006 reflection of a pyrolytic graphite (PG) plate was used to analyze the polarization of the scattered beam, with the 2θ angle at the Gd-L2 edge of PG fixed to 88 degrees. Our single crystal with large (110) planes was set in a cryostat equipped with a vertical 8 T superconducting magnet, so that the magnetic field was applied parallel to the c-axis.

The expression for the magnetic part of the scattering amplitude in REXS is written as34

$$f_{{\rm{res}}} = C_0{\mathbf{\varepsilon }}_f^ \ast \cdot {\mathbf{\varepsilon }}_i + {\rm{i}}C_1\left( {{\mathbf{\varepsilon }}_f^ \ast \times {\mathbf{\varepsilon }}_i} \right) \cdot {\mathbf{m}} + C_2{\mathbf{\varepsilon }}_f^\dagger O \, {\mathbf{\varepsilon }}_i$$
(3)

with constants Ci, local magnetization m, and εi and εf the initial and final polarization of the x-ray beam, respectively. The scattering intensity is Ires = fres ∙ fres*. Only the second term depends explicitly on m. In our experiment, the incident beam had the linearly polarized electric field component in the scattering plane (π-polarization). Recalling that ε ∙ k = 0 for light, let us separate two components of the scattered radiation: (1) εi and εf are both in the scattering plane π. As the angle between ki and kf is defined as 2θ, it follows that Iππ ~ (sin(2θ) mz)2. (2) εi is in the π-plane, and εf is perpendicular to it. Their cross product is directly proportional to ki, and we write Iπσ ~ (ki ∙ m)2. The π-σ intensity thus depends very sensitively on the relative alignment of the incoming beam and the modulated magnetization m = m(q). A comparison of magnetic scattering at different reflections (Fig. 3) can be used to separate the two components (m⊥q,c, m//q) of the in-plane magnetic moment.

Lorentz-transmission electron microscopy (L-TEM)

A (001) plate of Gd3Ru4Al12 was cut from a bulk sample with dimensions 2 × 2 × 0.1 mm3. The plate was further thinned by an Ar-ion milling process subsequent to mechanical grinding. Real-space observations were performed by L-TEM using a multifunctional transmission electron microscope (JEM2800, JEOL) equipped with a double-tilt helium cooling holder (ULTDT). We measured the transmittance of the incident electron beam with accelerating voltage of 200 kV. The temperature of the thin plate was carefully controlled from 7 to 50 K.

L-TEM is a useful method to observe magnetic order in thin plates with thickness below several hundred nanometers due to the interactions between the incident electron beam and magnetic moments29. Data in low field (μ0Hint = 0.59 T, μ0H = 0.70 T, Fig. 4a, b) were obtained at defocusing length ld = −3 μm, where only magnetic contrast is visible. The instrument cannot be operated at these same settings in high magnetic field (μ0Hint = 1.53 T corresponding to μ0H = 1.95 T, Fig. 4c–h): Instead, measurements were carried out at small ld = −193 nm, so that both the crystal lattice and the magnetic order leave their fingerprints in the data. These high-field experiments conveniently show that the magnetic q-vector is aligned with the crystallographic a-axis in the thin plate sample, in contrast to the x-ray and neutron scattering results on bulk samples, where q // a*. We found that further analysis of the data using the transport-of-intensity equation (TIE) approach was unreliable in the present case, where the size of the spin vortices is comparable to the real-space resolution of the L-TEM experiment and the magnetic contrast is relatively weak.

Small-angle neutron scattering (SANS) experiments

For our neutron scattering study, we used the time-of-flight (TOF) type small-and-wide-angle neutron scattering instrument TAIKAN(BL15) at the Material and Life Science Experimental Facility (MLF) of J-PARC in Japan. Pulsed neutrons from a spallation source with a distribution of wavelengths centered around λn,max = 2.5–3 Å are incident on the sample; some of the neutrons are scattered elastically from the crystal and deflected so that kf ≠ ki (ki, kf being the wave vector of the incoming and outgoing neutron, respectively). The time delay between the spallation event and the arrival at the detector bank allows to extract the wavelength λn = 2πkin−1 of each neutron, while its position on the detector screen corresponds to the scattering angle 2θ between ki and kf. From λn and 2θ, the components Qx, Qy, and Qz of the scattering vector Q = ki−kf can be calculated.

Figure 5a shows a two-dimensional Qx-Qy map of scattering intensity (a Qz = 0 cut of the full TOF dataset) measured using the small-angle neutron scattering (SANS) detector bank of TAIKAN in zero field and T = 3.58 K. These data, obtained before the resonant x-ray scattering experiment discussed in the main text, provided the first evidence for incommensurate magnetic order in Gd3Ru4Al12. The magnetic modulation vector was found to be along the a* crystallographic direction and those other directions equivalent to a* by hexagonal symmetry. The six-fold symmetric pattern on the detector arises due to scattering from multiple domains with different alignment of q, as we show in Supplementary Note 2. Due to the relatively large absolute value of q = 0.272 r.l.u. as compared to, for example, MnSi2, our sample has to be tilted slightly away from ki // c in order to observe the SANS reflections corresponding to incommensurate magnetic order (i.e., if ki // c, Q is not perfectly aligned within the a-b plane). In such a tilted configuration, only the three reflections with Qx > 0 are accessible out of the total six peaks allowed by hexagonal symmetry. However, an additional reflection at Qx < 0 in Fig. 5a arises due to a second crystallographic domain in the 160Gd enriched sample. We estimate that the minority domain has about 25% volume fraction of the total crystal, with its c-axis about 6° misaligned from the c-axis of the major domain.