Magnetic structure of monatomic Fe chains on Re(0001): emergence of chiral multi-spin interactions

TL;DR

This study uses first-principles calculations to explore how the magnetic structure of monatomic Fe chains on Re(0001) evolves with chain length, revealing the emergence of chiral multi-spin interactions influencing spin spiral states.

Contribution

The paper introduces a classical spin model incorporating chiral four-spin interactions derived from ab initio calculations, explaining the observed magnetic behaviors.

Findings

Transition from antiferromagnetic to spin spiral states with increasing chain length

Opposite chirality of spin spirals predicted by different methods

Chiral four-spin interactions are essential to accurately model the magnetic structure

Abstract

We present results of first-principles calculations of the magnetic properties of Fe chains deposited on the Re(0001) surface. By increasing the length of the chain, a transition is found from an almost collinear antiferromagnetic state for a 5-atom-long chain to a spin spiral state with the rotational plane slightly tilted from the surface of the substrate for the 15-atom-long chain. It is shown that a classical spin model derived from the ab initio calculations containing only two-spin interactions supports opposite chirality of the spin spiral compared to a direct optimization of the spin configuration within the ab initio method. The differences between the results of the two methods can be understood by introducing chiral four-spin interactions in the spin model.

Figures

Tables

Table 1. Table 1: Components of the NN and NNN DM vectors, D → i , i + 1 subscript → 𝐷 𝑖 𝑖 1 \vec{D}_{i,i+1} and D → i , i + 2 subscript → 𝐷 𝑖 𝑖 2 \vec{D}_{i,i+2} , respectively, for the 15-atom-long chain, given in units of meV. The components of the DM vectors for i > 7 𝑖 7 i>7 can be obtained by the symmetry relations Eq. ( 14 ) and are illustrated in Fig.  3 .

	${\overrightarrow{D}}_{12}$	${\overrightarrow{D}}_{23}$	${\overrightarrow{D}}_{34}$	${\overrightarrow{D}}_{45}$	${\overrightarrow{D}}_{56}$	${\overrightarrow{D}}_{67}$	${\overrightarrow{D}}_{78}$
$x$	–0.60	0.64	0.23	0.42	0.18	–0.03	–0.03
$y$	–3.95	–0.71	–0.98	–1.24	–1.76	–1.62	–1.43
$z$	–1.24	2.96	1.61	1.80	1.34	1.48	1.49
	${\overrightarrow{D}}_{13}$	${\overrightarrow{D}}_{24}$	${\overrightarrow{D}}_{35}$	${\overrightarrow{D}}_{46}$	${\overrightarrow{D}}_{57}$	${\overrightarrow{D}}_{68}$	${\overrightarrow{D}}_{79}$
$x$	0.05	–0.25	–0.30	–0.04	–0.02	–0.07	0.00
$y$	–4.06	–4.47	–5.12	–5.13	–4.92	–4.96	–5.03
$z$	0.44	0.84	0.73	1.16	1.16	1.13	1.12

Table 2. Table 2: Chiral interaction vectors between pairs of spins obtained from the rotational scheme based on MFT calculations, 𝒟 r , i ​ j α superscript subscript 𝒟 r 𝑖 𝑗 𝛼 \mathcal{D}_{\text{r},ij}^{\alpha} , Eq. ( 18 ), and from the torque method, 𝒟 t , i ​ j α superscript subscript 𝒟 t 𝑖 𝑗 𝛼 \mathcal{D}_{\text{t},ij}^{\alpha} , defined in Eq. ( 41 ). Values are given in meV for NN and NNN interactions in the middle of the 15-atom-long chain.

$i$	$j$	$\alpha$	${𝒟}_{\text{r},ij}^{\alpha }$	${𝒟}_{\text{t},ij}^{\alpha }$
7	8	$x$	–0.243	–0.033
7	8	$y$	–3.107	–1.428
7	8	$z$	–2.306	1.491
7	9	$x$	–0.376	0.000
7	9	$y$	–1.083	–5.030
7	9	$z$	–1.319	1.117
8	9	$x$	–0.167	0.033
8	9	$y$	–3.154	–1.428
8	9	$z$	–2.313	1.491

Equations

$$\Omega \left(\left\{\vec{e}\right\}\right)={\Omega }_0+\sum _{i=1}^N{\vec{e}}_i{\underline{\underline{K}}}_i{\vec{e}}_i-\frac{1}{2}\sum _{\begin{array}{c}i,j=1\\ i\ne j\end{array}}^N{\vec{e}}_i{\underline{\underline{J}}}_{ij}{\vec{e}}_j,$$

$${\underline{\underline{J}}}_{ij}=J_{ij}^I\underline{\underline{I}}+{\underline{\underline{J}}}_{ij}^S+{\underline{\underline{J}}}_{ij}^A,$$

$$J_{ij}^I=\frac{1}{3}\mathrm{Tr}\left({\underline{\underline{J}}}_{ij}\right)$$

$${\underline{\underline{J}}}_{ij}^S=\frac{1}{2}\left({\underline{\underline{J}}}_{ij}+{\underline{\underline{J}}}_{ij}^T\right)-J_{ij}^I\underline{\underline{I}}$$

$${\underline{\underline{J}}}_{ij}^A=\frac{1}{2}\left({\underline{\underline{J}}}_{ij}-{\underline{\underline{J}}}_{ij}^T\right),$$

$${\vec{e}}_i{\underline{\underline{J}}}_{ij}^A{\vec{e}}_j={\vec{D}}_{ij}\left({\vec{e}}_i\times {\vec{e}}_j\right)$$

$$\Omega \left(\vec{e}\right)={\Omega }_0^{\prime }+\sum _{i=1}^N\vec{e}{\underline{\underline{A}}}_i\vec{e},$$

$${\Omega }_0^{\prime }={\Omega }_0-\frac{1}{2}\sum _{i\ne j}{J}_{ij}^I(-1{)}^{i+j},$$

$${\underline{\underline{A}}}_i={\underline{\underline{K}}}_i-\frac{1}{2}\sum _{j=1}^N{\underline{\underline{J}}}_{ij}^S(-1{)}^{i+j}.$$

$${\vec{k}}_i=\left(k_i^h-k_i^e\right){\vec{e}}_i^e.$$

$${\vec{e}}_k^{\prime }(t_{n+1})={\vec{e}}_k(t_n)-\lambda {\vec{e}}_k(t_n)\times \left({\vec{e}}_k(t_n)\times {\vec{B}}_k^{\text{eff}}(t_n)\right),$$

$${\vec{B}}_k^{\text{eff}}=\sum _j{\underline{\underline{J}}}_{kj}{\vec{e}}_j-2{\underline{\underline{K}}}_k{\vec{e}}_k$$

$$J^k(q)=\sum _{j=1}^k{J}_{8,8+j}^I\cos (jaq),q\in \left[-\frac{\pi }{a},\frac{\pi }{a}\right],$$

$$(D_{ij}^x,D_{ij}^y,D_{ij}^z)=(D_{\sigma (i),\sigma (j)}^x,-D_{\sigma (i),\sigma (j)}^y,-D_{\sigma (i),\sigma (j)}^z),$$

$${\vec{\chi }}_i={\vec{e}}_i\times {\vec{e}}_{i+1}.$$

$${\vec{e}}_i=$$

$$\sin \alpha \sin \left[i\left(\pi +\delta \right)+{\phi }_0\right]),$$

$$\Omega \left(C_1\right)-\Omega \left(C_2\right)-\Omega \left(C_3\right)+\Omega \left(C_4\right)=-2\left(J_{ij}^{zx}-J_{ij}^{xz}\right)$$

$$+2\left(J_{ij}^{xz}+J_{ij}^{zx}\right)\cos 2\theta +2\left(J_{ij}^{xx}+J_{ij}^{zz}\right)\sin 2\theta$$

$${\mathcal{D}}_{\text{r},ij}^y\equiv -\frac{1}{4}⟨\Omega \left(C_1\right)-\Omega \left(C_2\right)-\Omega \left(C_3\right)+\Omega \left(C_4\right)⟩,$$

$$\Omega =$$

$$-\frac{1}{2}\sum _{i,j}{\vec{D}}_{ijij}\left({\vec{e}}_i{\vec{e}}_j\right)\left({\vec{e}}_i\times {\vec{e}}_j\right)$$

$$-\sum _{i,j,k}{\vec{D}}_{ijjk}\left({\vec{e}}_i{\vec{e}}_j\right)\left({\vec{e}}_j\times {\vec{e}}_k\right)$$

$$-\frac{1}{4}\sum _{i,j,k,l}{\vec{D}}_{ijkl}\left({\vec{e}}_i{\vec{e}}_j\right)\left({\vec{e}}_k\times {\vec{e}}_l\right),$$

$${\vec{D}}_{ijij}=-{\vec{D}}_{jiji},$$

$${\vec{D}}_{ijkl}$$

$${\vec{D}}_{ijkl}=\left(\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}\right){\vec{D}}_{\sigma (i)\sigma (j)\sigma (k)\sigma (l)},$$

$${\mathcal{D}}_{\text{r},ij}^{\alpha }=$$

$$+\sum _{k\ne l\notin \{i,j\}}\left(\frac{1}{2}{D}_{klij}^{\alpha }+D_{iklj}^{\alpha }\right).$$

$${\mathcal{D}}_{\text{r},79}^x=\sum _{k\notin \{7,9\}}{D}_{7kk9}^x+\sum _{k\ne l\notin \{7,9\}}{D}_{7kl9}^x,$$

$${\mathcal{D}}_{\text{t},ij}^{\alpha }=$$

$$-\sum _{k\notin \{i,j\}}{D}_{kjji}^{\alpha }+\frac{1}{2}\sum _{k\ne l\notin \{i,j\}}{D}_{klij}^{\alpha };$$

$$\frac{\partial }{\partial {\beta }_{1i}}=-{\vec{e}}_{2i}\frac{\partial }{\partial {\vec{e}}_i}+{\vec{e}}_i\frac{\partial }{\partial {\vec{e}}_{2i}},$$

Full text

I Introduction

The investigation of clusters of magnetic atoms on nonmagnetic surfaces has recently opened several intriguing prospects for the storage and transfer of information on the nanometer scale. The reduced dimensionality of the clusters often leads to an enhancement of the magnetic anisotropy energy Gambardella et al. (2002), stabilizing the magnetic structure in one of two states connected by time-reversal symmetry. These two states can in turn be used for designing logic gates Khajetoorians et al. (2011). Besides the anisotropy, the interactions between the magnetic adatoms mediated by the substrate also crucially influence the magnetic state. One prominent type of these couplings is the Dzyaloshinskii–Moriya (DM) interaction Dzyaloshinsky (1958); Moriya (1960), the presence of which can be attributed to the spin–orbit coupling and the inversion-symmetry breaking caused by the surface. This interaction leads to the formation of noncollinear structures with a preferred chirality by which the information may be encoded Menzel et al. (2012). Linear chains of magnetic atoms on a superconducting surface also offer a possibility for realizing Majorana bound states Alicea (2012) as fundamental elements of topological quantum computing, the signatures of which have been investigated experimentally in chains both with collinear Nadj-Perge et al. (2014) and noncollinear Kim et al. (2018) magnetic ground states.

Determining the ground state for an interacting magnetic system based on ab initio electronic structure calculations remains a considerable challenge. Such computations can efficiently be performed by mapping the energy or grand potential of the system to a classical spin model. It was demonstrated in Ref. Liechtenstein et al. (1987) how Heisenberg exchange interactions between pairs of spins may be determined based on the derivatives of the energy, that is, the torque acting on the spin directions. The torque method has been generalized to tensorial two-spin interactions appearing in the presence of spin–orbit coupling Udvardi et al. (2003); Ebert and Mankovsky (2009), and it was validated for various systems over the last decade Antal et al. (2008); Vida et al. (2016); Simon et al. (2018). A fully real-space calculation of the interactions in a magnetic cluster is presented in Ref. Cardias et al. (2016). However, only considering two-spin interactions in the spin model is not sufficient for describing all types of magnetic order. It was demonstrated in various ultrathin film systems that isotropic four-spin interactions may stabilize up-up-down-down states Al-Zubi et al. (2011); Krönlein et al. (2018); Romming et al. (2018), conical spin spirals Yoshida et al. (2012); Zimmermann et al. (2014), or nanoskyrmion lattices Heinze et al. (2011).

The problem of finding a spin model which contains all types of magnetic interactions relevant in the system may be circumvented by updating the directions of the magnetic moments during the ab initio calculations. Because of the higher number of degrees of freedom the computational complexity increases dramatically, but such methods enable a more accurate determination of the magnetic ground state. It was proposed in Refs. Stocks et al. (1998); Újfalussy et al. (1999) that the constrained local moment method within density functional theory is applicable for performing first-principles spin dynamics simulations. Using this method, it was demonstrated in Ref. Újfalussy et al. (2004) that the reduction of the symmetry leads to a canted magnetic configuration in a finite Co chain along a step edge on the Pt(111) surface. An alternative procedure for updating the spin directions based on the Landau–Lifshitz–Gilbert equation Lan ; Gil was introduced in Ref. Rózsa et al. (2014), where the torques acting on the spins are determined directly from the electronic structure at each time step within a fixed electronic potential.

In the present paper the magnetic properties of monatomic Fe chains are investigated on the Re(0001) substrate. Spin-polarized scanning tunneling microscopy measurements performed for a 40-atom-long Fe chain on superconducting Re in Ref. Kim et al. (2018) revealed a spin spiral ground state with both in-plane and out-of-plane spin components and a period of approximately four lattice constants.

The paper is organized as follows. In Secs. II.1 and II.2 the details of the ab initio calculations are discussed, performed using the Vienna Ab-initio Simulation Package (vasp) Kresse and Furthmüller (1996) and the embedding technique within the KKR method Lazarovits et al. (2002), respectively. In Sec. II.3 the spin model including two-spin interactions is introduced. The results for the parameters entering the spin model are discussed in Sec. III.1. The possible magnetic ground states of the chains obtained from the spin model and from a direct optimization within the ab initio method are compared in Sec. III.2. The deviations between the different methods observed for the chirality of the spin structure for the 15-atom-long chain are resolved by taking into account four-spin chiral interactions introduced in Sec. III.3. Finally, the results are summarized in Sec. IV.

II Methods

II.1 VASP calculations

To model the geometries of Fe atomic chains on the Re(0001) surface, the equilibrium structure of an Fe adatom on Re(0001) has been first calculated by using the VASP method. The obtained structure corresponds to the total energy minimum after geometry optimization. In the calculation the generalized gradient approximation (GGA) within density functional theory (DFT) has been used with the exchange–correlation (XC) functional parametrized following the work of Perdew, Burke, and Ernzerhof (PBE) Perdew et al. (1996). The system has been modeled as a $7\times 7$ surface cell in a slab geometry consisting of four atomic layers of Re (in total $4\times 7\times 7=196$ Re atoms), and an Fe adatom in the hcp hollow position Kim et al. (2018). The chosen geometry ensures that the interactions between Fe atoms in repetitive supercells are negligible due to their large separation of $\sim 19.3$ Å  that corresponds to $7\,a_{Re}$, where $a_{Re}=2.761$ Å  is the in-plane lattice constant of Re. In the (0001) direction a 10-Å-thick vacuum region has been considered to avoid interaction between repetitive slabs. The Brillouin zone was sampled by the Gamma point only due to the large size of the supercell. The Re atoms in the bottom three layers of the slab have been fixed to their hcp bulk positions, and the vertical positions of all Re atoms in the topmost layer and the Fe adatom have been optimized by using a force convergence criterion of 0.01 eV/Å  acting on the individual atoms. The Fe adatom is found to have a spin magnetic moment of $2.61\,\mu_{\textrm{B}}$, and it pulls out its three nearest-neighbor (NN) Re atoms slightly from the top Re layer, arriving at a Fe-Re vertical distance of $1.85$ Å  with respect to these nearest neighbors. We also find that the top Re layer relaxes toward the substrate which leads to a vertical Re-Re distance of $2.16$ Å  between the above-mentioned three NN Re atoms of the Fe adatom and the Re atoms in the subsurface layer, which is smaller than the bulk Re interlayer distance of $2.228$ Å. These Fe-Re and Re-Re vertical distances were used in the subsequent KKR calculations for the Fe adatom and the atomic chains on Re(0001).

II.2 KKR calculations

We used the Green’s function embedding technique based on the KKR multiple scattering theory Lazarovits et al. (2002) to determine the electronic and magnetic properties of the Fe clusters. The Re(0001) surface has been modeled as an interface region between semi-infinite bulk Re and vacuum consisting of eight atomic layers of Re and four atomic layers of empty spheres (vacuum). The energy integrals were performed using 16 points along a semicircle contour in the upper complex semiplane and a sampling of up to 3282 $\vec{k}$ points in the Brillouin zone was used to calculate the Green’s function of the host. The Ceperley–Alder-type of exchange-correlation functionals Ceperley and Alder (1980) as parametrized by Perdew and Zunger Perdew and Zunger (1981) and an angular momentum cutoff of $l_{\textrm{max}}=2$ was considered in the KKR calculations, similarly to Ref. Lászlóffy et al. (2017). A single Fe adatom and chains consisting of five, 10, and 15 Fe atoms were calculated by embedding them in the first vacuum layer with the layer relaxations described in the previous section.

Three different methods have been used to investigate the magnetic properties of the systems. First, the relativistic torque method Udvardi et al. (2003); Ebert and Mankovsky (2009) was applied to determine parameters of a classical spin model restricted to two-spin interactions. The energies of magnetic configurations within this description were compared by atomistic spin model simulations. The spin model is discussed in Sec. II.3, while the method for fitting the parameters is given in Appendix A. Second, the energies of selected spin configurations were compared, such as collinear states with different magnetic orientations or spin spirals with different periods. In the spirit of the magnetic force theorem (MFT) Liechtenstein et al. (1987), these energy differences between magnetic configurations are calculated with fixed electronic potentials based on the band energy, which is obtained by using Lloyd’s formula Lloyd (1967). Third, the ground state of the Fe chain was also determined completely within the ab initio formalism, by updating the spin directions based on the torque acting on them and also performing self-consistent calculations in the obtained spin configurations. This method enables finding local energy minima, ideally the ground state, in the whole configuration space of the spin directions Balogh et al. (2012). It should be noted that in the second and third methods there is no restriction on the possible types of magnetic interactions apart from those enforced by the symmetry of the system.

First, we performed self-consistent calculations for an Fe adatom in hcp position on the top of the Re(0001) substrate with the embedded cluster KKR technique. We considered clusters of different sizes and concluded that the spin magnetic moment of Fe, $m_{\rm Fe}=2.46\,\mu_{B}$, changes by less than 1 % when increasing the size of the cluster from 13 lattice sites including three Re atoms and nine empty spheres in the first NN shell to 122 lattice sites including the first three neighbor shells around the Fe adatom. Note that the obtained spin moment of Fe is about 6% less than the value of $2.61\,\mu_{B}$ from the VASP calculations. The magnetic moment of the Fe atom induces a small ($<0.1\,\mu_{B}$) magnetic moment in the Re atoms directly below it, while the induced moments of farther Re atoms are negligible.

We determined the anisotropy energy of the adatom in the spirit of the MFT by calculating the energy difference $\Delta E$ between the cases where the Fe spin is pointing in-plane ($E_{\parallel}$) and normal to the plane ($E_{\perp}$). Due to the small value of the induced Re moments the exchange-correlation field was set to zero at the Re sites while calculating the energy differences, and we obtained that the single Fe adatom has easy-plane anisotropy with $\Delta E=E_{\perp}-E_{\parallel}=0.905\,\text{meV}$. We confirmed that taking into account the exchange-correlation field on the Re sites leads to a change within about 5 % in the magnetic anisotropy energy; therefore, in all calculations of the Fe chains in terms of the MFT we chose the above approach for simplicity.

We considered close-packed monatomic chains of five, 10, and 15 Fe atoms along the nearest-neighbor direction on the top of Re(0001). In the following this direction will be denoted by $x$, the in-plane direction perpendicular to $x$ by $y$, and the normal-to-plane direction by $z$. Based on our investigations for the adatom, we considered clusters containing the atomic positions in a NN environment relative to the Fe atoms, including 11, 21, and 31 Re atoms, as well as 25, 45, and 65 empty spheres for the chains of five, 10, and 15 Fe atoms, respectively.

For the five-atom-long Fe chain we first performed self-consistent calculations with ferromagnetic (FM) order and used the torque method to generate a spin model. As will be discussed in Sec. III.1, the NN isotropic couplings are strongly antiferromagnetic (AFM), implying that an alternating AFM order is considerably lower in energy than the FM state. In order to check the preference for the AFM state, we recalculated the potentials with the AFM order of the spins and, by using these potentials, we calculated the energy difference between the AFM and FM states within the MFT. Indeed, the AFM state was by 15.7 meV/Fe atom lower in energy than the FM state. Based on the above results, for all the Fe chains under consideration we used the self-consistent potentials obtained from alternating AFM configurations to generate the spin-model parameters.

II.3 Spin model

The adiabatic decoupling of the electronic and spin degrees of freedom and the rigid spin approximation Antropov et al. (1996) make it possible to characterize the energy of a magnetic system by a set of unit vectors $\left\{\vec{e}\right\}\equiv$ $\left\{\vec{e}_{1},\vec{e}_{2},\dots\vec{e}_{N}\right\}$ describing the directions of atomic magnetic moments, where $N$ is the number of magnetic atoms in the system. Since the metallic substrate acts as a particle reservoir for the clusters considered in the calculations, instead of the energy $E$ we will consider the grand potential $\Omega=E-\varepsilon_{\rm F}N_{\rm e}$ at zero temperature, with $\varepsilon_{\rm F}$ and $N_{\rm e}$ being the Fermi energy of the reservoir and the number of electrons in the cluster, respectively. Taking into account one-spin terms and two-spin magnetic interactions, the spin model can be written as

[TABLE]

where $\Omega_{0}$ is a constant, the $\underline{\underline{K}}_{i}$ are traceless and diagonal second-order single-ion anisotropy matrices, and the $\underline{\underline{J}}_{ij}$ are tensorial exchange interactions Udvardi et al. (2003). The matrices $\underline{\underline{J}}_{ij}$ can be decomposed into three parts,

[TABLE]

where

[TABLE]

is the isotropic exchange interaction,

[TABLE]

is the traceless symmetric part of the matrix, with $T$ denoting the transpose. This is known to contribute to the so-called two-ion magnetic anisotropy of the system. The antisymmetric part of the matrix,

[TABLE]

is related to the DM interaction Dzyaloshinsky (1958); Moriya (1960),

[TABLE]

with the DM vector $D_{ij}^{\alpha}=\frac{1}{2}\varepsilon_{\alpha\beta\gamma}J_{ij}^{\beta\gamma}$, $\varepsilon_{\alpha\beta\gamma}$ being the Levi–Civita symbol and $\alpha,\beta,\gamma$ denoting Cartesian components.

Following Ref. Lászlóffy et al. (2017), site-resolved easy-axis directions and anisotropy energies have been determined for the monatomic chains from the spin model, taking into account both single-ion and two-ion contributions. Since the nearest-neighbor isotropic interactions are found to be antiferromagnetic (see Sec. III.1), in order to characterize the magnetic anisotropy we consider alternating local moments $\vec{e}_{i}=(-1)^{i}\vec{e}$. The grand potential of the system can then be expressed as

[TABLE]

with

[TABLE]

and the effective anisotropy matrices

[TABLE]

The normalized eigenvectors of the symmetric matrices in Eq. (9), $\vec{e}^{\,e}_{i}$, $\vec{e}^{\,i}_{i}$, and $\vec{e}^{\,h}_{i}$ correspond in order to the easy, intermediate, and hard directions, with the respective energy eigenvalues $k^{e}_{i}\leq k^{i}_{i}\leq k^{h}_{i}$. For illustrating the site-specific easy directions together with the magnetic anisotropy energies, we will use the following vector:

[TABLE]

The ground state of the spin model was determined by zero-temperature Landau–Lifshitz–Gilbert (LLG) spin dynamics simulations where only the damping term was kept. This is described by the time integration step

[TABLE]

where

[TABLE]

is the effective magnetic field, and a small damping parameter, $\lambda\sim 10^{-4}\,\text{meV}^{-1}$, was chosen. The new spin vectors were normalized after each step to preserve the unit length of the vectors. The simulations were stopped when the spin components changed less than $10^{-6}$ in $5\cdot 10^{5}$ subsequent LLG steps. For each system ten runs with independently chosen random initial configurations were performed which all led to the same final state, providing a strong indication that this is the actual ground state of $\Omega\left(\left\{\vec{e}\right\}\right)$ instead of a local minimum.

III Results

III.1 Spin-model parameters for the Fe chains

In this section we discuss the parameters of the spin model containing two-spin interactions described in Sec. II.3, calculated in terms of the KKR method and the relativistic torque method detailed in Sec. II.2 and in Appendix A, respectively. The variation of NN and next-nearest-neighbor (NNN) isotropic interactions $J_{ij}^{I}$ from Eq. (3) along the chains can be seen in Fig. 1. For all chains, the NN isotropic interactions are the strongest, and their negative sign means AFM coupling. The isotropic interactions become more and more homogeneous at the middle of the chain as the length of the chain is increased. It is also apparent from Fig. 1 that the isotropic NN interactions are considerably smaller in magnitude at the edges of the chain than inside the chain for all chain lengths.

The ferromagnetic NNN interactions are more than one order of magnitude smaller than the NN ones as shown in Fig. 1. This also holds true for the interactions for farther neighbors as can be seen in Fig. 2(a), where calculated values for the middle spin $8$ in the 15-atom-long chain are displayed. However, if one summarizes the effect of farther interactions, it turns out that they play an important role in determining the ground state. This is demonstrated by introducing the Fourier transform of the isotropic interactions as

[TABLE]

where $k$ is the number of neighbors taken into account. If one assumes that the interactions in the middle of the chain will no longer be significantly modified as the chain length is increased, then $J^{k}(q)$ may be used as an approximation for the energy contribution of the isotropic interactions to homogeneous spin spiral states in infinitely long chains, where $q=0$ corresponds to the FM and $q=\frac{\pi}{a}$ to the AFM state. The most favorable state is given by the maximum of $J^{k}(q)$. It can be seen in Fig. 2(b) that up to $k=5$ this corresponds to the collinear AFM state. However, considering more shells ($k>5$) in the sum in Eq. (13) a spin spiral state becomes the most favorable, which can be regarded as a long-wavelength modulation of the AFM state. This is a consequence of the frustration of the isotropic interactions, which in the AFM state is indicated by the fact that the isotropic interaction between atoms at the distance of an odd multiple $j$ of the lattice constant becomes FM and at an even $j$ it becomes AFM, although the spins at odd and even $j$ positions should be antiparallel and parallel in the AFM state, respectively. It is shown in Fig. 2(a) that such kind of frustration occurs for $j\geq 4$, explaining how the spin spiral state is formed as the number of shells is increased.

The NN and NNN DM vectors in the chains, see Eq. (6), are drawn in Fig. 3. Similarly to the isotropic exchange interactions, the DM vectors at the middle of the chain become stabilized as the chain length is increased. The DM vectors between the three atoms at the edges of the chains significantly differ from those in the middle of the chains, while these DM vectors are similar between the five-, 10-, and 15-atom-long chains. The obtained DM interactions satisfy the symmetry rules of Moriya Moriya (1960) with respect to the only crystal symmetry of the system, namely the mirroring at the $yz$ plane intersecting the middle of the chain. Since the DM vector transforms as an axial vector, this symmetry implies

[TABLE]

where $\sigma(i)=N+1-i$ is the mirror image of site $i$ in the chain of length $N$. Note that in Fig. 3 the vectors are displayed for $i<j$ and $\vec{D}_{ij}=-\vec{D}_{ji}$ by definition.

In the case of the 15-atom-long chain, the numerical values for the components of the NN and NNN DM vectors are given in Table 1. It can be seen that the NNN DM vectors are the largest in magnitude and they are almost parallel to the $-y$ direction. Although the $y$ components of the NN DM vectors have the same sign, these are actually competing with the NNN vectors due to the short-range AFM order (see Sec. III.2), similarly how the alternating isotropic interactions in Fig. 2(a) are competing with the NN interaction. The $z$ components of the NN and NNN DM vectors are mostly positive; they only change sign for the NN atoms at the edge of chain ($D^{z}_{12}<0$). In general, $D^{z}_{ij}$ is larger for the NNs than for the NNNs, but at the middle of the chain they become roughly similar in size. The $x$ components of the DM vectors are very small and they should disappear in the limit of infinitely long chains due to the mirror symmetry with respect to the $yz$ plane. For the 15-atom-long chain, the mirror symmetry implies $D^{x}_{79}=0$.

The site-resolved anisotropy vectors defined in Eq. (10) are visualized in Fig. 4, where the arrows point along the easy directions and their magnitude is proportional to the energy difference between hard and easy axes at the given site. At all sites the easy axis is almost parallel to the $y$ direction, while the hard axis is roughly along the $z$ direction. With increasing chain length the magnitude of the $\vec{k}_{i}$ vectors gets quite homogeneous with the maxima of 7.79, 9.03, and 8.92 meV at the middle of the five-, 10-, and 15-atom-long chains, respectively. However, since the site-resolved magnetic anisotropy energy drops at the edge of the chains, the average length of the anisotropy vectors is 6.25, 7.60, and 7.91 meV for the three chains in order. These values are close to the average anisotropy energies, $(\Omega(\hat{z})-\Omega(\hat{y}))/N$ calculated from Eq. (7), 5.79, 7.37, and 7.76 meV/Fe, respectively, since as noted the local easy and hard axes are very close to the $y$ and $z$ axes for all the Fe atoms in the chains.

III.2 Ground state

Here the ground states of the magnetic clusters will be discussed, focusing on the chirality of spin rotation inside the structures. We will use the convention that the rotation of the spins in the $xz$ plane is locally right handed at site $i$ if the projection of $\vec{e}_{i+1}$ on this plane may be obtained from the projection of $\vec{e}_{i}$ via a right-handed rotation by an angle smaller than $180^{\circ}$ around the positive $y$ axis, meaning that the angle between the projections is smaller than $180^{\circ}$ when rotating from the positive $z$ towards the positive $x$ direction. The opposite chirality is called left handed, in agreement with previous definitions of the chirality for spin spirals close to the ferromagnetic state in Ref. Heide et al. (2008). Analogously, we call the rotation in the $xy$ plane locally right handed if the projections of the spins rotate from the positive $x$ towards the positive $y$ direction.

First we determined the ground states of the chains from the spin model following the method described in Sec. II.3. The obtained configurations can be seen in Figs. 5-5. The five-atom-long chain prefers AFM ordering due to the strong AFM coupling between the NN spins. The $z$ component of the DM vectors is responsible for a slight noncollinearity in the AFM state. The negative $D_{12}^{z}$ value in Fig. 3 causes a slight left-handed rotation between spins 1 and 2, while the positive $D_{23}^{z}$ causes a right-handed rotation between spins 2 and 3. The ground state is tilted away from the $xy$ plane through a rotation around the $x$ axis by about 12.8∘, so the spins with positive $y$ component now have positive $z$ component, too. The tilting is caused by the competition of the easy $y$ axis anisotropy (see Fig. 4) and the large $y$ components of the DM vectors (similar to those in Table 1) preferring a rotation of the spins in the $xz$ plane.

The ground state of the 10-atom-long chain in Fig. 5 is again almost collinear AFM due to the strong NN AFM coupling between the spins. The $z$ components of the DM vectors shown in Fig. 3 cause a rotation of the spins around the $z$ axis along the whole chain, but no full period of a spin spiral state can be observed. The ground state is now tilted from the $xy$ plane through a rotation around the $x$ axis by $-6.2^{\circ}$.

In both the five- and 10-atom-long chains, the strong anisotropy confines the systems close to the $xy$ plane. The chirality is right handed over the middle three atoms in the five-atom-long chain and over the middle eight atoms in the 10-atom-long chain. Note that only looking at every second spin in the chain, this visually corresponds to a left-handed rotation, since the angle between the neighboring spins is close to $180^{\circ}$ because of the strong AFM NN couplings. This chirality is determined by the $z$ components of the DM vectors shown in Fig. 3. The direction of the tilting is defined by the chirality in the $xy$ plane on the one hand and the $y$ components of the DM vectors on the other hand, the latter influencing the rotational sense in the $xz$ plane. In both systems the NN and NNN DM vectors are characterized by large negative $y$ components, both of which would prefer a left-handed rotation if the angle between the NN spins would be small. However, since this angle is larger than $90^{\circ}$ in the present case, the apparent left-handed rotation between the NNN spins actually corresponds to a right-handed rotation between the NN spins, meaning that the NN and NNN DM vectors are competing in this AFM spin structure. In the five-atom-long chain the rotational plane is tilted around the $x$ axis by a positive angle, leading to a left-handed chirality in the $xz$ plane enforced by the NN DM vectors. For the longer chain length of $10$ atoms with the same chirality in the $xy$ plane the influence of the NNN DM vectors becomes stronger, leading to a tilting around the $x$ axis by a negative angle and a right-handed chirality in the $xz$ plane.

The ground state of the 15-atom-long chain shown in Fig. 5 can much better be characterized as a spin spiral state. As described in the context of Fig. 2, the frustrated isotropic interactions induce a spin spiral with a wave number of $q_{\text{max}}=0.88\frac{\pi}{a}$ (see the $k=7$ curve), which corresponds to a wavelength of about $17\,a$, slightly larger than the length of the chain. In Fig. 5 the spins from positions 4 to 11 visually form a full period, which can be understood as a $14\,a$ wavelength modulation of the AFM state. This shorter modulation period might easily be caused by the $z$ components of the NN DM vectors, which are not included in Fig. 2. The large negative $y$ components of the NNN DM vectors play an important role in tilting the rotational plane out from the $xy$ plane by $-26.1^{\circ}$ around the $x$ axis, for a simple explanation see Appendix C. A similar tilted spin spiral state was already attributed to the interplay of the DM interaction and easy-plane magnetic anisotropy in Ref. Simon et al. (2018). Similarly to the 10-atom-long chain, the DM vectors imply right-handed rotation both in the $xy$ plane and in the $xz$ plane over the whole chain. This ground state cannot satisfactorily be reconciled with the experimental observation of a four-atomic period in a chain of 40 Fe atoms reported in Ref. Kim et al. (2018), which would correspond to a spin spiral state where neighboring spins are perpendicular to each other.

The local rotational sense of the spins in a spin spiral discussed above can be quantitatively described by the site-dependent chirality vector defined as

[TABLE]

If only NN DM vectors were included in the model, $\vec{\chi}_{i}$ would be parallel to the direction of $\vec{D}_{i,i+1}$. The chirality is right-handed in the $xy$ plane and in the $xz$ plane for $\chi_{i}^{z}>0$ and $\chi_{i}^{y}>0$, respectively. In Fig. 6 the magnitude of the local chirality vectors $\left|\vec{\chi}_{i}\right|=\sin\varphi_{i}$, where $\varphi_{i}\in\left[0^{\circ},180^{\circ}\right]$ is the angle between the NN spins, as well as the $y$ and $z$ components of the chirality vectors are displayed for the 15-atom-long chain. The angle between the spins at the edges of the chain is almost $180^{\circ}$ and also the sign of $\chi_{i}^{y}$ is switched compared to the middle of the chain, which can be explained by the edge effects in the interaction parameters discussed in Sec. III.1. In the middle of the chain the local environment of the sites is similar, leading only to slight variations in $\varphi_{i}$. The $y$ component of the chirality vectors takes a value between 0.1 and 0.2 for the spins 2 to 13, indicating that the spins tilt away from the $xy$ plane as argued in Appendix C.

For comparison with the ground state obtained from the spin model, we calculated the energy of the 15-atom-long Fe chain in the homogeneous spin spiral configuration,

[TABLE]

in the spirit of the MFT by keeping the AFM potentials fixed. Here $\delta$ is the modulation angle of the AFM state, $\alpha$ is the tilting angle of the spiral from the $xy$ plane and $\varphi_{0}$ is a phase factor. The NN spin angle is $\varphi=\pi+\delta$, indicating right-handed rotation in the $xz$ and left-handed rotation in the $xy$ planes for $\delta>0^{\circ}$ and $0^{\circ}\leq\alpha\leq 90^{\circ}$, respectively. Both rotational senses switch under a sign change of $\delta$, and the rotation in the $xy$ plane proceeds in the opposite direction for $90^{\circ}\leq\alpha\leq 180^{\circ}$. We did not consider an additional angle variable which would differentiate between left- and right-handed rotational senses in the $yz$ plane, since in an infinitely long chain these two chiralities are equivalent due to the mirror symmetry with respect to the $yz$ plane. The phase factor $\varphi_{0}$ was determined by minimizing the anisotropy energy assuming a homogeneous magnetic anisotropy with $y$ easy axis, i.e., maximizing $\sum_{i=1}^{N}\sin^{2}\left[i\left(\pi+\delta\right)+\varphi_{0}\right]$ for every $\delta$.

The energy of the spin spirals normalized to one Fe atom is shown in Fig. 7 as a function of $\delta$ and $\alpha$, where the zero level corresponds to the state with the lowest energy. The obtained minimum is at $\delta=24^{\circ}$ and $\alpha=12^{\circ}$, with the corresponding spin configuration shown in Fig. 5. While the period of this spin spiral, $\lambda\simeq 15\,a$, shows good agreement with that obtained from the simulations based on the spin model shown in Fig. 5, the actual values of $\delta$ and $\alpha$ indicate a right-handed rotation in the $xz$ and a left-handed rotation in the $xy$ planes (see above), the latter being the opposite of the spin model simulation results. This can clearly be seen in Fig. 6(c), where $\chi^{z}_{i}$ for the ground state of the spin model and for the spin spiral with lowest energy are opposite in sign.

In the case of the 15-atom-long chain we also performed a fully ab initio spin dynamics energy minimization as described in Sec. II.2. The energetically most favorable configuration found by this method is shown in Fig. 5, which also resembles a flat spin spiral state apart from small deviations from the coplanar spin arrangement. This method predicts a right-handed rotation of the spins in the $xz$ and a left-handed rotation in the $xy$ planes, in agreement with the MFT calculations performed for the homogeneous spin spiral states. This approach results in a wavelength of $\lambda\simeq 10\,a$ ($\delta\simeq 36^{\circ}$) of the spin spiral modulation which is significantly shorter than those obtained from the previous two methods. On top of the AFM state, this visually leads to a $5\,a$ period of the spin structure, see spins 3, 8, and 13 in Fig. 5, which fits the experimental observation the most.

III.3 Four-spin chiral interactions

It was found in Sec. III.2 that the spin model Eq. (1) and the MFT calculation of homogeneous spin spirals yield opposite rotational senses of the spin components in the $xy$ plane for the 15-atom-long chain. This most likely indicates that multispin interactions play an important role in the present system, since these were not taken into account in the spin model, but are implicitly included in the MFT calculations. In order to estimate the magnitude of these interactions, it is worthwhile to calculate the chiral interactions from energy differences based on the MFT directly and compare them to the values of the spin model shown in Table 1. This method is illustrated for the $y$ component in Fig. 8, where four configurations in the $xz$ plane are shown. In the configuration denoted by $C_{1}$, during a global rotation of the spins around the $y$ axis by the angle $\theta$ most spins are rotated according to $\vec{e}=\left(\sin\theta,0,\cos\theta\right)$, while the perpendicular spin at site $j$ will follow $\vec{e}_{\bot}=\left(\cos\theta,0,-\sin\theta\right)$. The configurations $C_{2},C_{3},$ and $C_{4}$ are obtained by switching the signs of spins $i$ or $j$ as illustrated on the right side of Fig. 8.

Assuming the spin model containing only two-spin interactions, from the grand potentials associated to the configurations $C_{1}$-$C_{4}$ one obtains

[TABLE]

as a function of the rotation angle $\theta$. Note that due to the choice of configurations and the switching of the spin directions only the interaction between sites $i$ and $j$ remains in Eq. (17). Averaging Eq. (17) over the angle $\theta$ we define the $y$ component of the chiral interaction vector as

[TABLE]

which, by comparing with Eq. (6), corresponds to the $y$ component of the DM vector $\mathcal{D}^{y}_{\text{r},ij}=D_{ij}^{y}=\frac{1}{2}(J_{ij}^{zx}-J_{ij}^{xz})$, if only two-spin interactions are considered in the spin model. Here the index r denotes that the chiral interaction vector was obtained from the rotational scheme depicted in Fig. 8. Analogously, the quantities $\mathcal{D}^{x}_{\text{r},ij}$ and $\mathcal{D}^{z}_{\text{r},ij}$ corresponding to the other two components of the DM vector can be calculated by performing the rotation in the $yz$ and $xy$ planes, respectively.

The calculated values for $\mathcal{D}^{\alpha}_{\text{r},ij}$ ($\alpha=x,y,z$) are collected in Table 2 for NN and NNN spins in the middle of the 15-atom-long chain. In addition, the $\mathcal{D}_{\text{t},ij}^{\alpha}$ values obtained from the torque method, defined in Eq. (41) in Appendix A, are listed in the table, which within the spin model Eq. (1) should also coincide with $D_{ij}^{\alpha}$ (cf. Table 1). Due to the mirror symmetry with respect to the $yz$ plane going through atom $8$, the symmetry rules for the DM vectors, Eq. (14), imply $\left(\mathcal{D}_{\text{t},78}^{x},\mathcal{D}_{\text{t},78}^{y},\mathcal{D}_{\text{t},78}^{z}\right)=\left(-\mathcal{D}_{\text{t},89}^{x},\mathcal{D}_{\text{t},89}^{y},\mathcal{D}_{\text{t},89}^{z}\right)$ and $\mathcal{D}_{\text{t},79}^{x}=0$. Remarkably, these symmetry relations do not apply to the corresponding chiral interaction vectors obtained from the rotational method; in particular, $\mathcal{D}_{\text{r},79}^{x}$ does not vanish, but it has a comparable value to the other components. This means that the spin model parametrization of the band energy surface of the present system is not compatible by taking into account two-spin DM interactions only.

In order to explain the differences between $\mathcal{D}_{\text{r},ij}^{\alpha}$ and $\mathcal{D}_{\text{t},ij}^{\alpha}$ in Table 2, it is necessary to include multispin chiral interactions in the model description. Consider a grand potential of the form

[TABLE]

where the last three terms represent two-, three-, and four-site four-spin chiral interactions combining isotropic (scalar product) and DM (cross product) contributions. In the sums in Eq. (19), the $i$, $j$, $k$, and $l$ indices run over all lattice sites, and in each sum the different indices label different atoms. Equation (19) may be rewritten in an alternative way where the summations are performed over all pairs, three-site and four-site clusters; this is presented in Appendix B.

The two-site four-spin chiral interactions were recently investigated in Ref. Brinker et al. (2019). Following from the definition Eq. (19), they are antisymmetric in their indices,

[TABLE]

similarly to the two-spin Dzyaloshinsky–Moriya interaction in Eq. (6). For the three-site interactions we will consistently use the notation $\vec{D}_{ijjk}$ where the second and third site indices coincide, in which case there is no intrinsic symmetry relation connecting the coefficients in the sum in Eq. (19). By definition, the four-site interactions $\vec{D}_{ijkl}$ satisfy the symmetry relations

[TABLE]

therefore, a prefactor of $1/4$ is introduced in the last term of Eq. (19). Moreover, the mirror symmetry on the $yz$ plane in the center of the chain implies

[TABLE]

where $\sigma(i)$ was defined in the context of Eq. (14). Equation (23) is satisfied for two-, three-, and four-site chiral interactions.

If Eq. (18) is evaluated based on the model Eq. (19), one obtains

[TABLE]

Note that the two-site four-spin interaction $\vec{D}_{ijij}$ does not contribute to the grand potential of the configurations shown in Fig. 8, since it is only finite between pairs of spins which are not parallel and not perpendicular to each other. Considering the chiral interaction vectors in Table 2 it is possible to derive

[TABLE]

which generally does not vanish since the three-site and four-site chiral interactions are not antisymmetric with respect to their first and fourth indices. This indicates the presence of four-spin chiral interactions in the system on the order of $0.4\,\text{meV}$, which is not negligible compared to the total value of chiral interactions also containing two-site contributions displayed in Table 2. The relation $\left(\mathcal{D}_{\text{r},78}^{x},\mathcal{D}_{\text{r},78}^{y},\mathcal{D}_{\text{r},78}^{z}\right)\neq\left(-\mathcal{D}_{\text{r},89}^{x},\mathcal{D}_{\text{r},89}^{y},\mathcal{D}_{\text{r},89}^{z}\right)$, which breaks the symmetry rules for two-spin DM interactions, may similarly be explained by the presence of the four-spin chiral interactions.

Changing from the two-spin model in Eq. (1) to the model containing also four-spin interactions in Eq. (19) naturally modifies the interpretation of the chiral interaction energies obtained from the torque method,

[TABLE]

for details of the derivation see Appendix A. It should be noted that contrary to Eq. (24), Eq. (26) is antisymmetric with respect to the site indices $i$ and $j$, thus preserving the symmetries of the two-site DM vectors. Most importantly, the different treatment of the two-, three-, and four-site four-spin interactions between the rotational scheme in Eq. (24) and the torque method in Eq. (26) can explain why the $z$ components of the two kinds of chiral interaction vectors have different signs in Table 2. While the torque method supports positive $\chi_{i}^{z}$ or right-handed chirality in the $xy$ plane, the rotational scheme supports negative $\chi_{i}^{z}$ or left-handed chirality in the same plane. The chirality derived from the rotational method is in agreement with the lowest-energy spin spiral state obtained from MFT calculations and the ground state found by the ab initio energy minimization discussed in Sec. III.2.

IV Conclusion

Ab initio electronic structure calculations were performed to study the magnetic properties of Fe monatomic chains on the Re(0001) substrate. For all the considered chains strong antiferromagnetic couplings between the nearest-neighbor spins were observed, which in the case of the five-atom-long chain led to a nearly collinear antiferromagnetic ground state with the spins approximately aligned perpendicular to the chain and to the surface normal. As the length of the chain is increased, the frustration of the isotropic exchange interactions at farther neighbors transforms the ground state into a spin spiral state, with a single modulation period observable in the 15-atom-long chain. The experimental investigations in Ref. Kim et al. (2018) also concluded on a spin spiral ground state of 40-atom-long Fe chains, although with a smaller period.

For the 15-atom-long chain the spin components in the $xz$ plane were found to follow a right-handed rotation in the spin spiral. Regarding the rotation in the $xy$ plane, calculations based on a spin model only containing two-spin interactions yielded a right-handed rotation, while determining the optimal homogeneous planar spin spiral state using magnetic force theorem calculations and performing an optimization of the configuration directly within the ab initio scheme both indicated a left-handed chirality. This discrepancy was resolved by considering chiral multispin interactions in the spin model, the presence of which in the system is supported by calculations of specific rotated spin configurations sensitive to the chirality.

Since the experiments were carried out in the superconducting phase of the Re substrate, future ab initio calculations including the superconducting state of Re by solving the Bogoliubov–de Gennes equation Csire et al. (2015, 2018) may reveal the reasons behind the observed discrepancies between theory and experiment regarding the spin spiral period. Such first-principles calculations would also enable the investigation of the interplay between exotic magnetic states and topological superconductivity in clusters of magnetic atoms on a superconducting substrate.

Acknowledgements.

The authors would like to thank R. Wiesendanger for useful discussions. Financial support of the Nemzeti Kutatási Fejlesztési és Innovációs Hivatal (Hungary) under Projects No. K115575 and No. FK124100, the BME-Nanotechnology FIKP grant of Emberi Erőforrások Minisztériuma (BME FIKP-NAT), the Alexander von Humboldt Foundation, and the Deutsche Forschungsgemeinschaft via SFB668 are gratefully acknowledged.

Appendix A Torque method

Here we give a short summary of how the parameters of the classical spin model in Eq. (1) were determined by using the torque method; for more details see Ref. Rózsa (2016). As shown in Fig. 9, at each site $i$ the spin direction will be denoted by $\vec{e}_{i}$, and two orthogonal vectors $\vec{e}_{1i}$, $\vec{e}_{2i}$ are defined which form a right-handed basis together with $\vec{e}_{i}$, and around which $\vec{e}_{i}$ is rotated by the infinitesimal angles $\beta_{1i}$ and $\beta_{2i}$, respectively. The derivatives with respect to these angles may be expressed as

[TABLE]

Supposing the model of Eq. (1), the second derivatives of $\Omega$ with respect to the angles $\beta_{1i}$ and $\beta_{2i}$ can be calculated as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\alpha$ and $\beta$ label summations over Cartesian indices.

Within the torque method the derivatives on the left-hand sides of Eqs. (29)–(35) are obtained directly from first principles, using the expressions derived from Lloyd’s formula given in Ref. Udvardi et al. (2003). The diagonal elements of the anisotropy tensor $\underline{\underline{K}}_{i}$ were determined from Eqs. (33) and (34) via

[TABLE]

which simplifies to $K_{i}^{yy}-K_{i}^{zz}$ if $\vec{e}_{i}$ is pointing along the $x$ direction. This way it is only possible to determine the differences between the diagonal elements, but this information is sufficient since $\underline{\underline{K}}_{i}$ was defined as a traceless tensor because its trace only shifts the grand potential by a constant factor due to the normalization of the spins. The two independent components were computed by considering spin configurations where all spins are pointing along the $x$ and $y$ directions.

For calculating the exchange interaction tensor $\underline{\underline{J}}_{ij}$ one has to rely on Eqs. (29)–(32). By introducing the tensor products

[TABLE]

for $r,p=1,2$, and the notation for the derivatives

[TABLE]

with $\overline{r}=3-r$ the opposite angle index, Eqs. (29)–(32) can be summarized as

[TABLE]

where $\underline{\underline{J}}_{ij}$ is rewritten as a vector $\vec{J}_{ij}$ in the tensor product space. By introducing $\underline{\underline{V}}_{ri,pj}=\vec{v}_{ri,pj}\circ\vec{v}_{ri,pj}$, the orthogonal projection onto the subspace of $\vec{v}_{ri,pj}$, summing over Eqs. (29)–(32) one obtains the system of linear equations

[TABLE]

Since $\underline{\underline{V}}_{ri,pj}$ is a rank-1 matrix, the sum over four equations in Eq. (40) only enables the calculation of four components of $\vec{J}_{ij}$ as described in Ref. Udvardi et al. (2003). In order to determine the full tensor, Eq. (40) has to be summed up over calculations performed for at least three linearly independent directions of the $\vec{e}_{i}$ vectors, corresponding to a least-squares fitting procedure for $\vec{J}_{ij}$. Due to the $C_{3\textrm{v}}$ symmetry of the Re(0001) surface, during the calculations four ferromagnetic configurations of the $\vec{e}_{i}$ vectors were taken into account, including three NN directions parallel to the $x$ axis and at angles $60^{\circ}$ and $120^{\circ}$ with respect to this direction, as well as one along the out-of-plane $z$ direction.

The chiral interaction vectors $\vec{\mathcal{D}}_{\text{t},ij}$ introduced in Sec. III.3 were calculated in ferromagnetic configurations with all spins pointing along one of the $\alpha=x$, $y$, or $z$ directions, i.e. $\vec{e}_{i}=\vec{e}_{\alpha}$ in Fig 9. The components of $\vec{\mathcal{D}}_{\text{t},ij}$ are defined as

[TABLE]

which within the model Eq. (1) containing only two-spin interactions simplifies to the DM vector [cf. Eqs. (30) and (31)],

[TABLE]

On the other hand, inserting the spin model Eq. (19) instead of Eq. (1) into Eq. (41) and using the relations Eqs. (27) and (28) to calculate the second derivatives yields Eq. (26) in Sec. III.3,

[TABLE]

Appendix B Alternative notation for the four-spin chiral interactions

Equation (19) follows the convention where the summations are performed for each index separately over all lattice sites, treating the cases where some of the indices coincide individually. This is in agreement, e.g., with the notation for the two-site two-spin and four-spin chiral interactions used in Ref. Brinker et al. (2019). Another way of expressing the grand potential or the Hamiltonian is by performing the summations over all different plaquettes, i.e. pairs, triangles, and quadrilaterals in the case of the two-, three-, and four-site interactions, respectively. Such a convention is used, e.g., in Ref. Hoffmann and Blügel (2018) for the isotropic four-spin interactions.

Following this convention, Eq. (19) may be rewritten as

[TABLE]

Once again, the $i$, $j$, $k$, and $l$ indices all label different atoms. The sums over pairs contain half as many terms as the two-site summations in Eq. (19); however, this is resolved by taking into account the intrinsic symmetry relations $\underline{\underline{J}}_{ij}=\underline{\underline{J}}_{ji}^{T}$ and Eq. (20), which enable canceling the $1/2$ prefactor. As discussed in the main text, no such intrinsic symmetry relations exist for the three-site four-spin chiral interaction, which necessitates a summation over six different terms in Eq. (LABEL:eqn60), corresponding to the possible permutations of the site indices in a triangle. For the four-site interactions, the symmetry relations in Eq. (22) cancel with the prefactor $1/4$ in Eq. (19), simplifying the sum over $24$ permutations in a quadrilateral to six different terms in Eq. (LABEL:eqn60).

Note that Eq. (LABEL:eqn60) is the most general form of the grand potential containing four-spin chiral interactions. Crystal symmetries may reduce the number of independent coefficients for specific plaquettes. For example, the mirror symmetry defined in Eq. (23) implies that in the triangle with $i=\sigma(j)$ and $k=\sigma(k)$ there are only three inequivalent three-site four-spin chiral interactions instead of six in the general case.

Appendix C Tilted spin spirals

In this section we summarize why the spin spiral ground state becomes tilted away from the $xy$ plane due to the $y$ component of the DM vectors. We will consider a simplified spin model consisting of $N$ sites where the interactions are homogeneous along the chain, with $K^{yy}<0$ anisotropy accounting for the easy $y$ direction and $D^{y}$ DM interaction between NN spins. We will compare the energies of harmonic spin spiral configurations as defined in Eq. (16), simplifying the description to two parameters $\delta$ and $\alpha$. It is assumed that a spin spiral state with a specific $\delta$ value is formed by the frustration of the isotropic exchange interactions, as was discussed in Sec. III.2, and only the dependence of the grand potential of the spin model on the $\alpha$ parameter is considered. This is given by the expression

[TABLE]

where now $\Omega_{0}$ describes the contributions which do not depend on $\alpha$, such as the isotropic interactions. For $N\geq 3$, the average of $\sin^{2}\left[i\left(\pi+\delta\right)+\varphi_{0}\right]$ with respect to atomic indices equals $1/2$, and using a simple addition formula it is possible to arrive at

[TABLE]

One can obtain the value of $\alpha$ minimizing the grand potential by differentiation,

[TABLE]

for $\left|\left(N-1\right)\sin\delta D^{y}/NK^{yy}\right|\leq 1$, and

[TABLE]

otherwise. This indicates that for an arbitrarily small value of $D^{y}$, the most preferred state is tilted away from the $xy$ plane where $\alpha=0$. We note that including the $z$ component of the DM vectors influences the dependence of the grand potential on the $\alpha$ parameter, but does not change this qualitative conclusion.