Role of higher-order exchange interactions for skyrmion stability

Transition-metal interfaces and multilayers are a promising class of systems to realize nanometer-sized, stable magnetic skyrmions for future spintronic devices. For room temperature applications, it is crucial to understand the interactions which control the stability of isolated skyrmions. Typically, skyrmion properties are explained by the interplay of pair-wise exchange interactions, the Dzyaloshinskii-Moriya interaction and the magnetocrystalline anisotropy energy. Here, we demonstrate that higher-order exchange interactions – which have so far been neglected – can play a key role for the stability of skyrmions. We use an atomistic spin model parametrized from first-principles and compare three different ultrathin film systems. We consider all fourth-order exchange interactions and show that, in particular, the four-site four spin interaction has a large effect on the energy barrier preventing skyrmion and antiskyrmion collapse into the ferromagnetic state. Our work opens perspectives to stabilize topological spin structures even in the absence of Dzyaloshinskii-Moriya interaction.

M agnetic skyrmions-localized spin structures with a topological charge 1 -have raised high hopes for future magnetic memory and logic devices due to their nanoscale dimensions, stability, and ultra-low energy-driven motion [2][3][4][5][6] . Skyrmion lattices have been first observed in bulk magnets with a broken inversion symmetry in their crystal structure 7,8 . The discovery of a skyrmion lattice in a single atomic layer of Fe on the Ir(111) surface 9 has opened the door to a new class of systems: transition-metal interfaces and multilayers. Due to the possibility of varying film composition and structure, these systems allow to modify magnetic interactions and thereby the properties of skyrmions. At such transition-metal interfaces, individual magnetic skyrmions with diameters ranging from a few 100 nanometers down to a few nanometers have been realized as a metastable state in the field-polarized ferromagnetic background 8,10-16 as needed for applications.
A key challenge of skyrmion based data processing and storage technology is the robustness of information carriers, i.e., stability of the skyrmionic bits, against random thermal fluctuations at operating temperatures. At finite temperature, the magnetic moments of skyrmions are coupled to the environment, which induces fluctuations. Over time, a rare energy fluctuation can grow in excess of the barrier height and can prompt the skyrmion to overcome the barrier and collapse to the ferromagnetic background leading to a loss of topological charge. Therefore, an accurate assessment of barrier height is essential to determine the stability of skyrmions. To achieve data reading and writing capabilities of skyrmionic bits with high efficiency, a control over the barrier height is also necessary.
The existence of chiral magnetic skyrmions 17 is ascribed to a competition of the Heisenberg pair-wise exchange interaction, the Dzyaloshinskii-Moriya interaction (DMI) 18,19 , the magnetocrystalline anisotropy, and the dipole-dipole interactions. A prerequisite of DMI-which provides a unique rotational sense to skyrmions-is the concerted action of spin-orbit coupling and broken inversion symmetry, which can be achieved at interfaces of transition-metals 20 . The DMI further stabilizes metastable isolated skyrmions against annihilation into the ferromagnetic background 17,21 . Often the exchange interactions are treated in a micromagnetic or effective nearest-neighbor approximation. However, the exchange interactions are long-range in itinerant magnets such as 3d transition-metals. This can lead to a competition between exchange interactions from different shells of atoms resulting in an enhanced skyrmion stability 16,22 even in the absence of DMI 23 .
The itinerant character of 3d transition-metals limits the applicability of the Heisenberg model to describe their magnetic properties. Based on the spin-1/2 Hubbard model, it has been shown that the higher-order exchange interactions (HOI) beyond pair-wise Heisenberg exchange can arise such as the two-site four spin (biquadratic) or the four-site four spin interaction 24,25 . Such higher-order terms can lead to intriguing magnetic ground states due to a superposition of spin spiralsso-called multi-Q states-which have been predicted based on first-principles calculations 26 . The interplay of the four-site four spin interaction and DMI is the origin of the nanoskyrmion lattice of the Fe monolayer on Ir(111) 9 and the effect of the biquadratic interaction on skyrmion lattice formation has been studied systematically 27 . It has been further demonstrated that the HOI can compete with the DMI and stabilize novel magnetic ground states 28 . Based on a multi-band Hubbard model, a threesite four spin interaction has recently been proposed for systems with a spin beyond S = 1/2 in addition to the biquadratic and the four-site four spin interaction in fourth-order perturbation theory of the hopping parameter t with respect to the Coulomb energy U 29 . This term has been attributed to stabilize a double-Q or so-called up-up-down-down (uudd) state in an Fe monolayer on Rh(111) 30 . Despite the compelling experimental evidence of the relevance of HOI 9,28,30-32 , they have been neglected so far in the theoretical description of the properties of isolated magnetic skyrmions at transition-metal interfaces.
Here, we reveal the intriguing role played by the HOI for the stability of topologically non-trivial spin structures such as skyrmions and antiskyrmions at transition-metal interfaces. We use spin dynamics simulations based on an atomistic spin model with all parameters calculated via density functional theory (DFT). The energy barrier, preventing skyrmions and antiskyrmions to collapse into the ferromagnetic state, is obtained using the geodesic nudged elastic band (GNEB) method 33 . We consider three ultrathin film systems: (i) fcc-Pd/Fe bilayer on Rh(111) for which sub-10 nm skyrmions have been predicted at low magnetic field 34 , (ii) fcc-Pd/Fe bilayer on Ir(111), the most intensively studied ultrathin film system that hosts isolated skyrmions 10,11,22,[35][36][37][38][39][40][41][42] and an hcp-Fe/Rh bilayer on Re(0001) with an in-plane easy magnetization axis 43 .
Upon including the HOI, the stability of skyrmions and antiskyrmions in all of these films is greatly modified. Surprisingly, the effect of the biquadratic and the three-site four spin interaction concerning the energy barrier is to a good approximation already captured in the exchange constants obtained by mapping the DFT results to a spin model neglecting HOI. The four-site four spin interaction has a large effect on the saddle point and is responsible for the large change in energy barriers. We find a linear scaling of the barrier height with the four-site four spin interaction. The barrier is enhanced or reduced depending on its sign. Even small values of the four-site four spin interaction of 1-2 meV, typical for 3d transition metals, modify the energy barrier by 50-120 meV. This leads to a huge enhancement or reduction of the skyrmion or antiskyrmion lifetime. We further show that the HOI can stabilize topological spin structures in the absence of DMI.

Results
Atomistic spin model and DFT calculations. We describe the magnetic state of an ultrathin film by a set of classical magnetic moments {M i } localized on each atom site i of a hexagonal lattice and their dynamics is governed by the following Hamiltonian: where m i = M i /M i is a unit vector. The exchange constants (J ij ), the DMI vectors (D ij ), the magnetic moments (μ s ) and the uniaxial magnetocrystalline anisotropy energy (MAE) constant (K) were calculated based on DFT (see refs. 22,34,39,43 ). We neglect the dipole-dipole interaction since it is small in ultrathin films, which is of the order of 0.1 meV per atom, and it can be effectively included into the MAE 44,45 .
The last three terms are the biquadratic interaction (B ij ), the three-site four spin interaction (Y ijk ) and the four-site four spin interaction (K ijkl ), respectively. Since these terms arise from the fourth-order perturbation theory, we restrict ourselves to the nearest-neighbor approximation, i.e., up to the first term of these HOI. The corresponding constants are denoted as B 1 , Y 1 and K 1 . The evaluation scheme of the HOI on a hexagonal twodimensional (2D) lattice is illustrated in Fig. 1b.
First, we calculate the energy dispersion of homogeneous flat spin spirals via DFT without taking spin-orbit coupling (SOC) into account to determine the exchange constants. A spin spiral is characterized by a wave vector q in the 2D Brillouin zone (2DBZ, Fig. 1d) and the magnetic moments of an atom at lattice site R i is given by M i ¼ MðsinðqR i Þ; cosðqR i Þ; 0Þ with the size of the magnetic moment M. Fig. 1e shows the energy dispersion E(q) of spin spirals without SOC along two high symmetry directions obtained via DFT for an fcc-Pd/Fe bilayer on the Rh(111) surface 34 (Fig. 1a). At the high symmetry points of the 2DBZ, we find well-known magnetic states: the ferromagnetic (FM) state at the Γ point, the row-wise antiferromagnetic (AFM) state at the M point and the Néel state with angles of 120°between adjacent magnetic moments at the K point.
Clearly, the FM state is energetically lowest among these three states (Fig. 1e). Along both high symmetry directions, the 90°spin spirals ( Fig. 1c) are found at q ¼ ð1=2ÞΓM and at q ¼ ð3=4ÞΓK (Fig. 1d). The total energy of homogeneous spin spirals without SOC is fitted to functions obtained by expressing the Heisenberg model (first term of Eq. (1)) in reciprocal space to extract the exchange interaction parameters up to 11th nearest neighbors (Supplementary Table 1).
Since the functional form of the three-site four spin and the biquadratic interactions for homogeneous spin spirals resemble that of the first three exchange constants, we cannot separate the exchange and the higher-order constants by fits (see "Methods"). Therefore, we calculate the HOI constants from the energy difference between the spin spiral (single-Q) and multi-Q states without SOC and modify the exchange constants obtained from a fit to only the first term of Eq. (1).
The three selected multi-Q states are a superposition of spin spirals (neglecting SOC) corresponding to symmetry equivalent q vectors (Fig. 1c, d) and have a constant magnetic moment at every atomic site. Within the Heisenberg model of pair-wise interaction, the multi-Q and single-Q states are energetically degenerate. However, the HOI lift the degeneracy which provides a way to compute their strengths. In DFT calculations, all the interactions are implicitly included through the exchangecorrelation functional. Therefore, we can obtain the HOI constants from total energy calculations of multi-Q and single-Q states without SOC.
We consider two collinear states, the so-called uudd or doublerow wise antiferromagnetic states 46 and a three-dimensional noncollinear state, the so-called 3Q state 26 , to uniquely determine three higher-order exchange constants (for spin structures see insets of Fig. 1e). The biquadratic (B 1 ), the three-site four spin (Y 1 ) and the four-site four spin interaction (K 1 ) constants are computed from the energy differences between the multi-Q and the corresponding single-Q states (without SOC) by solving the equations 29 : The three multi-Q states are higher in energy compared to the corresponding spin spiral states without SOC for Pd/Fe/ Rh(111) (Fig. 1e) and far above the FM state (for energies see  (1)). The energies of the uudd state along ΓK (pink filled circle), the uudd state along the ΓM direction (blue filled circle) and the 3Q state (green filled circle) are denoted at the q value of the corresponding single-Q state. The spin structures of the uudd state along ΓK (pink), which is formed by a superposition of two 90°spin spirals at q = ± ð3=4ÞΓK, the uudd state along the ΓM direction (blue), which is formed by a superposition of two 90°spin spirals at q = ± ð1=2ÞΓM, and the 3Q state (green), which is a superposition of three spin spirals corresponding to three M points (green) in 2DBZ. Inset of e shows a zoom around the Γ point. Red circles denote DFT calculations without SOC and black circles calculations including SOC, i.e., the effect of DMI and MAE. The energetically lowest state is a spin spiral with a period of 4.8 nm.  Table 2). Nevertheless, they have a large effect on skyrmions in this film system as we show below.
The computed HOI constants modify the first three exchange constants, obtained from fits of the spin spiral energy dispersion neglecting HOI (see "Methods" for a derivation): where we denote the exchange parameters obtained from fits neglecting HOI as unprimed and the modified ones by considering the higher-order terms as primed. Note the special role played by the four-site four spin interaction which does not adjust any exchange parameter since its contribution to the energy dispersion of spin spirals is a constant value of −12K 1 independent of the spin spiral vector. Further the HOI only modify the first three exchange constants and, therefore, the other exchange constants used in atomistic spin dynamics simulations remain unchanged (Supplementary Table 1). The first three exchange and the higher-order exchange constants of Pd/Fe/Rh(111) are displayed in Table 1 along with those for Pd/Fe/Ir(111) and Fe/Rh/Re(0001) obtained by similar DFT calculations (see also Supplementary Table 2). We find that the Pd/Fe bilayer on Rh(111) and on Ir(111) behaves similarly in terms of exchange and higher-order exchange constants since Rh and Ir are isoelectronic 4d and 5d transition metals. In these two film systems, the signs of the nearest-neighbor exchange constant (J 1 ) and the second and third nearest neighbors are opposite which leads to exchange frustration 22,34 . The exchange interaction in Fe/Rh/Re(0001), in contrast, is dominated by the nearestneighbor exchange constant. Note that the sign of the biquadratic (B 1 ) and the four-site four spin constants (K 1 ) is negative in Fe/ Rh/Re(0001), while it is positive for the other two systems. As shown below, the sign of K 1 is essential for the skyrmion stability in these films.
SOC introduces two additional contributions: DMI and MAE (see "Methods" for details). The DMI lowers the energy of cylcoidal spin spirals with a clockwise rotational sense in the vicinity of the Γ point (see inset of Fig. 1e). The MAE shifts the energy of spin spirals by K/2 with respect to the FM state. The DMI constants and MAE are given in Supplementary Table 3. Note that, very recently, higher-order interactions arising from SOC have been proposed for transition metal systems 32,47,48 . However, experimental evidence for these interactions is missing and we exclude them from our current investigation.
Spin dynamics simulations. We use atomistic spin dynamics simulations (see "Methods") based on the Hamiltonian of Eq. (1) with all parameters obtained from DFT including DMI, MAE and total magnetic moments, as discussed in the previous section, to calculate the zero temperature phase diagram and the properties of isolated skyrmions and antiskyrmions in the three film systems. The results for Pd/Fe/Rh(111) are shown in Fig. 2  We first discuss the effect of the HOI on the zero temperature phase diagram (Figs. 2a, b). The FM state and the homogeneous spin spiral states remain almost unaffected by the higher-order terms. However, the skyrmion lattice loses a small amount of energy with respect to the homogeneous spin spirals which remains constant throughout the range of magnetic fields. This leads to an expansion of the spin spiral and FM phases at the expense of the skyrmion lattice phase, which is squeezed. Since isolated skyrmions can be stabilized in the FM (field-polarized) phase, it is of prime importance. The onset of the FM phase, characterized by the critical field B c , has shifted from 2.75 T to a lower value of 2.25 T due to the HOI. As expected from the magnetic interaction constants, a similar trend of the phase diagram is obtained for Pd/Fe/Ir(111) (Supplementary Fig. 1). However, since the higher-order constants are quite small for Fe/Rh/Re(0001) (cf. Table 1), the phase diagram is basically unchanged ( Supplementary Fig. 2).
In our spin dynamics simulations, we have created isolated skyrmions and antiskyrmions in the field-polarized background following the theoretical profile 49 and relaxed the spin structures with the full set of DFT parameters. The radius of skyrmions and antiskyrmions-defined as in ref. 49 increases with HOI for Pd/ Fe/Rh(111) on average by~35% and 30%, respectively (Fig. 2g). The skyrmion and antiskyrmion profiles at 5.2 T (see insets of Fig. 2g) can be fitted by the standard skyrmion profile. The antiskyrmion exhibits a steeper profile which reflects that it has a smaller radius than the skyrmion. Similar trends of the skyrmion and antiskyrmion radii are found for Pd/Fe/Ir(111) (Supplementary Fig. 1). Due to relatively small HOI, the skyrmion radii remain almost unchanged for low magnetic fields above B c for Fe/ Rh/Re(0001) (Supplementary Fig. 2 and Supplementary Note 1).
To study the stability of metastable isolated skyrmions and antiskyrmions, we calculate the minimum energy path (MEP) for the collapse of a single skyrmion or antiskyrmion into the FM background (see "Methods"). The point of maximum energy on this path, known as the saddle point, with respect to the initial state (skyrmion or antiskyrmion) is a measure of the barrier height. As seen in Fig. 2h, the HOI increase the energy barrier for skyrmion annihilation in Pd/Fe/Rh(111) by more than a factor of two at small magnetic fields above B c . For antiskyrmions, the barrier height is even increased by a factor of 5.
The energy barriers of skyrmions vary nonlinearly at small magnetic fields, and, thereafter, reduce almost linearly with increasing magnetic fields 22 . On an average, we notice an increase in barrier height of nearly 100 meV for skyrmions upon including the HOI. At low fields up to B − B c = 1.36 T, there is a transition from the normal radial collapse mechanism 22,34 without HOI to a chimera collapse mechanism 16,50 with HOI. Above B − B c = 1.36 T, the skyrmions merge into the FM background through the normal radial collapse without HOI, which remains unchanged after including HOI.
The barrier heights of antiskyrmions with HOI exhibit a similar variation with field as that of skyrmions. The energy barriers of antiskyrmions without HOI are extremely small (~20 meV) which implies that they are basically unstable even at cryogenic temperatures. However, after including HOI, the energy barriers become~100 meV at small fields (up to 1.5 T above B c ), which suggests that metastable antiskyrmions could be Table 1 Exchange constants including HOI. Exchange constants for i-th neighbor spins (J 0 i ), biquadratic (B 1 ), three-site four spin (Y 1 ) and four-site four spin interaction (K 1 ) constants for Pd/Fe/Rh(111), Pd/Fe/Ir(111) and Fe/Rh/Re (0001). The exchange constants are modified upon including the HOI according to Eqs. (5)- (7). This table shows only those exchange constants which are modified upon including the HOI. The full sets of exchange constants as well as the DMI, the MAE and the total magnetic moments used in the atomistic spin dynamics simulations are listed in Supplementary Tables 1 and 3 (see also Supplementary Data 1 and 2). All values are given in meV. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18473-x realized in experiments. The annihilation mechanism of antiskyrmions is via the radial collapse 22 , which is unaffected upon including HOI.
In Pd/Fe/Ir(111), the HOI increase the energy barriers of isolated skyrmions and antiskyrmions by similar values of~100 and~90 meV, respectively ( Supplementary Fig. 1). On the other hand, in Fe/ Rh/Re(0001), the stability of isolated skyrmions is reduced on an average by 70 meV upon including the HOI (Supplementary Fig. 2). This large barrier reduction shows that even small values of the higher-order constants (cf. Table 1) can have significant effects.
Analysis of collapse mechanisms. Now we focus on the question which of the HOI is responsible for the large changes of the energy barriers for skyrmion or antiskyrmion collapse. We consider both collapse mechanisms for skyrmions, i.e., the chimera collapse at low fields and the radial collapse at higher fields and the radial collapse mechanism for antiskyrmions (cf. Fig. 2h). In Fig. 3, the energy decomposition of three representative minimum energy paths are displayed at selected magnetic fields for Pd/Fe/Rh(111) with and without taking HOI into account (for the other two systems see Supplementary Figs. 3 and 4).
The total energy rises along the MEP as one moves from the initial (skyrmion) state to the saddle point and descend thereafter to the final (ferromagnetic) state (Fig. 3a). In the simulation neglecting the HOI, we find that the energy barrier is dominated by the energy contribution from the DMI, which favors the skyrmion state. Due to exchange frustration, there is also a small energy contribution to the barrier from the exchange energy. Naturally, the energy due to the Zeeman term and the magnetocrystalline anisotropy decrease in the ferromagnetic state.
Upon including the HOI (Fig. 3b), we find the large increase of the energy barrier as discussed in the previous section. In addition, the annihilation mechanism changes at this magnetic field from the radial collapse without HOI (Fig. 3a) to the chimera collapse mechanism with HOI (cf. Fig. 3j, k which show the spin structures in the vicinity of the saddle point of the two types of annihilation mechanisms). Interestingly, the chimera collapse mechanism has been previously discussed in ultrathin films with very strong exchange frustration 16,50,51 , which suggests that HOI act in a similar way. The energy decomposition shows that the DMI contribution is of similar magnitude at the saddle point of the path with and without HOI (Fig. 3c). However, one cannot compare the exchange interactions before and after the HOI are included in the simulations, since the higher-order terms modify the exchange constants according to Eqs. (5)-(7). Therefore, we add the contributions due to the exchange, the three-site four spin  interaction and the biquadratic terms which we denote as combined exchange. The comparison of Fig. 3a, b shows that the exchange and the combined exchange behave qualitatively quite similar along the path-which is also true for the other two collapse mechanisms (Fig. 3d, e and Fig. 3g, h). The absolute energy change from exchange to combined exchange at the saddle point is relatively small (Fig. 3c, f, i).
The four-site four spin interaction acts in a qualitatively different way compared to all other terms. For all considered paths (Fig. 3b, e, h), it gains in energy slowly as one approaches the saddle point, in the vicinity of the saddle point it becomes very steep, reaches a maximum at the saddle point and drops quickly thereafter. The energy contributions at the saddle point (Fig. 3c, f, i) show that it provides by far the largest difference between the simulations with and without HOI, irrespective of the collapse mechanism or the initial state. The energy contribution from the three-site four spin and the biquadratic interactions decreases along all collapse processes and the energy drop escalates after the saddle point (see insets of Fig. 3b, e, h). The difference in energy profile of the DMI and MAE is an associated effect of the HOI caused by the changes in relative spin angles during the collapse process. Fig. 3c, f, i show that the combined exchange can provide a tiny contribution to the energy barrier depending on the collapse mechanism, while the DMI and Zeeman terms assert only a little weight if not compensated by each other. Therefore, the four-site four spin interaction mainly controls the change of the barrier height.
The spin structure in the vicinity of the saddle point is shown in Figs. 3j, k for the chimera and the radial collapse of the isolated skyrmion including the effect of HOI. We see that the radial collapse of an isolated skyrmion is very similar to that found by neglecting HOI in ref. 22 . However, at the saddle point, there are four spins pointing towards each other while previously a three-spin structure was reported. The unusual saddle point including HOI is obtained throughout the studied field range and for annihilation of skyrmions in Pd/Fe/Ir(111). However, for the skyrmion collapse in Fe/Rh/Re(0001), a threespin structure at the saddle point similar to that of ref. 22 occurs. The chimera skyrmion collapse and the radial antiskyrmion collapse are similar to that found in simulations neglecting HOI 16,22,50 .
We have also performed atomistic spin simulations for Pd/Fe bilayers on Rh(111) and on Ir(111) without and with higherorder interactions of the MEP for the escape mechanism ( Supplementary Fig. 5) introduced previously by Bessarab et al. 38 . The saddle point along this MEP resembles a slightly deformed skyrmion ( Supplementary Fig. 5) and is distinctively different from the saddle point of the collapse mechanisms with large spin rotations on the atomic scale. Since the latter property is decisive for the large energy contribution of the four-site four spin interaction, there is almost no influence of HOI on the energy barrier of the escape mechanism.
Analysis of the four-site four spin interaction. To understand the prominent effect of the four-site four spin interaction on the energy barrier, we present its site-resolved energy at the saddle point with respect to the initial state (skyrmions or antiskyrmions) for the three MEPs of Fig. 3b, e, h in Fig. 4a-c. We notice that a group of only 14 spins around the core provide contributions to the four-site four spin interaction, while the surrounding spins do not add any significant value. This finding is independent of whether we consider the saddle point of the chimera collapse (Fig. 4a), the radial skyrmion collapse (Fig. 4b) or the radial collapse of the antiskyrmion (Fig. 4c). Similar observations are made for the other film systems ( Supplementary  Figs. 6-9). Therefore, the four-site four spin interaction at the saddle point exhibits a general behavior irrespective of the type of collapse mechanism or the initial spin configuration.
In order to explain this localized energy gain at the saddle point, we use a simplified model in which we consider only the site with the largest contribution at the origin. To evaluate the four-site four spin interaction, we need to consider at least the 6 nearest neighbors and the 6 next-nearest neighbors of the central site [cf. Fig. 4d-f]. To simplify the discussion, we slightly symmetrize the spin structure. For the radial skyrmion and antiskyrmion collapse, the 12 neighboring spins are nearly all inplane (Fig. 4e, f). We neglect any out-of-plane component as shown in Fig. 4h, i to calculate the contributions from the 12 diamonds for the four-site four spin interaction.
For the saddle point of the radial skyrmion collapse (Fig. 4h), we find three distinct types of diamonds which contribute to the four-site four spin interaction. We find a pair of diamonds with values +K 1 and −K 1 , two pairs of diamonds with values þ 1 2 K 1 and À 1 2 K 1 , which cancel out mutually. Out of the six remaining diamonds, there are three groups each containing two diamonds with values À ffiffi 3 p 2 K 1 , −K 1 and þ 1 2 K 1 , which results in a total energy at the saddle point of E ISk SP ¼ À2:73K 1 . For the saddle point of the radial antiskyrmion collapse (Fig. 4i), we identify three types of diamonds with the same magnitude as for the skyrmion saddle point (Fig. 4h). However, two uncompensated diamonds with þ 1 2 K 1 and À ffiffi 3 p 2 K 1 and two diamonds with − K 1 each lead to a total energy of E IASk SP ¼ À2:37K 1 .
The spin structure at the saddle point of the chimera collapse (Fig. 4d) is more complex. There are non-negligible out-ofplane components of the spins surrounding the central spin which we take into account in the symmetrization (Fig. 4g). As a consequence, we find six distinct types of diamonds (Fig. 4g). Similar to the other two saddle points, there is a mutual cancellation of many terms which leads to a total energy contribution of the four-site four spin interaction of E chimera SP ¼ À2:37K 1 . Note that the values of these three energies taking the exact spin structure at the saddle points are E ISk SP ¼ À2:0K 1 , E IASk SP ¼ À2:2K 1 and E chimera SP ¼ À2:14K 1 , which are very close to those obtained using the simplified spin structures.
We obtain similar values for Pd/Fe/Ir(111) at the saddle points corresponding to the skyrmion and antiskyrmion initial states and at the chimera saddle point, the value is only slightly different E chimera SP ¼ À2:58K 1 (Supplementary Figs. 6-8). For Fe/Rh/Re (0001), we find E ISk SP ¼ À ffiffi ffi 3 p K 1 (Supplementary Fig. 9). The energy contribution per site of the four-site four spin interaction for the ferromagnetic state or any flat spin spiral is −12K 1 , which is also relatively close to the energy in the skyrmion state (cf. Fig. 3b and Supplementary Figs. 3b, e and 4b). Therefore, we obtain an energy difference of E SP − E FM ≈ 10K 1 for the two symmetric central sites of the saddle point. The surrounding sites provide smaller contributions, however, they still scale with K 1 . In total, we find an energy contribution of the four-site four spin interaction to the energy barrier of roughly 40K 1 . Due to the linear dependence on K 1 , it is also clear that the sign of the four-site four spin interaction determines whether there is an energy gain (K 1 > 0) or loss (K 1 < 0) at the saddle point as observed for the two types of systems: Pd/Fe on Rh(111) and on Ir(111) vs. Fe/Rh/Re(0001) (cf. Table 1).

Discussion
Our simplified model states that the barrier height E SP − E ISk/IASk varies linearly with the magnitude of the four-site four spin constant K 1 . To verify this prediction, we have carried out spin dynamics simulations for three ultrathin film systems at a given magnetic field by changing only the four-site four spin interaction constant while leaving all other interactions the same. Note that the four-site four spin interaction does not affect the energy dispersion of spin spirals which is essential for the equilibrium properties of skyrmions and antiskyrmions. Figure 5 shows thatas expected from our model-the barrier heights for skyrmions and antiskyrmions exhibit a linear scaling with the four-site four spin constant. Only for Pd/Fe/Rh(111), we find a slight deviations from the linear dependence. The scaling constant α, defined as the ratio of the change in energy barrier to the change in four-site four spin constant, is the same for skyrmions and antiskyrmions consistent with our discussion of the energy contributions at the saddle point. We find a value of the scaling constant α of~40-60 depending on the system. Figure 3 implies that the HOI can stabilize skyrmions and antiskyrmions by themselves. To test this notion, we have performed spin dynamics simulations by completely switching off the DMI, i.e., setting the DMI from the DFT calculation to zero, while keeping all other magnetic interactions as before. As shown in Fig. 6a, we find stable skyrmions and antiskyrmions for Pd/Fe bilayers on Rh(111) with the same radius which is reduced to 1.4 nm just above B c compared to the case with DMI (cf. Fig. 2g). Due to vanishing DMI, clockwise and anticlockwise rotating skyrmions are degenerate. Large energy barriers of up to 90 meV at B c are obtained (Fig. 6b) due to the four-site four spin interaction, which are identical for skyrmions and antiskyrmions as the DMI is zero.
The decomposition of the energy contributions along the MEP (Fig. 6c), which is the same for skyrmions and antiskyrmions, shows that the sum of all exchange interactions, the biquadratic, and the three-site four spin interaction (combined exchange) results only in a minimal energy barrier. In contrast, the four-site four spin interaction provides the energy barrier and exhibits its characteristic peak-like curve along the MEP. Although the initial skyrmion or antiskyrmion state exhibits a reduced diameter, the energy barrier due to the four-site four spin interaction is almost the same as in the case with DMI (cf. Fig. 3e, h). The saddle point configuration of the skyrmion (Fig. 6d) and antiskyrmion (Fig. 6e) show the characteristic spin structure with four spins at angles of nearly 90°with respect to each other which was observed including DMI (cf. Fig. 3e, f).
As shown in Supplementary Fig. 10, we obtain similar results for Pd/Fe/Ir(111) upon setting the DMI to zero. This suggests that it is possible to stabilize both types of topological states with arbitrary rotational sense due to HOI at inversion symmetric transition-metal interfaces. Note that it was previously proposed that strong frustration of exchange interactions could stabilize skyrmions and antiskyrmions without DMI. However, a specific ratio between different pair-wise exchange interactions is required 23 .
The lifetime τ of skyrmions or antiskyrmions is given by the Arrhenius law τ = τ 0 exp(ΔE/k B T), where ΔE is the energy barrier and τ 0 is the prefactor. Typically, the lifetime is dominated by the energy barrier due to the exponential term. Recently, it has been reported that due to entropy the prefactor can vary drastically with external parameters, e.g., the magnetic field 52,53 . For bulk magnetic materials a change by 30 orders of magnitude 52 and for ultrathin films a variation of up to seven orders of magnitude have been found 53 . However, the effect depends on details of the magnetic interactions. For fcc-Pd/Fe bilayers on Ir(111) and Rh(111) 34 considered here, there is almost no change of the prefactor and the lifetime is governed by the energy barrier.
An enhancement of the barrier height, ΔE, by 100 meV, as observed for an isolated skyrmion in Pd/Fe/Rh(111), leads to an increase of skyrmion lifetime by orders of magnitude because of the exponential factor. For example, at a temperature of T = 10 K, at which the spin-polarized scanning tunneling microscopy experiments on such ultrathin films are typically performed 10,11 , we find an enhancement by 50 orders of magnitude, at 100 K, it is still a factor of about 10 5 , and even at room temperature, it is a factor of about 50. One can also discuss the effect of the HOI in terms of the temperature-dependent phase diagram of Pd/Fe/Rh (111) 34 . Without the HOI, skyrmions can be stable for an hour up to temperatures of 25 K 34 , which is increased to a temperature of about 50 K upon including HOI. For antiskyrmions, the change from a barrier of below 20 to~100 meV leads to an enhancement of their lifetime, which should allow their experimental discovery at least at cryogenic temperatures. We have demonstrated that higher-order interactions beyond pair-wise Heisenberg exchange can play a key role for the stability of skyrmions or antiskyrmions at transition-metal interfaces. While the biquadratic and three-site four spin interaction contribute in a similar fashion as pair-wise exchange interactions to the minimum energy path, we find a qualitative difference for the four-site four spin interaction. Due to the cyclic hopping on four sites it acts on the atomic scale. Therefore, it affects the saddle point of the collapse path strongly while its energy contribution to the initial (skyrmion) and final (ferromagnetic) state is rather similar, in particular, for relatively large skyrmions. This leads to a characteristic peak-like shape of its energy contribution along the collapse path which is otherwise only obtained from a concerted interplay of interactions. On the other hand, DMI provides only an energy difference between the skyrmion and ferromagnetic state. In this respect, the four-site four spin interactions plays a unique role among all considered interactions.
Depending on the sign of the four-site four spin interaction, the energy barrier preventing the collapse of a metastable topological spin structure can be greatly enhanced or reduced. Even for the small values of HOI typical for 3d transition metals, we find large changes of the energy barriers and therefore giant effects on the lifetime, which means that these interactions cannot be neglected. The energy barriers are so much enhanced due to HOI that isolated skyrmions and antiskyrmions can be stable in the absence of Dzyaloshinskii-Moriya interaction. Our study opens up another avenue to stabilize topological spin structures at transition-metal interfaces.

Methods
First-principles calculations. The computational details for calculating the exchange, the DMI, the MAE constants and the magnetic moments are shown in ref. 34 for Pd/Fe/Rh(111), in ref. 22 for Pd/Fe/Ir(111) and in ref. 43 for Fe/Rh/Re (0001). Here, we have evaluated the higher-order exchange constants for all three systems. The electronic structure was calculated using a spin-polarized DFT code based on the projected augmented wave (PAW) scheme 54 as implemented in the Vienna ab initio simulation package (VASP) 55 . It ranks among the best available DFT codes in terms of accuracy and efficiency 56 . We use the same structural parameters as mentioned in the above references. We have used two atomic overlayers on top of nine substrate layers to mimic the surfaces. To maintain consistency with spin spiral calculations, we have chosen local density approximation (LDA) for the exchange and correlation part of potential 57 . A high energy cut-off of 400 eV was used to precisely calculation the energy of the multi-Q states. The 2DBZ was sampled by a Monkhorst-Pack 58 mesh of 22 × 28 × 1 k-points for the uudd state in the ΓK direction, of 14 × 44 × 1 k-points for the uudd state in the ΓM direction and of 15 × 15 × 1 k-points for the 3Q state at the M point. The total energy calculations for the multi-Q states were performed without considering SOC and the convergence criteria were set to 10 −6 eV for all calculations.
Fitting function for HOI. The spin spiral is the exact solution of the classical Heisenberg model for a periodic lattice. The spin spiral, which is characterized by a wave vector q in the 2DBZ and the magnetic moments of an atom at lattice site R i , is given by 59 , where R q and I q are two vectors that span the xy-plane. They obey the following relation, where M is the magnitude of M i and without loss of generality, we set its norm to unity. The spins of a spin spiral rotate around the z-axis in the xy-plane as one moves from one lattice site to another in the direction of q. Using the above two equations, the scalar product of a pair of spins can be written as 59 ,