Prospecting chiral multi-site interactions in prototypical magnetic systems

Atomistic spin models have found enormous success in addressing the properties of magnetic materials, grounded on the identification of the relevant underlying magnetic interactions. The huge development in the field of magnetic skyrmions and other noncollinear magnetic structures is largely due to our understanding of the chiral Dzyaloshinskii-Moriya interaction. Recently, various works have proposed new types of chiral interactions, with seemingly different forms, but the big picture is still missing. Here, we present a systematic construction of a generalized spin model containing isotropic and chiral multi-site interactions. These are motivated by a microscopic model that incorporates local spin moments and the spin-orbit interaction, and their symmetry properties are established. We show that the chiral interactions arise solely from the spin-orbit interaction and that the multi-site interactions do not have to follow Moriya's rules, unlike the Dzyaloshinskii-Moriya and chiral biquadratic interactions. The chiral multi-site interactions do not vanish due to inversion symmetry, and comply with a generalized Moriya rule: If all sites connected by the interaction lie in the same mirror plane, the chiral interaction vector must be perpendicular to this plane. We then illustrate our theoretical considerations with density functional theory calculations for prototypical magnetic systems. These are triangular trimers built out of Cr, Mn, Fe and Co adatoms on the Re(0001), Pt(111) and Au(111) surfaces, for which $C_\mathrm{3v}$ symmetry applies, and Cr and Fe square tetramers on Pt(001) with $C_\mathrm{4v}$ symmetry. The multi-site interactions are substantial in magnitude and cannot be neglected when comparing the energy of different magnetic structures. Finally, we discuss the recent literature in light of our findings, and clarify several unclear or confusing points.


I. INTRODUCTION
Atomistic spin models provide the foundation to understand the properties of magnetic materials: complex magnetic ground state structures, elementary excitations (spin waves), solitons whether topologically trivial or non-trivial (domain walls and magnetic skyrmions, respectively), thermal effects and real-time dynamics 1 .In comparison to the full quantum-mechanical description, this type of model aims at a coarse-grained, low-energy description of a given material, by assuming that magnetism is well-described by assigning rigid magnetic moments (spins) to specific spatial positions (sites), and specifying how these spins interact among each other (magnetic interactions).Uncovering a new type of magnetic interaction often leads to novel magnetic states (such as quantum spin liquids 2 ), or even to new fields of research in magnetism (e.g.skyrmionics 3 ).It is then essential to have a systematic catalogue of the possible magnetic interactions, preferably coupled with a theory that can make material-specific predictions and help guide or interpret experimental efforts.
Magnetic interactions can be broadly divided into isotropic and anisotropic interactions.Isotropic interactions depend only on the relative angles between the spins, such as the original exchange interaction discovered by Heisenberg 4 , with detailed microscopic understanding provided for instance by Anderson's theory of superexchange 5 .The anisotropic interactions depend on how the spins are aligned with real-space directions (for instance, the bonds between magnetic sites or cer-tain crystal directions), and arise from relativistic effects, namely spin-orbit coupling (SOC).Interactions which are symmetric under exchange of the spin components include the single-ion anisotropy and the two-site symmetric anisotropic exchange or compass anisotropy 6,7 , which together lead to the magnetocrystalline anisotropy.An extreme case of symmetric anisotropic exchange is the Kitaev interaction 8,9 , that can stabilized exotic quantum spin liquids 2 .A different kind of two-site anisotropic interaction, which is antisymmetric upon exchange of the spin components, is the Dzyaloshinskii-Moriya interaction (DMI) 10,11 .The DMI lifts the chiral degeneracy of the magnetic structure, as it favors one sense of rotation, leading to spiral magnetic ground states, domain walls and magnetic skyrmions of well-defined chirality or handedness.All these interactions can also be classified by the number of spin components that they couple (2-spin interactions, sometimes more for the single-ion anisotropy) and by how many sites are coupled (1-site interaction for the single-ion anisotropy, the others being 2-site interactions).
The 2-spin interactions tend to favor relatively simple magnetic structures, such as ferromagnetic, antiferromagnetic, or spin spiral ground states -these are called single-Q states as their periodicity can be described by a single wavevector Q.The isotropic 4-spin interactions can couple single-Q states which are degenerate when considering only 2-spin interactions, and stabilize complex superpositions of such states.The general mechanism and the possible connection to long-range interactions can be understood resorting to Kondo lattice models [24][25][26][27] .For example, the up-up-down-down state is a double-Q state that was recently uncovered in magnetic monolayers with an hexagonal lattice [28][29][30] , with the 4-spin 3-site interaction playing the key role.The nanoskyrmion lattice in an Fe monolayer on Ir(111) 31 can be seen as another type of double-Q state stabilized by a combination of various interactions, including isotropic 4-spin ones.The prediction of a triple-Q state in this type of monolayers stabilized by the ring exchange 32 has also finally been experimentally realized 33 .Isotropic 4-spin interactions also explain why some bulk materials host short-period magnetic skyrmion lattices and other multiple-Q states [34][35][36][37][38] .Another class of materials where the isotropic 4-spin interactions play various roles are high-temperature superconductors, likely stabilizing the bicollinear antiferromagnetic ground state of FeTe 39,40 and modifying the spin wave spectrum of the parent compound La 2 CuO 4 41 .As a final example, solid 3 He is perhaps the most famous system where multi-site isotropic interactions are essential to understand its magnetic properties 42 , with up to 6-spin 6-site interactions quantitatively determined in an hexagonal monolayer 43 .
In contrast to the large body of knowledge concerning isotropic multi-site interactions, not much attention has been paid to their anisotropic counterparts.The anisotropic 2-spin interactions have been derived from an extended Hubbard model 44 , but to our knowledge no attempt has been made to reach their 4-spin counterparts, so their possible forms remain unclear.One can still proceed in various ways, which become very powerful if combined with a model-independent method of evaluating the energy of a magnetic structure, for instance resorting to density functional theory (DFT) calculations.Firstly, the energy of different magnetic structures can be systematically mapped to a spin cluster expansion [45][46][47][48][49] , of which the four-state mapping method is a simplified form able to determine all types of 2-spin interactions (see for instance Ref. 50 for an application to the Kitaev interaction).In Ref. 51, we used the spin-cluster expansion in combination with an intuitive microscopic model to catalogue all possible anisotropic 2-spin and 4-spin interactions in magnetic dimers on surfaces with strong SOC, uncovering the chiral biquadratic interaction (CBI) and also 3-and 4-site chiral interactions.Phenomenolog-ical considerations can also be used to identify allowed forms for the interactions consistent with the symmetry of a target material.This led to the discovery of chiral 4-spin 3-site interactions in MnGe 20 (for which a derivation based on multiple scattering theory was also provided), and motivated their existence in an Fe chain on Re(0001) 52 , while the magnetism of Ca 3 Ru 2 O 7 was rationalized by invoking higher-order Lifshitz invariants in connection to a Ginzburg-Landau theory 53 .Lastly, the energy can also be expanded in a Taylor series in small deviations from a reference magnetic structure.This was developed into the very successful infinitesimal rotation method for 2-spin interactions [54][55][56] , with an extension to 4-spin interactions recently proposed 20,57 , and a growing body of work addressing its application for noncollinear magnetic structures 22,23,58,59 .
In this work, we present a comprehensive study of isotropic and chiral multi-site interactions in prototypical magnetic systems.First we specify the form of our spin model, containing besides 1-site and 2-site interactions (including the DMI and the CBI) also isotropic and chiral 3-site and 4-site interactions.The employed forms of the interactions are fully justified by our microscopic model 51 (see Appendix B of that reference), and we supply simple heuristic arguments for their derivation.We then discuss the symmetry properties of the chiral 3-site and 4-site interactions, showing in detail that they are not bound by Moriya's rules, in contrast to the DMI and the CBI.After presenting our computational approach, we proceed to investigate these interactions in several magnetic systems.We chose clusters that are common atomic motifs in many periodic magnetic materials, namely homoatomic trimers on fcc(111) and hcp(0001) surfaces (C 3v symmetry), and tetramers on the fcc(001) surface (C 4v ), which illustrate the magnitude and symmetry properties of the considered interactions.Finally, we discuss several recent works in light of our findings, and present our conclusions.

II. SPIN MODEL WITH ISOTROPIC AND CHIRAL MULTI-SITE INTERACTIONS
Consider a magnetic material with well-localized spin moments on atomic sites labelled {1, . . ., N }.The energy of a given spin configuration can be described by a function E(S 1 , . . ., S N ; B) of the orientations of each spin moment, represented by classical unit vectors, |S i | = 1, and of the external magnetic field B. This energy function can contain several terms, describing different types of magnetic interactions, which we catalogue by how many spin components are involved (pspin interactions), and by how many sites are involved (q-site interactions).Time-reversal symmetry demands E(−S 1 , . . ., −S N ; −B) = E(S 1 , . . ., S N ; B).This implies that the interactions which are independent of the external magnetic field must contain an even number of spin moments (p mod 2 = 0).We now describe the in-teractions that we consider in this paper.All the presented forms have already been derived from a microscopic model in our previous work, including the multisite forms (see Appendix B of Ref. 51 for details).

A. Spin model with 3-and 4-site interactions
The magnetic interactions can be grouped in the following way, going from 1-site up to 4-site interactions: Here i, j, k, l ∈ {1, . . ., N } with N the number of magnetic atoms.The sums over sites are unrestricted except for the exclusion of repeated sites, i = j = k = l, which is indicated by the primes.We now discuss the interactions that will be computed for the prototypical magnetic systems.

B. 1-site interactions
The 1-site contribution to the magnetic energy is (α, β ∈ {x, y, z}): We have a 1-spin interaction with the external magnetic field (in energy units), and the 2-spin interaction describes the lowest-order contribution to the on-site magnetic anisotropy energy.As S α 1 S β 1 is a symmetric rank-2 tensor, only the symmetric part of K αβ 1 contributes to the sum, so it can be written as Its eigenvectors specify the easy, intermediate and hard local anisotropy axes, in increasing order of the corresponding energy eigenvalues.The non-vanishing elements of this matrix are determined by the local symmetry of the environment of the magnetic atom.Furthermore, the condition |S 1 | = 1 removes one parameter, so five independent parameters remain.The choice of free parameters used in this work is explained in Appendix A.

C. 2-site interactions
For the 2-site interactions we consider all possible 2spin interactions plus two kinds of 4-spin interactions: Here J 12 is the conventional isotropic Heisenberg exchange interaction 4 , D 12 is the vector defining the chiral Dzyaloshinskii-Moriya interaction 10,11 which can only be present in systems without inversion symmetry, and ∆J αβ 12 = ∆J βα 12 defines the symmetric exchange anisotropy 60 .Like the on-site magnetic anisotropy, the eigenvectors of the ∆J 12 matrix can be used to specify the easy, intermediate and hard axes of a given pair of magnetic atoms (see Appendix A).For very anisotropic systems this leads to the Kitaev bond-dependent symmetric exchange anisotropy 9,50 .The 2-site interactions are augmented with the isotropic biquadratic interaction 2 and with the recently-uncovered chiral biquadratic interaction

D. Heuristic arguments for the form of the chiral multi-site interactions
Instead of repeating the derivations already presented in Ref. 51, here we present heuristic arguments for the construction of chiral multi-site interactions starting from known isotropic multi-site interactions.The microscopic model discussed in that work presents a systematic expansion of the grand potential in terms of the spindependent parts of the electronic hamiltonian.These are assumed to be a local exchange coupling between the electrons and the local moments ∝ σ • S i , and the local SOC ∝ σ • L (its spatial location can be left unspecified), with σ the vector of Pauli matrices describing the electron spin, S i a unit vector describing the orientation of the spin moment at site i, and L the local orbital angular momentum operator.Isotropic interactions are represented by closed loops that contain an even number of magnetic sites and are connected by spin-independent Green functions, while chiral interactions additionally contain one spin-orbit site.The specifics of the diagrammatic expansion of the grand potential are not needed if one is only interested in the form of the interactions, and it is based this point of view that we now present our heuristic arguments.What is essential is in which order the magnetic and SOC sites appear, and the algebra of Pauli matrices.The final goal is to identify the elementary forms of the interactions which are needed to construct the atomistic spin model according to Eq. (1).
As shown in Table I, the isotropic 2-site interactions arise from the loops connecting 1 and 2 once (2-spin, Diagrams for the heuristic derivation of the form of the magnetic interactions.The numbered sites represent local exchange couplings of the type σ•Si, and the interaction is given by a closed loop connecting the sites with an even number of lines.The form of the isotropic interaction can then be obtained by a spin trace over the ordered product of sites, according to the order in which the loop is traveled.The form of the chiral interactions (first-order in SOC) can be obtained from the diagrams for the isotropic interactions by inserting σ • L in-between a pair of sites (dashed line).As a set of sites can be connected in various ways, this will lead to various forms for the interactions, which can be related to each other via symmetry operations, if applicable.bilinear) and twice (4-spin, biquadratic).The form of these interactions can be obtained by tracing a product of local exchange couplings, according to the order that the sites are travelled in the loop.For the isotropic bilinear interaction one then writes and for the biquadratic one The second term is a constant and can be disregarded.For these interactions it clearly does not matter if the loop starts from site 1 or from site 2. The corresponding chiral 2-site interactions can be generated by inserting spin-orbit coupling in-between two magnetic sites.For the bilinear interaction, inserting SOC between 1 and 2 gives (the connection containing SOC is indicated by the dashed line) This generates the DMI.Similarly, starting from the biquadratic interaction and inserting SOC between 1 and 2 produces the chiral biquadratic coupling S 2 ×S 1 (S 1 • S 2 ), with similar properties to the DMI.The chiral interaction vectors are then governed by the SOC, and it can also be shown that they comply with the symmetry operations relating 1 and 2 (see Ref. 51).In the following we will present the simplified forms of the isotropic and chiral 3-site and 4-site interactions, along with the heuristic arguments that justify them.

E. 3-site interactions
As with the 4-spin 2-site interactions, here we will also restrict our attention to isotropic and chiral interactions.Due to time-reversal symmetry, a 3-site interaction must contain an even number of spin moments, 4-spin being the minimum, which implies that at least one of them appears repeatedly.We adopt the convention of Laszloffy et al. 52 that the second site is the one that will appear repeated, see Table I.The isotropic 3-site interaction is thus expected to have the form The second term is of the form of the isotropic 2-site interaction and so it can be dropped, leaving the first term as our prototype for isotropic 4-spin 3-site interactions.
If we insert SOC between 1 and 2, we obtain The second term is of the form of the DMI and so it can be dropped, leaving the first term as our prototype for chiral 4-spin 3-site interactions.There is no relation between (1, 2, 3) and (3, 2, 1), due to the cross product in the first term.This can be understood from our heuristic argument: SOC can influence different pairs of sites, which can naturally lead to different chiral interaction vectors, depending on the symmetry of the system.If one had instead inserted SOC between 1 and 3 one would find a similar form, but with the dot product now between 1 and 2 and the cross product between 3 and 2. Inserting SOC in other places does not bring other forms for the spin coupling, up to minus signs.We then write the 3-site interactions as Note the convention that the second site is the one that appears repeated.The other possible combination of dot and cross products is covered by E 321 , which is included in Eq. ( 1).The isotropic interaction has the general symmetry B 123 = B 321 , which justifies the prefactor of 1/2 in Eq. ( 1).An alternative way of expressing the chiral interactions is using the symmetric and antisymmmetric combinations: Following our heuristic argument, the form of the isotropic 4-site interactions can obtained from a 4-site loop as shown in Table I, which translates to This is the known form of the ring exchange, including the minus sign.We take just the first term as our prototype for isotropic 4-spin 4-site interactions, as the site summations in Eq. (1) will reproduce the remaining possibilities for combining pairs of sites with dot products, and the interaction coefficients will cover the symmetry (see e.g.Ref. 21).
The form of the chiral interaction can be obtained as before, by inserting SOC between 1 and 2: This can be obtained from the form of the ring exchange by replacing a dot product by a cross product in every term, and leaving the other dot product, with an additional minus sign if the second dot product is replaced.
The perhaps unexpected complexity of this interaction can be understood from the corresponding diagram.In contrast to the chiral 3-site interaction, where SOC only affects one of two bubbles in the diagram, here SOC affects the entire loop, even if it is inserted between a specific pair of sites, and thus generates all possible kinds of pairwise chiral couplings between the four spins.As for the chiral 3-site interaction, inserting SOC between a different pair of sites leads to a different form with an independent chiral interaction vector.Once again, by exploiting the summation over sites in Eq. ( 1) it is sufficient to take the first term as our prototype for the chiral 4-spin 4-site interaction.
We thus express the contribution to the magnetic energy from the isotropic and chiral 4-site interactions as Note the convention in these interactions that the sites are paired as (1,2) and (3,4), and that the cross product applies to the first pair.If we repeat one site and write (1,2,2,3) this form reduces to the 3-site one.
We also have the general symmetries for the interaction coefficients , which justify the prefactor of 1/4 assigned to the 4-site interactions in Eq. ( 1).

III. SYMMETRIES OF MULTI-SITE INTERACTIONS
Following Neumann's principle, the magnetic energy function must respect the point group symmetry.If G is a symmetry operation of the point group, this means E(GS 1 , . . ., GS N ; GB) = E(S 1 , . . ., S N ; B).In the absence of an external magnetic field, all interactions must then be invariant under all symmetry operations of the crystallographic point group.The action of the symmetry operation can be separated as G = PO, where P maps the atomic sites to each other, and O transforms the orientations of the magnetic moments (rotation R, mirroring M, inversion I).The matrices O are orthogonal, O −1 = O T , and in particular for the mirror symmetries M = M T .For a magnetic interaction connecting q-sites (i 1 , . . ., i q ), a symmetry operation that maps these sites into themselves, P(i 1 , . . ., i q ) = (i 1 , . . ., i q ), can place constraints on the interaction coefficients (e.g.Moriya's rules).If this is not the case, i.e.P(i 1 , . . ., i q ) = (i 1 , . . ., i q ), then we only find relations between q-site interactions connecting different sets of sites.
The basic building blocks of the interactions that we discuss in this work are either dot products S i • S j or cross products S i × S j of the spin orientations, which are combined in various ways for the different types of interactions.A symmetry operation G acts on these building blocks as follows.The relation between atomic sites implied by the symmetry operation is expressed by the replacement P(i, j) = (k, l).The dot product transforms as and the last equality follows from the fact that the spatial transformations O leave the angle between vectors unchanged.If the interaction consists solely of dot products of spin orientations then the symmetry operations establish relations between the interaction coefficients, e.g.
For the cross product there is a subtlety: As the cross product is an axial vector, it transforms for proper rotations with det O = +1 and for improper rotations (inversion, mirroring) with det O = −1.For the chiral interactions, the cross product is combined with the chiral interaction vector, for instance D ij for the DMI, which then transfers the result of the spatial symmetry from the cross product to this vector, ) up to reordering we find constraints on the allowed components of the axial vector, which leads to Moriya's rules 11 .For instance, {S i , S j } → {S j , S i } if an inversion center is present between i and j, which leads to D ij = −D ij and confirms the vanishing of the DMI (and CBI) in this case.If (k, l) = (i, j), we find relations between chiral vectors connecting different sites, and this shows that the vectors must be related by a simple change in orientation, as the spatial symmetry O leaves the length of vectors invariant.The properties of the different multi-site interactions then follow from combining these principles with each type of interaction and the symmetry of the considered system.
For the benefit of the reader, we briefly recall the Moriya rules 11 for the DMI vector acting on the bond between sites i and j: 1.If there is an inversion center in the middle of the bond, D ij = 0.
2. If the bond is bisected by a mirror plane, D ij must lie in this plane.
3. If the bond is contained in a mirror plane, D ij must be perpendicular to this plane.
4. If a twofold rotation axis passes through the middle of the bond, D ij must be perpendicular to this axis.
5. If the bond lies on an n-fold rotation axis, D ij must be along this axis.
They follow from the general symmetry principles.
In the following, we consider C 3v and C 4v , which are the point groups of the magnetic trimers and tetramers for which we will present results for the magnetic interactions.The corresponding magnetic structures and symmetry operations are illustrated in Fig. 1 and the latter listed in Appendix C.

A. Symmetries for a trimer
For the trimer (C 3v symmetry, see Fig. 1a), the three mirror symmetries correspond to The isotropic interactions are only affected by a permutation of site labels, which implies J 12 = J 13 = J 23 and likewise for the isotropic biquadratic interaction B ij .Similarly, there is only one independent parameter for the isotropic 3-site interactions, B 123 .
To illustrate a symmetry operation that maps a set of sites onto itself, consider the DMI between atoms 1 and 2 and the mirror symmetry M 3 .The cross product transforms as The DMI thus transforms as which implies that M 3 D 12 = D 12 for the energy to remain invariant, or D 12 = (0, D y 12 , D z 12 ).This shows that D 12 lies in the mirror plane M 3 , as expected from Moriya's rules.The mirror symmetries M 1 and M 2 map the atom pair (1, 2) into (1, 3) and into (3, 2), respectively, so they exemplify symmetries that relate different sites.From these symmetries we get  4)) has the same transformation properties as the DMI vector.We thus have C 12 = (0, C y 12 , C z 12 ), and all the other vectors can be generated from this one with the same symmetry operations used for the DMI vectors.If the symmetry is increased from C 3v to D 3h then the mirror symmetry M z (z → −z) also applies, which would lead to −M z D 12 = D 12 and so D 12 = (0, 0, D z 12 ), again in accordance with Moriya's rules.
Similar considerations allow us to establish the symmetry properties of the chiral 3-site interaction vectors C ijk (see Eq. ( 10)).Take the interaction connecting atoms 1 and 3 through atom 2 as an example (see Fig. 1a).If SOC mediates the interaction between atoms 1 and 2 the interactions corresponds to Eq. ( 9) with the coefficient C 123 .In contrast to the isotropic multi-site interactions, as well as the DMI, there is no general relation between C 123 and C 321 , as justified by the microscopic model that assigns SOC to a specific bond.In addition, the C 3v point group symmetry does not map the form S 1 × S 2 (S 2 • S 3 ) onto itself, which would be necessary in order to find symmetry constraints for C 123 .Thus, the interaction vector C 123 for the trimer is a general 3-component vec-tor.To show how a constraint on this interaction can emerge, consider increasing hypothetically the symmetry to D 3h .The additional mirror symmetry M z leaves the site labels invariant and enforces C 123 = −M z C 123 , or C 123 = (0, 0, C z 123 ).Nonetheless, the C 3v symmetry of the trimer leads to a simplification of the chiral 3-site interaction since it can be used to relate the different interaction vectors to each other.Starting from Eq. ( 10) we can group the interactions in two subsets corresponding to a cyclic permutation of the sites, S  We now illustrate the convention of Eq. ( 1) for the 3site interactions with the case of the trimer, for which summing over all triples results in We repeat that all chiral interaction vectors can be related to the general vector C 123 using the C 3v symmetry.

B. Symmetries for a tetramer
For the tetramer (C 4v symmetry, see Fig. 1b), the four mirror symmetries correspond to the following mappings: There are two types of interactions, those that connect atoms only along the edges of the tetramer (nearestneighbors) and those that include connections across the diagonals (next-nearest-neighbors).For the 2-site interactions connecting (i, j) we can list 8 ordered pairs along the edges and 4 ordered pairs along the diagonals, with (j, i) being related to (i, j) by construction (e.g.
Once again, the CBI vector C ij has the same properties as the DMI vector.
The 3-site interactions are specified by triples (i, j, k), for which 8 connect only sites along the edges (nearestneighbors) and 16 include a diagonal connection (nextnearest-neighbors).For the isotropic 3-site interactions given in Eq. (10) we have B ijk = B kji (the repeated site is unchanged) and two independent parameters: B 123 if both dot products are along the edges, and B 124 if one of the dot products is along the diagonal.
The chiral 3-site interactions given in Eq. ( 10) can be separated into three groups.The diagrams from our heuristic arguments provide visual insight into this: The dashed line indicates the location of the cross product between spins arising from SOC.The C 4v symmetry places no constraints on a single chiral interaction vector C ijk , but establishes groups of vectors which are related to each other by symmetry.
The diagram given in Eq. ( 29 The 4-site interactions are specified by quadruples (i, j, k, l), for which 8 connect only sites along the edges (nearest-neighbors) and 16 include a diagonal connection (next-nearest-neighbors).For the isotropic 4-site interactions given in Eq. (10) we have B ijkl = B jikl = B ijlk = B jilk and two independent parameters: B 1234 if both dot products are along the edges, and B 1324 if both dot products are along the diagonal.
For the chiral 4-site interactions given in Eq. ( 10) we have C ijkl = C ijlk = −C jikl = −C jilk , which reduces the amount of independent terms by a factor of four.The six independent terms can be further grouped according to whether the cross product is along an edge of the tetramer (four terms) or along a diagonal (two terms), each generated by a reference interaction vector.In contrast to the chiral 3-site interactions, the reference chiral vectors have additional constraints enforced by the mirror symmetries.The chiral interaction vectors for the group of four terms are generated from C 1234 by rotational symmetry: 0 .This is analogous to the Moriya rule for the DMI vector: if the bond containing the cross product lies in a mirror plane, the chiral vector must be perpendicular to this mirror plane.If we again consider increasing hypothetically the symmetry, this time to D 4h , the additional mirror symmetry M z leaves the site labels invariant and enforces ).Thus, the chiral 4-site interaction vector C ijkl obeys symmetry rules similar to the well-known Moriya rules for the DMI vector, if the symmetry operation maps the pairs (i, j) and (k, l) onto themselves (up to reordering).
We now illustrate the convention of Eq. ( 1) for the 4site interactions with the case of the tetramer, for which summing over all quadruples results in IV. GLOBAL MAPPING FROM DFT CALCULATIONS TO THE SPIN MODEL The energy of a chosen magnetic structure for a given material can be obtained from first principles, for instance from a DFT calculation.Suppose that the relevant coarse-grained variables which are necessary to map to the spin model are already defined.These are the localized spin moments M i that exist on some subset of all the atomic sites of the magnetic material.The spin moments must be relatively rigid, so that their magnitude is not strongly dependent on the magnetic structure.It is then meaningful to separate the magnitude from the orientation, M i = M i S i (recall |S i | = 1), and to identify the orientations of the magnetic moments with those in the target classical spin model.One can then view the DFT total energy as a functional of the magnetic structure specified by those orientations, E DFT [S 1 , . . ., S N ].In practice, an arbitrary magnetic structure is not a stationary solution of the DFT total energy, which must then be stabilized by the addition of constraints [61][62][63] , satisfying the condition when the derivatives are evaluated for the chosen magnetic structure {S 1 , . . ., S N }.The self-consistently obtained constraining magnetic field B i is thus equal to the exact total energy derivative with respect to the orientation of the spin moment on site i, and this is the key observation that is used to map the DFT calculations to the classical spin model given in Eq. ( 1), The dependence of the DFT total energy on the spin orientations can be arbitrarily complex, and so the mapping must be chosen in a way that allows for systematic improvement and a controllable error between the energy obtained from DFT and the one computed from the parametrized spin model, for a selected set of magnetic structures.We have to generate a number of magnetic structures that is enough to obtain all the magnetic interaction coefficients from fitting the self-consistent constraining magnetic fields from the respective DFT calculations to the corresponding derivatives of the spin model.To this end, each spin moment is set to one of 14 predefined orientations (6 along the cartesian axes plus 8 along the diagonals of each octant).The magnetic structures are constructed by the tensor product of all possible orientations of each spin moment, leading to 14 N magnetic structures.For the systems we consider, we would have 2744 magnetic structures for the trimer (N = 3) and 38416 magnetic structures for the tetramer (N = 4).Enforcing time-reversal symmetry together with the corresponding point group symmetry (C 3v for the trimer and C 4v for the tetramer) achieves a reduction to 252 and 2513 inequivalent magnetic structures, respectively.As the number of configurations is still fairly large for the tetramer, we resort to randomly sampling 250 magnetic structures from this subset of inequivalent magnetic structures.The constraining magnetic fields for each magnetic structure are then self-consistently obtained from the corresponding DFT calculation.
The quality of the fit is quantified by the mean-averageerror (mae) between the constraining fields and the corresponding derivatives of the fitted spin model, summed over the used magnetic structures: Here the sum is over the used N s magnetic structures, and for each magnetic structure {S} s we compute for all the N a atoms in the trimer or tetramer the absolute error between cDFT and the fitted spin model.

V. COMPUTATIONAL DETAILS
We employ the all-electron Korringa-Kohn-Rostoker Green function method in full potential with spinorbit coupling added to the scalar relativistic approximation 64,65 .
Exchange and correlation effects are treated in the local spin density approximation (LSDA) as parametrized by Vosko, Wilk and Nusair 66 .The pure surfaces are modeled by 22 layers with two vacuum regions corresponding to four inter-layer distances each using the experimental lattice constants, which are for the considered fcc structures a Pt = 3.924 Å and a Au = 4.078 Å, and for the considered hcp structure a Re = 2.761 Å and c Re = 4.456 Å.The scattering wave functions are expanded up to an angular momentum cutoff of max = 3 and a k-mesh of 150 × 150 is used.The nanostructures are embedded in real space using a nearest-neighbor cluster.
The constrained DFT We considered compact fcc-top-stacked trimers for the Pt(111) and the Au(111) surfaces, and hcp-top-stacked trimers for the Re(0001) surface.On the Pt(001) surface two compact tetramers are used.We define the average relaxation of a nanostructure as r = 1−d/d0, where d is the average vertical distance between the atoms comprising the nanostructure and the atoms of the surface layer, and d0 is the bulk vertical interlayer distance.These are computed with Quantum Espresso, and inform the value used in the KKR calculations as shown.We also show the spin magnetic moments per atom M obtained from both types of calculations.
calculations for the magnetic structures required for the mapping to the spin model are performed as explained in the previous section, following the procedure introduced for magnetic dimers in our previous work, Ref. 51.
To account for structural relaxations we use the Quantum Espresso package 67,68 .The surfaces are modelled by five layers surrounded by a vacuum region of the same thickness.The nanostructures are placed on 4 × 4 supercells of the surfaces and a k-mesh of 2 × 2 × 1 is used.Exchange and correlations effects are treated in the generalized gradient approximation using the PBEsol functional 69 , and we used ultrasoft pseudopotentials from the pslibrary 70 with an energy cutoff of 100 Ry.The nanostructure as well as the first surface layer are allowed to relax.For the setup of the KKR geometry we use only the vertical relaxation of the nanostructure, since the relaxation of the first surface layer turns out to be negligible in agreement with previous studies 71 .The results of the structural relaxations are summarized in Table II.

VI. RESULTS
We consider Cr, Mn, Fe and Co trimers on the Pt(111), Re(0001) and Au(111) surfaces, and additionally Cr and Fe tetramers on the Pt(001) surface.The constrained magnetic configurations described in Section IV were used to fit the interaction parameters for the 1-site, 2-(p-spin, q-site) interactions: (  site, 3-site and 4-site interactions.The quality of the fits is quantified by the mean-average-error, Eq. ( 36).We show in Table III how the fit improves (or not) by adding more types of interactions to the spin model.A significant drop in the mean-average-error when adding a new class of interactions to the fit shows that these make an important contribution to the energy.However, no significant improvement in the fitting error can also arise due to the weakness of the additionally-fitted interactions.Overall, all of the systems could be well-fitted with our procedure.

A. Trimers with C3v symmetry
A trimer of identical magnetic atoms arranged as a compact triangle on an fcc(111) or hcp(0001) surface has C 3v symmetry.The chosen coordinate system, labelling of the magnetic atoms and mirror planes is shown in Fig. 1a.These are called top-stacked trimers, as they enclose a surface atom 60 .The isotropic and chiral interactions for the trimer are specified by Eqs. ( 4) and (10), with the detailed form for 3-site interactions given in Eq. ( 24).The reference interaction parameters are given in Table IV, from which the full spin model can be parametrized by applying the C 3v -symmetry operations.The symmetric anisotropy parameters (see Eqs. (2) and Eq. ( 4) and Appendix A) are not of our primary interest, but can be found in Table VI.
We first discuss possible scenarios considering the isotropic interactions alone.The corresponding magnetic structures are sketched in Fig. 3.The dominant interaction is almost always the isotropic 2-spin interaction J 12 .J 12 < 0 favors a ferromagnetic alignment, while J 12 > 0 favors the planar Néel state for the trimer (the two forms N + and N − shown in Fig. 3  .We give the reference parameters for all isotropic and chiral 2-site and 3-site interactions, Eqs. ( 4) and (10).The full set of interactions can be obtained from the shown ones by applying the C3v-symmetry operations to the trimer, see Section III A.  IV shows that for most trimers B 12 < 0 and B 123 > 0, so their combined action prefers an up-updown state.Its energy difference to the ferromagnetic state is E uud − E F = −4 (J 12 + B 123 ).We find this energy difference to be about 7 meV for the Mn trimer on Re(0001), but for all other trimers with J 12 < 0 the ferromagnetic state is much more stable than the other states.If J 12 > 0 then the relevant energy difference is to the planar Néel state, E uud − E N = 1 4 (2J 12 + 9B 12 − 7B 123 ), showing that the isotropic 4-spin interactions strongly penalize the planar Néel state.For instance, this energy difference is just 10 meV for Cr 3 on Pt(111), and becomes negative (−2 meV) when this trimer is placed on Re(0001), which would make the up-up-down state the ground state if only isotropic interactions were at play.
In order to easily compare the isotropic and chiral con-tributions to the magnetic energy of the trimer, we now consider a family of magnetic structures with 3-fold rotational symmetry.We parametrize the spins as with Considering only isotropic and chiral interactions, for these structures the total magnetic energy per trimer atom is (c.f.Eq. ( 24)) The isotropic interactions combine into an effective interaction J 12 (θ), while the chiral interactions form the combined chiral interaction vector D 12 (θ).Here C 123 = (0, C y 123 , C z 123 ), as for these magnetic structures C x 123 does not contribute.The chiral part of the energy can be split into an out-of-plane contribution, and an in-plane contribution This shows that the z-component distinguishes between the two planar Néel states shown in Fig. 3, favoring one or the other depending on its sign.The y-component only results in an energy gain for the structures with an anticlockwise rotation of the spins, which is maximized for φ = 0 • if D y 12 (θ) sin(2θ) > 0 and for φ = 180 • otherwise.A structure that would favor a clockwise rotation of the spins around the triangle might still be able to gain energy from the in-plane component of the effective chiral vector, if the opening angles are not fixed to be the same for all pairs of spins.
For these types of magnetic structures, the isotropic biquadratic and 3-site interactions contribute to the energy as B 12 + B 123 .As seen from Table IV, these interactions have opposite signs for most trimers, cancelling out almost completely for Cr 3 on Pt(111) -see Fig. 4(a), Mn 3 on Re(0001) and Fe 3 on Au(111).When J 12 is the dominant interaction, we can approximate α(θ) ≈ 1 − 3 2 θ 2 for J 12 < 0 (canted ferromagnetic), from which we obtain that the isotropic energy increases with the tilt angle with the coefficient J 12 + 2(B 12 + B 123 ), while for J 12 > 0 we have instead α(90 • + θ) ≈ − 1 2 + 3 2 θ 2 (canted planar Néel), with the corresponding coefficient being J 12 − (B 12 + B 123 ).Thus, the change in the energy from tilting the spins is influenced by the isotropic 4-spin interactions in different ways depending on the chosen reference configuration, and can be an indirect way of establishing their relative importance.Overall, we see that it is important to include the different types of isotropic interactions, and not just the 2-site biquadratic interactions, when extending the spin model.
The chiral interactions lead to a canting of the magnetic structure with a well-defined chirality.This is controlled by the effective chiral vector, which has contributions from the DMI, the CBI and the chiral 3-site interactions.As with the isotropic interactions, the latter two interactions can either cooperate and compete, and as they are vector interactions this can happen in different ways for the different vector components.The CBI and the chiral 3-site interaction for Cr 3 on Pt(111) have comparable magnitudes but are almost perpendicular to each other, combining into a vector C 12 +C 123 which is almost half of the DMI in magnitude.However, the additional dependence on α(θ) of the CBI and the chiral 3-site interaction conspires to cancel the enhancement of the DMI near θ = 90 • , instead leading to its weakening, as shown in Fig. 4(b,c).The different contributions to the energy are shown in Fig. 4(d).The strong antiferromagnetic J 12 favors the planar Néel state, which becomes canted due to the usual DMI -the two minima near θ = 90 • correspond to two canted N + states, with the spins tilting either towards the center of the trimer or away from it.The chiral 4-spin interactions weaken the canting of the Néel state.Finally, the contribution to the energy from the symmetric anisotropic interactions (not included in Eq. ( 38)) modifies the canting angle, due to the tilted easy-axis of the single-site magnetic anisotropy.
While the previously discussed rotational symmetric magnetic structures showed the importance of the Moriya-like components of the chiral 3-site interaction, given by J12 alone.The analytic form of the energy due to the isotropic and chiral interactions is given in Eq. (38).For each θ we consider the value of φ that minimizes the energy.The last curve includes the contribution from the symmetric anisotropic interactions.The energy minima are indicated by the grey stars.where all the mutual angles are the same, which implies (S 2 − S 1 )•S 3 = 0, as we discussed in relation to Eq. (38).However, if one started from an up-up-down state (see Fig. 3) then (S 2 − S 1 ) • S 3 = −2, the non-Moriya interaction becomes active and would lead to a canting of this magnetic structure.One way of comparing the different chiral interactions acting on the trimer is to find what is the magnetic structure that is favored by each type of interaction on its own.The out-of-plane DMI (D z 12 ) favors the planar Néel state, selecting either N + or N − depending on its sign (see Fig. 3).The in-plane DMI (D y 12 ) favors a canted form of the N + state, with all spins having the same polar angle of θ = 45 • and with the spins tilting either towards the center of the trimer or away from it.The chiral biquadratic and the Moriya components of the chiral 3-site interactions favor similar states as the usual DMI.The out-of-plane components (C z 12 and C z 123 ) also favor either N + or N − , while the in-plane components (C y 12 and C y 123 ) favors a canted form of the N + state, but with a different value of the polar angle, θ = 25.5 • .In contrast, the non-Moriya contribution favors a completely different state, shown in Fig. 5.One of the spins stays normal to the plane of the trimer, while the other two cant symmetrically towards each other (φ 1 = 173.0• and φ 2 = 7.0 • ) with a large opening angle with respect to the out-of-plane spin (θ 1 = θ 2 = 78.6 • ).

B. Tetramers with C4v symmetry
A tetramer of identical magnetic atoms arranged as a compact square on an fcc(001) surface has C 4v symmetry.The chosen coordinate system, labelling of the magnetic atoms and mirror planes are shown in Fig. 1b, and the symmetry matrices are described in Appendix C. The isotropic and chiral interactions for the tetramer are specified by Eqs. ( 4), ( 10) and ( 16), with the detailed form for 4-site interactions given in Eq. (32).The reference interaction parameters are given in Table V, from which the full spin model can be parametrized by applying the C 4vsymmetry operations.The symmetric anisotropy parameters (see Eqs. (2) and Eq. ( 4) and Appendix A) are not of our primary interest, but can be found in Table VII.
We first discuss possible scenarios considering the isotropic interactions alone.The dominant interaction is the nearest-neighbor interaction J 12 .J 12 < 0 favors a ferromagnetic alignment, while J 12 > 0 favors the collinear up-down-up-down antiferromagnetic state (going around the tetramer).As the tetramer allows for next-nearest-neighbor interactions, J 13 could in principle compete with J 12 and stabilize an up-up-downdown state, but this requires J 13 ≥ |J 12 |, a condition far from being satisfied in our tetramers.The isotropic biquadratic and 4-site interactions contribute to the energy in exactly the same way for these three collinear states, E = 4 (B 12 + B 1234 ) + 2 (B 13 + B 1324 ), and so do not favor any of them.The isotropic 3-site interactions do distinguish between the three collinear states: B 123 distinguishes the up-up-down-down state from the other two, while B 124 distinguishes all three of them.However, the magnitude of these isotropic 4-spin interactions for our tetramers is not strong enough to modify the collinear state favored by the nearest-neighbor interaction J 12 .
In order to easily compare the isotropic and chiral contributions to the magnetic energy of the tetramer, we now consider a family of magnetic structures with fourfold rotational symmetry.We parametrize the spins as in Eq. ( 37), but now setting As before, θ = 0 • or 180 • corresponds to the ferromagnetic state.However, θ = 90 • is not the up-downup-down state, but a planar noncollinear state with all spins perpendicular to their nearest-neighbors.We define α , which together characterize the openings of nearest-neighbor spins.We also define α and an in-plane contribution As for the trimer, only one chirality can gain energy from the in-plane components of the effective DMI.We now focus on the Fe tetramer, for which J 12 = −25 meV and J 24 = −2.4meV favor a ferromagnetic structure, and so an opening due to chiral interactions is expected.For a small opening, we can set α(θ) ≈ 1 − θ 2 and α (θ) ≈ 1 − 2θ 2 , so that the second derivative of the isotropic energy has the coefficients When all interactions are taken together, the ground state of the Fe tetramer on Pt(001) is found to be a ferromagnetic structure almost perpendicular to the plane of the surface, with a symmetric canting of all spins (s = +) away from the center of the tetramer, with θ = 8 • .The main origin of this canting are the chiral interactions across the diagonals of the tetramer.
The Cr tetramer has strong antiferromagnetic interactions favoring the collinear antiferromagnetic up-downup-down state.By setting θ 1 = θ 3 = θ and θ 2 = θ 4 = 180 • − θ, one can still use Eq. ( 43) but with α → −α and taking care of the change in handedness of some cross products, which leads to The out-of-plane contribution from the chiral interactions is unchanged.The strongest interactions are J 12 = 29 meV and J 24 = 12 meV.For a small opening, we can set α(θ) ≈ −1 + θ 2 and α (θ) ≈ 1 − 2θ 2 , so that the second derivative of the isotropic energy has the coefficients  001) is found to be a slightly canted up-down-up-down state almost collinear with either the x-or y-directions, due to the symmetric anisotropic interactions, and the spins tilt away from the xy-plane by ∆θ = ±2.3• .This magnetic structure is illustrated in Fig. 6.

VII. RELATION TO OTHER WORKS
We now relate our findings to other works addressing magnetic interactions in very disparate systems, highlighting common ground and clarifying several aspects concerning the multi-site interactions.
First we would like to mention that Bornemann et al. 72 presented an extensive survey of the magnetic properties of diverse clusters of Fe, Co and Ni on Ir(111), Pt(111) and Au(111), using the infinitesimal rotation method based on the ferromagnetic state as reference.However, direct comparison of our data in Table IV with their reported values is difficult due to their neglect of structural relaxations and various other computational differences, so we will not attempt this here.We have previously calculated the magnetic exchange interactions for Fe trimers on Pt(111), in connection to scanning tunneling microscopy experiments 60, 73 .In Ref. 73 we employed the infinitesimal rotation method for the 2-spin interactions with a correction due to the spin polarizability of Pt, obtaining J 12 = −54 meV, D y 12 = 1.3 meV, and D z 12 = −0.7 meV.As discussed in Sec.VI A, we should compare J 12 from the infinitesimal rotation method with J 12 + 2(B 12 + B 123 ) = −44 meV, and the DMI to D y 12 = 4.9 meV and D z 12 = −0.4meV, as defined in Eq. (38).We believe that the discrepancies can be explained by the different treatment of the potential: Ref. 73 employed the atomic sphere approximation while our present work makes no shape approximation to the potential.
To the best of our knowledge, there is only one previous work addressing isotropic 4-spin interactions in magnetic clusters, Ref. 47, using a different computational geometry and approximations, which once again cautions against quantitative comparisons.Ref. 47  Fe chains on Re(0001) were experimentally found to have a short-period spin-spiral ground state 74,75 .While exploring the magnetic properties of this system with DFT calculations, Lászlóffy et al. 52 found an inconsistency between the noncollinear spin structure obtained via magnetic force theorem calculations and the one obtained with an atomistic spin model containing only 2spin interactions.While trying to understand the origin of this inconsistency, they introduced on phenomenological grounds chiral multi-site interactions.Our systematic procedure for constructing chiral multi-site interactions fully justifies the phenomenological forms proposed by these authors.We note, however, that we did not recover the aforementioned inconsistency when revisiting the same system without any shape approximation to the potential, instead finding a spin structure which is consistent with the experimental one 75 .
Grytsiuk et al. 20 investigated theoretically the complex magnetism of MnGe 34,36 .Given that the Mn atoms in MnGe are coordinated in triangular plaquettes, the authors proposed so-called topological-chiral interactions which are built upon the scalar spin chirality of the three spins forming one such triangle, χ 123 = S 1 • (S 2 × S 3 ).The chiral-chiral interaction is a 6-spin 3-site interaction κ CCI 123 (χ 123 ) 2 that does not require SOC, while the spin chiral interaction is a 4-spin 3-site interaction built of terms such as C SCI 123 • S 1 χ 123 and is driven by SOC.DFT calculations showed that both types of interactions are quite strong in MnGe.Given our proposed classification of multi-site interactions into isotropic (no SOC required) and chiral (SOC required) interactions, the question arises as to how these topological-chiral interactions fit this classification.In fact, we show in Appendix B that the chiral-chiral interaction can be expressed solely using dot products of the involved spins, so it is of the generic form of the isotropic interactions, and that the spin-chiral interaction can be rewritten using combinations of dot products and cross products, thus being covered by the chiral 4-spin interactions that we discuss in the present work.However, interactions built solely out of the scalar spin chirality cannot reproduce all the different types of interactions that we uncovered in the present work (see Appendix B), and for a complete spin model our systematic forms for the interactions should be used.
Cardias et al. 59 recently submitted a preprint titled 'Dzyaloshinskii-Moriya interaction in absence of spinorbit coupling'.Our microscopic model and heuristic arguments concerning the forms of the magnetic interactions unequivocally show that without SOC the interactions are isotropic, no cross products of spins appear, and without cross products there are no DMI-like chiral interactions.We are convinced that these puzzling findings can be explained using the generalized atomistic spin model that we discussed in this work, by identifying the type of isotropic multi-spin and/or multi-site interactions that lead to the obtained angular dependence of the energy or its first derivative.
As we mentioned in the Introduction, several works have proposed ways of extending the infinitesimal rotation method to 4-spin interactions 20,22,23,57 .In particular, Ref. 57 discusses both isotropic and chiral multi-site interactions, and derives the corresponding expressions for their calculation in a DFT context.They showed that the chiral 4-spin 3-site interactions do not vanish for centrosymmetric systems with the example of bcc Fe, in full agreement with our symmetry analysis of these interactions.The authors also introduce a chiral 3-spin 3-site interaction which is defined through the scalar spin chirality χ 123 .This kind of interaction is not time-reversalinvariant unless the interaction coefficient also changes sign, and so we believe that the corresponding calculations need to be reinterpreted in a way that complies with time-reversal symmetry.

VIII. CONCLUSIONS
In this work, we presented a comprehensive framework for isotropic and chiral multi-site interactions, along with systematic calculations of these interactions for several magnetic trimers and tetramers.
First, we imposed the general requirement of timereversal invariance of the magnetic energy on the possible form of the interactions, and gave simple heuristic arguments that can be used to obtain the form of the multi-site interactions.We thus arrived at a generalized atomistic spin model containing 2-spin and 4-spin interactions that couple up to four distinct magnetic sites.Next we demanded that our interactions comply with the crystallographic point group symmetry, with the concrete examples of C 3v (trimers) and C 4v (tetramers).Contrary to the 1-site and 2-site interactions, those based on three or four sites are much less constrained by the point group symmetry operations.For instance, while the 2site Dzyaloshinskii-Moriya interaction vanishes if there is an inversion center in the middle of those two sites, this is not the case for the chiral 3-and 4-site interactions.We also found that the respective chiral interaction vectors can have components which are forbidden for the DMI due to Moriyas rules.The chiral multi-site interactions do comply with a generalized Moriya rule: If all sites connected by the interaction lie in the same mirror plane, the chiral interaction vector must be perpendicular to this plane.
After outlining our global mapping scheme from DFT calculations to a target spin model, we presented our results on a series of homoatomic trimers and tetramers on several surfaces with strong spin-orbit coupling.While in most-cases the dominant interaction is the familiar isotropic Heisenberg exchange, this is not so for the trimers on the Re(0001) surface.For the trimers, we found that the isotropic biquadratic and 3-site interactions tend to counteract each other, while the chiral biquadratic and 3-site interactions more easily combine due to their vector nature, supporting or hindering the DMI depending on the magnetic structure and on the type of atoms forming the trimer.We also discussed the magnetic structure favored by the non-Moriya component of the chiral 3-site interactions.For the tetramers on Pt(001), the isotropic 4-spin interactions were found to cooperate and have a strong contribution for the Cr case, while the chiral 4-spin interactions dominate over the DMI, playing a leading role in the canting of the ground state magnetic structures of the tetramers.Lastly, we briefly addressed recent proposals for multi-site interactions and placed them into the context of our work.
We believe that our work is a timely contribution to the growing research activity on materials with complex magnetic structures, such as multiple-Q states [28][29][30][31][32][33] , for which the role of the chiral multi-site interactions is yet to be explored.Our detailed exposition of the symmetry properties of the multi-site interactions, as well as concrete examples for their enumeration, should clarify the bookkeeping which is essential to properly account for all possible ways of combining a set of sites with a given type of magnetic interaction.Our example systems, trimers and tetramers, are ubiquitous building blocks (triangles and squares, respectively) of many lattices, which will help transfer our findings to extended systems, from a single layer to bulk magnets.This defines an ellipsoid with principal axes of length given by the eigenvalues and orientation given by the eigenvectors.As an example, consider the mirror symmetry M = 1−2 n⊗n, where n is the unit normal to the mirror plane.If the anisotropic interaction is invariant under this mirror symmetry, MAM = A, then n must be one of the eigenvectors of A.
Consider now the symmetric anisotropic interactions.For α,β K αβ i S α i S β i we have an on-site anisotropy matrix, for which only the energy differences when the spin is aligned with each of the principal axes is meaningful.We can thus set one of its eigenvalues to zero, and we do this with the one invariant under the mirror plane, say λ 3 = 0 if n = u 3 , which leaves a finite two-dimensional subspace K i = λ 1 u 1 ⊗ u 1 + λ 2 u 2 ⊗ u 2 .This corresponds to the matrix K i having three independent parameters, which only reduce to two if the coordinate axes are aligned with the eigenvectors of this matrix.If i = j we have a symmetric anisotropic exchange matrix, and we cannot a priori set any eigenvalue to zero.We can however rewrite α,β J αβ ij S α i S β j = J ij S i • S j + α,β ∆J αβ ij S α i S β j .Now the role of the anisotropy is isolated in the matrix ∆J ij , and we can choose to set to zero the eigenvalue that is invariant under the mirror symmetry, if the symmetry applies.
First we discuss the trimer with C 3v symmetry, see Fig. 1a.We take atom 3 as reference, for which M 3 K 3 M 3 = K 3 can be used to impose the form  VII.
It is the sum of a constant, an isotropic 6-spin 3-site interaction, and three isotropic biquadratic interactions.The isotropic 6-spin interaction written as sums over triples of dot products of spins was derived from a Hubbard model in Ref. 19

(see last row of Table II in that work).
We now address the spin-chiral interaction 20 (SCI) and its reduction to combinations of dot products and cross products.In three dimensions, any vector v can be written using three linearly-independent vectors {a, b, c} as The general chiral 4-site form given in Eq. ( 16) has four symmetries while the 4-site SCI has the six symmetries of the scalar spin chirality built from 234, so a mapping is not possible.The 3-site restriction of this formula (the one actually discussed in Ref.The first term contributes to the 2-spin DMI between sites 2 and 3, and the remaining two terms fall into the form of the 4-spin 3-site chiral interaction given in Eq. (10).In fact, this corresponds to an antisymmetrized form of the chiral 3-site interaction, C − 312 • S 3 × S 1 (S 1 • S 2 ) − S 2 × S 1 (S 1 • S 3 ) , so it cannot capture the symmetrized remainder of the chiral 3-site interaction.
Appendix C: Symmetry operations for C3v and C4v For C 3v symmetry we give two rotation and three mirror matrices.The 120 • rotations around the z-axis are R + being anticlockwise and R − clockwise.The matrices for the mirror planes in Fig. 1a are and M 3 = diag(−1, 1, 1) is a diagonal matrix.
For C 4v symmetry we define three rotation matrices and four mirror planes.The reference rotation matrix is the anticlockwise rotation by 90 • around the z-axis, from which all rotation matrices can be defined by a suitable power or transpose.We have two mirror planes along the cartesian axes, which follows the same Moriya rules as the 2-spin DMI.These were found to be the dominant 4-spin contributions to the 2site interactions in our previous study, Ref. 51.If the relativistic spin-orbit interaction can be neglected then only the isotropic Heisenberg and biquadratic interactions remain.The interaction coefficients have the general symmetries J 12 = J 21 , D 12 = −D 21 , ∆J αβ 12 = ∆J αβ 21 , B 12 = B 21 and C 12 = −C 21 , which justify the prefactor of 1/2 assigned to the 2-site interactions in Eq. (1).

FIG. 1 .
FIG. 1. Illustration of the symmetries of the considered nanostructures.The magnetic atoms (numbered) are illustrated by red spheres and the surface atoms by grey spheres.Mirror symmetries are indicated by dashed lines while the arrows indicate rotational symmetries.a) Compact trimer on a hexagonal surface with C3v symmetry.b) Compact tetramer on a square lattice with C4v symmetry.

or M 1 D
12 = D 31 and M 2 D 12 = D 23 .The rotations lead to D 23 = R + D 12 and D 31 = R − D 12 .As we demonstrated in Ref. 51, the CBI vector (see Eq. (

FIG. 2 .
FIG.2.Relations between the chiral 3-site interaction vectors under C3v symmetry.For comparison, the DMI vectors are along the mirror planes intersecting each edge of the triangle.

1 =
{C 123 , C 231 , C 312 } and S 2 = {C 321 , C 132 , C 213 }.The rotational symmetries can be used to relate the vectors in the set S 1 to each other, C 231 = R + C 123 and C 312 = R − C 123 .The vectors in the set S 2 transform among themselves in the same way, and the mirror symmetries connect the two sets.For instance, the mirror M 3 imposes which leads to C 213 = −M 3 C 123 , and similarly C 321 = −M 1 C 231 and C 132 = −M 2 C 312 .The pattern formed by the six chiral vectors is shown in Fig. 2.

1 √ 2 ,
neighbors) and J 13 = J 24 (next-nearestneighbors), and similarly for the isotropic biquadratic interaction B ij .The DMI vectors along the edges of the tetramer have the same form as for the trimer, as they contain a mirror plane perpendicular to the edge connecting each pair, M x D 12 = D 12 so D 12 = (0, D y 12 , D z 12 ).The pattern of DMI vectors around the edges of the tetramer can then be simply obtained by rotation starting from D 12 as reference: D 23 = RD 12 , D 34 = RD 23 and D 41 = RD 34 .The DMI vector across the diagonals have to comply with two mirror planes, M + and M − .We have D 24 = −M − D 24 and D 24 = +M + D 24 so D 24 = D 24 1 √ 2 , 0 , and the other DMI vector is obtained from ) has connections with both cross and dot products along the edges (nearestneighbors).It can be drawn in 8 symmetry-related ways, and the corresponding interaction vectors are all related to C 123 .The rotations give directly C 234 = RC 123 , C 341 = RC 234 , and C 412 = RC 341 .The remaining four chiral vectors can be related to C 123 using a mirror, C 214 = −M x C 123 , and rotations, C 321 = RC 214 , C 432 = RC 321 and C 143 = RC 432 .The connections including a diagonal have to be distinguished by whether the cross product occurs on an edge or on a diagonal.The first case is represented by the diagram in Eq. (30), which can be drawn in 8 symmetry-related ways, with all chiral vectors being related to C 124 .The rotations give C 231 = RC 124 , C 342 = RC 231 , and C 413 = RC 342 .Applying a mirror we find C 213 = −M x C 124 , and with rotations we get C 324 = RC 213 , C 431 = RC 324 and C 142 = RC 431 .The second case is represented by the diagram in Eq. (31), which can be drawn in 8 symmetry-related ways, with all chiral vectors being related to C 421 .The rotations give C 132 = RC 421 , C 243 = RC 132 , and C 314 = RC 243 .Applying a mirror we find C 423 = −M − C 421 , and with rotations we get C 134 = RC 423 , C 241 = RC 134 and C 312 = RC 241 .This completes the list of the chiral 3-site vectors.
2341 and C 4123 = RC 3412 .The reference chiral vector is constrained by the mirror plane M x combined with the general symmetry of the interaction, which yields −M x C 1234 = C 2143 = −C 1234 , and so C 1234 = (0, C y 1234 , C z 1234 ).This is analogous to the Moriya rule for the DMI vector: if the bond containing the cross product is bisected by a mirror plane, the chiral vector must lie in this mirror plane.The chiral interaction vectors for the group of two terms are specified by C 1324 and C 2431 = RC 1324 .Now the M + symmetry imposes a stronger constraint on the reference vector, −M + C 1324 = C 1342 = C 1324 , from which follows

FIG. 3 .
FIG. 3.Basic magnetic structures for a trimer.F: Ferromagnetic state.uud: up-up-dow state.N+: planar Néel state with anticlockwise rotation of the spins.N−: planar Néel state with clockwise rotation of the spins.

αB 12 JJ
FIG. 4. Contributions to the energy of C3v-symmetric magnetic structures of Cr3 on Pt(111).(a) Effective isotropic interaction J12(θ) = J12 + α(θ)(B12 + B123).(b) y-component and (c) z-component of the effective chiral interaction vector D12(θ) = D12 + α(θ)(C12 + C 123 ).(d) Total magnetic energy per trimer atom relative to the energy of the planar Néel state (θ = 90 • )given by J12 alone.The analytic form of the energy due to the isotropic and chiral interactions is given in Eq.(38).For each θ we consider the value of φ that minimizes the energy.The last curve includes the contribution from the symmetric anisotropic interactions.The energy minima are indicated by the grey stars.
reports for the isotropic 4-spin interactions of Cr 3 on Au(111) B 12 = −4.4meV vs. our B 12 = −5.1 meV, and B 123 = 7.1 meV vs. our B 123 = 8.1 meV.Both values are in fair agreement with ours, also concerning the opposite sign of these interactions.They also report J 12 = 145 meV vs. our J 12 = 88 meV, and |D 12 | = 1.8 meV vs. our |D 12 | = 4.3 meV, both quite different from ours, so the numerical agreement at the level of the 4-spin interactions might be fortuitous.

TABLE III .
Mean-average-error (mae), Eq.(36), in units of [meV] of the spin model fits to the constrained DFT calculations for the different trimers and tetramers that we studied.

TABLE IV .
have the same energy).The isotropic biquadratic interactions B 12 and the 3-site ones B 123 tend to have comparable Magnetic interactions in different compact top-stacked trimers in units of[meV]