Introduction

Frustration in magnetic systems has attracted great interest due to the ability to support diverse unconventional states, such as spin-glass behavior 1, incommensurate magnetic order2,3, and noncoplanar magnetic patterns4,5. Magnetic frustration can arise from the geometry of the lattice or competing exchange interactions. The antiferromagnetic (AFM) Heisenberg model on the two-dimensional (2D) hexagonal lattice is a paradigm for frustrated magnets. With only AFM nearest-neighbor (n.n.) interactions, the ground state of this system is the three-sublattice Néel structure consisting of coplanar spins forming ± 120° angles between nearest neighbors, which is commensurate to the underlying lattice. Including higher-order multi-spin interactions can stabilize more complex magnetic states such as superposition of spin spirals6,7,8,9 or atomic-scale spin lattices10,11. Those multi-Q states are intriguing, as they often exhibit nonzero vector chirality, Si × Sj, and scalar chirality, Si ⋅ (Sj × Sk). These chiralities lead to unconventional phenomena, such as the spin Hall effect12,13, the topological Hall effect14,15,16,17,18, the noncollinear Hall effect 19, nonreciprocal transport20,21, complex magnetoresistances22,23, and the formation of topological orbital moments24,25,26,27, as well as topological high-order interactions 28. A notable example of a multi-Q state is the triple-Q state6, characterized by the superposition of three symmetry-equivalent spin spirals, where the relative angle between all nearest-neighbor spins is the tetrahedron angle of 109.47°29. The triple-Q state in 2D was first predicted in an Mn monolayer on the Cu(111) substrate and then experimentally observed in an hcp-stacked Mn monolayer on Re(0001)7. Subsequently, the triple-Q state has also been observed in Pd/Mn and Rh/Mn bilayers on Re(0001) through SP-STM studies30.

In contrast to the 2D hexagonal lattice, the classic ground state of Heisenberg kagome antiferromagnets can exhibit macroscopic degeneracy, arising from infinitely many spiral states with the same energy. If the system does become magnetically ordered, the ground state is then often selected by the order-by-disorder mechanism31,32, driven either by quantum or thermal fluctuations. Recently, triple-Q states, which can be described as skyrmion lattice phases, have been discovered in the rare-earth compound Gd3Ru4Al1233,34, where Gd forms stacked breathing kagome lattice layers. For quantum spins on the kagome lattice, exotic ground states such as spin liquids or valence-bond crystal states may be stable35,36,37,38. Mn is a good candidate for investigating magnetic frustration as it often forms lattices composed of triangles and tends to have AFM interactions29,39,40. In this work, we explore the magnetic breathing kagome lattice formed by a Mn monolayer on a heavy-metal substrate, capped by h-BN, where the kagome structure arises via kagomerization—a spontaneous rearrangement of transition-metal atoms catalyzed by the h-BN overlayer, as uncovered in our prior study41. Our magnetic breathing kagome lattice structurally differs from conventional ones by maintaining near-ideal kagome geometry without significant bond-length alternation. The breathing behavior arises from magnetic modulation induced by the nearby h-BN layer, which modulates the electrostatic potential and creates strong bonds between Mn and N. The magnetic properties and magnetic interactions, including Heisenberg exchange interactions and Dzyaloshinskii-Moriya interactions (DMI), have been investigated. Compared to the triangular lattice of Mn on Pt substrate, the kagome structure of Mn induces strong frustration, as evidenced by the almost flat bands obtained for spin spiral energies not only when considering the n.n. interactions but also when including further-neighbor interactions. Due to the microscopic properties of this magnetic breathing kagome lattice, the spin spiral energy minimum is found on the Γ-K line, and a triple-Q state is uncovered through atomistic spin dynamics. The strong noncoplanarity of this spin texture is characterized via its spin structure factor (SSF) and by its non-vanishing topological charge density. Additionally, a more complex spin state, characterized by several Q-pockets in the spin structure factor, emerges due to the flat bands. Such spin textures enabled by the kagomerization mechanism may offer a platform for the investigation of novel topological phenomena.

Results

Magnetic interactions in breathing magnetic kagome lattice

We show in Fig. 1 the atomic structure of the considered kagome structure Pt/Mn/h-BN. As we found in ref. 41, this metastable kagome arrangement of Mn atoms is enabled by the h-BN overlayer and should be experimentally realizable under appropriate growth conditions. In the optimized structure of Pt/Mn/h-BN, all N atoms of the unit cell are on top of an Mn atom, except for one N that instead attracts one Pt atom from the substrate towards the hexagonal hole of the kagome lattice. More details on how the kagome lattice is formed can be found in our previous work41.

Fig. 1: Atomic structure of Pt/Mn/h-BN.
figure 1

a Top view with the hexagonal 2D Brillouin zone. a1 and a2 are the in-plane Bravais vectors, and a is the lattice parameter of the surface unit cell. b1 and b2 are the vectors of the corresponding reciprocal lattice. b Side view. In each unit cell, there are three Mn atoms, labeled as Mn(1), Mn(2), and Mn(3).

Full size image

To analyse the source of the magnetic frustration, we calculated the pairwise interactions, including Heisenberg exchange (J) and DM vectors (D), for the Pt/Mn/h-BN kagome and Pt/Mn 2D hexagonal structures. The magnetic interactions were calculated from the disordered local moment (DLM) state (see Methods), which serves as an excellent reference for the magnetic state due to its neutrality towards both ferromagnetic (FM) and AFM preferences. In order to assist the analysis, we also removed the h-BN overlayer from the Pt/Mn/h-BN kagome structure while retaining the kagome arrangement of the Mn atoms. We present the spatial distribution of the Heisenberg exchange interactions J of the three structures in Fig. 2a–c. Each circle represents an Mn atom and is colored as a function of J values with respect to the central Mn (gray circle), which is taken as a reference. The value of J is given by the color scale, with red (blue) for FM (AFM) coupling. Figure 2a depicts the results obtained for the Pt/Mn triangular structure with C3v symmetry. The n.n. distance is d1 = 2.81 Å, the spin moment of the Mn atoms is 4.20 μB, and the strong AFM interactions in the six nearest neighbors lead to the 120° Néel state. In the kagome Pt/Mn structure shown in Fig. 2b, the n.n. bonds are shorter and unequal, so that four nearest neighbors form one large and one small equilateral triangle, with side lengths of d1 = 2.46 Å and \(d^{\prime} =2.41\) Å, respectively. The spin moment is reduced to 3.66 μB, and the slightly shorter distance of \(d^{\prime}\) compared to d1 leads to a stronger AFM exchange interaction, with a difference of 7 meV. With the presence of h-BN, as illustrated in Fig. 2c, the magnetic moment of Mn decreases further to 2.15 μB, and the asymmetry in the exchange interactions increased to ~16 meV. Such a significant asymmetry in the nearest-neighbor J values motivates us to term it a magnetic breathing kagome lattice, with \(J^{\prime} /{J}_{1}=1.6\) measuring the strength of the breathing bond alternation 42. More details about the magnetic interactions can be found in Supplementary Fig. S1 and Table S1, as well as in the input file for the Spirit simulations (“pairs-5R.txt”) found in Supplementary data package.

Fig. 2: Heisenberg exchange interactions and their lattice Fourier transform.
figure 2

Spatial distribution of the Heisenberg exchange interactions for: a Pt/Mn with Mn in triangular lattice; b Pt/Mn with Mn in kagome lattice; c Pt/Mn/h-BN kagome structure. In a d1 = 2.46 Å (\(d^{\prime}\)= 2.41 Å) is the side length of the large (small) equilateral triangle in the kagome structure. The small white and black circles in panel c represent N and B atoms, respectively. d–f Eigenvalues of the Fourier-transformed magnetic interactions, including Heisenberg exchange interactions and DMI for: a Pt/Mn hexagonal structure; b Pt/Mn kagome structure; c Pt/Mn/h-BN kagome structure. The flat bands exhibit a narrow bandwidth of ~6 meV in panel f.

Full size image

In Fig. 2d–f, we present the eigenvalues of the Fourier-transformed magnetic interactions \({\mathcal{J}}({\bf{q}})\) as a function of reciprocal momentum vector q, including Heisenberg exchange and DMI with a converged real-space cut-off radius of 5a, where a = 4.87 Å is the lattice constant of Pt, as defined in the “Methods”. The eigenvector associated with the lowest energy eigenvalue is a proxy for the magnetic ground state, and to explore the energy landscape, we plot the eigenvalues along high-symmetry lines in the 2D Brillouin zone. The Heisenberg exchange interactions control the overall dispersion of the eigenvalues of \({\mathcal{J}}({\bf{q}})\), while DMI lifts their degeneracy. In Fig. 2d, the global minimum is found at the Γ point of the Pt/Mn triangular lattice and corresponds to the Néel state—AFM merons43 can be found in this Néel background (see Supplementary Fig. S2). In contrast to the triangular lattice, the eigenvalues of \({\mathcal{J}}({\bf{q}})\) for the Pt/Mn kagome structure shown in Fig. 2e have a rather flat band at −100 meV, which is due to the geometric frustration built into the kagome lattice. Such flat bands can lead to the quasi-degeneracy of multiple magnetic states, increasing the importance of quantum and thermal fluctuations. The asymmetry of the n.n. Heisenberg exchange interactions shown in Fig. 2b lead to the gap opening at the K point seen in Fig. 2e. The presence of h-BN increases the size of the gap at the K point due to the larger asymmetry in the n.n. J (Fig. 2c), and modifies the dispersion of the flat bands with the energy minimum at q = (0.195, 0.195) (in reciprocal lattice units), as shown in Fig. 2f.

To further understand the properties of the flat band and the origin of the energy minimum we have repeated this analysis considering different real-space cutoff radii in the calculation of \({\mathcal{J}}({\bf{q}})\), see Supplementary Fig. S3. We found that the flat bands persist even when considering only the highly asymmetric n.n. interactions, and that the energy minimum is essentially determined by the long-range Heisenberg exchange interactions, with the DMI making a small correction. The RKKY effect, describing an indirect exchange interaction between localized magnetic moments mediated by conduction electrons, manifests in our system as a long-range, RKKY-like magnetic coupling due to the complex Fermi surface and itinerant nature of the magnetism. We conclude that unconventional phenomena related to flat bands can also be expected in this type of magnetic breathing kagome lattice.

Noncoplanar triple-Q state

The information obtained from the \({\mathcal{J}}({\bf{q}})\) was used to inform and interpret atomistic spin dynamics simulations to fully characterize the magnetic ground state. We also took into account the magnetic anisotropy of the system, which is calculated to be −0.14 meV, indicating a weak easy-plane magnetic anisotropy. We found a highly noncollinear state in Fig. 3a with a magnetic unit cell comprising 36 sites (dashed rhombus). To gain more insights, we plot next to the magnetic structure the intensity of the spin structure factor \(S(q)=\sqrt{{\bf{S}}{({\bf{q}})}^{* }\cdot {\bf{S}}({\bf{q}})}\) (see Methods). The spin structure factor analysis can reveal if the simulated noncollinear state is composed by a few superimposed spin spirals (ordered state) or by very many (disordered state). The six dark blue spots in the 2D Brillouin zone show that this is a triple-Q state with the dominant q at the K/2 points. However, the periodicity revealed by the spin structure factor is different from what we obtained in the analysis of the eigenvalues of \({\mathcal{J}}({\bf{q}})\) in Fig. 2e. The eigenvectors corresponding to the lowest eigenvalue of \({\mathcal{J}}({\bf{q}})\) lead to single-Q states which do not necessarily have unit spin length in real space, and so cannot be directly compared to the results of the atomistic spin dynamics. In addition, the in-plane anisotropy of the Pt/Mn/h-BN structure also favors the stabilization of the commensurate spin texture observed in Fig.3a, and this was not included in the \({\mathcal{J}}({\bf{q}})\) picture. The noncollinear state in Fig. 3a can be decomposed into three sublattices, labeled as Mn(1), Mn(2), and Mn(3). The spin structure factors for each sublattice (Fig. 3b–d) reveal that each sublattice has a different dominant q, given by the green spots circled in black. The magnetic arrangement of each sublattice can be approximately understood as a distorted spin spiral state, and we can obtain the main rotation axis about which the magnetic moments spiral by computing S(q)* × S(q) from the dominant q contribution to the spin structure factor. These rotation axes are tilted from the z-axis by about 25° and are rotated 120° around the z-axis relative to each other. Repeating the simulations starting from different initial random states leads to seemingly different magnetic configurations of this noncollinear 36-site magnetic state, but they can be verified to be equivalent by having the same properties of the spin structure factor just described.

Fig. 3: The triple-Q state and its spin structure factor.
figure 3

a The triple-Q state for the Pt/Mn/h-BN kagome structure. The spins are colored as shown in the colorbar: red is +z, blue is −z, and gray is in-plane. The dashed rhombus indicates the magnetic unit cell. The intensity of its spin structure factor is shown as a colormap on the 2D Brillouin zone, where the darker the color, the higher the intensity. The triple-Q state can be decomposed into three sublattices: b sublattice of Mn(1); c sublattice of Mn(2); d sublattice of Mn(3). Insets in b–d are the spin structure factor of each sublattice. The main q direction of each sublattice is circled in black, marked as K1/2, K2/2, or K3/2. The atomistic spin dynamics simulations were performed on a 30 × 30 supercell in the zero temperature limit.

Full size image

The discovered noncollinear triple-Q state is noncoplanar, as shown in the three-dimensional image in Fig. 3a, so it will imprint an emergent magnetic field on the conduction electrons, which defines the topological charge (see Method section and Fig. S4). A finite integer charge confers topological properties to the underlying spin-texture, enhancing its stability by giving rise to a protecting topological barrier while inducing the topological Hall effect. Within the 30 × 30 simulation lattice, the topological charge reaches an integer value of ~−300. It is interesting to remark that, unlike the topological charge of skyrmions, which is governed by the DMI, the non-zero topological number of our triple-Q state is primarily contributed by the Heisenberg exchange interactions. We would expect a rich set of nontrivial Hall effects in this topologically nontrivial spin texture as an implication of its strong local noncoplanarity.

Q-pockets and thermal perturbations

However, when we used a larger simulation lattice with 120 × 120 unit cell, we could also observe a more complex spin state with an energy 0.18 meV/spin higher than the triple-Q state, as shown in Fig. 4a. We can discern several magnetic domains, and the spin structure factor for this magnetic state with a slightly higher energy than the ground state is shown in Fig. 4b, where several Q-pockets can be observed. Those Q-pockets are defined by a distribution of favorable Q-states, which leads to a superposition of Q-states with different periodicities residing in each pocket. Such features are reminiscent of the Q-space behavior of self-induced spin glass in the magnetic ground state of elemental crystalline Nd44, which also arises from competing Heisenberg exchange interactions. The small energy cost for creating such a superposition of different Q-states led us to investigate how the magnetic state of the system changes under thermal perturbations, with the results summarized in Supplementary Fig. S5. The main observation can be rationalized by considering the flat bands in Fig. 2f, which have a bandwidth of about 6 meV and corresponds to a temperature of 70 K. Therefore, when the temperature in the atomistic spin dynamics simulations is increased past that threshold, the peaks of the spin structure factor melt away and we find a complex disordered magnetic state. This demonstrates that the magnetic interactions embodied in the flat bands are indeed central in determining not only the ground state but also the thermal properties of this system.

Fig. 4: Image of the magnetic state and the spin structure factor with a 120 × 120 simulation lattice.
figure 4

a Magnetic domains with an energy 0.18 meV/spin higher than the triple-Q state. b The spin structure factor of the magnetic state on the Pt/Mn/h-BN kagome structure, relaxed within a 120 × 120 simulation lattice at 0 K.

Full size image

Discussion

We have previously proposed a method to form a monolayer 3d transition metal kagome lattice using h-BN as an aid and focused on its implications for ferromagnetic layers 41. Inspired by this approach, here we investigated the frustrated magnetism of AFM Mn on a metastable kagome lattice realized in such a sandwiched structure, with Pt(111) as the substrate and h-BN as the overlayer. We adopted a two-pronged approach based on extracting all the relevant magnetic interactions from first-principles calculations combined with an exploration of the magnetic states through atomistic spin dynamics. Our analysis shows that the geometrical frustration in the Mn kagome lattice leads to rather flat bands in the Fourier-transformed magnetic interactions, which signify the almost degeneracy of multiple magnetic states. h-BN breaks the symmetry of the ideal kagome lattice without significantly affecting the magnetic flat bands, resulting in a magnetic breathing kagome lattice.

Analysis of the Fourier transform of the magnetic interactions indicated a tendency for a noncollinear magnetic ground state, and atomistic spin dynamics simulations revealed that the ground state is indeed a specific type of triple-Q state that has not been previously reported, to the best of our knowledge. We have fully characterized the properties of this triple-Q state by considering its real and reciprocal space structure, in particular, exploiting the information encoded in its spin structure factor. Additionally, we found that each of the three sublattices can be characterized as a distorted spin spiral state with a dominant q vector, with the dominant q-vectors for different sublattices making 120° angles relative to each other, and their combination leading to the observed triple-Q state.

Triple-Q states on the breathing kagome lattice have been previously discussed, but we expand on them in significant ways as we explain by comparing to selected works. Reference 42 presents numerical simulations for a simplified spin model, which provides guidance on the possible magnetic phases and their relation to the model parameters, but not discussing concrete material realizations of the proposed physics. This work discusses a type of triple-Q state based on a superposition of the q-vectors corresponding to the three inequivalent M points of the Brillouin zone. Comparing to our DFT results for Pt/Mn/h-BN, a qualitative difference arises from the range of the considered exchange interactions. Reference 42 considered a model with only three distinct exchange interaction parameters, while in our case, the complex Fermi surface and itinerant nature of the magnetism lead to a long-range, RKKY-like set of magnetic exchange interactions. The different range of the magnetic interactions also explains why the star of q-vectors involved in the triple-Q state for Pt/Mn/h-BN is near the K/2 points, while in ref. 42 it is the M/2 points.

An example for a real magnetic material is the rare-earth compound Gd3Ru4Al1233,34, which is a three-dimensional system hosting a skyrmion lattice phase characterized by q-vectors of about M/2. The experimental evidence for the proposed magnetic state is very convincing, although the respective topological charge has not been reported. There is a qualitative difference between the breathing kagome lattice of Gd3Ru4Al12 and ours. Letting d and \(d^{\prime}\) be the two different nearest-neighbor bond lengths (\(d^{\prime} \le d\)), the breathing kagome lattice for Gd3Ru4Al12 is mostly structural (short and long bonds), with \(d^{\prime} /d=0.73\). In contrast, the breathing kagome lattice that we found for Pt/Mn/h-BN is very close to the ideal kagome structure with \(d^{\prime} /d=0.98\). The strong spatial anisotropy of the magnetic interactions, which is characteristic of the breathing kagome lattice, is actually imprinted by the h-BN layer on the Mn lattice—it both modulates the electrostatic potential and creates strong bonds between Mn and N. Hence, besides the different dimensionality of the two systems, also the specific mechanisms by which the magnetic kagome lattice physics is generated are quite different.

As a final example, we consider Co1/3TaS245, for which the magnetic lattice is composed of stacked triangular layers, and in this case, the triple-Q state needs additional magnetic interactions in order to become the ground state—the authors proposed a biquadratic coupling as a possible mechanism. We found that for Pt/Mn/h-BN the triple-Q state could be stabilized with only isotropic Heisenberg exchange interactions, as also proposed in ref. 42. Comparing Co1/3TaS2 with Pt/Mn/h-BN, we conclude that the stabilization mechanisms for the triple-Q states are qualitatively distinct, besides their different periodicities and overall noncollinear spin arrangements.

The most striking feature of this triple-Q state is its very strong local noncoplanarity, which we quantified by mapping its topological charge density. Such a strong noncoplanarity will have a dramatic impact on the motion of conduction electrons through this magnetic state, and we anticipate highly nonlinear Hall effects as a consequence. The relative simplicity of fabrication of such heterostructures in laboratories that are currently investigating 2D materials makes the experimental realization of our proposed system a very concrete possibility, and its unconventional magnetic and transport properties can be accessed by various techniques, most notably spin-polarized scanning tunneling microscopy.

Concerning the almost flat bands obtained for spin spirals, we also found more complex spin states with a higher energy than the triple-Q state, where the spin structure factor exhibits characteristics reminiscent of self-induced spin glass behavior. In this context, it is also interesting to make a comparison with the results reported for Nd in ref. 46 (see their Fig. 2B). That paper introduces the concept of self-induced spin glassiness, which is generated when the interactions are long-ranged and competing but there is no chemical disorder in the system, as in more familiar spin glass systems such as Cu-Mn alloys. We believe that this concept is relevant to the results of our simulations, as the Q-pockets shown in Fig. 4b correspond to the magnetic structure shown in Fig. 4a, which can be described as patches formed by domains of different triple-Q states separated by narrow domain walls. The broadening of the spin structure factor peaks centered at the K/2 positions is then mostly contributed by the existence of the domain walls and the resulting local distortions to the ideal magnetic structure of the triple-Q states. A competing explanation is provided by the concept of short-range order, as discussed, for example, in ref. 47. This is illustrated by increasing the temperature of the system for the 30 × 30 simulation cell. At low temperature, we have the single domain triple-Q state, but for T = 70 K, Fig. S4 shows the typical signature of short-range order around the originally sharp peaks at K/2. These are not elongated as the Q-pockets in Fig. 4b and indeed resemble more the fuzzy structures reported in Fig. 4c of ref. 47. So, although there are similarities between the Q-pocket picture and the short-range order picture, we believe that they are describing different aspects with regard to our simulations.

In conclusion, the atomic structure enabled by the kagomeization mechanism provides a novel platform for investigating frustration-related properties and complex noncollinear magnetic states. Our findings offer valuable insights into the interplay between geometry and magnetism, paving the way for advancements in exploring novel spin textures.

Methods

First-principles calculations

Optimization of the considered structures was performed using the density functional theory (DFT) approach implemented in the Quantum Espresso computational package48, including the van der Waals correction (DFT-D349). We employed the projector augmented wave pseudopotentials from the pslibrary50, with the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE)51 as the exchange and correlation functional. Convergence tests led to a plane-wave energy cut-off of 90 Ry and an 18 × 18 × 1 k-point grid. We optimize the lattice constant of the HM substrates and use those values for the heterostructures with TM layers and with or without h-BN.

The magnetic properties and magnetic interactions were computed using the all-electron full-potential scalar-relativistic Korringa-Kohn-Rostoker (KKR) Green function method52,53,54, including spin-orbit coupling self-consistently as implemented in the JuKKR computational package. The angular momentum expansion of the Green function was truncated at ℓmax = 3 with a k-mesh of 18 × 18 × 1 points. The energy integrations were performed, including a Fermi-Dirac smearing of 502.78 K, and the Perdew–Wang 91 generalized gradient approximation was employed53. The Heisenberg exchange interactions and DM vectors were extracted using the infinitesimal rotation method55,56 with a finer k-mesh of 80 × 80 × 1 from the DLM57 state. The DLM theory is built upon the identification of robust local magnetic moments at certain sites {n}, with orientations prescribed by unit vectors \(\{{\hat{{\bf{e}}}}_{n}\}\) under a specified probability distribution. The approach is based on the assumption of a time-scale separation between the evolution of \(\{{\hat{{\bf{e}}}}_{n}\}\) in comparison to a rapidly adapting underlying electronic structure57. Following this adiabatic approximation, DFT calculations constrained to different magnetic configurations \(\{{\hat{{\bf{e}}}}_{n}\}\) can be performed to describe different states of magnetic order arising at finite temperatures. We extracted the magnetic interactions from fully disordered local moments whose orientations average to zero (\(\langle {\hat{{\bf{e}}}}_{n}\rangle =0\)), which is the high-temperature paramagnetic state. Adopting the DLM state as the reference state instead of the usual ferromagnetic or simple collinear AFM states has the benefit of not biasing the obtained magnetic interactions by choosing a magnetic reference state that might be quite different from the actual ground state.

Magnetic interactions and atomistic spin dynamics

We consider a classical extended Heisenberg Hamiltonian including Heisenberg exchange coupling (J), DMI (D), the magnetic anisotropy energy (K), and Zeeman term (B). All parameters were obtained from first-principles calculations. The energy of the reference unit cell in this spin-lattice model reads as follows (∣S∣ = 1):

$$\begin{array}{l}E=-\sum _{n\mu \nu }{J}_{0\mu ,n\nu }\,{{\bf{S}}}_{0\mu }\cdot {{\bf{S}}}_{n\nu }-\sum _{n\mu \nu }{{\bf{D}}}_{0\mu ,n\nu }\cdot \left({{\bf{S}}}_{0\mu }\times {{\bf{S}}}_{n\nu }\right)-\sum _{\mu }{\bf{B}}\cdot {{\bf{S}}}_{0\mu }\\\qquad\,-\sum _{\mu }{K}_{\mu }{\left({S}_{0\mu }^{z}\right)}^{2}\end{array}\ .$$
(1)

where n labels unit cells with 0 being the one including the origin. μ and ν label different magnetic sites within a unit cell. The magnetic properties can be characterized by analysing the Fourier-transformed magnetic interactions, which in reciprocal space give access to the magnetic ground state and the energy of spin spiral magnetic states:

$${J}_{\mu \nu }({\bf{q}})=\sum _{n}{J}_{0\mu ,n\nu }{e}^{-{\rm{i}}{\bf{q}}\cdot ({{\bf{R}}}_{0n}+{\bf{R}}\mu \nu )}$$
(2)
$${D}_{\mu \nu ,\alpha }({\bf{q}})=\sum _{n}{D}_{0\mu ,n\nu ,\alpha }{e}^{-{\rm{i}}{\bf{q}}\cdot ({{\bf{R}}}_{0n}+{\bf{R}}\mu \nu )}$$
(3)

where α = {x, y, z}, R0n is a vector connecting unit cell 0 and n, while Rμν is a vector connecting atoms μ and ν in the same unit cell. The contribution to the energy from the two-site interactions can be expressed in reciprocal space with the following quantity,

$${\mathcal{J}}({\bf{q}})=\left(\begin{array}{lll}J({\bf{q}})&{D}_{z}({\bf{q}})&-{D}_{y}({\bf{q}})\\ -{D}_{z}({\bf{q}})&J({\bf{q}})&{D}_{x}({\bf{q}})\\ {D}_{y}({\bf{q}})&-{D}_{x}({\bf{q}})&J({\bf{q}})\end{array}\right),$$
(4)

which is a 9 × 9 matrix (three magnetic sublattices × three spatial dimensions). The eigenvectors of \({\mathcal{J}}({\bf{q}})\) provide a set of independent possible magnetic states, and the eigenvalues are their respective energies. There is an important caveat of the energy analysis based on the eigenvalues of \({\mathcal{J}}({\bf{q}})\), which is essentially the Luttinger-Tisza approach: it only rigorously applies to a Bravais lattice, and not to a multi-sublattice structure such as the breathing kagome structure (see in ref. 58). This is because the eigenvectors that are obtained in the multi-sublattice case cannot in general be used to form magnetic states that satisfy the condition that the spins on all sites must have unit length, while this is possible and trivial for a Bravais lattice. Our experience with the \({\mathcal{J}}({\bf{q}})\) approach is that it is a good way to obtain a quick estimate of the relevant magnetic length scale, but only atomistic spin dynamics can provide the correct length scale and magnetic state subject to the unit spin constraint in the general case.

A definite answer regarding complex magnetic states is afforded through atomistic spin dynamic simulations using the Landau-Lifshitz-equation (LLG) as implemented in the Spirit code and performed using the Hamiltonian from Eq. (1). We employed periodic boundary conditions to model the extended two-dimensional system and used the simulated annealing method to navigate the intricate energy landscape and arrive at the ground state. We started from a random spin state at 1000 K which we let equilibrate, then cooled the system in steps by reducing the temperature to half of its previous value and equilibrating again, until we reached below 1 K. During the annealing process, a 1 T external field was applied along the z-axis to facilitate the formation of a single domain. The obtained magnetic structures from the annealing process were then reoptimized without the Zeeman field and found to be stable. The spin structure factor is defined by \({\bf{S}}({\bf{q}})=\frac{1}{{N}_{c}}{\sum }_{i}{{\bf{S}}}_{i}\,{e}^{-{\rm{i}}{\bf{q}}\cdot {{\bf{R}}}_{i}}\) with Nc the numbers of spins in the simulation cell, and Si is the magnetic moment orientation at position Ri.

There are two simple commensurate q-vectors close to q* = (0.195, 0.195), which is where the minimum of \({\mathcal{J}}({\bf{q}})\) is found—these are q< = (1/6, 1/6) = K/2 and q> = (1/5, 1/5). This motivated our choice for the smaller simulation cell size of 30 × 30, which can accommodate a single-q state of either q< or q>, and led to the triple-Q state formed by the star of q<. For the larger simulation cell size of 120 × 120, we could also accommodate q> = 24/120 = 1/5 = 0.2 ≈ 0.195. However, despite making several attempts using this larger simulation cell size with varying simulated annealing protocols, we always obtained a more or less disordered multi-domain state with the dominant periodicity characterized by the star of q<, as shown in Fig. 4b. Our interpretation is that the single-Q state is indeed energetically unfavorable when compared to the triple-Q state, and that the glassiness reflects both the multi-domain triple-Q state and the small energy cost for introducing spin deviations to the triple-Q states of the type belonging to the almost flat band of eigenvalues of \({\mathcal{J}}({\bf{q}})\) seen in Fig. 2f. This also explains the strong thermal fluctuations that we found in our atomistic spin dynamics, exemplified in Fig. S4, and confirms the need for a careful simulated annealing protocol.

Topological charge

We know that magnetic skyrmions are characterized by an integer topological number corresponding to the total solid angle \({n}_{{\rm{sk}}}=\frac{1}{4\pi }{\sum }_{i,j,k}{\Omega }_{ijk}\), which can be related to the topological Hall effect. The solid angle for the three spins Ωijk is related with the scalar chirality χijk = Si ⋅ (Sj × Sk) via \({\Omega }_{ijk}=2{\tan }^{-1}(\frac{{\chi }_{ijk}}{| {{\bf{S}}}_{i}| | {{\bf{S}}}_{j}| | {{\bf{S}}}_{k}| +{\sum }_{{\rm{cyclic}}}{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}| {{\bf{S}}}_{k}| })\)59. We use the topological charge density defined by summing the solid angle for all triangles neighboring a given spin to characterize the degree of noncoplanarity of the triple-Q state on the magnetic breathing kagome lattice. In Supplementary Fig. S4b, we provide a schematic representation of the calculation method. We use blue, green, and orange to represent spins in three sublattices. The colored rectangles (blue, green, and orange) indicate the cells used for calculating the solid angle subtended by the central spin (matching the shadow’s color) and its surrounding spins. The cell with colored shading encompasses six triangles, numbered in counterclockwise order. We show the calculated topological charge density of each cell in Supplementary Fig. S4c, the value of each cell is calculated by \(n=\frac{1}{4\pi }{\sum }_{i,j,k}{\Omega }_{ijk}\), where i, j, and k represent the three spins within each of the six triangles. We note that in Supplementary Fig. S4c, there are dark blue and dark red points with values greater than 1, indicating that these cells exhibit very strong noncoplanarity. The topological charge is defined as the sum over the solid angles assigned to all the sites in the 30 × 30 simulation lattice, and we obtain an integer value of −300. However, this value can change depending on the details of the magnetic configuration, which we believe signals a breakdown in the lattice version of the topological charge calculation; this concept is mathematically well-defined only for continuum vector fields.