Navigation on the energy surface of the noncollinear Alexander-Anderson model Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

NAVIGATION ON THE ENERGY SURFACE OF THE NONCOLLINEAR ALEXANDER-ANDERSON MODEL

P. F. Bessarab1'2, A. Skorodumov3, V. M. Uzdin2,3, H. Jonsson4,5

1 Royal Institute of Technology KTH, Stockholm, Sweden 2 St. Petersburg State University, St. Petersburg, Russia

3 ITMO University, St. Petersburg, Russia

4 University of Iceland, Reykjavik, Iceland

5 Aalto University, Espoo, Finland

bessarab@kth.se, skorodumov@vingrad.ru, v_uzdin@mail.ru, hj@hi.is PACS 05.20.Dd, 75.10.-b

Implementation of the multiple impurity, noncollinear Alexander-Anderson model is described in detail and an analytical expression given for the force which determines the orientation of the magnetic momenta as well as a corresponding magnetic force theorem. Applications to trimers of Cr, Mn and Fe adsorbed on a metal surface are described, including the energy surface as a function of the the angles specifying the orientation of the magnetic momenta and minimum energy paths for transitions between stable states, which necessarily involve noncollinear ordering. A simple model for the interaction of a magnetic STM tip with a Cr dimer on a surface is briefly described. A finite range approximation is also formulated, which simplifies the self-consistency calculations and results in linear scaling of the computational effort with the number of magnetic atoms in the system. The theoretical approach described here can be used to study magnetic systems with complex energy landscapes, including stable states and magnetic transitions in frustrated magnetic systems, over a range in length scale, from a few to several thousands of magnetic atoms.

Keywords: itinerant magnetism, Alexander-Anderson model, magnetic force theorem, noncollinear ordering, energy surface, minimum energy path.

Received: 7 November 2014 Revised: 15 November 2014

1. Introduction

Noncollinear magnetic states are of great importance, not only because the stable states of many systems are noncollinear, as has been demonstrated in recent experimental measurements of nano-scale islands [1] and thin layers [2,3], but also because magnetic systems driven out of equilibrium by some external perturbation such as a magnetic field, spin-polarized current or thermal fluctuation, as well as spin dynamics at finite temperature necessarily involve configurations where the magnetic momenta can be far from parallel to each other. In particular, a proper description of non-stationary, noncollinear magnetic states is important in studies of thermally activated magnetic transitions. Within harmonic transition state theory (TST) [4], minimum energy paths (MEPs) in a configuration space connecting stable states have to be found. Following an MEP means rotating the magnetic momenta in an optimal way so as to minimize the energy with respect to all degrees of freedom perpendicular to the path. Of particular importance for estimating the rate are the maxima along the path which correspond to first order saddle points on the energy surface. These define the activation energy for transitions between

stable magnetic states and give the temperature dependence of the transition rate. Being a path of maximal statistical weight, an MEP also gives a detailed description of the optimal transition mechanism. Even if the stable states of the system are collinear, the configurations along the MEP can be noncollinear, representing complex, non-uniform rotations of the spins [5-7].

The basic physical quantity from which most magnetic properties of an electronic system can be derived is the electron density, which becomes a 2 x 2 matrix in a noncollinear, spin-polarized theory:

where + and — denote the spin projections. The sum of the diagonal elements of p(r) gives the charge density while their difference gives the projection of the spin density on the quantization axis. The off-diagonal elements of the matrix (1) determine in quasi-classical framework the components of the magnetization density perpendicular to the quantization axis. In general, at any point in space, the expectation value of magnetization density can be calculated from

Here, a = (ax, ay, az) are the Pauli matrices and Trs means the trace operation in spin space. A continuous vector field of magnetization density can be calculated using, for example, density functional theory (DFT) [8-10], for many important magnetic materials including 3d-transition metals and rare-earth magnets. But, in most cases, the magnetization density is highly localized on atomic sites. An atomic moment approximation can then be used where the orientation of the magnetization within a region surrounding each atom is considered to be fixed and the magnetization density is assumed to be zero in the interstitial regions between the atomic sites. A magnetic configuration is then defined by a set of magnetic momentum vectors, one vector associated with each of the magnetic atoms.

A further approximation can be made based on the hierarchy of relaxation times. The relaxation of charge and magnitude of the magnetic moments is assumed to be much faster than the relaxation of the orientation of the magnetic vectors [11]. The fast degrees of freedom, charge and magnitude of magnetic moments, are assumed to adjust instantaneously to changes in the slow degrees of freedom, the direction of the magnetic moments as defined by polar and azimuthal angles 9 and 0. The magnetic properties of a system, including the total energy, are then completely characterized by the slow degrees of freedom only. This is analogous to the Born-Oppenheimer approximation in atomic systems where the fast degrees of freedom are associated with the electrons while the positions of the nuclei are the slowly varying degrees of freedom.

Within the atomic moment and adiabatic approximations, spin dynamics simulations as well as iterative searches for stable magnetic states and MEPs between them can be viewed as navigation on the systems energy surface as a function of the angles defining the orientation of the magnetic moments. Efficient navigation requires efficient evaluation of not only the energy but also the changes in energy with orientation, that is the 'forces', or 'torques' acting on the magnetic moments. The energy surface can in principle be characterized within the tight-binding approach [12] and DFT [13,14]. Such calculations are, however, complicated and computationally intensive. In order to describe nonstationary magnetic states, local constraining fields are needed to orient the magnetic vectors in the predefined orientations [15]. These local fields are not known a priori, so an additional iterative cycle needs to be added to the self-consistency procedure to find them. As a result, the calculation of an arbitrary, nonstationary, noncollinear state in a system of several non-equivalent magnetic moments is a challenging task within DFT and tight-binding methods.

(1)

m(r) = Trs [ap(r)].

(2)

In most cases, the evolution of extended magnetic systems is described using simple, phenomenological models, in particular Heisenberg-type models, where the total energy and gradients of the energy can be obtained analytically. However, in order to describe magnetic systems accurately enough, the model Hamiltonian may need to include several phenomenolog-ical terms. In addition to the usual magnetic exchange, anisotropy and dipole-dipole interaction, more elaborate interactions, such as biquadratic exchange and Dzyaloshinsky-Moriya interaction have been invoked to reproduce experimental observations accurately enough [2]. The magnitude of the magnetic moments and interaction parameters in such Heisenberg-type models are typically kept unchanged as the magnetic vectors rotate. This approach can be accurate enough for small deviations from collinear stable states, but is expected to fail for large rotation angles in itinerant electron systems [16] where the magnitude of the magnetic moments and coupling parameters depend on the relative orientation of the moments. While additional parameters and elaborate expressions for the dependence of the parameters on the orientation of the magnetic moments can, in principle, be used to make a Heisenberg-type model reproduce any magnetic system, the transferability of the parameter values obtained in this way may be quite limited. Models that better describe the underlying physics and rely on only a few, well defined parameters are then preferable.

The Alexander-Anderson (AA) model [17,18] describes the interaction of magnetic impurities in a system containing itinerant electrons. It includes two electronic bands: a quasilocal-ized band of d-electrons and a band of itinerant s(p)-electrons. A noncollinear extension of the AA model (NCAA) has been developed within mean-field approximation as well as an efficient implementation of the self-consistency calculations for an arbitrary number of non-equivalent magnetic impurities using the recursive Green function method [19,20] and analytical transformations of the density of states [21-23]. The NCAA model can be applied to large and complex magnetic systems where a self-consistent calculation of the number of d-electrons and magnitude of the magnetic moments is carried out for a given orientation. For an arbitrary, stationary or non-stationary orientation, only the number of d-electrons and magnitude of the magnetic moments are modified during the self-consistency calculations. The orientation of the magnetic vectors remains unaffected, i.e. in the NCAA model, changes in orientation are completely decoupled from the self-consistency procedure. This is different from DFT calculations, where the orientation of the magnetic moments at a non-stationary point is modified during self-consistency calculations unless local constraining fields fixing the orientation of the magnetic moments to the predefined direction are included [15].

In our recent paper [23], a magnetic force theorem for the NCAA model has been derived which makes it possible to calculate the energy gradient without repeated self-consistency calculations and also provides analytic expression for the force acting on the orientation of the magnetic moments. This theorem is particularly important for large scale simulation of dynamics, calcualtion of minimum energy paths, or, in general, navigation on the energy surface of a magnetic system.

The NCAA model, however, does not describe the band structure of an itinerant magnet in detail. Similar to the Heisenberg-type models, the NCAA model only focuses on the energetics due to magnetic interactions. But, unlike Heisenberg-type models, the NCAA model correctly accounts for the itinerant nature of 3d-transition metal systems: both the magnitude of the magnetic moments and the interatomic exchange parameters (analogous to exchange parameters in a Heisenberg Hamiltonian), are predicted to vary upon rotation of the magnetic momentum vectors. The NCAA model has, for example, been used successfully to describe

magnetism of 3d-metal surfaces and interfaces [22]. Moreover, a noncollinear ordering of magnetic moments in nanoclusters of 3d-metal atoms was obtained in calculations using the NCAA model [24-26], and this prediction was later confirmed by DFT calculations [27,28].

This article is organized as follows. In the following section, the NCAA model is described as well as the method used in the self-consistency calculations. In Sec. III, the magnetic force theorem is presented and an expression given for the magnetic force. In Sec. IV, the method is applied to magnetic trimers and to a simple model of a magnetic tip of a scanning tunneling microscope (STM) interacting with a dimer adsorbed on a metal surface. In Sec. V, the finite range approximation is discussed as well as an application to Fe-atoms adsorbed on a W(110) surface, as both a full monolayer and a large island. Section VI gives a summary.

2. Noncollinear Alexander-Anderson model

The AA model [18] extended to multiple-impurities and noncollinear magnetic ordering has been described elsewhere [21,24]. However, we present here a comprehensive description of the model and its implementation for completeness.

The electronic structure of a 3d transition metal system is approximated by an itinerant s(p)-electron band and five degenerate, quasi-localized d-orbitals. The Hamiltonian is as follows:

H = ^Gk^ka + ^£0iUia + ^ (vik4aCka + Vkiclad

k,a i,a k,i,a

t 1 (3)

+ 5] dldja + -Y, Uiniani-a,

i=j,a i,a

where only one of the five d-orbitals is considered explicitly. Here, dja(dia) and cka(cka) are creation (annihilation) operators for d-electrons localized on atom i and itinerant s(p)-electrons with the wave vector k, respectively; nia = d|adia, nka = ckacka are the corresponding occupation number operators. Greek indices denote spin projection (a,0 = ±). The energy of non-interacting s(p) electrons, ek, and d-electrons, e0, hybridization parameters, uik, hopping parameters, uij, and Coulomb repulsion between electrons with opposite spin projection, Ui, are spin independent. The last term in the Hamiltonian, Uiniani-a, describes the interaction between d-electrons localized on atom i.

This Hamiltonian is invariant with respect to the choice of quantization axis. In order to describe noncollinear magnetic states, a mean field approximation is invoked at each site i for the d-electron operators, d\a and dia, where the quantization axis, zi, is chosen to be along the local magnetic moment associated with atom i. The last term in Eqn.(3) is transformed according to

niani-a ~ nia(ni-a) + (nia)ni-a (nia)(ni-a). (4)

Operators with a tilde correspond to a local reference frame associated with the quantization axis zi and angular parentheses denote expectation values of operators. The last term in Eqn. (4) represents the double counting term. The mean-field Hamiltonian is rewritten in terms of dja and dia whose quantization axis is the laboratory z axis, the same for all sites i. Creation (annihilation) operators are transformed using the spin-2 rotation matrix [29]:

= ufdr) , (5)

d i- di-

where U is given by:

U = / exp (i^i/2) cos (0i/2) exp (-i&/2) sin (0i/2) \ (6)

exp (i&/2) sin (0i/2) exp (-i&/2) cos (0i/2) J ' ()

Here, the polar angle 0i and the azimuthal angle ^ define the direction of ith magnetic moment with respect to the laboratory quantization axis z. After performing the operations (4) and (5) for each site i, the mean-field Hamiltonian is given by:

Hmf = £ SkUka + £ £«nia + £ (^d^Ck« + U^C^d« )

k,« i,a k,i,a

+ £ Ui?di« j - - £ Ui (n2 - Mi2)

(7)

ij

i'j,a,,

where

£« = + U (Ni - a cos 0iMi), (8)

uj = U (¿«^ - 1) ¿j exp (-ai^i) sin 0iM + (1 - ¿j) ¿«^ Uj. (9)

Ui 2

The number of d-electrons, N^ and the magnitude of the magnetic moments, M^ in Eqs. (7)-(9) are defined in terms of the expectation value of occupation number operators:

Ni = <dl+d~i+> + (di-di-) = (ni+) + (Ui-), (10)

Mi = <dl+ di+> - <di-di-> = <Ui+> - (Ui-). (11)

In this model, the magnetic structure of a system consisting of P 3d-metal atoms supported on a metallic substrate is described by a set of values of N and Mj, i = 1,..., P, which need to be found from self-consistency iterations for any given orientation of the magnetic moments. The number of d-electrons, N, and the magnitude of the magnetic moments, M^ can be obtained from the Green's function, S(e - is) = [e - is - HMF]-1, s = +0. Equations for the matrix elements of the Green function are derived from

(E - Hmf) S(E) = I, (12)

where

E = e - is. (13)

This gives the following:

(E - ek) S«f (E) - £ UkiS«/(E) = 0, (14)

i

(E - e«) Sj(E) - £ u«YS j (E) - £ UikSa,(E) = ¿j¿^. (15)

l,Y k

By substituting the matrix elements S«f (E) derived from Eqn. (14) into Eqn. (15), an equation for the matrix elements of the d-electron Green function can be obtained as:

(e - e« - £ E^) Sj(E) - £ L + £ E^ ) Sj (E) V k E ev i=i \ k E ev (16)

+£ ui:7szj3(e) = ¿ij'

7

The magnetism of 3d-transition metal systems is mostly determined by the d-electrons. Therefore, only the Green function for d-electrons, G(E) = S(d) (E), will be considered explicitly

in what follows. According to Eqn. (16), G(E) is a resolvent of an effective, mean field Hamiltonian for the d-electrons, H = H^F, with parameters that are renormalized due to the

influence of the itinerant s(p)-electrons. This Hamiltonian is given by:

4

H = £ + £ Vf dLj - 1 £ U (N - M2) , (17)

where

Ef = E0 + y (Ni - a cos , (18)

Ui 2

Here, E0 is a renormalized energy of unperturbed d-states:

Vf = U (¿^ - 1) ¿j exP (—ai^i) sin 0iMi + (1 - ¿j) Vj. (19)

E0 = + Re£ ^, (20)

k

E - £k

which now acquires non-zero width due to the s(p)-d hybridization and the width parameter, r, is given by

r ^ Im £ . (21)

k

r is included in the argument of the Green function when computing the density of states [17]. The Vj are referred to as hopping parameters. They represent both a direct exchange between d-states localized on sites i and j and as well as a contribution from indirect d-s(p)-d coupling through the conduction band

Vj = Uj + £ EkUkj. (22)

k

Following [18], it is assumed that the on-site s(p)-d coupling is stronger than d-s(p)-d interaction of d-electrons at different sites and, hence, the imaginary part of Vj is neglected. E0, U, Vj and r are assumed to be constant parameters for the model. The choice of values for E0 and U depends mainly on the type of atom i, while the hopping parameters Vj also depend on the geometry of the system, in particular the distance between atoms i and j.

The number of d-electrons, N, and the magnitude of the magnetic moments,Mj, defined in Eqs. (10) and (11), are expressed in terms of the Green function in the local frame of reference, G(S), using the following standard relations:

o

Nj = 1 f de ImTr G* (e - ir), (23)

o

Mi = 1 j de ImTr

aGii(e - ir) . (24)

Here, Gn(e — ir) is a 2x2 matrix with elements G°f (e — ir) and is the z component of the Pauli spin matrix. It is assumed that the magnetic system is in contact with a bath of s(p)-electrons which fixes the Fermi level. The zero of energy is set to be at the Fermi level (£f = 0).

In order to be able to use Eqs. (23) and (24) in self-consistency calculations, the Green function has to be transformed from the local reference frame to a global one using the spin-2 rotation matrix U, Eq. (6):

(5(8) = UG(£)Ut. (25)

By substituting this relation into Eqs. (23) and (24), one obtains the following:

0

Ni = 1 J de Im [C++.+(e - ir) + C-- (e - ir)], (26)

—<x

0

Mi =1 [ de Im [C++(e - ir) - C--(e - ir)] cos 9.

(27)

7T

— oo

0

+ 1/ de Im [C+-(e - ir)e^ + C-+(e - ir)e-i<^] sin9..

—oo

Since the matrix elements of the Green function depend on N and M, Eqns. (26) and (27) should be solved self-consistently for a given orientation of the magnetic vectors starting from some initial estimate for the number of d-electrons and magnitude of the magnetic moments. If the system consists of more than one atom, the self-consistency procedure contains two nested loops. The inner loop involves finding the constrained solution of Eqns. (26) and (27) for each individual atom, while the number of d-electrons and magnitude of magnetic moments at all other atoms are kept fixed. Values of N and M found in the output of the inner loop are then used as the revised constraints for the next iteration of the outer loop. The iterations are repeated until N and M change by less than a predefined magnitude in an iteration. The implementation of the self-consistency procedure is illustrated in Fig. 1.

It is important to realize that this procedure is the same for both stationary and non-stationary orientations of the magnetic momentum vectors and that only the number of d-electrons and magnitude of the magnetic moments are modified during the self-consistency calculation. The orientations of the magnetic vectors remain unaffected, i.e. spin rotations are completely decoupled from the self-consistency procedure in the NCAA model. This is different from DFT calculations, where the orientations of magnetic moments at a non-stationary point are modified during self-consistency calculations, unless local constraining fields are introduced [15].

After self-consistency has been achieved, the total energy of d-electrons can be found

from:

0

E = 5 J dee ImTr C*(e - ir) - 5 ^ U (N*2 - M*2), (28)

where the factor of 5 is due to the five-fold degeneracy of the d-orbitals. Quantities marked with an asterisk correspond to self-consistent values.

In practice, the procedure described above implies that integration over the density of d-states needs to be carried out repeatedly (see Eqs. (26) and (27)). An efficient approach for this has been described in the literature [19-21]. It is reviewed in Appendices A and B for completeness. Initially, the recursive Green's function method is applied in order to represent the Green function in terms of a continued fraction (see §3 in [19] and §2 in [20]). Then, the

JV(0),M<°)

outer loop

inner loop inner loop

JD t ^^^ • • • NV

¿B <b ^

Fig. 1. An illustration of the self-consistency procedure. The magnetic system consists of P atoms whose magnetic moments are denoted by arrows. The direction of the magnetic vectors is fixed. Self-consistent values of the magnitude, M*, of the magnetic moments and the number of J-electrons, N*, are found using the nested loop scheme. In the inner loop, the constrained solution of Eqns. (26) and (27) is found for each atom separately, keeping the number of J-electrons and magnitude of the magnetic moments for all other atoms fixed. After the k-th iteration, the inner loop produces a set of N(k) and M(k) values that are then used as revised constraints in iteration (k+1) of the outer loop. Self-consistency is reached when all the N and Mj do not change upon further iterations.

continued fraction is expanded in a series of partial fractions [21]. A matrix element of the Green function then takes the form (see Appendix B)

GOT00 = £ (29)

k qk

where the numbers pk, qk depend on the orientation of the magnetic moments as well as on indices i, j, a and p. The density of states is expressed in terms of Lorentzian functions and the integration then carried out analytically. Namely, since the J-levels contain an imaginary part, ir, integrands in Eqs. (26) and (27) include the terms

Pk ^ Pk

ImV-= r

—' f — q — iF

fc f - qk- v (f - qk)2 + r2'

that can be integrated analytically to give

0

IdfT E (f - J)2 + F2 = Epk arcc0t (T

The total energy of the system is also expressed analytically in terms of parameters of the self-consistent Hamiltonian as:

0 2P ^ P

E=;/dff E (7-^- 5 E U (N*2 - M*2)

M V ßJ (30)

5 _2P_ r f* P /f*2 \ 1 P tt

= - E <• arccot f + 2 'M P2 + 0 - 5 E J (N*2 - M*2) ,

ß=1 L ^ ' 1 i=1

where e* are the eigenvalues of H*. Integral in Eq. (30) diverges at the lower limit. However, this infinite contribution is the same for all magnetic states and, therefore, can be omitted.

In the self-consistency procedure, the number of d-electrons and the magnitude of the magnetic moments are found for a given orientation of the magnetic momentum vectors. This assumes a hierarchy of relaxation times. The relaxation of the diagonal components of the spin density matrix, which give the number of d-electrons and magnitude of the magnetic moments, is assumed to be much faster than the relaxation of the off-diagonal components which give the orientation of the magnetic momentum vectors [11]. Thus, N and M, are treated as fast degrees of freedom that adjust instantaneously to changes in the orientation of the magnetic moments, the slow degrees of freedom. This is analogous to the Born-Oppenheimer approximation for electronic and nuclear degrees of freedom.

3. Magnetic force theorem

Gradients of the energy with respect to the orientation of the magnetic momentum vectors, i.e. magnetic forces, can be used to guide the search for minimum energy configurations of the magnetic moments, identify minimum energy paths, and to simulate dynamics of magnetic systems. A simple approach for evaluating the energy gradient is to approximate it by finite differences where the energy is evaluated self-consistently for slightly different orientations. However, this is inefficient for large systems as at least 2P + 1 self-consistency calculations need to be carried out for a system with P magnetic moments each time the force is evaluated.

Below, a method is presented for evaluating all components of the energy gradient without having to perform any additional self-consistency calculations. The method is based on a magnetic force theorem (MFT) for the NCAA model. The theorem states that for an arbitrary orientation of the magnetic vectors, stationary or non-stationary, the energy as a function of the fast degrees of freedom, N and M, reaches an extremum at the self-consistent values, N * and M *

dE

Vi :

dNi

= 0 and _ „ ^

N=N * dMi

M=M *

= 0. (31)

N=N * M=M *

This simplifies greatly the calculation of the gradient and speeds up exploration of the energy surface of the system.

Two lemmas for the Green function that are proved in Appendix C are needed to derive the MFT:

= -UI (Gi+ G--«) • (32)

and

dTrG(e) Ui d r(Gii( ) )) „

=Ti(e) - G— (e)) (33)

+ (G++-(e)e^ + G-+ (e)e—^) sin.

According to eqns. (28) and (32)

0

H=1/<fa Im dN t g(£ - ir) -5 U Ni

0

5 U

2

1 f de Im (G+ii(e - if) + G—— (e - if)) - Ni

where integration by parts has been invoked. According to Eqn. (26), the expression in square brackets is equal to zero when self-consistency has been reached.

The equation for the derivative of the energy with respect to Mj in Eq. (31) is proved in the same way. Using (28) and (33), we obtain the following:

0

H = U deeim dMTr G(e - ,r)+s " m

—TO

✓ 0

2 1 n

5-1 i de im [(G++(e - if) - G—— (e - if)) cos

—to

+ (G+— (e - ir)ej^i + G—+ (e - ir)e—sin 6j + Mj

Due to Eqn. (27), the expression in curly brackets vanishes when M has the self-consistent value, M*.

By using the MFT, the computational effort involved in the calculation of magnetic forces in NCAA is significantly reduced. According to the MFT, the energy change due to the infinitesimal rotation of magnetic vectors does not contain a contribution from the variation of the number of d-electrons and the magnitude of the magnetic moments, i.e. fast degrees of freedom. In particular, this means that in the finite difference scheme, the energy of the perturbed spin state obtained by a small rotation of magnetic vectors from a particular orientation can be approximated using the same values for the fast degrees of freedom as those calculated self-consistently for the unperturbed spin state and, therefore, no additional self-consistency calculations are needed.

A simple finite difference scheme is, however, problematic because it involves evaluating the difference between two numbers of similar magnitude. This procedure can lead to a significant loss of accuracy. Instead, first order perturbation theory can be used to derive an expression for the magnetic force in terms of the self-consistent values. According to the MFT, a derivative of the energy, E = E(A), with respect to any adiabatic parameter A (a slow degree of freedom) can be computed from the explicit A dependence only, without having to include implicit dependence

0

dE(A) = dEA) = 5 f dtt Im T 3G'- ^ A) (34)

dA SA W 3A ' ( )

— TO

Here, 5G*(e - ir; A)/5A can be found by using the resolvent identity

^ = G(e; A)^G(e; A), (35)

which, together with Eqn. (34), gives

0

( ) — [ dee im Tr

dA n

— TO

G*(e - ir; A)d ;A) G*(e - ir; A) dA

0 (36)

— I de im Tr n J

—TO

G'(e - iF; A)^

that is, the derivative of the total energy with respect to a parameter coincides with the expectation value of the derivative of the Hamiltonian with respect to that parameter. This is analogous to the Hellmann-Feynman theorem.

In practice, it is convenient to calculate the trace in Eqn. (36) using the basis in which H*(A) and G*(e; A) are diagonal

dH *(A)n 2P

Tr

G*(e; A)-

dA

t *

£ ■ (37)

where are the diagonal elements of dH*(A)/dA in the relevant basis. The integral in Eqn. (36) can then be evaluated analytically leading to

^ = 5 £ « (r ). (38)

With A = Oi or A = 0i and i = 1,..., P, this gives the gradient of the energy with respect to the angles defining the orientation of the magnetic moments.

The procedure for evaluating the energy gradient is as follows: First, derivatives of the self-consistent Hamiltonian, dH*(6, ^)/d9i and dH*(6, 0)/d0i, which are given explicitly by

dH*(6 S) \ af 1

-¿r^ = 1 jiSkiUiM* [a^f sin Oi + {Saf - 1) exp (-ai0j) cos Oj , (39)

ddj Jkj 2

(^Hd^) a = t (1 - } sin Oj (40)

are transformed to a basis where H* (6, 0) is diagonal. Then, their diagonal matrix elements, t*(Oi), are inserted into Eqn. (38) and the derivatives with respect to O and 0 evaluated.

The MFT significantly reduces the computational effort involved in calculations of the first derivatives of the energy with respect to O and 0. Second derivatives of the energy are also important, because they are used to calculate magnetic exchange interaction parameters, Jj, between magnetic moments i and j. However, the MFT is not valid for the second derivatives. One cannot neglect the change in the fast degrees of freedom when computing the second variation of the energy.

4. Applications

4.1. Trimers adsorbed on a metallic substrate

Trimers of 3J transition metal atoms supported on a metallic surface are good test systems for studying noncollinear magnetism at the atomic scale. Noncollinear magnetic ordering in Cr, Mn, Fe trimers was predicted theoretically using a model Hamiltonian approach [24-26] and was later obtained in ab initio calculations [27,28,30,31]. The angles formed between the magnetic moments in the stable magnetic states depend on the type of atoms in the trimer and the geometrical arrangement of the atoms. Without spin-orbit interaction, the configuration and spin spaces are uncoupled. For coplanar magnetic ordering, two angles completely determine the configuration of the three magnetic momentum vectors. Therefore, the energy surface for such systems can be visualized easily.

Within the NCAA model, the parameters and Ui are determined mainly by the chemical element whereas the hopping parameters Vjj depend also on the arrangement of the atoms in the trimer and the hybridization of the 3J-states with electronic states of the substrate. The parameter associated with the width of the J-states, r, may depend on the position of the

trimer on the surface as well as on the trimer geometry. The unit of energy is taken to be r and we give below the values chosen for the scaled energy parameters, U = U/r, E0 = E0/r and

Vj = Vj/r.

Magnetic ordering in the Cr, Mn and Fe trimers is expected to be quite different. To illustrate this, we choose the Coulomb integral, U = 13, and the hopping parameters, V12, V13 and V23, to be the same for all the trimers, while the value of E0 is chosen to ensure that the number of d-electrons per atom in each case is the same as in the corresponding solid (N « 5, 6 and 7, for Cr, Mn and Fe, respectively). The definitions of the two variables, the angles 62 and 63, are shown in Fig. 2.

Fig. 2. A noncollinear magnetic state of a trimer of atoms adsorbed on a surface. The definitions of the atom labeling, 1, 2 and 3, are shown as well as the two angles, $2 and 03, specifying the two degrees of freedom. All three distances between the atom pairs are different, leading to different values for the three hopping parameters, Vj.

The calculated energy surfaces for the three trimers, Cr3, Mn3 and Fe3, are shown in Fig. 3. The stable configurations of the magnetic momentum vectors, corresponding to minima on the energy surface, are shown, as well as MEPs for transitions between these states calculated using the climbing image nudged elastic band (NEB) method [32,33]. The MEP gives the mechanism for low energy collective excitations of the magnetic trimers and the maximum energy along the MEP gives the minimum energy needed for the transitions.

A Cr crystal has spin density wave magnetic structure with antiferromagnetic coupling of near neighbor spins. In an equilateral Cr trimer, the ground state is noncollinear with angles 2n/3 between the magnetic momentum vectors [24-27]. However, a collinear ordering of the magnetic moments is also possible for a different choice of distances between the Cr atoms and, thereby, the hopping parameters Vj. Such a case is shown in Fig. 3a, where V23 = 1.0, V12 = 0.6 and V13 = 0.9. Maxima on the energy surface correspond to ferromagnetic ordering of two of the magnetic moments, while the magnetic moment of the third atom is pointing in the opposite direction. All minima on the energy surface correspond to the same magnetic state. There are two different MEPs for transitions between these states, as shown in Fig. 3a. One of the MEPs corresponds to an almost coherent rotation of the magnetic moments of atoms 2 and 3, maintaining a nearly antiparallel ordering during the rotation. At the saddle point, both moments have reversed their direction. The activation energy for this transition is 5.4 ■ 10-2 in r-units. Along the second MEP, the magnetic moment of atom 2 rotates monotonously by an angle of 2n, while the magnetic moment of atom 3 only rotates slightly from the initial orientation and then returns. At the saddle point, the magnetic moment of atom 3 is aligned in the same direction as in the ground state but the moment of atom 2 has been reversed. The activation

—71 0 63 71 ^

Fig. 3. Energy surfaces of supported trimers: (a) Cr3, (b) Mn3 and (c) Fe3. The zero of energy (red) is the energy for the lowest energy configuration of the magnetic moments. White lines represent minimum energy paths and the dots indicate positions of images in the NEB calculations. Minima are marked with circles, while saddle points are indicated with crosses. The insets show the magnetic momentum vectors at the minima and at the saddle points.

energy is larger for this second mechanism, 7.8 ■ 10-2 Г, because of strong antiferromagnetic coupling between atoms 2 and 3. Note that all minima and saddle points on the energy surface correspond to collinear ordering of the magnetic moments.

For the Mn trimer, we obtain two distinct stable states, both corresponding to noncollinear magnetic ordering, as shown in Fig 3b. Mn atoms often form noncollinear structures [24,26,28]. The two stable states have the same energy, but the MEPs between them are quite different and give different activation energies. The saddle points correspond to collinear magnetic ordering. The activation energy is 1.9 ■ 10-2 Г for the saddle point at 02 = п, 03 = п, 1.7 ■ 10-2 Г for the saddle point at 02 = п, 0з = 2п and 0.7 ■ 10-2 Г for the saddle points at

02 = 0, 03 = п.

For the Fe trimer, all the minima on the energy surface correspond to the same ferromagnetic state (see Fig. 3c). The saddle points between those states also correspond to collinear magnetic ordering, but with one of the magnetic moments pointing in the opposite direction to the other two. The energy surface for the Fe trimer is similar to the one for the Cr trimer, but shifted by п along the 03 axis. There are also two possible MEPs between the states. The MEP with the lower activation energy, 2.7 ■ 10-2 Г, corresponds to nearly coherent rotation of the second and the third magnetic moments. The other MEP corresponds to relatively small deviation of 03 from zero. The barrier for this transition mechanism is higher, 4.5 ■ 10-2 Г.

4.2. STM tip interacting with adsorbed dimer

The energy surface of a magnetic cluster can be modified by applying a local external magnetic field. Another way to deform the surface is by hybridization with an additional atom placed close enough to the cluster. The latter effect was recently demonstrated experimentally by using magnetic atoms at the tip of a scanning-tunneling microscope (STM) [34]. The magnetic tip was used to modify the energy barrier between two magnetic states of a magnetic nanostructure on a non-magnetic substrate.

As a simple model of this experiment, we consider the isosceles Cr trimer, where atom 1 represents the STM tip and atoms 2 and 3 represent the dimer adsorbed on a non-magnetic substrate. The hopping parameters between the tip atom and each one of the two dimer atoms are set to be the same, V12 = V13 = V, and are varied to mimic changes in the distance between the tip and the dimer, while the third hopping parameter is kept fixed with a value of V23 = 1. For each value of V, the energy of the system is minimized with respect to the angles 02 and

03 defining the orientation of the magnetic moments in the dimer. Due to the symmetry of the system, 02 differs from 03 by an angle of п at the minimum energy configuration. The value obtained for one of the angles, 02, is shown as a function of V in Fig. 4. For small V, the magnetic moments of the dimer point in the opposite direction and are perpendicular to the magnetic moment of the tip. However, when V becomes larger than V23, the magnetic moments of the dimer tend to point in the same direction and opposite to the magnetic moment of the tip. The energy surface of the system for three different values of y is shown in the insets of Fig. 4.

5. Finite range approximation

In the self-consistency calculations, the computational effort of the tridiagonalization of the NCAA Hamiltonian which is needed to evaluate the matrix elements of the Green function scales as P3, where P is the number of magnetic atoms in the system. This operation has to be performed for each magnetic atom, so the computational effort in each iteration of the self-consistency procedure scales as P4. The number of iterations needed to reach self-consistency

0.7 0.9 1.1 1.3

v/r

Fig. 4. Angle between the magnetic moment of an STM tip and the magnetic moment of one of the Cr-atoms in a Cr2 addimer as a function of the corresponding hopping parameter. The three insets show the energy surface of the system for the values indicated by the three gray circles on the curve. White and black dots in the insets mark energy minima.

could also have some dependence on the system size, but is difficult to assess. The self-consistency calculations become too demanding when the number of non-equivalent magnetic atoms becomes large, on the order of thousands. One way of addressing this problem is to use parallel computing as the self-consistency algorithm can be parallelized efficiently. Another option is to use approximations in order to make the computational effort scale linearly with P.

One approach is to assume that the magnetic state of an atom is mostly affected by its neighbors, while the effect of distant atoms can be neglected in the self-consistency calculations. More specifically, when calculating the number of d-electrons and the magnitude of the magnetic moment of a particular atom, only neighbors within a sphere of a certain radius, rc, centered on the atom are included. The accuracy of this approximation can be controlled by varying the radius, rc. Typically, the calculated results converge quite fast as rc increases. This is demonstrated below in calculations for the electronic and magnetic structure of a monolayer island of Fe-atoms on a W(110) surface.

The position of the Fe-atoms in the monolayer are determined by the atomic structure of the W(110) substrate. Each Fe-atom has four nearest neighbors at a distance of \[3/2a, where a is a lattice constant of the W-crystal (see Fig. 5(a)). For the Fe monolayer, the number of atoms inside the sphere, Nc, changes stepwise as a function of rc as shown in Fig. 5(b). The position of the steps, r(z), Z = 0,1, 2,..., is defined by the distance between the neighbors.

Any value of rc in the range will result in the same self-consistent values for the

number of J-electrons and magnitude of the magnetic moments. Below, we will refer to Z as the number of shells of interacting neighbors and study how the choice of Z affects the results of the self-consistency calculations.

(a) • • • •

(b)

m • *

vW1) rc \

• 1 •

[110] \

• •

—► [001]

a

35 30 25 20 15 10 5 0,

[

r(!) c r(2) ' c r(3) ' c

J

I H

[ —' J

rc, A

Fig. 5. (a) Positions of Fe-atoms in a monolayer of Fe on a W(110) surface. (b) Number of Fe atoms inside a coordination sphere of radius, rc, in a Fe monolayer on a W(110) surface.

The parameters E0 and U in the NCAA model were chosen to have values that are typical for Fe: -12 and 13, respectively [7]. For simplicity, only nearest neighbor hopping parameters were included and the value V = 0.9 was chosen to reproduce DFT calculated magnetization of an Fe overlayer on a W(110) surface [35].

A self-consistent calculation was carried out for a relatively large Fe island containing 17 x 17 atomic rows. Observe that due to the structure of the (110) surface the total number of atoms in the island is not equal to the product of the number of atomic rows along [001] and [110] directions. Reference values of Nref and Mref at each Fe-atom were first calculated by including all 144 atoms explicitly. Then, approximate calculations were carried out for the same orientation of the magnetic vectors. The coordination sphere of a particular size was chosen for each atom in the island and corresponding self-consistent values of N and M evaluated and compared with the reference values. Fig. 6(a) shows the deviation in the number of J-electrons and the magnitude of the magnetic moments from the reference values as a function of the number of coordination spheres included. For both N and M, the error drops rapidly as Z increases.

The MFT is a rigorous statement only when all atoms in the system are included explicitly. When the coordination sphere approximation is used, the theorem is, strictly speaking, not satisfied and the self-consistency state is no longer a stationary point of the energy as a function of the number of J-electrons and magnitude of the magnetic moments. This is illustrated in Fig. 6(b) where the derivative of the energy with respect to N and M as a function of Z is shown. The magnitude of the derivative, however, vanishes as Z increases.

Thus, the coordination sphere approximation can provide reasonably accurate results while significantly reducing the computational effort. After the system becomes larger than the

(a)

0.12 0.1 0.08 0.06 0.04 0.02 0

♦ X

x * 1 * * * _

0.12

10.1 *

0.08 ^

<D

10.06 "7

o

6 8

10

0.04 0.02 0

c

(b)

0.5

0.4

fe; 0.3

CO

0.2

0.1

0.5

» * * «

8 10

c

Fig. 6. (a) Maximum deviation of the number of d-electrons (crosses, and left vertical axis) and the magnitude of the magnetic momentum vectors (diamonds, and right vertical axis) from the reference values, Nref and Mref, as a function of the number of coordination shells included in the self-consistency calculation. (b) Maximum derivative of the energy, E, with respect to N (crosses) and M (diamonds), as a function of the number of coordination shells included in the self-consistency calculation.

chosen number of coordination shells included in the self-consistency calculation, the computational effort increases only linearly with the system size. The level of approximation can readily be increased by increasing Z, i.e. including more coordination shells.

With this approximation, the NCAA model can be used to study large systems including thousands of non-equivalent atoms.

6. Summary

A detailed description of the NCAA model and its implementation is given in this article as well as an application to trimers of Cr, Mn and Fe adsorbed on a metal surface and a simple model for STM tip interaction with a dimer. While the Cr and Fe trimers are found to have collinear stable states, the former with antiferromagnetic ordering and the latter ferromagnetic, the Mn trimer has two different stable states, both non-collinear. In all three cases, the maximum energy configurations along the minimum energy path, i.e. first order saddle points on the energy surface, correspond to collinear ordering. These simple systems were chosen as illustrations of the method because the energy surface can be visualized easily. A finite range approximation is also presented which leads to linear scaling of the computational effort with the number of magnetic atoms in large systems. The theoretical approach described here can be used in studies of magnetism in 3d-metal systems, including simulations of spin dynamics and, in general, navigation on the energy surface characterizing such systems to find stable state and transition state for magnetic transitions.

Acknowledgements

This work was supported by the Government of Russian Federation (Grant No. 074-U01), RFBR Grants No. 14-02-00102, and No. 14-22-01113 ofi-m, the Icelandic Research Fund, University of Iceland research fund and the Nordic-Russian Training Network for Magnetic

Nanotechnology (NCM-RU10121). PB greatly acknowledges support from Goran Gustafsson Foundation.

APPENDIX

A. Recursion method and continued-fraction expansion of the Green function

An efficient algorithm for calculating the matrix elements of the Green function G (in this Appendix, a dependence of G on e is implied) using the recursion method involves finding a new orthonormal basis in which the Hamiltonian has a tridiagonal form, while the matrix element to be evaluated, G°f, remains unchanged.

We begin by discussing how diagonal matrix elements,

G

mm

| G |xM), can be

evaluated. Here, (|x1),..., |x2P)} is an initial basis set and indices v enumerate both atomic site and spin projection. The first vector of the new basis is chosen so that |y1) = |xM). Then, the diagonal matrix element of the Green function is given by:

= (yi| G |yi) = G11,

where tilde indicates a quantity in the new basis. Other vectors of the new basis {|y1) are found given that the Hamiltonian has a tridiagonal form:

(41)

\V2P)}

( «1 bi 0

H,

mv

bi a2 b2

0

b2 «3

0 0

b3

0 0 0 0 0 0 0 0 \0 0 0 0

0 0 0

0 0 0

\

. «2P—2 b2P—2 . b2P—2 «2P— 1

b

0 0 0

0

2P-1

0

(42)

b2P —1 «2P )

Not only basis vectors |yM) are found within the recursion scheme, but also elements {a1, a2,..., a2P} and {b1, b2,..., b2P-1} are successively calculated. From the equation:

H |y1) = a1 |y1) + b1 |y2) ,

the second basis vector |y2) as well as a1 and b1 can be found as follows:

a1 = (y1| H |y1) , |H |y1)- a1 |y1>H , H |y1) - a1 |y1)

b1 = \îfc)

b1

(43)

(44)

(45)

(46)

Then, the Hamiltonian acts upon the second basis vector and the next portion of elements of the tridiagonal matrix is found. This operation is reapplied until all elements of H in the new basis have been calculated. For example, on the k-th step (k < 2P), when a1,... ,ak-1, b1,..., bk-1 and |y1),..., |yk) are known, we have:

= (Vk1 H \yk) , bk = ||H \vk) - bk— 1 \Vk—1 ) - «k \Vk)

|Vk+1)

H \vk) - bk—1 \Vk—1) - «k \Vk) bk '

(47)

(48)

(49)

Gii is found in the following way. Let us introduce a notation:

Ak = <yi| G \yk). (50)

The following set of equations is then valid for Ak, which follows from the definition of the Green function, (e — H) G = I:

(e — ai) Ai — 61A2 = 1, (e — ak) Ak — bk-iAk-i — bkAk+i = 0, k =1,k = 2P, (e — a2p) A2P — b2p—iA2P-i = 0.

After straightforward algebra, we obtain the representation of Gii in the continued-fraction form:

e — ai

1

"b1

Gßß — Gii — Ai —-jr"2--(51)

e —

b2

b2P-1

e — a2P

The nondiagonal elements of the Green function can also be expanded in terms of continued fractions. Namely, if , ^ = v, is to be found, one has to perform tridiagonalization four times, where the starting vectors are:

|y?> — + |Xv »

1

72

|yî> — 7 (M +i |xv»

— — i |Xv>) .

In each basis, matrix element G11 — (y{

(

rn

ued fractions as described above. Continuec part of the nondiagonal element, GßV, is thus given by:

1 2

yz), Z = a, b, c, d, is expanded in terms of contin--fraction representation of the real and imaginary

Re Gßv — ^ ( GGai — G?i) , (52)

Im Gßv — - (3?i — G^) . (53)

B. Partial-fraction expansion of the Green function

We proceed with the continued-fraction representation of the matrix element of the Green function, Eqn. (51). Let us consider the last level of the continued fraction:

f (1)(e} = e - a2p- - . (54)

e — a2P

If b2P-1 = 0, then the equation f (1)(e) = 0 has two real roots, q(2) and q22), q(2) < q(2), which can be found either analytically or numerically. Then:

f<1»(e) = il^^^e - Of

e - qS1'

where q(1) = a2P. As a result, the next level of the continued fraction acquires the form:

b2p-2 (e - qS1')

f (2)(e) = e - a2P-2 -

e - qS2M e - q<2)

The ratio of polynomials here can be transformed into a sum of two partial fractions:

b2p-2 (e - q(1)) _ P12) + p22)

where

e - if) (e - q22)) e - q<2) e - qf'

q(2) q(1) q(2) q(1) p(2) = b 2 q1 - q1 p(2) = ,2 q2 - q1

p1 = b2P-2 (2) (2) , p2 = b2P-2 (2) (2) '

q() - q2) q2) - q()

Thus, the function f 2(e) is represented in the form analogous to Eqn. (54):

(2) (2)

f (2)(e) = e - a2P-2 - - "T", (55)

e - q1 e - q2

and the number of levels in the continued fraction is reduced by one.

The same technique is sequentially applied in order to transform the continued fraction into a sum of partial fractions. For example, at the k-th step we have:

k P(fc)

f (k)(e) = e - a2P_fc -j) • (56)

j=1e - qj)

Zeros of the function f (k)(e) are well separated (see Fig. 7) and can be found numerically without problems.

The function f (k+1)(e) is then represented as follows:

k+1 p(k+1) f (k+1)(e) = e - a2P-k-1 - £ •

j=1e - qj

Here q(k+1), • • •, q(.++11) are the roots of f(k) (e) = 0 and pjk+1) are given by:

P(k+1) = b2 p7 = b2J

k n i=1 (qf'+1) - qf)

k-1k+1 ( n i=1 ?+1) - q(k+1)

As a result, any continued fraction of the form (51) and thus an arbitrary element of the Green function can be expanded in a sum of partial fractions (see Eqn. (29)).

Fig. 7. Graphical solution of the equation f (k)(e) = 0. Exactly one root is located between poles q(k),..., qkk) of f (k)(e).

C. Proof of lemmas (32) and (33)

We proceed with the equation for the matrix elements of the Green function:

(e — Ef) Off (e) — £ G j(e) = ^. (57)

k,Y

This equation can be solved iteratively giving the following result for the diagonal elements of the Green function:

1 V"1 V Ya

Gaa( ) _ 1__| \ " _V ik V ki_ |

* ( )= e — Ef + ¿k (e — Ef)(e — E^)(e — Ef) + ...

V aYi V71Y2 V Yna (58)

i x iki kik2 kni

+ ¿Yn (e — Ef)(e — EY1) ... (e — E£)(e — Ef) + ...

k i ,...,kn

where summation over ki,k2,...,kn and 7i,72,...,7n runs over all atomic sites and spin projections, respectively. Each term in Eq. (58) can be considered as a sequence of vertices forming a path with endpoints at (a, i):

VfYl V?1?2... VJ"a

(a, i) ^ (yi, ki) ^ (72, k2) ^ ... ^ (Yn, k„) ^ (a, i) = ikl klk2 kni

(e — Ef )(e — EY1)... (e — Ef)"

(59)

Therefore, Gfia(e) can be thought of as the sum of all such closed paths. Each path describes some multiple scattering process, where a single scattering event ('hopping') can be either an on-site spin-flip or a non-spin-flip intersite transition. The sum of all possible paths of the same

length and endpoints fixed at (a, i) can be represented as:

E

ki ,...,kn

v ?71 VY?2 ... V7na iki kik2 k;i

(e - Ef)(e - EY1 )... (e - E£ )(e - Ef)

[Go(e) VGo(e)V... VGo(e)]^ STte n).

V

n+1 V-factors

(60)

Here, V is a hopping matrix (see Eq. (19)) which contains information about all scattering events in the system and G0(e) is the Green function for the non-interacting magnetic impurities. G0(e) is diagonal and its matrix elements are given by:

[Go(e)r =

(e - Ef )■

(61)

In order to prove lemmas (32)-(33), it is sufficient to derive corresponding identities for the sum of the paths, S(e; n), for an arbitrary path length, n. The derivative of Tr S(e; n) with respect to N is given by:

d d

Tr S(e; n) = — Tr [Go(e)VGo (e) V ... VGo(e)] = Tr

dNi

dN, + Tr

dGo(e) dNi

VGo(e)V... VGo(e)

Go(e)VdGo(e) V... VGo (e)

Only two matrix elements of

dNi dGo(e)

dNi

+ ... + Tr

Go(e)VGo(e)V...V

dGo (e)

dNi

(62)

are non-zero (see Eq. (61)). According to Eqs. (61) and

(18), they can be expressed in terms of derivatives with respect to e :

dGo(e) Ui dGo(e)

dNi

2 de

so

where is 2P x 2P matrix defined so that only the elements (++, ii) and (to unity, while the others are zero. Therefore, Eq. (62) takes the form:

(63)

ii) are equal

d U

m Tr S (e; n) = - f<Tr

dGo(e)^ o

de

SoVGo(e)V...VGo(e)

+ Tr + ... + Tr

Go(e)VdGo(e)SoV...VGo(e)

de

Go(e)VGo(e)V...V

(64)

dGo(e)^ o

The factor dGo(e) ^0 found in every term in the curly brackets can be transformed into:

de

^Ge^so = -Go(e)SoGo(e),

where use has been made of the identity

dGo(e) de

- Go(e),

(65)

(66)

1

as well as commutativity of G0(e) and S0. Using Eqs. (65) and (66) and invariance of a trace under cyclic permutations, each term in Eq. (64) is transformed to obtain the following:

d U

m Tr S (e; n) = — f<Tr

S0 ^ VGo(e)V...VGo(e) de

+ Tr

+ ... + Tr

U- d

E0Go(e)VdG0(e) V...VGo(e)

de

E0Go(e)VGo(e)V...V

U d

dGo(e) de

U d

— T^ Tr [s0Go(e)VGo(e)V... VGo(e)]

= — y^ Tr [SOS(e; n)] = — (S++(e; n) + S--(e; n))

2 de

(67)

which proves Lemma (32).

Lemma (33) is proved in a similar way. However, when computing the derivative of Tr S(e; n) with respect to Mi, one needs to observe that both G0(e) and V depend on Mi. More specifically:

dGo(e) Ui dGo(e)E

"dMT = Tcos Ez,

dM = — U sin Q(cos + sin ). The only nonzero matrix elements of EX, Ey and Ef are defined as follows

(SX)+- = 1

- 'H

(sy )+

++

— 1

(S )++ = 1

(SX)- + = 1,

(sy )- + = 1,

(S )-- = —1.

(68)

(69)

(70)

(71)

(72)

Note that if only one magnetic impurity is present in the system, Ef, Ey and Ef coincide with ordinary 2 x 2 Pauli matrices. This gives the following:

d U d

— Tr S(e; n) = U cos Qi— Tr [EfGo(e) VGo(e) V... VGo(e)] +

T T

U sinQi — (cosTr [EfGo(e)VGo(e)V ... VGo(e)] + sin&Tr [EyGo(e)VGo(e)V ... VGo(e)] UU d f

= u~ — < cos Qi Tr [EfS(e; n)] + sin Qi ^cos Tr [EfS(e; n — 1)] + sin Tr [EyS(e; n — 1)]

2 de

(S++(e; n) — S--(e; n)) cos + (S+ -(e; n — 1)e*< + S- +(e; n — 1)e-^) sin 0

(73)

which proves lemma (33).

For large (or infinite) systems, the calculation of Green functions is often performed approximately, using a finite range approximation (see Sec. 5) or by truncating the recursion method (see appendix A) at a certain step. In this case, the scattering matrix V may be different for different sites i and Lemmas (32)-(33) and, therefore, the MFT, are strictly not valid. This was demonstrated in Sec. 5 where direct calculations of the derivatives of energy with respect to N and Mi for the self-consistent values N* and M* were found not to equal zero because of the finite range approximation. However, if any subsystem of atoms is chosen and the same

set of hopping contributions are included for all the atoms of this subsystem, then the scattering matrix will be the same for all atoms in the system and the proof of Lemmas (32)-(33) will hold just as in the case when all atoms and hopping contributions are included. In particular, one can choose as a subsystem the atoms included in one closed loop in Eqn.(59) and the corresponding hopping contributions to form the scattering matrix V. Then, the MFT will be valid for this subsystem. Formally, it is a consequence of the fact that all closed paths in Eqn.(59) which make a contribution to the matrix element of the Green function G^, make the same contribution to any Gj if the path in Eqn.(59) goes through site j. If, however, we try to improve the calculation by taking into account additional atoms and hopping contributions that are different for different atoms, as in the case of the finite range approximation described in Sec. 5, then the MFT is violated.

References

[1] R. Wiesendanger. Spin mapping at the nanoscale and atomic scale. Rev. Mod. Phys., 81, 1495 (2009).

[2] S. Heinze, K. von Bergmann, et al. Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions. Nature Phys. 7, 713 (2011).

[3] M. Bode, O. Pietzsch, et al. Experimental Evidence for Intra-Atomic Noncollinear Magnetism at Thin Film Probe Tips. Phys. Rev. Lett. 86, 2142 (2001).

[4] P.F. Bessarab, V.M. Uzdin, H. Jonsson. Harmonic transition state theory of thermal spin transitions. Phys. Rev. B 85, 184409 (2012).

[5] H.B. Braun. Topological effects in nanomagnetism: from superparamagnetism to chiral quantum solitons. Advances in Physics 61, 1 (2012).

[6] P.F. Bessarab, V.M. Uzdin, H. Jonsson. Size and shape dependence of thermal spin transitions in nanoislands. Phys. Rev. Lett. 110, 020604 (2013).

[7] P.F. Bessarab, V.M. Uzdin, H. Jonsson. Potential Energy Surfaces and Rates of Spin Transitions. Z. Phys. Chem. 227, 1543 (2013).

[8] T. Oda, A. Pasquarello, R. Car. Fully Unconstrained Approach to Noncollinear Magnetism: Application to Small Fe Clusters. Phys. Rev. Lett. 80, 3622 (1998).

[9] O. Ivanov, V.P. Antropov. Molecular magnetism: Noncollinear ordering and spin dynamics. J. Appl. Phys. 85, 4821 (1999).

[10] L. Nordstrom, D.J. Singh. Noncollinear intra-atomic magnetism. Phys. Rev. Lett. 76, 4420 (1996).

[11] V.P. Antropov, M.I. Katsnelson, et al. Spin dynamics in magnets: equation of motion and finite temperature effects. Phys. Rev. B 54, 1019-1035 (1996).

[12] L.M. Small, V. Heine. A couple method for calculating interatomic interactions in itinerant electron magnetic systems. J. Phys. F 14, 3041 (1986).

[13] G. Bihlmayer, 'Density-functional theory of magnetism', in 'Handbook of magnetism and advanced magnetic materials', ed. H. Kronmuller, S. Parkin. Vol. 1, page 3.

[14] L.M. Sandratskii. Noncollinear magnetism in itinerant-electron systems: Theory and applications. Advances in Physics 47, 91 (1998).

[15] G.M. Stocks , B. Ujfalussy, et al. Towards a constrained local moment model for first principles spin dynamics. Philos. Mag. B 78, 665 (1998).

[16] S. Lounis, P.H. Dederichs. Mapping the magnetic exchange interactions from first principles: Anisotropy anomaly and application to Fe, Ni, and Co. Phys. Rev. B 82, 180404(R) (2010).

[17] P.W. Anderson. Localized Magnetic States in Metals. Phys. Rev. 124, 41 (1961).

[18] S. Alexander, P.W. Anderson. Interaction between localized states in metals. Phys. Rev. 133, A1594 (1964).

[19] R. Haydock, V. Heine, M.J. Kelly. Electronic structure based on the local atomic environment for tight-binding bands. J. Phys. C: Solid State Phys. 5, 2845 (1972).

[20] R. Haydock, V. Heine, M.J. Kelly. Electronic structure based on the local atomic environment for tight-binding bands: II. J. Phys. C: Solid State Phys. 8, 2591 (1975).

[21] V.M. Uzdin, N.S. Yartseva. Periodic Anderson model for the description of noncollinear magnetic structure in low-dimensional 3d-systems. Comput. Mater. Sci. 10, 211 (1998).

[22] S.V. Uzdin. On the calculation of the magnetic structure of surfaces, near surface layers, and interfaces of 3d metals. Physics of the Solid State 51, 1260 (2009).

[23] P.F. Bessarab, V.M. Uzdin, H. Jonsson. Calculations of magnetic states and minimum energy paths of transitions using a noncollinear extension of the Alexander-Anderson model and a magnetic force theorem. Phys. Rev. B 89, 214424 (2014).

[24] S. Uzdin, V. Uzdin, C. Demangeat. Magnetic trimer on non-magnetic substrate: from frustration towards non-collinearity. Europhys. Lett. 47, 556 (1999).

[25] S. Uzdin, V. Uzdin, C. Demangeat. Non-collinear structure of Cr trimer on the surface of non-magnetic metals. Comput. Mater. Sci. 17, 441 (2000).

[26] S. Uzdin, V. Uzdin, C. Demangeat. Non-collinear magnetism of Cr, Mn and Fe trimers supported on the non-magnetic metal surface. Surf. Sci. 482, 965 (2001).

[27] H.J. Gotsis, N. Kioussis, D.A. Papaconstantopoulos. Evolution of magnetism of Cr nanoclusters on Au(111): First-principles electronic structure calculations. Phys. Rev. B 73, 014436 (2006).

[28] A. Bergman, L. Nordstrom, et al. Magnetic interactions of Mn clusters supported on Cu. Phys. Rev. B 73, 174434 (2006).

[29] K. Hirai. Electronic Structure of Helical Spin Density Wave State in fcc Iron. J. Phys. Soc. of Jpn. 61, 2491 (1992).

[30] S. Lounis, P. Mavropoulos, P.H. Dederichs, S. Blügel. Noncollinear Korringa-Kohn-Rostoker Green function method: Application to 3d nanostructures on Ni (001). Phys. Rev. B 72, 224437 (2005).

[31] R. Robles, L. Nordstrom. Noncollinear magnetism of Cr clusters on Fe surfaces. Phys. Rev. B 74, 094403 (2006).

[32] H. Jonsson, G. Mills, K.W. Jacobsen. In Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B.J. Berne, G. Ciccotti, D. F. Coker (World Scientific, Singapore, 1998), p. 385.

[33] G. Henkelman, H. Jonsson. 'Improved Tangent Estimate in the NEB Method for Finding Minimum Energy Paths and Saddle Points'. J. Chem. Phys. 113, 9978 (2000).

[34] R. Schmidt, A. Schwarz, R. Wiesendanger. Ground State of Magnetic Dimers on Metal Surfaces. Phys. Rev. B 86 174402 (2012).

[35] A.T. Costa, R.B. Muniz, et al. Magnetism of an Fe monolayer on W(110). Phys. Rev. B 78 054439 (2008).