First-principles Dzyaloshinskii–Moriya interaction in a non-collinear framework

We have derived an expression of the Dzyaloshinskii–Moriya interaction (DMI), where all the three components of the DMI vector can be calculated independently, for a general, non-collinear magnetic configuration. The formalism is implemented in a real space—linear muffin-tin orbital—atomic sphere approximation (RS-LMTO-ASA) method. We have chosen the Cr triangular trimer on Au(111) and Mn triangular trimers on Ag(111) and Au(111) surfaces as numerical examples. The results show that the DMI (module and direction) is drastically different between collinear and non-collinear states. Based on the relation between the spin and charge currents flowing in the system and their coupling to the non-collinear magnetic configuration of the triangular trimer, we demonstrate that the DMI interaction can be significant, even in the absence of spin-orbit coupling. This is shown to emanate from the non-collinear magnetic structure, that can induce significant spin and charge currents even with spin-orbit coupling is ignored.

In this expression e i represents the direction of the magnetic moment at site i, and the sum is done over pairs of spins. The interatomic exchange interaction between spins is J ij , and theoretical calculations of it and the DMI, D ij , offers the possibility to predict important properties, e.g. magnon excitations 16 , which are important for fields such as spintronics and its technological applications. However, it was shown that these parameters are strongly affected by non-collinearity in case of non-Heisenberg systems 25 and, therefore, the properties need to be calculated under a non-collinear formalism. In fact, it has been shown that the calculation of magnon excitations in bcc Fe are closer to the experimental results if one takes into account the calculation of J ij in the non-collinear (NC) framework 26 . Furthermore, in more complex systems, DMI can play an important role in the magnon spectra calculation 23,27 . Also, as the systems get more and more complex, the magnetic structure needs to be taken into consideration for a more precise description of the static and dynamic properties. For that, specifically for www.nature.com/scientificreports/ the DMI, a NC description is needed. In this work, we show that the DMI is quite different, both regarding the strength and direction, for different magnetic configurations. For instance, for triangular trimers considered in this work, where there is a local minima at the Néel configuration (in-plane magnetic moments with an angle of 120 • between each other), the DMI is different in direction and strength compared with the DMI calculated considering a collinear configuration. Similar to what the improved J ij could provide in Ref. 26 , we believe that a general expression for the DMI can lead to new and better understanding of experimental nano-magnetism.
In the absence of SOC interaction, the interatomic exchange coupling can be mapped onto a scalar Heisenberg spin model derived in terms of infinitesimal spin rotations based on Green's function formalism that is known as the Liechtenstein-Katsnelson-Antropov-Gubanov (LKAG) formalism 28 . The method has been extended for NC atomistic spin configurations in Ref. 26 . Here, we present a calculation method for the DMI vectors when the SOC is present as a perturbation resulting in (only) anti-symmetric off-diagonal elements of the exchange tensor. This method has two advantages. Firstly, it may give an efficient and expedient estimation for the DMI vectors from a computational point of view. Secondly, this method is able to handle NC atomistic spin configuration by having explicit formulas for the three components for the DMI vectors in the NC framework 26 . We describe here the implementation of this method and illustrate its applicability on selected magnetic nano-sized objects. In order to reach the established goals, we used the linear muffin-tin orbital-atomic sphere approximation method (RS-LMTO-ASA) to calculate self-consistently the electronic structure for each studied structure.
We have considered Cr triangular trimer supported on the Au(111) surface (Cr 3 /Au(111)) and Mn triangular trimers supported on the Ag(111) (Mn 3 /Ag(111)) and Au(111) (Mn 3 /Au(111)) surfaces. As described in the "Methods" section, we have employed a linear muffin-tin orbital method in the atomic sphere approximation for the calculation of the electronic structure. Furthermore we considered a real-space version of this method, which enables calculations without periodic boundary conditions, such as clusters of atoms on a substrate. The fcc(111) substrate has been modeled by a slab of 4500 atoms with the experimental lattice parameter of Au (or Au). To simulate the vacuum, outside the surface, we added two layers of empty spheres above the Ag (or Au) surface, in order to provide a basis for the wave function in the vacuum and to treat charge transfers correctly. The calculations of the Mn (or Cr) adsorbed clusters have been performed by embedding the clusters as a perturbation on the previously self-consistently converged Au(111) or Ag(111) surfaces. The Mn (Cr) sites and the first-nearest-neighbor (NN) atoms (Ag or Au, and empty spheres) around the defect were recalculated self-consistently, while the electronic structure for atoms far from the Mn (Cr) cluster were kept unchanged. The collinear and the NC solutions presented here were obtained from fully relativistic calculations, where the spin-orbit interaction was treated at each variational step 29 . We have performed calculations without structural relaxation, where the Mn (Cr) atoms occupy the unrelaxed hollow positions, assuming the experimental lattice parameter of the Ag (or Au) substrate. The DMI parameters, D ij 's, have been computed using an implementation following the theory described in Ref. 26 and further details are shown in the next section. The obtained values of D ij are then used to analyse the effects of non-collinearity onto its strength and direction, aiming to show their differences and consequently their importance for further studies that use them.

Spin interactions
Let us consider a system that is described by the spin Hamiltonian given by Eq. (1). If a spin at site i is rotated by infinitesimally small angles and the rest of the spins are unperturbed then an energy variation term, δE i can be calculated. In the case when two spins in the spin system at site i and j are simultaneously rotated by infinitesimally small angles and the rest of the spins are unperturbed then the total energy variation is given as where the last term characterizes the interaction between the two spins while the other first two terms describe the interaction with the rest of spin system. We will focus on the interacting term δE ij which is given as The reason we do it is that the same scheme can be applied for the electronic grand potential variation, δ� , under the same perturbation. It will have also three terms when two spins are simultaneously rotated by infinitesimally small angles at the same time. The analog interacting term is denoted by δ� ij . For this terms we have derived an expression in terms of multiple scattering electronic structure theory, see Eq. (19) in the "Methods" section. Having Eq. (19), it is easy to identify −2J ij δe i δe j and −2D ij δe i × δe j in Eq. (3) with −2J * ij δe i δe j and −2D * ij δe i × δe j , respectively. In other words, the Heisenberg exchange parameter, J ij can be mapped onto an expression that contains information from the electronic structure, J * ij , defined by Eq. (20) in the "Methods" section. Expression J * ij is reduced to (Eq. (22)) in the collinear limit, which is commonly referred to as the LKAG formula. On the top of this, D ij of the spin-Hamiltonian can be mapped onto information of the electronic structure energy, represented by D * ij in Eq. (21) of the "Methods" section. This expression is derived for the general, NC, spinalignment. The corresponding expressions for the DMI in the collinear limit were obtained in Refs. 30,31 . We show in the "Methods" section that the presence of SOC interaction mainly contribute to −2D * ij δe i × δe j , see Eq. (37). We note that the derivation of how to calculate the J ij and D ij parameters from the electronic structure is based on the analysis of two-site energy variation, i.e, the mapping procedure is here made by the comparison of δE ij to δ� ij . It should be also noted that the exchange and DMI parameters in Ref. 32 have been derived from one-spin rotation when the application of the NC, generalized, LKAG sum rule makes possible to make a mapping between electronic structure properties and an effective spin-Hamiltonian. This approach was also shown in Ref. 33 . The detailed presentation of the connection of the two approaches will be the part of a further study.  Fig. 1) and, in this particular coordinate system, this DMI vector has two independent components, an in-plane and an out-of-plane, D 12 = D �n + D zẑ , with n = (x +ŷ)/ √ 2. Note that DMI for the Mn trimer on Ag(111) is comparable with the one of the Mn trimer on Au(111), despite the fact that Au has a larger spin-orbit coupling than Ag. One would therefore expect that an interaction www.nature.com/scientificreports/ that commonly is ascribed to pure spin-orbit effects, would be larger for the latter substrate. However, the Au 5d bands hybridize only weakly with the 3d states of Mn 34 . This weakens the influence of the large spin-orbit coupling of the Au atom. In contrast, the hybridization between Mn and Ag is stronger. Although Ag has a weaker spin-orbit strength compared to Au, this increase in hybridization results in a DM interaction that is similar for the two systems. If one considers the Heisenberg Hamiltonian for the triangular trimer, for systems where the J ij favours antiferromagnetism, which is the case for all of the systems studied here, the magnetic configuration that minimizes the equation is the known Néel configuration, where the magnetic moments are in-plane with an angle of 120 • between each other. In that case, the D z component is responsible to lift the degeneracy between different signs of the vector chirality χ = S 1 × S 2 + S 2 × S 3 + S 3 × S 1 in the Néel configuration 14 . It has been discussed in Ref. 26 that the Heisenberg picture can be broken and non-Heisenberg quantities can play an important role in case of non-collinearity, therefore changing the values of J ij and D ij . In the next section, we show that the calculated DMI's highly depend on their magnetic configuration reference by calculating the DMI's for their respective NC, minima.
Non-collinear DMI. At this stage of the calculation, the magnetic moments were allowed to relax in orientation, to find the lowest energy configuration. The co-planar 120 • Néel state was found to be the local minimum energy configuration. For this co-planar (CP) configuration, the Heisenberg exchange and DMI interaction were calculated, according to the discussion of the previous section and the "Methods" section. The results were significantly different, both in direction and strength, compared to the ones found for the collinear case. In the CP case, all three studied systems have the DMI direction almost exclusively perpendicular to the surface, with the in-plane components D x and D y being very weak compared with D z . The strengths of these DM interactions are quite intriguing, since they are approximately one order of magnitude larger than the strength of the collinear case. A similar behaviour was discussed in Ref. 32 for the compound Mn 3 Sn. It is important to mention that in Ref. 32 , the J ij 's and DMI terms were derived considering linear, single site variations of the magnetic moment, which is sufficient in the case of non-collinear arrangement. In this present work, the non-collinear arrangement is also central, but we have followed an approach of making two site spin rotations for the calculations of interaction parameters of the spin-Hamiltonian.
The comparison between collinear and CP case can be seen in Table 2 and the DMI directions are shown in Fig. 2. The strength of the DMI is seen to be significantly enhanced for the co-planar orientation. It is unlikely that this enhancement is due to spin-orbit coupling, that is expected to influence the electronic structure in a similar way, for any magnetic configuration. Our explanation of this enhancement, presented below, follows instead an analysis of spin currents, that are known to originate from either spin-orbit coupling or a co-planarly aligned magnetism 35 . Influence of spin and charge current on the DMI. It is known that non-collinearity can lead to spin and charge currents in a system 36 . Moreover, the relation between the DMI and spin-current was discussed in Refs. 37 and 38 . The spin-current that flows between a pair of atoms that the DMI couples, is the main reason for the difference between the collinear and NC DMI strengths. Also, since both the collinear and the CP Néel state have no scalar spin chirality, i.e. C 123 = S 1 · (S 2 × S 3 ) = 0 , the charge-current flowing in the system is weak and is exclusively due to the spin-orbit coupling 36 . These mechanisms have been recently discussed in Ref. 39 .
In order to analyze in detail the role of the spin and charge currents for the DMI, we first divide the Green's function (GF) into spin independent G 0 ij and spin dependent G µ ij parts. Here µ = x, y, z and G ij being 9 × 9 matrices in case of spd-orbitals, with orbital indexes written here as α, β when needed explicitly as G ijαβ . In this representation, one can describe the behaviour of the GF under a given symmetry operation T, that is an operator that transposes simultaneously the orbital and site indices. As discussed previously 32 , one can further decompose the GF as and (4) Table 2. Comparison of the strength of the DMI between the collinear and CP configuration, in meV. Note that (a), (b) and (c) denotes the systems shown in Fig. 1 (4) and (5), denotes whether that part of the GF is is even (0) or odd (1) under T. This property is valid for a real basis, in case of spherical harmonics. Further details of this derivation can be seen in Ref. 32 . The definition of the components of the GFs in Eqs. (4) and (5)  With these relationships, one can rewrite the DMI formula, Eq. (35) from the "Methods" section, as www.nature.com/scientificreports/ While G 00 and G µ0 are the result of the magnetic system, G 01 and G µ1 are connected with spontaneous currents that may result from the NC magnetic texture of the system. We will exemplify these mechanism later in the paper. With this formulation, one may divide the DMI in two parts: one related with the spin-current, called D S ∝ G µ1 ; and other one related with the charge-current D C ∝ G 01 . In this interpretation of the DMI, its origin is due to inter-atomic spin and charge currents. The origin and the mechanism behind the DMI has been intensively discussed in the literature 34,[36][37][38]40 . The DMI is known to be a direct effect from the broken inversion symmetry, caused by the SOC. In the case analyzed here, a triangle trimer, if the magnetic configuration is collinear, e.g. with the magnetic moments aligned parallel to the surface normal, the SOC induces a current that has direct influence on the DMI. If one turns off the SOC (see Computational details Section), of this collinear structure, both spin and charge currents vanish, due to the lack of both vector and scalar spin chiralities. However, there are other mechanisms which induce spin and charge currents into a system. The non-collinearity is the relevant mechanism of the presently analyzed systems, from which an additional current flows. Therefore, the currents that may flow in the systems analyzed here, has two different sources which will affect differently the DMI. The electronic structure method used here allows to scale the spin-orbit strength and it is therefore possible to analyse the relative importance of relativistic effects, and the influence of NC magnetic moments. Our results show that in the absence of SOC, the current driven by the non-collinearity is large, with the result that a significant DMI appears, even in the absence of SOC. This has recently also been discussed for Mn 3 Sn 32 . A different perspective of these ideas, based on multi-spin and/ or multi-scale interaction has recently been discussed by Brinker et al in Ref. 41 . Efforts have been made trying to map a many-body problem onto spin-lattice models attempting to take in consideration the multi-site spininteraction and their origins [42][43][44] . The model used in this paper only considers terms up to second order, meaning that high-order terms are folded into the terms here presented. We suggest then that non-collinear exchange interactions are incorporated into our DMI-like term.
In this work, we have performed calculations of the DMI for three different trimers with different magnetic configurations, as represented in Fig. 2d, namely: (i) a rotation of the magnetic moments from a ferromagnetic configuration to a CP Néel magnetic structure and back to a ferromagnetic configuration. This is achieved by rotating the out of plane component of each atomic moment from θ = 0 • to 180 • (data shown from θ = 60 • to 120 • ), while keeping the in-plane component of each magnetic moment at a 120 • with respect to the neighbouring Mn (or Cr) atom. We also performed calculations (ii) of the DMI while scaling of the SOC strength, for a configuration when the magnetic moments have an in-plane angle of 120 • between each other and an angle of 60 • with respect to the z-axis ( θ = 60 • ). Furthermore, we investigated the effects of a global spin rotation, R α , of the magnetic configuration around the y axis with an angle α (iii). Once the rotation is done, the new DMI was calculated and the vectors were rotated back to the original reference frame.
In Fig. 3, case (i), we show results from a calculation of the DMI for the triangular trimers of Cr on Au(111), Mn on Au(111) and Mn on Ag(111) at every 10 • , starting from an out-of plane angle of the moment ( θ ) of 60 • and ending at 120 • . A self consistent calculation for the electronic structure was done for every step, and the DMI was evaluated. The calculations were done including SOC (full line) and without SOC (dashed line), and it was found that in general the difference between the two calculations is minor. Our results also show that if the magnetic configuration is collinear, the DMI only exists in the presence of the SOC (data not shown in Fig. 3). However, in the NC magnetic configuration case, the DMI is clearly seen from the figure to be non-zero even in the absence of the SOC. This is due to the spin-and charge currents flowing in the system, that are driven by the non-collinearity. Note that Fig. 3 describes the strength of the three components of the Dzyaloshinskii-Moriya interaction. It is noteworthy that when θ = 90 • , i.e. the Néel magnetic structure, the scalar spin chirality C 123 is zero and the charge current flowing in the system is zero. This means that the charge-current dependent part of the DMI, D C , vanishes and all the contribution comes from the spin-current dependent D S , whose z components are allowed by symmetry. Finally, in the CP Néel state and in the lack of SOC, the x and y components of the spin-current dependent part of DMI are also zero.
In Fig. 4, case (ii), we show results from a converged calculation, with the magnetic moments having an angle of θ = 60 • with respect to the z-axis. In this plot we show results of the x-, y-and z-components of the DMI as function of a scaling parameter that was used to tune the strength of the SOC. In this plot, a value of 1 on the abscissa means 100% of the true, calculated SOC strength, and e.g. a value 0.5 corresponds to 50% of the true, calculated SOC strength. One may note from the figure that for the particular configurations considered, the influence of the SOC is fairly minimal. This demonstrates that the DMI can have a huge, non-relativistic source, that emanates from currents driven by the non-collinearity of the magnetism.
In Fig. 5, case (iii), the set-up of Fig. 4 is repeated with θ = 60 • , but now a global spin rotation, R α , of the magnetic moment around the y-axis is done. Here, α is the varied angle that is used to rotate moments around the y-axis (see Fig. 2d). The non-relativistic part of the DMI should be a constant and therefore, for zero SOC one expects D ij (α) · R αêi × R αêj = D ij (0) ·ê i ×ê j , for any value of α . In the limit of weak SOC, the quantity R −1 α D ij (α) ≈ D ij (0) should be fairly independent of α with any anisotropy directly connected to the SOC. It is possible from Fig. 5 to verify an almost constant value for the DMI with very small dependence on α . This corroborates that the DMI is present, in this NC scenario, mostly due to non-relativistic effects, a conclusion also reached in the analysis of Fig. 4.  www.nature.com/scientificreports/ reach this conclusion, the DMI was evaluated from a NC magnetic configuration, for all trimers. The results show a drastic difference, both for the direction and strength of the DMI, between a collinear and a NC configuration. In addition, we argue here that the DMI carries a dependence on the spin and charge currents flowing in the system, that can be induced by either spin-orbit effects or by a NC magnetic configuration. Thus, the fact that the non-collinearity induces a spin and charge current highly contributes to the final value and direction of the DMI. In the particular case of the triangular trimer and the NC magnetic configuration, we have shown that the DMI is mostly influenced by non-relativistic effects. In fact, the results presented here show that the non-relativistic contribution can be orders of magnitude larger than the contribution from spin-orbit coupling. Lastly, it is noteworthy that this paper describes a simple way of calculating all the three components of the DMI from a single calculation that works for any considered magnetic configuration.

Methods
The fundamental equation of a scalar relativistic multiple scattering theory formalism is given as  www.nature.com/scientificreports/ where τ ij stands for the scattering path operator (SPO), p i , denotes the inverse of single site scattering operator (ISO). The double underline stands for the fact that they are, in case of spd orbitals (single underline), 18 × 18 matrices in the spin and orbital (double underline). In Eq. (11) L = (l, m) stands for the angular momentum and magnetic quantum numbers, σ refers to the spin-index, G 0 ij is the free (or bare) electron Green's function and indices i and j refer to the considered lattice sites. G 0 ij is calculated from the Hamiltonian of the free particle, hence it is spin-independent. We introduce a general notation for the single site scattering operator in a NC framework as where the unit vector e i refers to the magnetic spin moment at site i (as it was already defined in the Introduction), σ stands for the vector formed by Pauli-matrices, I 2 is the unit matrix in spin space, t 0 i denotes the nonmagnetic (charge) part, and t i stands for the magnetic (spin) part of the single site scattering operator, ε is the energy variable.
For the ISO, one can introduce a similar notation as for the t i in Eqn. (12) as follows, Later we will need to deal with the variation of the ISO under a small rotation that can be written as where δe i stands for the deviation of a spin moment after an infinitesimal rotation at site i. Finally, the SPO has a structure as where T 0 ij denotes the charge while T ij = T x ij , T y ij , T z ij stands for the spin part of the SPO. We have defined the quantities required to calculate the pairwise total energy (grand potential-) variation in a NC framework, i.e., when a spin at site i and an other spin at site j are rotated by an infinitesimal angle simultaneously. As it has been shown in Ref. 26 , the expression for changes of the grand potential when spins at site i and j are rotated with infinitesimal angles can be written where both the ISO and the SPO depend on the energy, ε . This expression is commonly evaluated using the Lloyd formula 45 . For this purpose, we introduce the matrix and the matrix Inserting Eqs. (14) and (15) into Eqn. (16), and using Eq. (17), one can formally obtain that where expression and have been introduced. The first term in Eq. (19) is denoted δ� H ij (for Heisenberg like contributions) while the third term is denoted δ� DM ij (for DMI like contributions). Here, the new superscript * was introduced in order to differentiate the terms derived from the multiple scattering theory from the ones presented in Eq. (1). However, a practical calculation of exchange parameters used in an effective spin-Hamiltonian, amount to identifying the connection between the parameters in Eq. (1) and those in Eqs. (20) and (21). www.nature.com/scientificreports/ We analyse in this paragraph the case when the atomic moments are in a collinear arrangement (and the spin orbit-interaction is absent). One can then write that T ij = 0, 0, T z ij , and the component of the SPO for the up and down spin channels can be defined as T It can be shown that D * µ ij = 0 in this case and the second term in Eq. (19) gives only higher (fourth) order contribution to in terms of the angle variation. Only the first term is left in which J * ij is simplified to The expression defined in Eq. (22) is widely used for ab-initio calculations of Heisenberg exchange, and is commonly referred to as the LKAG formula.
In the rest of the section, we analyse for a general NC magnetic configuration what contributions are added to the variation of , when the SOC is considered as a perturbation, i.e. only the terms that being of first order in the SOC parameter, ξ , are kept. One can then write that where the perturbed quantities are denoted by the symbol ′ . Our task now is to determine the δp ′ -s and τ ′ -s. First, we determined the perturbed single scattering matrix with the help of perturbed Green function. Next, we introduce the perturbed G ′ ij where G ′ ij = G 0 ij + G ij . The matrix G 0 ij has the structure of G 0 ij I 2 , while G ij can be written as ξŴ ij σ , where its vector component Ŵ µ ij ( µ = x, y and z) is obtained as and where L µ is a component of the angular momentum operator. This implies that Ŵ ij transforms under T as Keeping the leading terms, we get that t ′ i = t i + ξ t i Ŵ ii σ t i which can alternatively be written as see Eqs. (12)(13). This implies that i.e., we have to calculate the simpler expression that contains the perturbed SPO, however, the ISO is present in a non-perturbed form.
Next we need to determine the perturbed SPO. It can be shown that where It can furthermore be demonstrated that �τ ij has the same structure as τ ij , i.e., Using Eqs. (24) and (30), it can be shown by using Eq. (25) that where α runs over the 0, x, y and z. This implies that  (33) and (34), it can be shown that A ′αα ij = A αα ij . In addition, by using the same equations, it can be also written that A ′µν ij + A ′νµ ij = 2A µν ij . This implies that Eq. (36) will be reduced to the expression From this expression we note that primed quantities, that signify that they contain linear contributions of the SOC, are only found for the DMI interaction. Hence, we conclude from perturbation theory that spin-orbit interaction influences only the DMI interaction.
Summarizing the "Methods" section, we have shown that the exchange parameter J ij defined in the spin model, Eq. (1), can mapped onto J * ij given by Eq. (20) for a general, NC, spin arrangement. Note that the mapping is correct if the second term in Eq. (19) is not relevant. This is the case for instance with symmetry resolved exchange parameter in the T 2g channel in bcc Fe bulk 25 . We can also note that the J ij parameter in Eq. (1) for collinear magnets is the LKAG parameter given by Eq. (22). In addition, we have shown that the DMI vector D ij can be calculated from the electronic structure by the D * ij formula given by Eq. (21). This formula is general, and holds for any kind of NC spin configuration. We have also shown, see Eq. (37) that the SOC interaction contributes only to the DMI interaction in leading order.
Computational details. The calculations were performed using the ab initio RS-LMTO-ASA method 46-52 , which is suitable to describe the physics of isolated clusters supported on surfaces in an efficient way, since it is real-space based and does not depend on translational symmetry. Moreover, the RS-LMTO-ASA method has been generalized to describe NC magnetism 46,49,53 , and is based on the Haydock recursion method 54 . Our Hamiltonians are constructed within an RS-LMTO-ASA formalism 55 , and therefore, all calculations presented here are fully self-consistent, and the spin densities were treated within the local spin-density approximation (LSDA) 56 . The continued fraction, that occurs in recursion method, have been terminated with the Beer-Pettifor terminator 57 after 21 recursion levels. The electronic structure calculations were also performed using the ASA and the full potential approximation (latter not shown) and no distinguished differences were seen that would change the conclusions of the paper. The spin-orbit coupling interaction is treated by adding a H SOC = ξ L · S to the Kohn-Sham Hamiltonian, with ξ calculated self-consistently. When calculations without SOC are mentioned in the text, it means that we take ξ → 0 after the self-consistent calculation is done.