Nanometric skyrmion lattice from anisotropic exchange interactions in a centrosymmetric host

Skyrmion formation in centrosymmetric magnets without Dzyaloshinskii–Moriya interactions was originally predicted from unbiased numerical techniques. However, no attempt has yet been made, by comparison to a real material, to determine the salient interaction terms and model parameters driving spin-vortex formation. We identify a Hamiltonian with anisotropic exchange couplings, local ion anisotropy, and four-spin interactions, which is generally applicable to this class of compounds. In the representative system Gd3Ru4Al12, anisotropic exchange drives a fragile balance between helical, skyrmion lattice (SkL), and transverse conical (cycloidal) orders. The model is severely constrained by the experimentally observed collapse of the SkL with a small in-plane magnetic field. For the zero-field helical state, we further anticipate that spins can be easily rotated out of the spiral plane by a tilted magnetic field or applied current.


Introduction
The modern push to realize complex magnetism with strong coupling to the electronic degrees of freedom is underpinned, to a large extent, by ever more powerful computational techniques. Among these are unbiased numerical simulations, such as Monte-Carlo methods, for magnetic ground states, excitation spectra, and spin dynamics starting from an (effective) spin-spin Hamiltonian. Monte-Carlo studies originally predicted that centrosymmetric triangular lattice magnets, be they insulators with (superexchange) interactions [1,2] or metals with effective Ruderman-Kittel-Kasuya-Yosida (RKKY) couplings [3][4][5], can host nanometer-sized magnetic spin vortices with net topological charge (i.e. skyrmions), and a multitude of even more complex states [3,6,7]. The theoretical work has introduced a new playground for the application of concepts from topology in condensed matter, but only with the observation of skyrmion lattices in the hexagonal intermetallics Gd 2 PdSi 3 [8] and Gd 3 Ru 4 Al 12 [9] has a more systematic comparison of theory and experiment become possible.
In this article, we demonstrate via a combination of theory and experiments that moderately weak anisotropic exchange from dipolar and spin-orbit interactions is sufficient to realize the SkL state in centrosymmetric magnets. Anisotropic exchange hence represents a new, third avenue toward the SkL in bulk crystals with inversion center, beyond frustrated exchange and higher-order RKKY interactions. Next to the present focus on intermetallic systems, this notion may have implications also for spin-vortex formation in insulating materials, e.g. on the magnetic diamond lattice of MnSc 2 S 4 [10,11]. In noncentrosymmetric magnets and at interfaces, spin-orbit coupling (SOC) provides the crucial underpinning of spin-spiral and skyrmion formation via the Dzyaloshinskii-Moriya interaction [12][13][14][15]. Likewise, we here show that SOC can play an important role in centrosymmetric magnets with SkL, as a driver of anisotropic exchange. Our work on weakly anisotropic exchange also creates a conceptual link (c) The distorted Kagome motif of rare earth sites is approximated as a triangular lattice of trimer plaquettes (magenta triangles). Three possible directions, equivalent by symmetry, of the magnetic ordering vector q ν (ν = 1, 2, 3) are indicated in the figure. Specifically for q 1 , we illustrate the anisotropy of parallel and perpendicular magnetic exchange interactions by black arrows. Data points in (b) adopted from reference [9].
between SkL formation and the movement to realize new quantum-disordered states in frustrated magnets via strongly bond-anisotropic exchange [16,17].
The focus of our study is Gd 3 Ru 4 Al 12 , a representative of the new class of centrosymmetric skyrmion hosts with coupled local moments and itinerant electrons. The Gd 3+ magnetic moments in this structure are arranged in quasi-layered Kagome nets, which are distorted by alternate stretching and compression of Kagome bond distances [18][19][20]-termed breathing Kagome lattice [figure 1(a)]. The magnetic phase diagram for field along the crystallographic c-axis (perpendicular to the Kagome plane) harbors five distinct regimes [9,21]: helical order (HE), a transverse conical state (TC), fan-like order (F), the SkL, and the as-yet uncharacterized pocket labeled by the roman numeral V [figure 1(b)]. In our analysis, a key guiding factor for the numerical simulation is the stability range of the SkL as compared to TC.

Model Hamiltonian
Motivated by the structural feature of rare earth triangles, and by previous considerations of the specific heat [22], we approximately treat Gd 3 Ru 4 Al 12 as a two-dimensional triangular lattice of strongly coupled superspin trimers [small triangles in the projection plane of figure 1(a)]. Figure 1(b) shows that the zero-field magnetic ground state of hexagonal Gd 3 Ru 4 Al 12 is a spiral, where the magnetic modulation vector q is aligned within the basal plane. In fact, three q ν with ν = 1, 2, 3 must be degenerate by symmetry, and previous work has shown that the q ν point along the a * axis and equivalent directions [9,21]. Backed by the knowledge that the directorsq ν are identical in all phases experimentally studied so far, we choose a minimal Hamiltonian consistent with the six-fold symmetry of the lattice and suitable for discussion of the magnetic phase transitions [23]. Only the dominant Fourier components S α q ν of the local spin S i at lattice site i are carried over [24], so that with Here, J is the energy scale for the dominant two-spin interaction, while (K/N), A, and h are dimensionless parameters corresponding to the strength of four-spin interaction, local ion anisotropy, and external magnetic field, respectively (normalized by J in all three cases). N is the number of spins in the system, and the direction of the applied magnetic field is labeled asĥ. The matrices Γ αβ q ν are also dimensionless and characterize the anisotropic exchange interactions. For example for ν = 1, is a vector component of the modulated magnetic moment at q ν with ν = 1, 2. The simulated annealing calculation is carried out [28] for parameters I ani = 0.01, A = −0.007 at low temperature T = 0.01 · J. Inset of (a), the three-fold degeneracy of the q ν in the hexagonal plane. Colored bands at the top define parameter regimes for phases HE, transverse conical (TC), SkL, F (fan-like), and field-aligned ferromagnet. In (c-e) magnetic textures are shown by arrows (magnetic moment in the hexagonal plane) and the color map, where yellow and blue correspond to positive and negative c-component of the moment, respectively. In (f, g) we depict both a three-dimensional sketch of the orders in phases HE and TC, and an illustration in terms of spiral planes and cones (lower side).
we have Γ q 1 = diag (I 0 − I ani , I 0 + I ani , I z ); the other Γ q ν are obtained from Γ q 1 by a simple rotation operation of ±120 • . In contrast to models of magnetic interactions in real space [25], the present (momentum-space) approach for the two-spin and four-spin terms is not sensitive to the choice of boundary conditions. For simplicity of the numerical treatment, the calculations were carried out for q ν rotated by 90 • from those in figure 1(c), i.e. for q ν along the triangular lattice's bond direction.
Anisotropic interactions of the I ani -type can be caused by relativistic SOC [23], but we note that dipolar interactions at the nearest neighbor level are also included in this term. In Discussion, we comment on the prospect of fully ab initio calculations of I ani in this class of compounds. More broadly speaking, Hamiltonians analogous to equation (1) are believed to be a good starting point as well for the description of phase transitions in centrosymmetric metallic magnets of different symmetries, such as tetragonal or cubic, if the q ν are strongly pinned to preferred crystal axes. Moreover, the formalism can be adapted to the description of non-centrosymmetric materials [26] by inclusion of antisymmetric terms Γ αβ q ν = Γ βα q ν , although this implies a larger number of adjustable parameters. When allowing for moderate (%-scale) changes of the modulus q ν = |q ν | at phase boundaries-including transitions from orders incommensurate with the underlying crystal lattice, to commensurate ones-the energy landscape is modified but weakly due to the typically very broad magnetic susceptibility (Lindhard function) χ(q) in metallic magnets with large, localized magnetic moments [27].
Without loss of generality, I z = I 0 = 1 is used in the following, as small differences between I 0 and I z can be absorbed into A, where large (small) I z corresponds to positive (negative) A. Previous magnetization measurements indicated mild easy-plane anisotropy A < 0 of local ions in Gd 3 Ru 4 Al 12 [9,22]. As A < 0 favors the TC state already in zero magnetic field, I ani > 0 is crucial to obtain phase HE with spiral vectors q ν //a * at h = 0. In other words, weaker interactions parallel to the propagation vector q ν are necessary to obtain phase HE in the present compound [ figure 1(c)].
Detailed supporting calculations show that excess I ani strongly favors multi-q phases, promoting the SkL over the TC state at all temperatures and eventually causing an instability of HE toward a multi-q state [28]. To describe Gd 3 Ru 4 Al 12 , we focus on 0 < I ani < 0.1 and set K = 0, the latter choice being revisited below. Within this range of parameters, the phase boundaries atĥ / / c are rather robust. For example, the simulations demonstrate that in absence of single ion anisotropy (A = 0), the overall properties of the c-axis phase diagram are nearly unchanged as compared to (I ani , A) = (0.01, −0.007). As we will show now, a good match of experiment and theory is obtained for I ani = 0.01 and A = −0.007, especially when aiming to also describe the case ofĥ tilted away from the c-axis. Figure 2(a) shows the components S α (q ν ) of the calculated magnetic structure factor, which is proportional to the squared Fourier component of the magnetic moment m α (q ν ) 2 . We use S α (q ν ) and m α (q ν ) 2 interchangeably in the following. Here, the indices α = , ⊥, z relate to a Cartesian frame of reference with unit vectors e ν =q ν , e zν =ẑ, and e ⊥ν perpendicular to both. For convenience, we also define m 2 xy = m 2 + m 2 ⊥ . For the chosen model parameters, the zero-field ground state corresponds to a HE with q = q 1 [figures 2(d) and (g)] with very weak admixtures of fan-like components at q 2 , q 3 [28]. Whenĥ is applied parallel to the c-axis, the system transitions into a TC state (only m xy = 0 for a single q ν ), then to the SkL (m xy = 0, m z = 0 with equal weight for all three q ν ), back into the TC phase, and finally into a high-field order which can be either transverse conical or fan-like, depending on the exchange parameters [28]. Sharp changes in S α q ν indicate first order phase transitions, which bound the SkL state and also characterize the spin-flop transition between HE and TC. Experimentally, these transitions were found to be accompanied by strong hysteresis. The strong history dependence of the experimental observations in phase TC, cf figure 1(b) and reference [9], in fact represents a challenge for the comparison to the theoretical result [figures 2(a) and (b)]. We discuss agreement of the two phase diagrams in more detail in the Supplementary Information https://stacks.iop.org/NJP/23/023039/mmedia [28].

Experiment: destruction of SkL with in-plane magnetic field
The general characteristics of the calculated phase diagram forĥ // c are unchanged for a rather broad range of Hamiltonian parameters as long as I ani > 0 > A and I ani + A 0 [28]. To further constrain the two free parameters A and I ani , we experimentally studied the effect of tilting the magnetic field away from c by an angle θ, c being the direction of six-fold symmetry ( figure 3). We hence establish that in Gd 3 Ru 4 Al 12 , the critical angle θ c of the SkL phase is only 10-15 • . In the experiment, we use the Hall conductivity σ xy as a sensor for the SkL phase. This observable is the transverse element of the conductivity tensor σ αβ relating electric field E α and resulting charge current via J α = σ αβ E β . The Hall conductivity acquires an additional contribution in the SkL phase, the topological Hall effect (THE), which effectively measures the winding number of a single skyrmion [29][30][31]. In our measurements of angle-dependent σ xy , the THE emerges as a bell-shaped anomaly on top of a smooth background.
Outside the SkL phase, the curves of σ xy (θ) collapse nicely onto a cosine-shaped profile. This behavior is consistent with the understanding that the oscillating background, on top of which the topological Hall signal σ T xy ≈ 250 Ω −1 cm −1 of the SkL develops, is due to spin-orbit coupling and the Karplus-Luttinger type anomalous Hall effect σ KL xy [32,33]. For σ KL xy , temperature dependence is not expected as long as the conduction electron's spin polarization remains unchanged [34]. Moreover, σ KL xy is proportional to the c-component of the net magnetization, consistent with the cosine-law observed here [32].
The experimental σ xy data is reported as a function of the applied magnetic field μ 0 H in units of Tesla, where μ 0 is the vacuum permeability. We corrected for the demagnetization effect by determining the internal field μ 0 H int = μ 0 (H − NM), where M is the magnetization per unit volume, measured independently, and 0 N 1 is the demagnetization factor calculated in elliptical approximation [28]. The present experimental observations contrast previous findings on a far more stable SkL with larger θ c ∼ 45 • in the related hexagonal compound Gd 2 PdSi 3 , with a triangular net of rare earth moments [8]. In Gd 3 Ru 4 Al 12 , the delicate interplay between TC and SkL for H//c, and also the low value of θ c , are both consequences of a fine balance between the parameters A and I ani .

Modeling the phase diagram in tilted field
Let the in-plane component ofĥ be aligned along q 1 . We use the simulated annealing framework to first discuss the character of magnetic order above the critical angle θ c . Figures 4(a) and (b) show how only a single modulation direction, q 1 , survives under these conditions, with spins arranged in a conical fashion around the axisn. In the simulation,n//ĥ in this regime, so that the conical axis rotates smoothly as θ is changed. The resulting dependences m 2 ∼ cos 2 θ and m 2 z ∼ sin 2 θ are indicated by black solid lines in figure 4(a).
We are now in a position to explore the parameter range of I ani and A suitable to describe the experimental situation in Gd 3 Ru 4 Al 12 , when using the trimer approximation (figure 5). The calculations demonstrate that for each I ani > 0 there is a critical A c which separates the stability regime of two zero-field ground states: for A > A c , a helical spiral is the preferred spin configuration, but it succumbs to the transverse conical state when A < A c . Figure 5 shows the A c values as dashed vertical lines for I ani = 0.01, 0.02, 0.03. The key point here is that only for small I ani = 0.01-0.02, the critical angle θ c of the SkL can be continuously tuned to zero before the collapse of the zero-field HE state. Following Fig. 5, we conclude that the experimental situation in Gd 3 Ru 4 Al 12 is well described using the trimer approximation and 0 < I ani < 0.02 with correspondingly small, finely tuned A.
We return to a question which was postponed earlier: what is the role of the four-spin term K in equation (1), or rather: could K = 0 explain the experiments for Gd 3 Ru 4 Al 12 ? K is well known to strongly favor multi-q order such as the SkL [4], and it comes as no surprise that K > 0 enhances the critical angle θ c for a variety of combinations of (A, I ani ). For example when I ani = 0.01 and K = 0.03, the SkL is stable even whenĥ is fully aligned along the in-plane direction x. As this is inconsistent with the observed small critical angle, we maintain that K = 0 is suitable for the present breathing Kagome compound.

Discussion
Many previous numerical studies have relied on Heisenberg Hamiltonians which yield, in a certain parameter range, a SkL phase when aided by thermal fluctuations or easy-axis anisotropy A > 0 [1,2,35]. These models are not well suited to the present material, where the SkL is realized despite A < 0. Likewise, RKKY Hamiltonians with biquadratic interactions (K = 0) [3,4] are inconsistent, in our simulations, with the very small critical angle θ c observed experimentally for Gd 3 Ru 4 Al 12 . Hence, a new approach, emphasizing the role of anisotropic exchange I ani > 0 together with easy-plane anisotropy of rare earth ions, is used to model the stability ranges of helical, transverse conical, and SkL phases on the triangular lattice of superspin trimers, corresponding to the trimerized limit of the breathing Kagome network in Gd 3 Ru 4 Al 12 [21,22]. SOC [23] and next-neighbor dipolar interactions both contribute to I ani . Note that easy-plane single-ion anisotropy was previously thought to be detrimental to skyrmion formation in magnets with inversion center [2,35]. However, A < 0 can stabilize Neel skyrmions in noncentrosymmetric polar materials [36,37].
The prediction of I ani for itinerant electron systems through fully ab initio electronic structure calculations remains an open challenge in this field. I ani is expected to be strongly dependent on the shape of Fermi surface sheets, as well as on the atomic spin-orbit coupling parameter λ SOC of the orbitals constituting the conduction bands.
It is worth reiterating that the conical axisn in Gd 3 Ru 4 Al 12 may be easily tilted away from being parallel to q (figures 4 and 5). This experimental condition has potential applications in the field of reading and writing skyrmions by very weak external perturbations. Our study also indicates that depinning of the SkL using an applied DC electrical current should be achievable in this material class.
A recent focus of research on spiral magnets is the detection of the emergent electric field resulting from spin dynamics driven by electrical currents, i.e. the realization of an emergent (or quantum) inductor [38]. In fact, the emergent inductance signal was recently observed by some of us in Gd 3 Ru 4 Al 12 [39], exploring the zero-field HE phase-discussed above-as a prototypical spiral order of large local moments coupled to the Fermi sea. The present model indicates that the spiral plane in Gd 3 Ru 4 Al 12 is easily depinned and rotated, enabled by a balance of anisotropic exchange and local-ion anisotropy. This property facilitates spin dynamics, helping to make the emergent electric field observable in transport experiments. Beyond shedding light on the issue of skyrmion formation, we hope that our work can guide the ongoing search for room-temperature material hosts of the emergent inductance effect.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.