Multiple-scattering approach for multi-spin chiral magnetic interactions: application to the one- and two-dimensional Rashba electron gas

Non-collinear spin-textures, such as the one-dimensional (1D) domain walls [1], 2D magnetic skyrmions [2] and 3D hedgehog states [3, 4] or magnetic hopfions [5] are at the heart of investigations exploring their potential realization as magnetic bits for future information technology. Likewise, such complex spin-states are key players in quantum computing by being paramount in the physics of non-Abelian Majorana bound states [6–15] induced by magnetic nanostructures at the vicinity of superconductors and in several proposals for topological superconductivity [16–20].

The occurrence of complex magnetism is driven by the existence of various competing magnetic exchange interactions (MEI). The initially proposed interactions are bilinear in nature such as the isotropic MEI, ${J}_{ij}^{\text{iso}}\enspace \left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}\right)$ [21, 22], involving two spins while the Dzyaloshinskii–Moriya interaction (DMI), D_ij ⋅ (S_i × S_j), requires the presence of spin–orbit coupling (SOC) and broken inversion symmetry [23, 24]. The direction of the DM vector, D, settles the sense of rotation of the magnetic moments, defining thereby and usually the chirality of the magnetic structures. Besides the isotropic MEI and DMI, a compass term, J^ani [25, 26], can be finite, which define the tensor of bilinear (two-site two-spin) MEI. ${J}^{\text{ani}}\in {\mathbb{R}}^{3{\times}3}$ represents a (trace-less symmetric) matrix in the space of spin-moment orientation giving rise to a non-local magnetic anisotropy energy and behaves like the latter quantity.

The importance of going beyond the bilinear magnetic interactions was recognized several decades ago. Solid ³He is a good example, where multi-spin isotropic interactions are a key ingredient to grasp its magnetic behavior [27, 28]. It is established that if not properly taken into account, such higher-order interactions contribute to effective bilinear ones (see e.g. references [29, 30]). In fact, these interactions can involve an unequal number of magnetic moments and number of sites [31–33] but without breaking time reversal symmetry, the number of magnetic moments should be even. For instance, the biquadratic interaction, connecting four spin positioned on two sites, can be isotropic [34–36], ${\mathcal{B}}_{ij}{\left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}\right)}^{2}$, or chiral as recently discovered [37], ${\boldsymbol{\mathcal{C}}}_{ij}\cdot \left({\mathbf{S}}_{i}{\times}{\mathbf{S}}_{j}\right)\left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}\right)$. The latter follows Moriya's rules similarly to the DMI. Generalization to four-spin interactions can be performed at the isotropic [38–40] or the chiral level [37, 41–43], which are found to be large with non-negligible impact on the magnetic behavior of materials [37, 42, 43]. The microscopic origin of the chiral four-spin interactions was given in reference [37], where the focus was on terms proportional to ${\boldsymbol{\mathcal{C}}}_{ijkl}\cdot \left({\mathbf{S}}_{i}{\times}{\mathbf{S}}_{j}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{l}\right)$, while those of the form $\left({\boldsymbol{\mathcal{C}}}_{ijkl}\cdot {\mathbf{S}}_{i}\right){\chi }_{jkl}$, with χ_jkl = S_j ⋅ (S_k × S_l) being the scalar three-spin chirality, were addressed and derived in reference [42] after a systematic multiple-scattering expansion in the spirit of the infinitesimal rotation method of the magnetic moments. The scalar three-spin chirality is finite whenever the three magnetic moments are non-coplanar. In such non-coplanar spin-configurations, an orbital moment, chiral in nature, can emerge even without SOC [44–47]. This implies that the chiral four-spin interaction can be interpreted as an interaction between the spin and chiral orbital moments. One notes that subsequently a derivation of some of the four-spin terms has been proposed in [48]. As shown in reference [43], the various forms of the chiral four-spin interactions can be re-expressed in terms of sums involving ${\boldsymbol{\mathcal{C}}}_{ijkl}\cdot \left({\mathbf{S}}_{i}{\times}{\mathbf{S}}_{j}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{l}\right)$. Also it is demonstrated that these interactions do not have to follow Moriya's rules of the DMI. Surprisingly, even isotropic six-spin interactions, which can be written in terms of a product of scalar three-spin chiralities χ_ijkχ_lmn, can be considerably large as shown by Grytsiuk et al [42] for B20-MnGe.

The recent developments triggered by the investigations of multi-spin interactions motivate the present study. Their potential impact on the magnetism of materials and their magnetodynamical properties remains to be explored. The focus here is on the long-range behavior of such interactions, which is textbook knowledge for the isotropic bilinear Ruderman–Kittel–Kasuya–Yosida (RKKY) interactions [49–51]. The latter have key implications in the interlayer exchange coupling of multilayers [1, 52–54] and when taken to their isotropic four-spin form, they can trigger the superposition of multiple Q-states [55–58] and the stabilization of 2D magnetic skyrmions [59] or possibly 3D magnetic textures [60–64]. In the context of nuclear spins embedded in a Luttinger liquid, the RKKY interactions emerging from weak hyperfine interaction between the nuclear spin and the conduction electron spin can trigger complex helimagnetic phases. A strong renormalization of the nuclear Overhauser field is predicted accompanied with a universal reduction in the electric conductance and a potential gaping of the electronic spectrum [65, 66].

Note that long-range interactions are intimately connected to conventional and complex magnetic Friedel oscillations [67–69], which do not only affect the isotropic MEI, the DMI [70–76] but also the magnetic anisotropy energy [77, 78], the magnetic stability and spin-fluctuations of nanostructures [79]. The goal of this article is to establish and analyse the long-range behavior of chiral multi-spin interactions within the Rashba model [80, 81] used to describe the 1D and 2D electron gas subjected to the Rashba SOC. For that a multiple-scattering approach based on Green functions partly presented in reference [42] is utilized after reviewing it. The multi-spin interactions are quantified on the basis of parameters describing the Rashba 1D and 2D surface states characterizing Au surfaces [82–84]. I consider localized magnetic atomic spins, which can be manipulated and probed on surfaces utilizing scanning tunneling microscopy/spectroscopy (STM/STS). This tool became ubiquitous to quantify the isotropic [85–90] and DM [73] long-range bilinear interactions. In some cases, puzzling magnetic responses were probed in exotic man-crafted structures [89] indicating that there is still a lot to unravel in this field of research. STM/STS would be an ideal tool to probe the multi-spin interactions addressed in this work.

The extraction of the MEI characterizing a magnetic material can be made via the celebrated infinitesimal rotation method of the magnetic moments [91–93], which defines a mapping procedure between the energy obtained from electronic structure calculations and the energy of an extended Heisenberg model. The magnetic force theorem [94, 95] permits to use the band energy instead of the total energy to evaluate the impact of rotating magnetic moments. The latter modifies the atomic potential by δV, which alters the density of states δn and the corresponding energy

$$\begin{equation}\delta E\enspace ={\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \left({\epsilon}-{{\epsilon}}_{\mathrm{F}}\right)\delta n\left({\epsilon}\right)=-{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \delta N\left({\epsilon}\right),\end{equation} \tag{ 1 }$$

where $\delta N\left({\epsilon}\right){\;\;=\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \delta n\left({\epsilon}\right)$ represents the change of the integrated local density of states up to the Fermi energy _F.

Since the density of states can be obtained from the imaginary part of the one-electron retarded Green function G() = ( − H)⁻¹ corresponding to the Hamiltonian H of the system:

$$\begin{equation}n\left({\epsilon}\right)=-\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}{\int }_{{{\Omega}}_{i}}\mathrm{d}\mathbf{r}\enspace {\mathrm{T}\mathrm{r}}_{s}\sum _{i}\enspace {G}_{ii}\left(\mathbf{r},\mathbf{r},{\epsilon}\right),\end{equation} \tag{ 2 }$$

where the integration is performed over the volume Ω_i defining the atomic site i and the trace is taken over the spin index s, one can use multiple-scattering theory as implemented in the Korringa–Kohn–Rostoker (KKR) Green-function method [96] to evaluate δn, δN and therefore δE. Because of the perturbation δV, the new green function g() = ( − H − δV), can be calculated from the original green function G using the Dyson equation: g() = G() + G()δVg(). The change in the Green function

$$\begin{equation}\delta g\left({\epsilon}\right)=g\left({\epsilon}\right)-G\left({\epsilon}\right)=G\left({\epsilon}\right)\delta VG\left({\epsilon}\right){\left(1-G\left({\epsilon}\right)\delta V\right)}^{-1}\end{equation} \tag{ 3 }$$

combined with the relation $\frac{\mathrm{d}G}{\mathrm{d}E}=-GG$ permits to find $\delta N\left({\epsilon}\right)=-\frac{1}{\pi }\enspace \mathrm{T}\mathrm{r}\enspace \mathrm{ln}\left[1-G\left({\epsilon}\right)\delta V\right]$, where the trace is taken over the position r, site i and spin index. Consequently, after perturbation, the modified energy reads

$$\begin{equation}\delta E=\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \mathrm{Tr}\enspace \mathrm{ln}\left[1-G\left({\epsilon}\right)\delta V\right],\end{equation} \tag{ 4 }$$

known as the Lloyd's formula [97, 98], leading to a power series expansion

$$\begin{equation}\delta E=-\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \sum _{p}\enspace \frac{1}{p}\enspace \mathrm{Tr}\enspace {\left[G\left({\epsilon}\right)\delta V\right]}^{p}\end{equation} \tag{ 5 }$$

given in a matrix notation. More details on the derivation can be found for example in references [29, 91–93].

In practice, the evaluation of the inter-site interactions requires the off-site part of the Green function instead of the full one. The off-site elements can be expanded as

$$\begin{equation}{G}_{ij,s{s}^{\prime }}\left(\mathbf{r},{\mathbf{r}}^{\prime };{\epsilon}\right)=\sum _{L{L}^{\prime }}\enspace {R}_{iLs}^{{\times}}\left(\boldsymbol{r};{\epsilon}\right){G}_{ij,L{L}^{\prime },s{s}^{\prime }}^{\text{str}}\left({\epsilon}\right){R}_{j{L}^{\prime }{s}^{\prime }}\left({\boldsymbol{r}}^{\prime };{\epsilon}\right)\end{equation} \tag{ 6 }$$

within multiple scattering theory as implemented in KKR. G^str is the structural Green function, which in contrast to the full Green function has no radial dependence but is site-dependent. It is also a matrix in angular orbital momentum L = (l, m) and spin representation, L ⊕ s. R and R^× are respectively the left- and right-hand wave functions resulting from the potential V without the perturbation. They are scattering solutions of either the Dirac or of the Schrödinger/scalar-relativistic equations augmented with the SOC. This allows one to utilize the single-site scattering t-matrix instead of the radial-dependent atomic potential V when evaluating the modified energy:

$$\begin{equation}\delta E=-\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \sum _{p}\enspace \frac{1}{p}\enspace \mathrm{T}\mathrm{r}{\left[{G}^{\text{str}}\left({\epsilon}\right)\delta t\left({\epsilon}\right)\right]}^{p},\end{equation} \tag{ 7 }$$

where the trace is taken this time over the sites i, orbital momentum and spin indices. Similarly to the structural Green function G^str, δt is a matrix in angular momentum and spin representation, L ⊕ s. Equation (7) is in practice more convenient than equation (5) since the latter hinges on cumbersome radially dependent quantities and radial integrations.

Adopting the rigid spin approximation implies a change of the magnetic part of the single-site t-matrix, which can be expressed as

$$\begin{equation}\delta {t}_{i}\left({\epsilon}\right)={t}_{i}^{\sigma }\left({\epsilon}\right)\enspace \delta {\mathbf{S}}_{i}\cdot \boldsymbol{\sigma }.\end{equation} \tag{ 8 }$$

Similarly to the proposal made in reference [93], ${t}_{i}^{\mathrm{s}}\left({\epsilon}\right)$ describes the scattering at the magnetic part of the atomic potential or of the exchange and correlation potential, B_xc, if one works in the framework of density functional theory: ${t}_{\begin{subarray}{c}i,L{L}^{\prime },s{s}^{\prime }\end{subarray}}^{\sigma }\left({\epsilon}\right){\enspace =\int }_{{{\Omega}}_{i}}\mathrm{d}\boldsymbol{r}{R}_{iLs}^{{\times}}\left(\boldsymbol{r};{\epsilon}\right){B}_{\mathrm{x}\mathrm{c}}\left(\boldsymbol{r}\right){R}_{i{L}^{\prime }{s}^{\prime }}\left(\boldsymbol{r};{\epsilon}\right)\enspace$.

The Green function can be decomposed into two Green functions, A that is non-magnetic and diagonal in spin space and B that contains the magnetic part and a contribution induced by the SOC:

$$\begin{equation}{G}_{ij}^{\text{str}}={A}_{ij}\enspace {\sigma }_{0}+{\mathbf{B}}_{ij}\cdot \boldsymbol{\sigma },\end{equation} \tag{ 9 }$$

where σ is the vector of Pauli matrices and σ₀ is the identity matrix. If inversion symmetry is not broken, A_ij is symmetric with respect to site exchange and likewise for ${\mathbf{B}}_{ij}^{\sigma }$ if SOC is not present. ${\mathbf{B}}_{ij}^{\sigma }$ could have a magnetic contribution, which should behave like A_ij, and a contribution from SOC. Within the Rashba model, where for example broken-inversion symmetry and SOC are incorporated, ${\mathbf{B}}_{ij}^{\sigma }=-{\mathbf{B}}_{ji}^{\sigma }$. Note that the decomposition of the Green functions as proposed by Cardias et al [99] is useful for the investigation of multi-spin excitations.

2.1. Two-spin interaction terms

Considering infinitesimal rotations, I proceed to a Taylor expansion of the logarithm in equation (7) and extract systematically high-order terms ordered according to powers of δt. As aforementioned terms with odd powers of δt cancel since they are not compatible with the time-reversal symmetry requirement of the total energy. The standard second order term

$$\begin{equation}\delta {E}_{2-\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=-\frac{1}{2\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \sum _{ij}\enspace {G}_{ij}^{\text{str}}\left({\epsilon}\right)\delta {t}_{j}\left({\epsilon}\right){G}_{ji}^{\text{str}}\left({\epsilon}\right)\delta {t}_{i}\left({\epsilon}\right)\end{equation} \tag{ 10 }$$

$$\begin{equation}=-\sum _{ij}\left[{J}_{ij}^{\text{iso}}\delta {\mathbf{S}}_{i}\cdot \delta {\mathbf{S}}_{j}+\delta {\mathbf{S}}_{i}\cdot \underline{\underline{J}}{}_{ij}^{\text{ani}}\cdot \delta {\mathbf{S}}_{j}+{\mathbf{D}}_{ij}\cdot \left(\delta {\mathbf{S}}_{i}{\times}\delta {\mathbf{S}}_{j}\right)\right].\end{equation} \tag{ 11 }$$

The integrand of equation (10) can be rewritten as

$$\begin{equation}\mathrm{Tr}\enspace {G}_{ij}^{\text{str}}\left({\epsilon}\right)\delta {t}_{j}\left({\epsilon}\right){G}_{ji}^{\text{str}}\left({\epsilon}\right)\delta {t}_{i}\left({\epsilon}\right)=\mathrm{Tr}\enspace \left[{A}_{ij}\enspace {\sigma }_{0}+{\mathbf{B}}_{ij}\cdot \boldsymbol{\sigma }\right]\enspace {t}_{j}^{\sigma }\left({\epsilon}\right)\enspace \delta {\mathbf{S}}_{j}\cdot \boldsymbol{\sigma }\enspace \left[{A}_{ji}\enspace {\sigma }_{0}+{\mathbf{B}}_{ji}\cdot \boldsymbol{\sigma }\right]\enspace {t}_{i}^{\sigma }\left({\epsilon}\right)\enspace \delta {\mathbf{S}}_{i}\cdot \boldsymbol{\sigma }.\end{equation} \tag{ 12 }$$

Here one can use known relations involving products with Pauli matrices: (a ⋅ σ)(b ⋅ σ) = (a ⋅ b) σ₀ + i(a × b) ⋅ σ to recover the interaction energy between two magnetic moments from which the bilinear MEI can be extracted:

$$\begin{equation}{J}_{ij}^{\text{iso}}=\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}{\mathcal{J}}_{ij}^{\text{iso}}\left({\epsilon}\right),\end{equation} \tag{ 13 }$$

$$\begin{equation}{\mathbf{D}}_{ij}=\frac{1}{\pi }\enspace \mathrm{R}\mathrm{e}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace {\boldsymbol{\mathcal{D}}}_{ij}\left({\epsilon}\right)\enspace \enspace \enspace ,\end{equation} \tag{ 14 }$$

$$\begin{equation}\delta {\mathbf{S}}_{i}\cdot \underline{\underline{J}}{}_{ij}^{\text{ani}}\cdot \delta {\mathbf{S}}_{j}=\frac{1}{\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\delta {\mathbf{S}}_{i}\cdot \underline{\underline{\mathcal{J}}}{}_{ij}^{\text{ani}}\enspace \enspace \left({\epsilon}\right)\cdot \delta {\mathbf{S}}_{j}.\end{equation} \tag{ 15 }$$

The integrands are energy-dependent products of matrices of the same size than the inter-site Green functions and are given by the following forms:

$$\begin{equation}{\mathcal{J}}_{ij}^{\text{iso}}\left({\epsilon}\right)={A}_{ij}\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right){A}_{ji}\left({\epsilon}\right){t}_{i}^{\sigma }\left({\epsilon}\right)-\left[{\boldsymbol{B}}_{ij}^{\sigma }\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right)\right]\cdot \left[{\boldsymbol{B}}_{ji}^{\sigma }\left({\epsilon}\right){t}_{i}^{\sigma }\left({\epsilon}\right)\right],\end{equation} \tag{ 16 }$$

$$\begin{equation}{\boldsymbol{\mathcal{D}}}_{ij}\left({\epsilon}\right)={A}_{ij}\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right){\boldsymbol{B}}_{ji}^{\sigma }\left({\epsilon}\right){t}^{\sigma }i\left({\epsilon}\right)-{\boldsymbol{B}}_{ij}^{\sigma }\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right){A}_{ji}\left({\epsilon}\right){t}_{i}^{\sigma }\left({\epsilon}\right)\enspace \enspace \end{equation} \tag{ 17 }$$

$$\begin{equation}\delta {\mathbf{S}}_{i}\cdot \underline{\underline{\mathcal{J}}}{}_{ij}^{\text{ani}}\enspace \enspace \left({\epsilon}\right)\cdot \delta {\mathbf{S}}_{j}=\left[{\boldsymbol{B}}_{ij}^{\sigma }\left({\epsilon}\right)\cdot \delta {\mathbf{S}}_{j}\right]{t}_{j}^{\sigma }\left({\epsilon}\right)\left[{\boldsymbol{B}}_{ji}^{\sigma }\left({\epsilon}\right)\cdot \delta {\mathbf{S}}_{i}\right]{t}_{i}^{\sigma }\left({\epsilon}\right)+\left[{\boldsymbol{B}}_{ij}^{\sigma }\left({\epsilon}\right)\cdot \delta {\mathbf{S}}_{i}\right]{t}_{j}^{\sigma }\left({\epsilon}\right)\left[{\boldsymbol{B}}_{ji}^{\sigma }\left({\epsilon}\right)\cdot \delta {\mathbf{S}}_{j}\right]{t}_{i}^{\sigma }\left({\epsilon}\right),\end{equation} \tag{ 18 }$$

where the decomposition of the Green functions introduced in equation (9) leads to forms different from but still fully equivalent to the usually presented ones [92, 93]. Indeed, the elements of the tensor of MEI are conventionally obtained from equation (10) and then combined to calculate the isotropic MEI, DMI and the compass term. Interestingly, this decomposition permits to identify two contributions to J: one is independent of magnetism and SOC while the other one depends on both with at least a quadratic dependence on SOC but being rotationally invariant in spin space owing to the inner product of the different B. D is, however, linear in B and consequently in SOC. The direction of the DMI vector is dictated by B and for lattices with inversion symmetry, I naturally recover that the DMI vanishes because of the symmetry of A and B with respect to site exchange.

2.2. Four-spin interaction terms

The fourth-order term in the Taylor expansion of the energy,

$$\begin{equation}\delta {E}_{4-\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=-\frac{1}{4\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \sum _{ijkl}\enspace {G}_{ij}^{\text{str}}\left({\epsilon}\right)\delta {t}_{j}\left({\epsilon}\right){G}_{ji}^{\text{str}}\left({\epsilon}\right)\delta {t}_{k}\left({\epsilon}\right){G}_{kl}^{\text{str}}\left({\epsilon}\right)\delta {t}_{l}\left({\epsilon}\right){G}_{li}^{\text{str}}\left({\epsilon}\right)\delta {t}_{i}\left({\epsilon}\right),\end{equation} \tag{ 19 }$$

involves plaquettes of four magnetic moments. One notes that, in general, the four-spin contracted-site interaction terms can be recovered from the aforementioned forms by equaling some of the indices, for example by replacing the fourth site l with j to get the three-site terms, while replacing both (k, l) with (i, j) leads to the biquadratic terms. Note that this can lead to vanishing of some of the derived terms and the recovery of bilinear-like terms.

Isotropic interaction. Analogous to the two-spin interaction, E_4‐spin gives rise to the conventional isotropic four-spin interactions proportional to (S_i ⋅ S_j)(S_k ⋅ S_l) including quadratic and fourth order contributions of the SOC. The isotropic change of the energy reads

$$\begin{equation}\delta {E}_{ijkl}^{\text{iso}}=-{\mathcal{B}}_{ijkl}\left[\left(\delta {\mathbf{S}}_{i}\cdot \delta {\mathbf{S}}_{j}\right)\left(\delta {\mathbf{S}}_{k}\cdot \delta {\mathbf{S}}_{l}\right)-\left(\delta {\mathbf{S}}_{i}{\times}\delta {\mathbf{S}}_{j}\right)\left(\delta {\mathbf{S}}_{k}{\times}\delta {\mathbf{S}}_{l}\right)\right]\end{equation} \tag{ 20 }$$

with which one recovers the usual form of the isotropic four-spin energy

$$\begin{equation}{E}_{ijkl}^{\text{iso}}=-{\mathcal{B}}_{ijkl}\left[\left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{l}\right)-\left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{k}\right)\left(\delta {\mathbf{S}}_{j}\cdot {\mathbf{S}}_{l}\right)+\left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{l}\right)\left({\mathbf{S}}_{j}\cdot {\mathbf{S}}_{k}\right)\right],\end{equation} \tag{ 21 }$$

with

$$\begin{equation}{\mathcal{B}}_{ijkl}=\frac{1}{2\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\left\{{A}_{ij}\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right){A}_{jk}\left({\epsilon}\right){t}_{k}^{\sigma }\left({\epsilon}\right){A}_{kl}\left({\epsilon}\right){t}_{l}^{\sigma }\left({\epsilon}\right){A}_{li}\left({\epsilon}\right){t}_{i}^{\sigma }\left({\epsilon}\right)\right\},\end{equation} \tag{ 22 }$$

in zero-order with respect to B.

Four-spin chiral interactions. The fourth-order term in the Taylor expansion of the energy gives also rise to terms linear in the SOC. The latter are the recently derived [37] and postulated [41] four-spin vector-chiral interactions proportional to (S_i ⋅ S_j)(S_k × S_l) and to the scalar spin-chirality S_i ⋅ (S_j × S_k) [42]. Besides these two terms, one additional term proportional to the three-spin vector chirality S_i × (S_j × S_k) shows up. We note that the different terms can be rewritten as function of each others. The product of the four Green-function elements in equation (19) and the requested linearity in B leads finally to four terms of the type BAAA, where B can be placed at four different positions relative to A. For instance, BAAA can be written in various ways depending on how the products involving the Pauli matrices are grouped:

$$\begin{equation}\delta {E}_{ijkl}^{\mathbf{B}AAA}=-{\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}\cdot \left[\left(\delta {\mathbf{S}}_{j}\cdot \delta {\mathbf{S}}_{k}\right)\left(\delta {\mathbf{S}}_{l}{\times}\delta {\mathbf{S}}_{i}\right)+\left(\delta {\mathbf{S}}_{j}{\times}\delta {\mathbf{S}}_{k}\right)\left(\delta {\mathbf{S}}_{l}\cdot \delta {\mathbf{S}}_{i}\right)+\left(\delta {\mathbf{S}}_{j}{\times}\delta {\mathbf{S}}_{k}\right){\times}\left(\delta {\mathbf{S}}_{l}{\times}\delta {\mathbf{S}}_{i}\right)\right]\end{equation} \tag{ 23 }$$

or

$$\begin{align}\hfill \delta {E}_{ijkl}^{\mathbf{B}AAA}& =-\left({\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}\cdot \delta {\mathbf{S}}_{j}\right)\left[\left(\delta {\mathbf{S}}_{k}{\times}\delta {\mathbf{S}}_{l}\right)\cdot \delta {\mathbf{S}}_{i}\right]-\left({\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}{\times}\delta {\mathbf{S}}_{j}\right)\left[\left(\delta {\mathbf{S}}_{k}\cdot \delta {\mathbf{S}}_{l}\right)\cdot \delta {\mathbf{S}}_{i}\right]\hfill \\ \hfill & \quad +\left({\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}{\times}\delta {\mathbf{S}}_{j}\right)\cdot \left[\left(\delta {\mathbf{S}}_{k}{\times}\delta {\mathbf{S}}_{l}\right){\times}\delta {\mathbf{S}}_{i}\right],\hfill \end{align} \tag{ 24 }$$

with

$$\begin{equation}{\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}=\frac{1}{2\pi }\enspace \mathrm{R}\mathrm{e}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace {\mathbf{B}}_{ij}\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right){A}_{jk}\left({\epsilon}\right){t}_{k}^{\sigma }\left({\epsilon}\right){A}_{kl}\left({\epsilon}\right){t}_{l}^{\sigma }\left({\epsilon}\right){A}_{li}\left({\epsilon}\right){t}_{i}^{\sigma }\left({\epsilon}\right),\end{equation} \tag{ 25 }$$

where the bars on the indices indicate the sites connected by B. It is interesting to note that the integrand showing up in the previous equations can be interpreted as products involving the bilinear terms ${\mathcal{J}}^{\text{iso}}$ and $\boldsymbol{\mathcal{D}}$ (see equations (16)–(18)).

In equation (25), one can identify the three-spin scalar chirality, χ_ijk = S_i ⋅ (S_j × S_k), and if I introduce what I call a three-spin vector chirality ℵ_ijk = S_i × (S_j × S_k), the four-spin chiral energy linear in B reads

$$\begin{align}\hfill {E}_{ijkl}^{\mathbf{B}}& =-{\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}\cdot \left[{\mathbf{S}}_{j}\enspace {\chi }_{ikl}-\left({\mathbf{S}}_{i}{\times}{\mathbf{S}}_{j}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{l}\right)+{\mathbf{S}}_{j}{\times}{\boldsymbol{\aleph }}_{ikl}\right]-{\boldsymbol{\mathcal{C}}}_{i\overline{j}\overline{k}l}\cdot \left[{\mathbf{S}}_{k}\enspace {\chi }_{jli}-\left({\mathbf{S}}_{j}{\times}{\mathbf{S}}_{k}\right) \left.\right]\left({\mathbf{S}}_{l}\cdot {\mathbf{S}}_{i}\right)+{\mathbf{S}}_{k}{\times}{\boldsymbol{\aleph }}_{jli}\right]\hfill \\ \hfill & \quad -\enspace {\boldsymbol{\mathcal{C}}}_{ij\overline{k}\overline{l}}\cdot \left[{\mathbf{S}}_{l}\enspace {\chi }_{kij}-\left({\mathbf{S}}_{k}{\times}{\mathbf{S}}_{l}\right)\left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}\right)+{\mathbf{S}}_{l}{\times}{\boldsymbol{\aleph }}_{kij}\right]-{\boldsymbol{\mathcal{C}}}_{\bar{i}jk\overline{l}}\cdot \left[{\mathbf{S}}_{i}\enspace {\chi }_{ljk}-\left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\left({\mathbf{S}}_{j}\cdot {\mathbf{S}}_{k}\right)+{\mathbf{S}}_{i}{\times}{\boldsymbol{\aleph }}_{ljk}\right].\hfill \end{align} \tag{ 26 }$$

This equation indicates that the ground state favored by such a four-spin chiral energy is the result of a complex competition of various terms: the three-spin scalar chirality, the three-spin vector chirality and a mix of cross and dot products. For instance, the three-spin scalar chirality pertaining to an equilateral triangular plaquette reaches its largest value for a polar angle $\mathrm{cos}\left(\theta \right)=1/\sqrt{3}$ while the largest length of the three-spin vector chirality is found for $\mathrm{cos}\left(\theta \right)=\sqrt{7}/3$. By focusing on the first line of equation (26), one notices that the inner product ${\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}\cdot {\boldsymbol{S}}_{j}$ favors a parallel alignment between the four-spin chiral interaction vector and the magnetic moment j. The cross products S_i × S_j, however, drives the spin-configuration towards a state where the two perpendicular magnetic moments S_j and S_i live in a plane perpendicular to ${\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}$. The third term favors the chiral vector and the three magnetic moments at sites j, k, l to lie in the same plane perpendicular to S_i with S_k and S_l perpendicular to each other while ${\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}$ and S_i would be parallel. Although, I present the four-spin energy as a sum of three types of terms, as shown in reference [43] one can rewrite all of them as function of the latter one.

While E^B contains the first-order contributions in terms of B to E_4‐spin, second-order contributions give rise to terms correcting the dot-product and cross-product of the magnetic moments:

$$\begin{equation}{E}_{ijkl}^{\mathbf{B}\mathbf{B}}={E}_{ijkl}^{\mathbf{B}\mathbf{B}AA}+{E}_{ijkl}^{A\mathbf{B}\mathbf{B}A}+{E}_{ijkl}^{AA\mathbf{B}\mathbf{B}}+{E}_{ijkl}^{A\mathbf{B}A\mathbf{B}}+{E}_{ijkl}^{\mathbf{B}A\mathbf{B}A}+{E}_{ijkl}^{\mathbf{B}AA\mathbf{B}}\end{equation} \tag{ 27 }$$

As it can be noticed in the following example, where I provide the contribution arising from BBAA, one finds products of bilinear terms similar to ${\mathcal{J}}^{\text{ani}}$ multiplying J^iso or $\boldsymbol{\mathcal{D}}$ as defined in equations (16)–(18):

$$\begin{align}\hfill {E}_{ijkl}^{\mathbf{B}\mathbf{B}AA}& =-\frac{1}{2\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\left\{\right. \left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot {\mathbf{S}}_{j}\right)\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot {\mathbf{S}}_{k}\right)\left({\mathbf{S}}_{l}\cdot {\mathbf{S}}_{i}\right)\hfill \\ \hfill & \quad -\left[\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }{\times}{\mathbf{S}}_{j}\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right]\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot {\mathbf{S}}_{k}\right)-\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot {\mathbf{S}}_{j}\right)\left[\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }{\times}{\mathbf{S}}_{k}\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right]\hfill \\ \hfill & \quad +\left[\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }{\times}{\mathbf{S}}_{j}\right){\times}\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }{\times}{\mathbf{S}}_{k}\right)\right]\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right) \left.\right\}{A}_{kl}{t}_{l}^{\sigma }{A}_{li}{t}_{i}^{\sigma }.\hfill \end{align} \tag{ 28 }$$

The first rhs term looks like the anisotropic bilinear term (equation (18)) weighted by a dot product between two magnetic moments of the four-spin plaquette. The second and the third terms can be rewritten in a form similar to the first one:

$$\begin{equation}\left[\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }{\times}{\mathbf{S}}_{j}\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right]\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot {\mathbf{S}}_{k}\right)=\left[\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot {\mathbf{S}}_{l}\right)\left({\mathbf{S}}_{j}\cdot {\mathbf{S}}_{i}\right)-\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot {\mathbf{S}}_{i}\right)\left({\mathbf{S}}_{j}\cdot {\mathbf{S}}_{l}\right)\right]\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot {\mathbf{S}}_{k}\right)\end{equation} \tag{ 29 }$$

and

$$\begin{equation}\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot {\mathbf{S}}_{j}\right)\left[\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }{\times}{\mathbf{S}}_{k}\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right]=\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot {\mathbf{S}}_{j}\right)\left[\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot {\mathbf{S}}_{l}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{i}\right)-\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot {\mathbf{S}}_{i}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{l}\right)\right]\end{equation} \tag{ 30 }$$

while the fourth term can be written as

$$\begin{align}\hfill \left[\left({\mathbf{B}}_{ij}{t}_{j}^{\sigma }{\times}{\mathbf{S}}_{j}\right){\times}\left({\mathbf{B}}_{jk}{t}_{k}^{\sigma }{\times}{\mathbf{S}}_{k}\right)\right]\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)& =\left[{\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot \left({\mathbf{S}}_{j}{\times}{\mathbf{S}}_{l}\right)\right]\left[{\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot \left({\mathbf{S}}_{k}{\times}{\mathbf{S}}_{i}\right)\right]\hfill \\ \hfill & \quad -\left[{\mathbf{B}}_{ij}{t}_{j}^{\sigma }\cdot \left({\mathbf{S}}_{j}{\times}{\mathbf{S}}_{i}\right)\right]\left[{\mathbf{B}}_{jk}{t}_{k}^{\sigma }\cdot \left({\mathbf{S}}_{k}{\times}{\mathbf{S}}_{l}\right)\right].\hfill \end{align} \tag{ 31 }$$

2.3. Six-spin interaction term

The sixth order term obtained from equation (7) involves a plaquette of six magnetic moments. One expects the largest contribution to be the isotropic, rotationally invariant one, as this term can be finite without spin–orbit interaction. Although being isotropic, one finds contributions involving products of three-spin scalar chiralities, or of three-spin vector chiralities:

$$\begin{align}\hfill \delta {E}_{ijklmn}^{\text{iso}}& =-{\varkappa }_{ijklmn}^{6-\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\left[\left(\delta {\mathbf{S}}_{i}\cdot \delta {\mathbf{S}}_{j}\right)\left(\delta {\mathbf{S}}_{k}\cdot \delta {\mathbf{S}}_{l}\right)\left(\delta {\mathbf{S}}_{m}\cdot \delta {\mathbf{S}}_{n}\right)-\left(\delta {\mathbf{S}}_{i}\cdot \delta {\mathbf{S}}_{j}\right)\left(\delta {\mathbf{S}}_{k}{\times}\delta {\mathbf{S}}_{l}\right)\cdot \left(\delta {\mathbf{S}}_{m}{\times}\delta {\mathbf{S}}_{n}\right)\right.\hfill \\ \hfill & \left.\quad -\left(\delta {\mathbf{S}}_{i}{\times}\delta {\mathbf{S}}_{j}\right)\left(\delta {\mathbf{S}}_{k}\cdot \delta {\mathbf{S}}_{l}\right)\cdot \left(\delta {\mathbf{S}}_{m}{\times}\delta {\mathbf{S}}_{n}\right)-\left(\delta {\mathbf{S}}_{i}{\times}\delta {\mathbf{S}}_{j}\right)\cdot \left(\delta {\mathbf{S}}_{k}{\times}\delta {\mathbf{S}}_{l}\right)\cdot \left(\delta {\mathbf{S}}_{m}\cdot \delta {\mathbf{S}}_{n}\right)\right],\hfill \end{align} \tag{ 32 }$$

and thus

$$\begin{align}\hfill {E}_{ijklmn}^{\text{iso}}& =-{\varkappa }_{ijklmn}^{6-\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}\left[\right. \left({\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}\right)\left({\mathbf{S}}_{k}\cdot {\mathbf{S}}_{l}\right)\left({\mathbf{S}}_{m}\cdot {\mathbf{S}}_{n}\right)-{\chi }_{imn}{\chi }_{jkl}-{\chi }_{kmn}{\chi }_{lij}-{\chi }_{mkl}{\chi }_{nij}\hfill \\ \hfill & \quad -\enspace {\boldsymbol{\aleph }}_{ikl}\cdot {\boldsymbol{\aleph }}_{jmn}-{\boldsymbol{\aleph }}_{kij}\cdot {\boldsymbol{\aleph }}_{lmn}-{\boldsymbol{\aleph }}_{mij}\cdot {\boldsymbol{\aleph }}_{nkl} \left.\right],\hfill \end{align} \tag{ 33 }$$

with χ_ijk = S_i ⋅ (S_j × S_k) and ℵ_ijk = S_i × (S_j × S_k). In general, the effective six-spin interaction is given by

$$\begin{equation}{\varkappa }_{ijklmn}^{6-\mathrm{s}\mathrm{p}\mathrm{i}\mathrm{n}}=\frac{1}{3\pi }\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace {A}_{ij}\left({\epsilon}\right){t}_{j}^{\sigma }\left({\epsilon}\right){A}_{jk}\left({\epsilon}\right){t}_{k}^{\sigma }\left({\epsilon}\right){A}_{kl}\left({\epsilon}\right){t}_{l}^{\sigma }\left({\epsilon}\right){A}_{lm}\left({\epsilon}\right){t}_{m}^{\sigma }\left({\epsilon}\right){A}_{mn}\left({\epsilon}\right){t}_{n}^{\sigma }\left({\epsilon}\right){A}_{ni}\left({\epsilon}\right){t}_{i}^{\sigma }\left({\epsilon}\right).\end{equation} \tag{ 34 }$$

Here I discuss the derived multi-spin interactions considering localized spins embedded in 1D and 2D electron gas, with SOC added in the form of Rashba coupling [80, 81]. The conduction electrons are then either confined along a given direction (1D) or in the xy-plane (2D). The effective electric field emerging from SOC points along the z-axis. The Hamiltonian of both systems consists of the kinetic energy of the free electrons to which a linear term in momentum $\overrightarrow {p}$ is added:

$$\begin{equation}\boldsymbol{H}=\frac{{\hslash }^{2}}{2{m}^{{\ast}}}{\mathbf{\nabla }}^{2}\enspace {\boldsymbol{\sigma }}_{0}-\alpha \left(-\mathrm{i}\hslash \mathbf{\nabla }{\times}\overline{z}\right)\cdot \boldsymbol{\sigma },\end{equation} \tag{ 35 }$$

where α is the so-called Rashba parameter representing the strength of the spin–orbit interaction, is the unit vector along the z-axis, m* is the effective mass of the electron, σ₀ is the identity matrix and σ is the vector of Pauli spin matrices.

Similarly to reference [72], the potential of each of the localized spins is described via the s − d interaction J^sd, which has units of energy multiplying length to the power of dimensionality, such that the change of the magnetic part of the potential upon rotation of the spin moment is given by

$$\begin{equation}\delta V\left(\boldsymbol{r}\right)=\sum _{i}{J}_{i}^{sd}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\right)\enspace {\mathbf{S}}_{i}\cdot \sigma .\end{equation} \tag{ 36 }$$

In the following, I use the s − d potential instead of the energy-dependent t-matrix in order to get analytical forms of the multi-spin interactions similarly to what is done in the RKKY approximation. This means that I will use the expansion given by equation (5) based on the full Green function instead of equation (7), where the structural Green function is utilized.

3.1. Rashba model—one dimension

The Rashba Green function corresponding to the previous Hamiltonian for the 1DEG reads

$$\begin{equation}G\left({\boldsymbol{r}}_{ij};\varepsilon +\mathrm{i}{0}^{+}\right)={A}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)+{B}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)\enspace {\hat{\beta }}_{ij}\cdot \boldsymbol{\sigma },\end{equation} \tag{ 37 }$$

where r_ij = r_ij(cos β_ij, sin β_ij) is the vector connecting sites i and j, while ${\hat{\beta }}_{ij}=\left(\mathrm{sin}{\beta }_{ij},-\mathrm{cos}{\beta }_{ij}\right)$ is the unit vector perpendicular to r_ij. Defining $q=\sqrt{\frac{2m}{{\hslash }^{2}}\varepsilon +{k}_{\mathrm{R}}^{2}}$ with ${k}_{\mathrm{R}}=\frac{{m}^{{\ast}}}{{\hslash }^{2}}\alpha$, the components of the Green functions are [72]:

$$\begin{equation}{A}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)=-\mathrm{i}\frac{{m}^{{\ast}}}{{\hslash }^{2}}\frac{{\mathrm{e}}^{\mathrm{i}q\vert {r}_{ij}\vert }}{q+\mathrm{i}{0}^{+}}\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{ij}\right)\quad \text{and}\quad {B}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)=\frac{{m}^{{\ast}}}{{\hslash }^{2}}\frac{{\mathrm{e}}^{\mathrm{i}q\vert {r}_{ij}\vert }}{q+\mathrm{i}{0}^{+}}\enspace \mathrm{sin}\left({k}_{\mathrm{R}}{r}_{ij}\right).\end{equation} \tag{ 38 }$$

Now, I can proceed with the multi-spin interactions by plugging in A and B into the MEI derived in section 2. I start with the basic bilinear types [72, 74].

3.1.1. Two-spin interactions

Within the 1DEG Rashba model, utilizing equations (13)–(18), I recover the bilinear magnetic interactions that were already derived in reference [72]:

$$\begin{equation}{J}_{ij}^{\text{iso}}={\mathcal{F}}_{ij}^{\left(1\text{D}\right)}\enspace \mathrm{cos}\left(2{k}_{\mathrm{R}}{r}_{ij}\right),\end{equation} \tag{ 39 }$$

$$\begin{equation}\delta {\mathbf{S}}_{i}\cdot \underline{\underline{J}}{}_{ij}^{\text{ani}}\cdot \delta {\mathbf{S}}_{j}={\mathcal{F}}_{ij}^{\left(1\text{D}\right)}\left(1-\mathrm{cos}\left(2{k}_{\mathrm{R}}{r}_{ij}\right)\right)\left[{\hat{\beta }}_{ij}\cdot \delta {\mathbf{S}}_{j}\right]\left[{\hat{\beta }}_{ij}\cdot \delta {\mathbf{S}}_{i}\right],\end{equation} \tag{ 40 }$$

$$\begin{equation}{\mathbf{D}}_{ij}={\mathcal{F}}_{ij}^{\left(1\text{D}\right)}\enspace \mathrm{sin}\left(2{k}_{\mathrm{R}}{r}_{ij}\right){\hat{\beta }}_{ij}.\end{equation} \tag{ 41 }$$

Within the RKKY approximation, the direction of the DMI is along ${\hat{\beta }}_{ij}$, i.e. perpendicular to the bond connecting the two sites i and j. Naturally, without SOC, D and ${\underline{\underline{J}}}^{\text{ani}}$ vanish. For small k_Rr_ij, the DMI is linear with the Rashba parameter α while the compass-term interaction shows a quadratic dependence similarly to the SOC correction to the isotropic MEI. These dependencies change at large k_Rr_ij, where the various bilinear interactions can be of the same order of magnitude with a decay dictated by the range function ${\mathcal{F}}_{ij}^{\left(1\text{D}\right)}$, which is in fact the isotropic bilinear MEI before application of the SOC [100, 101]:

$$\begin{align}\hfill {\mathcal{F}}_{ij}^{\left(1\text{D}\right)}& =-{J}_{i}^{sd}{J}_{j}^{sd}\frac{{m}^{{\ast}2}}{\pi {\hslash }^{4}}\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\enspace \frac{{\mathrm{e}}^{\mathrm{i}2q{r}_{ij}}}{{\left(q+\mathrm{i}{0}^{+}\right)}^{2}}\hfill \\ \hfill & =-{J}_{i}^{sd}{J}_{j}^{sd}\frac{{m}^{{\ast}}}{\pi {\hslash }^{2}}\left[\text{Si}\left(2{q}_{\mathrm{F}}{r}_{ij}\right)-\frac{\pi }{2}\right],\hfill \end{align} \tag{ 42 }$$

where Si() is the sine integral function. Owing to the asymptotic behavior of the sine function: $Si\left(x\right)\simeq \frac{\pi }{2}-\frac{\mathrm{cos}\left(x\right)}{x}\left(1-\frac{2!}{{x}^{2}}+\frac{4!}{{x}^{4}}-\frac{6!}{{x}^{6}}+\dots \right)-\frac{\mathrm{sin}\left(x\right)}{x}\left(\frac{1}{x}-\frac{3!}{{x}^{3}}+\frac{5!}{{x}^{5}}-\frac{7!}{{x}^{7}}+\dots \right)$, the range function behaves like:

$$\begin{equation}{\mathcal{F}}_{ij}^{\left(1\text{D}\right)}\simeq {J}_{i}^{sd}{J}_{j}^{sd}\frac{{m}^{{\ast}}}{\pi {\hslash }^{2}}\left[\frac{\mathrm{cos}\left(2{q}_{\mathrm{F}}{r}_{ij}\right)}{2{q}_{\mathrm{F}}{r}_{ij}}\left(1-\frac{1}{2{q}_{\mathrm{F}}^{2}{r}_{ij}^{2}}\right)+\frac{\mathrm{sin}\left(2{q}_{\mathrm{F}}{r}_{ij}\right)}{4{q}_{\mathrm{F}}^{2}{r}_{ij}^{2}}\right].\end{equation} \tag{ 43 }$$

Thus, the bilinear interactions are characterized by various decays, ${r}_{ij}^{-1}$ being the lowest one, but interference effects are expected with the term proportional to ${r}_{ij}^{-2}$. The wave length of the oscillations is expected to be complex depending on the interatomic distances. It is given by 2λ_F = π/q_F at small q_Fr before interference effects kick in, which originate from either or both the 1D nature of the electron gas (equation (43)) and from SOC (equations (39)–(41)).

3.1.2. Four-spin interactions

The isotropic four-spin interaction. As derived in equation (22), this interaction involves a product of four cosine functions that depend on k_R:

$$\begin{equation}{\mathcal{B}}_{ijkl}={\mathcal{F}}_{ijkl}^{\left(1\text{D}\right)}\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{ij}\right)\mathrm{cos}\left({k}_{\mathrm{R}}{r}_{jk}\right)\mathrm{cos}\left({k}_{\mathrm{R}}{r}_{kl}\right)\mathrm{cos}\left({k}_{\mathrm{R}}{r}_{li}\right).\end{equation} \tag{ 44 }$$

For small k_Rr, the isotropic four-spin interaction experiences corrections with even power of the Rashba coupling parameter α. The lowest dependence being quadratic in this regime, this implies a SOC contribution to the magnetic energy that is similar to the isotropic bilinear contribution. Also, and in analogy to the latter term, a range function takes care of the distance-dependent behavior of the four-spin interaction:

$$\begin{align}\hfill {\mathcal{F}}_{ijkl}^{\left(1\text{D}\right)}& ={J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{m}^{{\ast}4}}{2\pi {\hslash }^{8}}\enspace \mathrm{I}\mathrm{m}\enspace \mathrm{T}\mathrm{r}{\int }_{-\infty }^{{{\epsilon}}_{\mathrm{F}}} \mathrm{d}{\epsilon}\frac{{\mathrm{e}}^{\mathrm{i}q{L}_{ijkl}}}{{\left(q+\mathrm{i}{0}^{+}\right)}^{4}}\hfill \\ \hfill & =-{J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{m}^{{\ast}3}}{2\pi {\hslash }^{6}}\left[\frac{1}{2}\frac{\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}^{2}}+\frac{{L}_{ijkl}}{2}\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}}+\frac{{L}_{ijkl}^{2}}{2}\left(\text{Si}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)-\frac{\pi }{2}\right)\right].\hfill \end{align} \tag{ 45 }$$

Here I define L_ijkl = r_ij + r_jk + r_kl + r_li as the perimeter of the 1D plaquette (ijkl). Using once more the asymptotic form of the sine integral function, the long-range behavior of ${\mathcal{F}}_{ijkl}^{\left(1D\right)}$ simplifies to:

$$\begin{equation}{\mathcal{F}}_{ijkl}^{\left(1\text{D}\right)}\simeq -{J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{m}^{{\ast}3}}{2\pi {\hslash }^{6}}\left[\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}^{3}{L}_{ijkl}}+3\frac{\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}^{4}{L}_{ijkl}^{2}}\right],\end{equation} \tag{ 46 }$$

which indicates that the decay of the interaction is similar to the bilinear ones and is related to the distance spanned by the scattered electrons involved in the processes giving rise to the MEI. Depending on q_F and on the magnitude of the s − d interaction, the isotropic four-spin interaction can be of the same order of magnitude than the bilinear MEI. Compared to the latter, the wavelength of the oscillations decreases for the four-spin interaction. Assuming magnetic moments equidistant by R, the range function becomes:

$$\begin{align}\hfill {\mathcal{F}}_{ijkl}^{\left(1\text{D}\right)}& \simeq -{J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{m}^{{\ast}3}}{2\pi {\hslash }^{6}}\left[\frac{1}{2}\frac{\mathrm{sin}\left(4{q}_{\mathrm{F}}R\right)}{{q}_{\mathrm{F}}^{2}}+2R\frac{\mathrm{cos}\left(4{q}_{\mathrm{F}}R\right)}{{q}_{\mathrm{F}}}+8{R}^{2}\left(\text{Si}\left(4{q}_{\mathrm{F}}R\right)-\frac{\pi }{2}\right)\right]\hfill \\ \hfill & \quad -\enspace {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{m}^{{\ast}3}}{2\pi {\hslash }^{6}}\left[\frac{\mathrm{cos}\left(4{q}_{\mathrm{F}}R\right)}{4{q}_{\mathrm{F}}^{3}R}+3\frac{\mathrm{sin}\left(4{q}_{\mathrm{F}}R\right)}{16{q}_{\mathrm{F}}^{4}{R}^{2}}\right],\hfill \end{align} \tag{ 47 }$$

The four-spin chiral interaction. The chiral interaction that is linear in B has the following form

$$\begin{equation}{\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}=-{\mathcal{F}}_{ijkl}^{\left(1\text{D}\right)}\enspace \mathrm{sin}\left({k}_{\mathrm{R}}{r}_{ij}\right)\mathrm{cos}\left({k}_{\mathrm{R}}{r}_{jk}\right)\mathrm{cos}\left({k}_{\mathrm{R}}{r}_{kl}\right)\mathrm{cos}\left({k}_{\mathrm{R}}{r}_{li}\right){\hat{\beta }}_{ij},\end{equation} \tag{ 48 }$$

whose direction is perpendicular to the bond connecting the sites mediating the Green function B similarly to the DMI vector. Note that this is particular to the RKKY approximation based on the Rashba model since in general, the direction of ${\boldsymbol{\mathcal{C}}}_{\bar{i}\overline{j}kl}$ can be more complex as demonstrated in reference [43]. In comparison to the isotropic four-spin interaction, one of the cosines is replaced by a sine function. For small k_Rr, similarly to the DMI a linear dependence with respect to the Rashba coupling parameter α is expected for relatively short distances. Therefore, increasing SOC would permit to increase the amplitude of such chiral four-spin interaction. Since it is quadratic in B, the anisotropic chiral four-spin interaction involves two sine functions depending on α and therefore the short-distance behavior is quadratic with respect to SOC:

$$\begin{equation}{\mathcal{C}}_{\bar{i}\overline{j}\overline{k}l}^{\text{ani}}=-{\mathcal{F}}_{ijkl}^{\left(1\text{D}\right)}\enspace \mathrm{sin}\left({k}_{\mathrm{R}}{r}_{ij}\right)\enspace \mathrm{sin}\left({k}_{\mathrm{R}}{r}_{jk}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{kl}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{li}\right).\end{equation} \tag{ 49 }$$

This parameterizes the energy contributions such as the one given in equation (28), which in the Rashba model reads:

$$\begin{align}\hfill {E}_{ijkl}^{\mathbf{B}\mathbf{B}AA}& =-{\mathcal{C}}_{\bar{i}\overline{j}\overline{k}l}^{\text{ani}}\left\{\left({\hat{\beta }}_{ij}{t}_{j}^{s}\cdot {\mathbf{S}}_{j}\right)\left({\hat{\beta }}_{jk}{t}_{k}^{s}\cdot {\mathbf{S}}_{k}\right)\left({\mathbf{S}}_{l}\cdot {\mathbf{S}}_{i}\right)\right.\hfill \\ \hfill & \quad -\left(\left({\hat{\beta }}_{ij}{t}_{j}^{s}{\times}{\mathbf{S}}_{j}\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right)\left({\hat{\beta }}_{jk}{t}_{k}^{s}\cdot {\mathbf{S}}_{k}\right)-\left({\hat{\beta }}_{ij}{t}_{j}^{s}\cdot {\mathbf{S}}_{j}\right)\left(\left({\hat{\beta }}_{jk}{t}_{k}^{s}{\times}{\mathbf{S}}_{k}\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right)\hfill \\ \hfill & \quad +\left.\left(\left({\hat{\beta }}_{ij}{t}_{j}^{s}{\times}{\mathbf{S}}_{j}\right){\times}\left({\hat{\beta }}_{jk}{t}_{k}^{s}{\times}{\mathbf{S}}_{k}\right)\right)\cdot \left({\mathbf{S}}_{l}{\times}{\mathbf{S}}_{i}\right)\right\}{A}_{kl}{t}_{l}^{s}{A}_{li}{t}_{i}^{s}.\hfill \end{align} \tag{ 50 }$$

3.1.3. Six-spin interactions

The isotropic six-spin interaction is given by

$$\begin{equation}{\varkappa }_{ijklmn}^{6-\text{spin}}={\mathcal{F}}_{ijklmn}^{\left(1\text{D}\right)}\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{ij}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{jk}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{kl}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{lm}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{mn}\right)\enspace \mathrm{cos}\left({k}_{\mathrm{R}}{r}_{ni}\right),\end{equation} \tag{ 51 }$$

where

$$\begin{align}\hfill {\mathcal{F}}_{ijklmn}^{\left(1\text{D}\right)}& ={J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{3\pi {\hslash }^{10}}\left[\frac{1}{4}\frac{\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}^{4}}+\frac{{L}_{ijklmn}}{12}\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}^{3}}-\frac{{L}_{ijklmn}^{2}}{12}\right.\hfill \\ \hfill & \quad \left.{\times}\left(\frac{1}{2}\frac{\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}^{2}}+\frac{{L}_{ijklmn}}{2}\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}}+\frac{{L}_{ijklmn}^{2}}{2}\left(\text{Si}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)-\frac{\pi }{2}\right)\right)\right],\hfill \end{align} \tag{ 52 }$$

which in the asymptotic long-range regime simplifies to

$$\begin{equation}{\mathcal{F}}_{ijklmn}^{\left(1\text{D}\right)}\simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{3\pi {\hslash }^{10}}\left(\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}^{5}{L}_{ijklmn}}+\frac{5\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}^{6}{L}_{ijklmn}^{2}}\right).\end{equation} \tag{ 53 }$$

For completeness, I give the form of the range function for equidistant atoms:

$$\begin{equation}{\mathcal{F}}_{ijklmn}^{\left(1\text{D}\right)}\simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{3\pi {\hslash }^{10}}\left(\frac{\mathrm{cos}\left(6{q}_{\mathrm{F}}R\right)}{6{q}_{\mathrm{F}}^{5}R}+\frac{5\mathrm{sin}\left(6{q}_{\mathrm{F}}R\right)}{36{q}_{\mathrm{F}}^{6}{R}^{2}}\right).\end{equation} \tag{ 54 }$$

Interestingly, the decay function is similar to that of the bilinear and four-spin interactions highlighting their potential relevance at large distances.

To summarize the 1D case, the multi-spin interactions similarly decay like L⁻¹. The modification of q_F can considerably affect the amplitude of the interactions owing to their ${q}_{\mathrm{F}}^{N-1}$ dependence, N being the number of interacting magnetic moments defining the investigated plaquette. Moreover, the chiral multi-spin interactions are expected to be at short distances linear, at least, with the Rashba SOC.

3.2. Rashba model—two dimension

In the 2D case, A and B defining the Green function are given by linear combinations of Hankel functions of zero and first order, respectively:

$$\begin{align}\hfill {A}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)& =-\frac{\mathrm{i}{m}^{{\ast}}}{4{\hslash }^{2}}\left[\left(1+\frac{{k}_{\mathrm{R}}}{q}\right)\enspace {H}_{0}\left[\left(q+{k}_{\mathrm{R}}\right){r}_{ij}+\mathrm{i}{0}^{+}\right]+\left(1-\frac{{k}_{\mathrm{R}}}{q}\right)\enspace {H}_{0}\left[\left(q-{k}_{\mathrm{R}}\right){r}_{ij}+\mathrm{i}{0}^{+}\right]\right]\hfill \\ \hfill {B}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)& =+\frac{{m}^{{\ast}}}{4{\hslash }^{2}}\left[\left(1+\frac{{k}_{\mathrm{R}}}{q}\right)\enspace {H}_{1}\left[\left(q+{k}_{\mathrm{R}}\right){r}_{ij}+\mathrm{i}{0}^{+}\right]-\left(1-\frac{{k}_{\mathrm{R}}}{q}\right)\enspace {H}_{1}\left[\left(q-{k}_{\mathrm{R}}\right){r}_{ij}+\mathrm{i}{0}^{+}\right]\right].\hfill \end{align} \tag{ 55 }$$

Similarly to what was done for the 1D case, I proceed to the derivation of the multi-spin MEI in the RKKY approximation. The functional dependence with respect to the magnetic moments is not affected and all the physics is primarily encoded within the range functions, ${\mathcal{F}}_{ijkl}^{\left(\text{2D}\right)}$, i.e. the MEI without the Rashba SOC, which are listed below in the asymptotic regimes, i.e. when qr ≫ 1 and k_R ≪ q. There, the Hankel functions are ${H}_{0}\left(x\right)\simeq \sqrt{\frac{2}{\pi x}}\enspace {\mathrm{e}}^{\mathrm{i}\left(x-\frac{\pi }{4}\right)}$ and ${H}_{1}\left(x\right)\simeq \sqrt{\frac{2}{\pi x}}\enspace {\mathrm{e}}^{\mathrm{i}\left(x-\frac{3\pi }{4}\right)}$ for large x, which allow to write

$$\begin{equation}{A}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)\simeq -\frac{\mathrm{i}{m}^{{\ast}}}{{\hslash }^{2}}\frac{1}{\sqrt{2\pi q{r}_{ij}}}\enspace {\mathrm{e}}^{\mathrm{i}\left(q{r}_{ij}-\frac{\pi }{4}\right)}\enspace \mathrm{cos}{k}_{\mathrm{R}}{r}_{ij}\quad \text{and}\quad {B}_{ij}\left(\varepsilon +\mathrm{i}{0}^{+}\right)\simeq \frac{{m}^{{\ast}}}{{\hslash }^{2}}\frac{1}{\sqrt{2\pi q{r}_{ij}}}\enspace {\mathrm{e}}^{\mathrm{i}\left(q{r}_{ij}-\frac{\pi }{4}\right)}\enspace \mathrm{sin}{k}_{\mathrm{R}}{r}_{ij},\end{equation} \tag{ 56 }$$

Two-spin interactions. Here the two-spin range function

$$\begin{equation}{\mathcal{F}}_{ij}^{\left(\text{2D}\right)}\simeq {J}_{i}^{sd}{J}_{j}^{sd}\frac{{m}^{{\ast}}}{4{\pi }^{2}{\hslash }^{2}}\enspace \frac{\mathrm{sin}\left(2{q}_{\mathrm{F}}{r}_{ij}\right)}{{r}_{ij}^{2}}\end{equation} \tag{ 57 }$$

gives the usual decay function proportional to ${r}_{ij}^{2}$ for the 2D electron gas.

Four-spin interactions. In the asymptotic regime the four-spin range function is similar to the one characterizing the two-spin 1D interactions

$$\begin{align}\hfill {\mathcal{F}}_{ijkl}^{\left(\text{2D}\right)}& \simeq -{J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{{m}^{{\ast}}}^{3}}{8{\pi }^{3}{\hslash }^{6}}\frac{1}{\sqrt{{P}_{ijkl}}}\left[\text{Si}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)-\frac{\pi }{2}\right]\hfill \\ \hfill & \simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{{m}^{{\ast}}}^{3}}{8{\pi }^{3}{\hslash }^{6}}\frac{1}{\sqrt{{P}_{ijkl}}}\left[\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}{L}_{ijkl}}\left(1-\frac{1}{{q}_{\mathrm{F}}^{2}{L}_{ijkl}^{2}}\right)+\frac{\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}^{2}{L}_{ijkl}^{2}}\right],\hfill \end{align} \tag{ 58 }$$

with the already defined perimeter L_ijkl = r_ij + r_jk + r_kl + r_li and the product of interatomic distances P_ijkl = r_ijr_jkr_klr_li. The slowest decay function of the four-spin interactions is given by ${\left(\sqrt{P}L\right)}^{-1}$, which simplifies to R⁻³ when assuming a plaquette of equidistant atoms (r_ij = r_jk = r_kl = r_li = R), in which case the range function becomes

$$\begin{align}\hfill {\mathcal{F}}_{ijkl}& \simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{{m}^{{\ast}}}^{3}}{8{\pi }^{3}{\hslash }^{6}}\frac{1}{{R}^{2}}\left[\text{Si}\left(4{q}_{\mathrm{F}}R\right)-\frac{\pi }{2}\right]\hfill \\ \hfill & \simeq -{J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}\frac{{{m}^{{\ast}}}^{3}}{8{\pi }^{3}{\hslash }^{6}}\frac{1}{{R}^{2}}\left[\frac{\mathrm{cos}\left(4{q}_{\mathrm{F}}R\right)}{4{q}_{\mathrm{F}}R}\left(1-\frac{1}{16{q}_{\mathrm{F}}^{2}{R}^{2}}\right)+\frac{\mathrm{sin}\left(4{q}_{\mathrm{F}}R\right)}{16{q}_{\mathrm{F}}^{2}{R}^{2}}\right].\hfill \end{align} \tag{ 59 }$$

Thus and in contrast to the 1D-case, the four-spin interactions are expected to decay faster than the bilinear ones. A similar observation can be made for the following six-spin interactions.

Six-spin interactions. Here, the range function is approximately given by

$$\begin{align}\hfill {\mathcal{F}}_{ijklmn}^{\left(\text{2D}\right)}& \simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{24{\pi }^{4}{\hslash }^{10}}\frac{1}{\sqrt{{P}_{ijklmn}}}\left(\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)}{{q}_{\mathrm{F}}}+{L}_{ijklmn}\left[\text{Si}\left({q}_{\mathrm{F}}{L}_{ijklmn}\right)-\frac{\pi }{2}\right]\right)\hfill \\ \hfill & \simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{24{\pi }^{4}{\hslash }^{10}}\frac{1}{{L}_{ijkl}\sqrt{{P}_{ijklmn}}}\left[2\frac{\mathrm{cos}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}^{3}}-\frac{\mathrm{sin}\left({q}_{\mathrm{F}}{L}_{ijkl}\right)}{{q}_{\mathrm{F}}^{2}}\right]\hfill \end{align} \tag{ 60 }$$

and if the atoms in the plaquette are equidistant as they would be in an hexagonal lattice (r_ij = r_jk = r_kl = r_lm = r_mn = r_ni = R), the range function becomes

$$\begin{align}\hfill {\mathcal{F}}_{ijklmn}& \simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{24{\pi }^{4}{\hslash }^{10}}\frac{1}{{R}^{3}}\left(\frac{\mathrm{cos}\left(6{q}_{\mathrm{F}}R\right)}{{q}_{\mathrm{F}}}+6R\left[\text{Si}\left(6{q}_{\mathrm{F}}R\right)-\frac{\pi }{2}\right]\right)\hfill \\ \hfill & \simeq {J}_{i}^{sd}{J}_{j}^{sd}{J}_{k}^{sd}{J}_{l}^{sd}{J}_{m}^{sd}{J}_{n}^{sd}\frac{{{m}^{{\ast}}}^{5}}{144{\pi }^{4}{\hslash }^{10}}\frac{1}{{R}^{4}}\left[2\frac{\mathrm{cos}\left(6{q}_{\mathrm{F}}R\right)}{{q}_{\mathrm{F}}^{3}}-\frac{\mathrm{sin}\left(6{q}_{\mathrm{F}}R\right)}{{q}_{\mathrm{F}}^{2}}\right].\hfill \end{align} \tag{ 61 }$$

One notices that decay function pertaining to the 2D interactions is given by ${\left({P}^{\frac{1}{2}}L\right)}^{-1}$ which is different from the 1D case because of the additional factor ${P}^{\frac{1}{2}}$ that simplifies to ${R}^{\frac{N}{2}}$ for N equidistant magnetic moments. This ignites a strong difference in the long-range behavior since the power-law decay associated to the 1D electon gas is independent from the number of atomic spins, while it increases linearly in the 2D case. Similarly to the 1D interactions, the higher-order 2D magnetic interactions are expected to change when modifying q_F but the power-law dependence is different, $\approx {q}_{\mathrm{F}}^{N-2}$ for 2D instead of ${q}_{\mathrm{F}}^{N-1}$. This points to the possibility of manipulating the impact and magnitude of the magnetic interactions by controlling either the position of the magnetic moments or the Fermi energy.

To summarize, one can state that in general the decay of the multi-spin MEI can be cast into the following formula with which the behavior of the two-spin interactions is recovered ($\approx {q}_{\mathrm{F}}^{d-2}{R}^{-d}$):

$$\begin{equation}\text{Multi}-\text{spin}\;\text{MEI}\approx {\left\{{q}_{\mathrm{F}}^{\left(N-d\right)}{P}^{\frac{1}{2}\left(d-1\right)}L\right\}}^{-1}\approx {\left\{{q}_{\mathrm{F}}^{\left(N-d\right)}{R}^{\left[1+\frac{N}{2}\left(d-1\right)\right]}N\right\}}^{-1},\end{equation} \tag{ 62 }$$

with d = 1 or 2 depending on the dimension of the electron gas mediating the interactions. Noteworthy is the impact of SOC on amplitude of the chiral MEI. It is expected to be at least linear with α independently from the nature of the electron gas.

3.3. Numerical evaluation of multi-spin interactions

I consider the case of Au surfaces characterized by a surface state confined in 2D or 1D with the effective mass m* = 0.26m_e, the Fermi wave vector q_F = 0.167 Å⁻¹ and k_R = 0.0135 Å⁻¹ [68, 74]. I address the particular case of equidistant magnetic moments and therefore the MEI are discussed and analysed in terms of the interatomic distance R. It is important to keep in mind, that the contribution of the multi-spin interactions to the total energy increases in a factorial fashion when increasing the number of magnetic moments of the plaquettes. Therefore, although the coefficients defining the MEI might be smaller than the bilinear ones, their overall impact on the energy can be prominent. The parameter J^sd is chosen such to give the right order of magnitude of the first minimum of the bilinear RKKY interaction, which corresponds to an antiferromagnetic coupling of two magnetic moments. To allow a simple comparison between the 1D and 2D cases, the first largest antiferromagnetic coupling is assumed to be the same and equal to 1 meV, which is about the value measured for the RKKY interactions between Co adatoms on Pt(111) surface [87].

I start by analysing the isotropic magnetic interactions without the Rashba SOC shown in figures 1(a) and (b). Note that for the bilinear interactions, a positive (negative) value corresponds to a ferromagnetic (antiferromagnetic) coupling. As expected, the number of oscillations increases with the number of sites involved in the interactions since the wavelength is given by $2\frac{\pi }{N{q}_{\mathrm{F}}}=\frac{{\lambda }_{\mathrm{F}}}{N}$, N being the number of spins involved in the interaction while R is the interatomic distance, which for simplicity is assumed to be identical for all pairs of magnetic moments. This gives rise to situations where the higher order MEI are larger in magnitude than the bilinear ones. As expected, the 2D MEI decay more strongly than the 1D ones. This advocates for the investigation of the long-range multi-spin interactions in systems with reduced dimensionality.

Zoom In Zoom Out Reset image size

As illustrated in figures 1(c) and (d) (see also figures 2(d)–(f) for a better resolution of the six-spin interaction), the inclusion of SOC reduces the amplitude of the isotropic multi-spin interactions with the emergence of a beating effect leading to a vanishing of the interactions accompanied with a phase switch of their oscillations, which is well known for the two-spin interactions [74]. This occurs at a length scale defined by the Rashba wave vector, k_R, i.e. around 6 nm.

Zoom In Zoom Out Reset image size

In figures 2(a)–(d), the bilinear MEI are shown for the spins embedded respectively in the 1D and 2D Rashba electron gas. Similarly to the isotropic MEI, the amplitude of DMI and J^ani decreases by increasing the dimension of the electron gas. Making the moments further apart eventually favors the enhancement of the SO-driven interactions, which contrast strongly with the evanescence of the isotropic interactions around 6 nm.

The four-spin interactions are displayed in figures 2(b)–(e). These are initially, i.e. for relatively short interatomic distances, the largest in 2D but they experience a faster weakening than those of the 1D Rashba gas. Note that in contrast to the bilinear counterparts, the chiral four-spin interactions oscillate in opposite phase with respect to that of the isotropic ones.

It is instructive to compare the behavior of the chiral bilinear with that of the equivalent four-spin interactions (see figure 3). Although the amplitude of the four-spin interactions is weaker than that of the bilinear MEI, one notices that at short distances, their sign is opposite and owing to their different wavelengths, there are regions where the chirality behavior is settled by the four-spin interaction. For instance, this occurs at ≈1 nm in the 2D case, where the isotropic MEI is negligible (see figures 1(b)–(d)).

Zoom In Zoom Out Reset image size

In this article, I present a theoretical framework for the evaluation of multi-spin interactions utilizing multiple-scattering theory. This fits methodologies based on the calculation of Green functions, such as the KKR Green functions. I discuss the bilinear, four- and six-spin interactions with a particular focus on isotropic and chiral terms. Then I use this theory for the evaluation of multi-spin interactions for localized spins embedded in one- and two-dimensional electrons described by the Rashba model. Utilizing the RKKY approximation, the latter model offers the possibility of extracting analytically the long-range/asymptotic behavior of isotropic and chiral multi-spin interactions.

Within this approach, there is for each couple of sites i and j a four-spin chiral vector perpendicular to the bond connecting the two sites similarly to their DMI vector. The reported study shows that the strong contrast between the 1D and 2D bilinear MEI survives for higher-order interactions. I recover the power-law decay pertaining to the two-spin MEI, $\approx {q}_{\mathrm{F}}^{d-2}{R}^{-d}$, which I generalize to the N-spin case as ${\left\{{q}_{\mathrm{F}}^{\left(N-d\right)}{P}^{\frac{1}{2}\left(d-1\right)}L\right\}}^{-1}\approx {\left\{{q}_{\mathrm{F}}^{\left(N-d\right)}{R}^{1+\frac{N}{2}\left(d-1\right)}N\right\}}^{-1}$, where d is the dimension of the electron gas mediating the interactions, L the perimeter of the plaquette of N spins while P is the product of interatomic distances and q_F the Fermi wave vector. The two-, four- and six-spin 1D MEI experience a similar same decay, i.e. same power law, with respect to the spatial separation of the spins. In contrast, the more spins involved in the 2D MEI, the stronger is the decay. Moreover, the dependence with respect to q_F provides a path of engineering the magnitude of the higher-order MEI by tuning the electronic occupation.

Regarding relativistic effects, obviously the chiral multi-spin interactions cancel out without SOC. The smallest dependence with respect to the latter is linear at short distances. At atomic separations of the order of SOC length scale, beating effects reduce the amplitude of the isotropic MEI giving the opportunity for the chiral interactions to have a strong impact on the magnetic behavior of the investigated systems.

Numerical results were presented on the basis of Rashba parameters mimicking surface states residing on Au surfaces. This permits a visual comparison of the various interactions. Their distinct oscillatory behavior offers the possibility of exploring rich magnetic phase diagrams, which can be strongly altered when modifying the interatomic distances. This is probably possible with atomic manipulation based on scanning tunneling microscopy. Additional tuning parameters are the strength of SOC and the electronic occupation, which can be changed by gating or by modifying the nature of the substrate [102–104]. In the future, it would be interesting to extend the current study to the 3D case, either by using a 3D Rashba model [105] or by extending the Fert–Levy model [71] to the multi-spin interactions.

I acknowledge fruitful discussions with S Brinker, M dos Santos Dias, J Bouaziz, S Grytsiuk and S Blügel. Funding is provided from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC-consolidator Grant No. 681405 DYNASORE).