Cyrille Barreteau
email: cyrille.barreteau@cea.fr

Daniel Spanjaard
email: daniel.spanjaard@u-psud.fr

Marie-Catherine Desjonquères

An efficient magnetic tight-binding method for transition metals and alloys Un modèle de liaisons fortes magnétique pour les métaux de transition et leurs alliages

Keywords:

Liaisons fortes Magnétisme Modèle de Stoner Couplage spin-orbite Anisotropie magnéto-cristalline Hartree

r é s u m é

Nous présentons un modèle de liaisons fortes paramétré et auto-cohérent utilisant une base d'orbitales atomiques s, p, et d pour décrire les électrons de valence des métaux de transition. Les paramètres du modèle sont déterminés à partir d'un ajustement non linéaire sur des résultats de calculs ab initio d'éléments purs en volume. Notre procédure ne nécessite aucun paramètre ni ajustement supplémentaire pour l'étendre aux systèmes avec plusieurs atomes de natures chimiques différentes. Nous avons généralisé notre modèle aux matériaux présentant une polarisation de spin et orbitale à l'aide de termes de Stoner et de couplage spin-orbite. Nous traitons aussi bien le magnétisme colinéaire que non colinéaire ainsi que les spirales de spin. Enfin nous montrons comment prendre en compte l'interaction électron-électron intra-atomique dans l'approximation de Hartree-Fock. Cela introduit une dépendance orbitale des interactions qui peut s'avérer importante dans les systèmes de basse dimensionalité et pour décrire correctement l'anisotropie magnéto-

Introduction

Even though Density Functional Theory has become an extremely efficient method widely used in many areas of physics, chemistry, and material science, the tight-binding (TB) description of the electronic structure remains very popular, since it provides a physically transparent interpretation in terms of orbital hybridization and bond formation. In addition, its moderate computational cost permits to handle rather large and complex systems, and its straightforward implementation allows many generalizations and applications. In addition, in recent years, with the increasing interest in electronic transport and the explosion of studies in graphene nanostructures, there has been a renewal of interest for TB calculations.

Historically the TB method was introduced by Slater and Koster [START_REF] Slater | Simplified LCAO method for the periodic potential problem[END_REF]. It was originally thought as a semi-empirical model to describe the electronic structure of solids with a reasonably small number of parameters that can provide reliable semiquantitative results when these parameters are determined from a fit on first-principles calculations. Jacques Friedel in the 1960s was one of the pioneers in the application of TB to the physics of transition metals [START_REF] Friedel | Metallic alloys[END_REF][START_REF] Friedel | Sur la structure électronique et les propriétés magnétiques des métaux et alliages de transition[END_REF]. These models were essentially based on a physical description of the band structure, but no real arguments about the total energy were developed. Qualitative explanations of the trends in the variation of the total energy when some parameters are varied (number of d electrons, concentration) were proposed, but only based on the band contribution to the total energy. Jacques Friedel was particularly talented in developing simple models with a simplified schematic description of the electronic density of states (such as the rectangular band model [START_REF] Friedel | On the possible impact of quantum mechanics on physical metallurgy[END_REF]) that could nevertheless describe surprisingly well many physical properties of materials. Later on, physicists started to add a phenomenological repulsive pair-potential to the band energy [START_REF] Ducastelle | Modules élastiques des métaux de transition[END_REF]. It was however not very clear what was "hidden" behind this phenomenological term. Over the years, TB methods have acquired a more solid fundamental basis. In particular, with the work of Harris and Foulkes [START_REF] Harris | Simplified method for calculating the energy of weakly interacting fragments[END_REF][START_REF] Foulkes | Tight-binding models and density-functional theory[END_REF], it was shown how a tight-binding formalism can be derived from Density Functional Theory. Nowadays, there exist many electronic structure codes based on various versions of TB with different degrees of approximations [START_REF] Mehl | Applications of a tight-binding total-energy method for transition and noble metals: elastic constants, vacancies, and surfaces of monatomic metals[END_REF][START_REF] Sankey | Ab initio multicenter tight-binding model for molecular-dynamics simulations and other applications in covalent systems[END_REF][START_REF] Koskinen | Density-functional tight-binding for beginners[END_REF].

Very early, TB in combination with a Stoner-like model [START_REF] Stoner | Collective electron ferromagnetism[END_REF] was recognized as an adequate tool to describe magnetism in transition metals, and Jacques Friedel was indeed very active in this field [START_REF] Friedel | On the nature of the magnetic couplings in transitional metals[END_REF][START_REF] Friedel | On the role of d-d electron correlations in the cohesion and ferromagnetism of transition metals[END_REF]. Indeed, magnetism in a crystal is intimately related to its band structure. TB has been applied to a large variety of magnetic systems in various crystallographic structures, dimensionalities (from bulk to clusters), ordered alloys or presenting some kind of disorder [START_REF] Brouers | Van Der Rest, Local environment and magnetic properties in transitional binary alloys. I. (theory) -Abstract -IOPscience[END_REF]. It is not the goal of this paper to provide an exhaustive presentation of this wealth of research in magnetism. We will rather concentrate on the presentation of a TB model that we have developed over the years and that is able to describe accurately and efficiently a wide range of magnetic phenomena and materials. It is an empirical TB method with parameters fitted on ab-initio data. We will show how, with a limited number of simple and well controlled approximations, we have been able to generalize our model to alloys and include non-collinear magnetism, spin-spirals as well as spin-orbit coupling.

The paper is organized as follows: we will present the general concepts of the tight-binding description of electronic structure and its implementation in an s, p and d atomic orbital basis set for non-magnetic materials (Section 2). We will pay particular attention to describe properly features that are often not discussed thoroughly in publications: non-orthogonality of the TB basis set, self-consistent treatment, proper definition of local quantities, etc. Then, in Section 3, we will show how, using a simple Stoner model, spin-polarization can be included, first for collinear magnetism (Section 3.1), then for non-collinear configurations (Section 3.2). Spin-orbit coupling and magneto-crystalline anisotropy will also be discussed in detail. Section 3 will be ended by a discussion of more elaborated Hartree-Fock like Hamiltonians that can play an important role in low-dimensional or anisotropic systems. Finally we will draw conclusions in Section 4.

The tight-binding method

The tight-binding approximation is a kind of counterpart of the free-electron jellium model of solids. Indeed, while plane waves are the natural functions to describe the delocalization of electrons in a jellium, the Bloch functions are expanded over a set of atomic-like functions in the TB approach. Therefore, the TB model will in principle apply better to systems where the valence electrons originate from well-localized atomic functions that do not overlap too much between neighboring atoms. However, its domain of applicability has proved to be larger than expected. For a review of the different aspects of the tight-binding method, the reader can refer to the following books or review articles [START_REF] Desjonquères | Concepts in Surface Physics[END_REF][START_REF] Ducastelle | Order and Phase Stability in Alloys[END_REF][START_REF] Paxton | An introduction to the tight binding approximation -implementation by diagonalisation[END_REF][START_REF] Pettifor | Bonding and Structure of Molecules and Solids[END_REF][START_REF] Sutton | The tight-binding bond model[END_REF][START_REF] Sutton | Electronic Structure of Materials[END_REF]. The specificity of our approach consists in the derivation of a systematic procedure to determine the parameters of our TB model, which is able to describe the electronic, magnetic and energetic properties of transition metals and their alloys with an accuracy close to first-principle calculations.

Generalities

The electron potential in a crystal

A good first approximation for the potential seen by an electron in a crystal (or assembly of atoms) is to write it as a sum of isolated neutral atomic-like potentials,

V (r) = n V at (r -R n ) (1)
where the summation runs over the lattice vectors R n of the crystal. V at (r) is an atomic-like potential of radial symmetry.

Note that this decomposition is not unique and moreover the potential is not necessarily a truly atomic potential, but could be "modified" to take into account the confinement felt by the electrons in the crystal. In addition, since we are only interested in describing the valence electrons, this potential is rather a "pseudo" potential valid in a given range of energy for describing the interaction of valence electrons with the ions. Due to the lattice periodicity, the crystal potential is periodic. It can also be generalized to the case of several atoms per unit cell. We will then write:

V (r) = i,n V at i (r -R n -τ i ) (2)
The summation over i is running over the N at atoms in the unit cell. To avoid lengthy notations, we will usually write the potential in the following condensed manner by making use of an operator formalism V = i,n

V at i,n . The one-electron Hamiltonian operator is therefore written:

Ĥ = T + i,n V at i,n (3)
T being the kinetic energy operator.

Linear combination of atomic orbitals

Once the potential of the crystal has been written as a sum of atomic-like potentials, we consider the corresponding atomic wave functions φ at iλ , solutions of the Schrödinger equation for atom i:

( T + V at i )φ at iλ = ε at iλ φ at iλ ( 4
)
As mentioned previously, there is some flexibility in the way potential V at i is defined. In particular, for practical reasons, it is wished that the spatial extension of the associated atomic wave functions is of a smaller range than that of the true atomic wave functions. In addition, since the valence wave functions usually include more than one type of atomic orbital, the index λ (= 1 . . . N orb ) refers to the symmetry of the angular part of the real orbitals s, p

φ at iλ (r) = R iλ (r) λ(r) (5)
where R iλ (r) and λ(r) are respectively the radial and angular parts of the atomic orbital (r being the angular coordinates of the vector r with respect to some given x, y, z axes). In the following, for convenience, we will rather use the Dirac notations |inλ to denote the ket associated with the atomic orbital φ at

iλ (r -R n -τ i ) = r|inλ of atom i in the cell n.
The central approximation of the tight-binding model relies on the assumption that the restricted Hilbert space spanned by atomic-like orbitals is sufficient to describe the wave functions solutions of the Schrödinger equation (at least in a restricted energy range) in the crystal. The wave functions are therefore written as combinations of atomic-like orbitals:

| = inλ C inλ |inλ (6)
Note that even though the atomic orbitals form an approximate basis of the restricted "valence" space, this basis has no reason to be orthonormal. Indeed, the scalar product of two atomic wave functions located on different atomic sites defines the so-called overlap integrals:

S inλ, jmμ = inλ| jmμ [START_REF] Foulkes | Tight-binding models and density-functional theory[END_REF] For convenience, the overlap integrals have often been neglected to build up simplified models where it is assumed that the non-orthogonality can be taken into account through a renormalization of the other tight-binding parameters. These models can be very efficient and quite precise when one is interested more in the energetics rather than in the electronic structure, since energetics is mainly governed by d electrons, for which orthogonality is a good approximation. However, this approximation is rather crude if one wants a good description of the band structure, which needs the introduction of s and p electrons in the basis set since they are delocalized.

One-, two-, and three-center integrals

Let us now write the various matrix elements of the Hamiltonian in this atomic basis set. Making use of Eq. ( 4), the diagonal or on-site terms can be written:

ε inλ = inλ| Ĥ|inμ = ε at iλ δ λ,μ + inλ| jm =in V at jm |inμ (8)
The second term of pure electrostatic character is the so-called crystal-field integral α iλμ . The additional contribution α iλ obtained for λ = μ only shifts the value of the atomic level. Off diagonal term will be ignored in the following.

In the same spirit, the off-diagonal matrix elements of the Hamiltonian obtained for (in) = ( jm) can be split into two types of terms:

inλ| Ĥ| jmμ = ε at jμ S inλ, jmμ + inλ| V at in | jmμ + inλ| kp =in kp = jm V at kp | jmμ ( 9
)
where atom of index i = 0 is taken as the origin.

The first two terms [START_REF]The first two terms can also be written ε at jμ S inλ, jmμ + inλ|[END_REF] involve two-center integrals (the center of the potential is the same as one of the two wave functions) and are the only ones that will be retained since their contribution is much larger than the last three-center integral term. In the following, the off-diagonal elements of the Hamiltonian (involving two-center integrals only) will be denoted β inλ, jmμ . These are the so-called hopping integrals, which are crucial ingredients in the tight-binding model since they measure the ability of an electron to "jump" from one atom to the other and decay rapidly with the distance between neighboring sites. Therefore, a usual approximation consists in cutting-off the interaction above a certain threshold radius R c . This is a very important point since the rather short-ranged extension of the hopping integrals makes the TB Hamiltonian matrix sparse and allows the use of specific efficient algorithms. In addition, due to the spherical symmetry of the atomic potentials, the hopping integrals between two sites connected by a vector with an arbitrary orientation can be written as a linear combination of a set of hopping integrals β γ = ssσ , spσ , sdσ , ppσ , ppπ , pdσ , pdπ , ddσ , ddπ , ddδ defined in the case where the orientation is along the z-axis. This set of Slater-Koster integrals are the main ingredients in the TB model, and the β γ (R) are decaying functions with distance.

Bloch states and Hamiltonian matrix in the reciprocal space

In a periodic crystal, the Hamiltonian can be diagonalized using the Bloch states written in the form:

| iλ (k) = 1 √ N n e ik•(R n +τ i ) |inλ ( 10
)
where k is a vector lying in the reciprocal space and N the number of primitive unit cells in the crystal (where we use the Born-von Karman periodic boundary conditions). The eigenfunctions | (α) (k) of the TB Hamiltonian are written as a combination of Bloch states:

| (α) (k) = iλ C (α) iλ (k)| iλ (k) (11)
The

C (α ) inλ and C (α)
iλ (k) expansion coefficients are thus related by the expression:

C (α ) inλ = 1 √ N C (α) iλ (k) e ik•(R n +τ i ) ( 12
)
Each eigenstate α is now labelled by a band index α (running from 1 to N at N orb ) and a wave vector k. The Hamiltonian matrix of size (N at N orb × N at N orb ) can be written as:

H iλ, jμ (k) = iλ (k)|H| jμ (k) = m e ik•(R m +τ j -τ i ) i0λ|H| jmμ ( 13
)
A similar expression is derived for the overlap matrix

S iλ, jμ (k) = iλ (k)| jμ (k) = m e ik•(R m +τ j -τ i ) i0λ| jmμ (14)

Solving the TB Schrödinger equation

In a periodic crystal, writing the Schrödinger equation with a wave function given by Eq. ( 11) leads to a matrix equation of the form:

H (k)C α (k) = ε α (k)S(k)C α (k) (15)
H (k) and S (k) are the k-dependent Hamiltonian and overlap matrices while C α (k) is the vector built from the coefficients C α iλ (k). Due to the presence of the definite and positive S (k) matrix, Eq. ( 15) is called a generalized eigenvalue problem that needs to be solved for each k vector in the irreducible Brillouin zone.

The eigenvalues ε α (k) form the so-called band-structure of the crystal when the k vector is spanning the Brillouin zone. This generalized eigenvalue equation simplifies into a standard one when the overlap integrals are neglected. The presence of the overlap matrix slightly complicates the formalism, especially when one needs to define local quantities as will be seen in the following sections.

The local quantities in the TB method 2.2.1. Quantum mechanics in a non-orthogonal basis set

When using a non-orthogonal basis in quantum mechanics, many formulae have to be modified. This is the case of the so-called closure relation. For the sake of simplicity (and generality), let us drop the i λ and n indexes and consider a non-orthogonal basis set |a , and its corresponding overlap matrix S a,b = a|b . It is very convenient to introduce the dual basis |ã defined as follows [START_REF] Ballentine | Recursion, non-orthogonal basis vectors, and the computation of electronic properties[END_REF][START_REF] Artacho | Nonorthogonal basis sets in quantum mechanics: representations and second quantization[END_REF]:

|ã = b (S -1 ) a,b |b (16)
By definition, ã|b = a| b = δ a,b , and it is possible to define a generalized closure relation:

a |ã a| = a |a ã| = Îd ( 17
)
where Îd is the identity operator. It is then straightforward to prove that the trace of an operator  takes the form

Tr( Â) = a ã| Â|a = a a| Â|ã (18)
A very useful operator in quantum mechanics is the density operator, which is defined in a condensed way as ρ = f ( Ĥ) where f is the Fermi function. ρ takes a more transparent expression in the eigenfunction basis set |α of the Hamiltonian.

ρ = α |α f (ε α ) α| (19) f α = f (ε α )
being the occupation factor of the eigenstate of eigenvalue ε α . In the case of a periodic system, parameter α denotes the combination of the discrete index α previously introduced and of the continuous k-vector. Therefore using the previous relation of the trace, it comes out that the trace of the density operator that defines the total number of electrons in the system can be written as: N e = Tr( ρ) = Tr(ρ S ) [START_REF] Sutton | Electronic Structure of Materials[END_REF] The elements of the density matrix ρ are given by:

ρ a,b = α C α a C α b f α = c + b c a ( 21
)
The C α a = ã|α are the expansion coefficients of the eigenfunction |α in the basis set |a and c + b , c a the creation and annihilation operators in the atomic orbitals |b and |a , respectively. Finally, a natural way of defining a charge projected on orbital |a is to take the real part of the matrix element N a = Re(ρ S ) a,a , which can also be written in a more symmetric way as 1 2 ( ã| ρ|a + a| ρ|ã ). This compact expression is in fact nothing else than the usual Mulliken charge of a given orbital a:

N a = Re α C α a C α a f α , C α b = b S a,b C α b ( 22
)
Obviously, with this definition, the summation over all the orbitals |a allows us to recover the total charge, i.e., N e = a n a . If the overlap integrals are ignored in the expression of the local charge given by Eq. ( 22) (i.e. by replacing C α a by C α a ), then the conservation of the charge is no longer preserved and we call the corresponding quantity the "net" charge. These relations can readily be transcribed to define the local charge of a given orbital symmetry N iλ in the context of the TB formalism.

Density of states

The density of states is a fundamental quantity that provides a crucial information about the level distribution of the Hamiltonian spectrum. The total density of states is defined as:

n(E) = α δ(E -ε α ) = Tr(δ(E -Ĥ)) (23)
It simply "counts" the number of states between E and E + dE divided by dE. By definition, an integration up to the Fermi level gives the total number of electrons in the system. In a similar way, as done for the Mulliken charge, a local density of states projected on a given orbital can be defined:

n a (E) = Re α C α a C α a δ(E -E α ) (24)
The local charge is then recovered by an integration weighted by the Fermi function:

N a = n a (E) f (E) dE (25)

Band energy

As will be seen in the following, the band energy that can be obtained, at zero temperature, by summing all eigenvalues below the Fermi level is an essential component of the total energy of the system. At non-zero temperature, it is easily generalized and can be formulated in various ways:

E band = Tr( ρ Ĥ) = α ε α f α = En(E) f (E) dE (26)
In a similar manner, as for the charge or the projected DOS, one can define a local contribution to the band energy:

E band a = En a (E) f (E) dE (27)
In some cases, it is more convenient to work with the grand-potential , which is defined at zero temperature as [START_REF] Liechtenstein | Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys[END_REF] which can be generalized at finite temperature T [START_REF] Ducastelle | Order and Phase Stability in Alloys[END_REF]. The local contribution from site a is obviously obtained in the same way as the local band energy.

E band -N e E F : = E F (E -E F )n(E) dE

An spd tight-binding Hamiltonian for transition metals and alloys

In the following, we will present the TB model that we have developed over the years and applied to many different systems. The non-magnetic part of this tight-binding Hamiltonian is similar to that of Mehl and Papaconstantopoulos. We use a minimal non orthogonal basis set containing s, p, and d orbitals centered at each site i, |i, λ

(λ = s, p x , p y , p z , d xy , d yz , d zx , d x 2 -y 2 , d z 2 )
. Let us first consider the onsite matrix elements in the case of a system with a single chemical element. In a similar way to Mehl and Papaconstantopoulos [START_REF] Mehl | Applications of a tight-binding total-energy method for transition and noble metals: elastic constants, vacancies, and surfaces of monatomic metals[END_REF], we assume that at each interatomic distance the reference energy is chosen in such a manner that the total energy E tot is obtained by summing up the occupied one electron energy levels: [START_REF] Li | Magnetocrystalline anisotropy energy of Fe(001) and Fe(110) slabs and nanoclusters: a detailed local analysis within a tight-binding model[END_REF] Note that in this condition, the total energy can be calculated without introducing an empirical repulsive potential as currently done in pure d band calculations. Indeed, the repulsive contribution to the total energy is accounted for by the onsite terms, which, consequently, must depend on the local environment (number of neighbors and bond lengths). Following Ref. [START_REF] Mehl | Applications of a tight-binding total-energy method for transition and noble metals: elastic constants, vacancies, and surfaces of monatomic metals[END_REF], we write:

E tot = α ε α f α
iλ = a λ + b λ ρ 1/3 i + c λ ρ 2/3 i + d λ ρ 4/3 i + e λ ρ 2 i ( 30
)
with:

ρ i = j =i exp[-2 R ij ]F c (R ij ) (31)
where F c (R ij ) is a cut-off function and is a parameter. Let us consider now the hopping and overlap integrals. In order to reproduce more closely their variation with distance on a large scale, we have found useful to use a law somewhat more complicated than a simple exponential, as done in [START_REF] Mehl | Applications of a tight-binding total-energy method for transition and noble metals: elastic constants, vacancies, and surfaces of monatomic metals[END_REF]:

β γ (R) = (p γ + f γ R + g γ R 2 ) exp[-h 2 γ R]F c (R) ( 32
)
where γ indicates the type of interaction (e.g., ssσ , spσ , etc.) and F c (R) is a cut-off function:

F c (R) = 1 1 + exp[(R -R 0 )/l] (33)
In practice, this function fixes the number of neighbors that are taken into account. In the calculations we took R 0 equal to 14 Bohr and l equal to 0.5 Bohr, the interactions are strictly cut-off for distances above R c = R 0 + 5l = 16.5 Bohr. Let us recall that in this parameterization, the band energy is equal to the total energy. Limiting ourselves to s, p, and d orbitals, we have introduced 16 parameters in the expression of the on-site terms and four for each expression of the ten hopping and ten overlap integrals. Thus, to completely specify our model, we need to determine at most 96 parameters. These parameters can be obtained by a fit to ab-initio calculations of the non-magnetic band structure and total energy for several distances and simple crystallographic structures.

If we consider now a metallic alloy made of two chemical elements A and B, the following procedure has been carried out with success. A fit for both chemical elements is performed separately, but with the same value of . Then the intra-atomic terms of a given atom in the system will only depend on the nature of the considered atom by the coefficients

a λ , b λ , c λ , d λ and e λ .

The values of the hopping and overlap integrals β

A-A γ (β B-B γ
) between two identical A (B) atoms are the same as in the pure element, while for a mixed couple of neighbors A -B, β A-B γ is obtained by using the common approximation that consists in taking the arithmetic average of the corresponding homonuclear quantities:

β A-B γ = 1 2 (β A-A γ + β B-B γ ) ( 34
)
and therefore does not necessitate additional parameters. Finally, we apply a local charge neutrality procedure such that the local Mulliken charge of a given atom i is equal to the number of valence electrons of the specie occupying site i. Its practical implementation will be explained in the next section.

The advantage of our approach is that it does not require any further fitting to ab-initio data for the binary systems. This procedure works very well for the electronic and magnetic properties of alloys; however, if one is interested in the energetics of alloys (mixing energies etc.), it is sometimes necessary to slightly rescale the value of the hetero-nuclear hopping and overlap integrals with respect to the arithmetic average. This is the case for the FeCr system, which necessitates a small increase in the Fe-Cr hopping and overlap parameters to reproduce accurately quantities such as the mixing or interface energies [START_REF] Soulairol | A magnetic tight binding model for FeCr[END_REF].

Self-consistent corrections 2.4.1. Local charge neutrality and the notion of penalization

When dealing with inhomogeneous systems using this type of semi-empirical TB Hamiltonian, one is usually faced with the problem of charge transfer. Indeed, if such a model is applied for example to the case of a surface, or a cluster where atoms have very different geometrical environments, unphysically large charge transfers are obtained. In metallic systems in which strong screening effects prevent the accumulation of charges, it is important to avoid such effects. A straightforward manner to solve this problem is to use a so-called penalization technique that consists in adding a positive term to the total energy, which is zero if the charge neutrality is fulfilled, and very large otherwise. Let us start from the expression of the total energy of the system written in terms of the expansion (real) coefficients C α inλ on the atomic basis, which for simplicity is assumed to be orthogonal (the generalization to the case of a non-orthogonal basis set is presented at the end of this section). The total energy of a system described by an Hamiltonian H 0 before penalization is written as:

E tot,0 = α occ inλ, jmμ C α inλ H 0 inλ, jmμ C α jmμ ( 35
)
A minimization of this function with respect to C α inλ under the normalization constraint

inλ (C α inλ ) 2 = 1 leads to the usual eigenvalue Schrödinger equation H 0 C α = ε 0 α C α .
Let us now add a penalization term of quadratic form to the total energy:

E pen = in U i 2 (N in -N 0 i ) 2 (36)
where U i is a large positive quantity and N 0 i the charge that one wants to impose on site i. Minimizing the total energy E tot,0 to which the penalization term is added leads to a similar eigenvalue equation where the diagonal matrix elements of the Hamiltonian have been modified:

H inλ, jmμ = H 0 inλ, jmμ + U i (N in -N 0 i )δ inλ, jmμ (37)
In a periodic crystal, N in = N i , ∀n, and the application of the Bloch theorem leads to the corresponding modified Hamiltonian in k space:

H iλ, jμ (k) = H 0 iλ, jμ (k) + U i (N i -N 0 i )δ iλ, jμ (38)
This equation can be generalized to the case of a non-orthogonal basis set. Provided that the charge is defined in the Mulliken fashion, the modified Hamiltonian reads:

H iλ, jμ (k) = H 0 iλ, jμ (k) + 1 2 (U i (N i -N 0 i ) + U j (N j -N 0 j ))S iλ, jμ (k) (39)

Self-consistency in the tight-binding method

From Eq. ( 37) it is clear that this equation must be solved self-consistently, since the Hamiltonian now depends on the local charges. A straightforward procedure would be to start from an initial "input" charge distribution, and once the corresponding Hamiltonian is diagonalized, the "output" charges are used as new inputs and the process is iterated until the input and output charges differ by less than a given threshold number. However, it is well known that such a method generally fails to converge. The simplest way to get rid of this problem is to introduce only a portion of the output charge by a linear mixing strategy. Unfortunately, this is very uneffective and some more refined mixing techniques are necessary, such as Broyden mixing [START_REF] Johnson | Modified Broyden's method for accelerating convergence in self-consistent calculations[END_REF], based on a kind of Newton method that greatly improves the convergence. Nevertheless, convergence in highly anisotropic systems can still be challenging, especially when magnetism is introduced.

Double-counting corrections

Once the self-consistency has been achieved, total energy can be expressed in terms of the C α iλ (k) coefficients; however, it is usually more practical to use the band energy of the modified Hamiltonian, which can be written as:

E band = αkocc ε α (k) = αkocc C α iλ (k)H 0 iλ, jμ (k)C α jμ (k) + i U i (N i -N 0 i )N i ( 40
)
Therefore it comes out that the total (penalized) energy is not solely given by the band energy, but should be corrected by a so-called double counting term. Then total energy reads:

E tot = αkocc ε α (k) - i U i 2 (N 2 i -(N 0 i ) 2 ) ( 41
)
This expression also holds for a non-orthogonal TB model as long as the charges are defined in Mulliken's manner.

Force Theorem

In quantum mechanics, the absolute energies of a system are not meaningful, since they depend on the energy reference. Therefore, most of the time, what is really needed is the difference of total energies between cases (or systems) that involve changes in the interacting potential. When the modification of the potential is large, a self-consistent calculation is necessary to compare the total energies of two systems. However, there are many situations where these changes are relatively small, and in these cases a very useful approximation, commonly called the Force Theorem (FT), can be applied. FT is also known as Andersen force theorem, since it was first introduced by Andersen [START_REF] Mackintosh | Electrons at the Fermi Surface[END_REF] (see also a detailed discussion in Ref. [START_REF] Heine | Solid State Phys[END_REF] by V. Heine). A magnetic version of this theorem is presented in Ref. [START_REF] Liechtenstein | Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys[END_REF] where it is applied to calculate effective exchange interactions in metallic ferromagnets. We will now briefly recall its principles, but a more detailed derivation can be found in Ref. [START_REF] Li | Magnetocrystalline anisotropy energy of Fe(001) and Fe(110) slabs and nanoclusters: a detailed local analysis within a tight-binding model[END_REF] in the context of a tight-binding scheme. Let us consider a perturbative external potential δV ext , which in a tight-binding scheme is often a perturbation of the on-site levels of the Hamiltonian, but could be any other small perturbation. Due to self-consistent effects, the external potential slightly modifies the local charges, which consequently produces an additional induced potential δV ind . However, it can be easily shown that the variation of total energy brought by this induced perturbative potential is exactly compensated (to first order) by its corresponding double counting term. Therefore, the change in total energy is equal to the change of band energy induced by the external potential only: [START_REF] Lashell | Spin splitting of an Au(111) surface state band observed with angle resolved photoelectron spectroscopy[END_REF] calculated by ignoring the self-consistent corrections. This means that the eigenvalues of the perturbed system are obtained after a single diagonalization of the "perturbed" TB Hamiltonian (including the external potential only and not the induced potential). The Force Theorem, besides saving much computational time, is also often more precise than the brute force approach, which consists in performing full self-consistent calculations, therefore requiring an extreme precision.

E tot ≈ δ αocc f α ε α
Let us also add that the first-order variation of the free energy F at fixed number of electrons is equal the variation of the grand-potential at fixed chemical potential [START_REF] Ducastelle | Order and Phase Stability in Alloys[END_REF]. Therefore the Force Theorem can be equivalently applied to the grand-potential.

Magnetism in the tight-binding method

Collinear magnetism

Until now, the spin of the electron has been ignored and the two spin orbitals were degenerated since no spin-dependent potential was present to remove the degeneracy between "up" and "down" spins. One of the simplest, yet amazingly effi-cient, way to introduce magnetism in a TB Hamiltonian is based on a model initially proposed by Stoner [START_REF] Stoner | Collective electron ferromagnetism[END_REF], which will be described below. The wave functions as well as the pseudo-atomic orbitals forming the TB basis |inλσ now have two spin components, and σ = ±1 for up and down spins, respectively. The spin magnetization expressed in unit of Bohr magnetons is then defined as the difference between the number of majority (up) and minority (down) electrons, M = N ↑ -N ↓ .

Stoner model

To remove the degeneracy between the two spins, it is necessary to introduce a potential that generates a band splitting between "up" and "down" spins. In the TB model, a natural way of creating this spin splitting is to consider a local (on-site) potential that shifts down (up) the majority (minority) spin bands. Moreover, since this potential intrinsicallly originates from electron-electron interactions, exchange splitting should depend itself on the magnetization of the system and vice versa. The spin-dependent Stoner potential V Stoner σ σ takes the simple form in an orthogonal TB basis with only one type of orbitals (d for the transition metals):

V Stoner σ σ = - 1 2 inλ |inλσ (σ I i M i δ σ σ ) inλσ | ( 43
)
where I i are the Stoner parameters, and M i the spin magnetization of atom i (i.e. the difference between the number of electrons with up and down spins summed over all the orbitals of same character). The Stoner parameter is usually expressed in energy units, while the magnetization is in Bohr magnetons. I i only depends on the chemical nature of the atom sitting at site i. This local term produces a shift exch i = I i M i , called exchange splitting.

An alternative (but equivalent) manner to introduce the Stoner model is to start from the total energy of the nonmagnetic system E tot,0 to which a negative component is added, -

1 4 i I i M 2 i .
Then following the same procedure as for charge penalization, it comes out that the original non-magnetic Hamiltonian should be corrected by a local potential that has exactly the form of the Stoner potential. Then, taking into account the double counting terms, the total energy of the magnetic system can be written:

E tot = ασ ,occ ε α,σ + i I i 4 M 2 i ( 44
)
where ε α,σ are the eigenvalues of the spin-dependent Hamiltonian.

Stoner criterion

Let us consider the case of a system where all atoms are equivalent (typically a crystal with one atom per unit cell).

The eigenvalues of the Hamiltonian can be split into up and down spin energies, and the density of states of spin σ , n σ (E)

is simply obtained by a rigid shift of the density of states per spin of the non-magnetic system n 0 (E) (without exchange interaction):

n σ (E) = n 0 (E + σ 2 I M) ( 45
)
From this unique property of the density of states, a criterion for the existence (or not) of a magnetic solution can be derived. Indeed, the equation defining the Fermi level E F from the total number of electrons N and the one defining the total spin magnetization M are sufficient to completely define the system (i.e. determine the magnetization) as long as the original non-magnetic density of states is known:

N = E F n 0 (E + I 2 M) dE + E F n 0 (E - I 2 M) dE (46a) M = E F n 0 (E + I 2 M) dE - E F n 0 (E - I 2 M) dE (46b)
Eq. (46a) defines the Fermi level provided that the total number of electrons is fixed. It uniquely defines E F (M). Substituting E F (M) for E F in Eq. (46b) leads to an equation of the type M = F (M) where:

F (M) = E F (M) (n 0 (E + I 2 M) -n 0 (E - I 2 M)) dE (47)
F (M) is an odd function that necessarily saturates to a maximum value. A non-zero solution for M exists only when the slope of F (M) at M = 0 is larger than one that is equivalent to F (0) ≥ 1. It is then easily verified that this is equivalent to the inequality:

In 0 (E F ) ≥ 1 (48)
This is the famous Stoner criterion [START_REF] Stoner | Collective electron ferromagnetism[END_REF], which can be derived in many different manners. For example, it can also be obtained by minimizing total energy.

Table 1

Stoner parameter in eV for various elements obtained from a procedure of fine tuning of the m(a alat ) curve for a given crystallographic structure as explained in the text. In the case of Fe we used the bcc structure, and for Co, Ni and Pt the fcc structure. For Pt the equilibrium structure is non-magnetic but this element undergoes a ferromagnetic transition at a given expanded lattice. In the case of Cr we used the bcc antiferromagnetic phase and fitted the atomic magnetization projected on one the atom of the lattice.

Element

H = H TB + V LCN + V Stoner (49)
where H TB is the non-magnetic (diagonal in spin) Hamiltonian matrix described in Section 2.3, V LCN is the local charge neutrality correction whose matrix elements are given by:

V LCN inλσ , jmμσ = 1 2 (U i (N i -N 0 i ) + U j (N j -N 0 j ))S inλ, jmμ δ σ σ ( 50
)
This contribution does not depend on the spin either since the local charges N i , N j include a summation over both spin components. V Stoner is the Stoner potential:

V Stoner inλσ , jmμσ = - σ 2 (I i,λ M i,d )δ inλ, jmμ δ σ σ (51)
Here we have considered Stoner parameters that depend on the atomic nature of the atom centered at site i, but also on the orbital character (s, p, or d). M i,d is the spin magnetization of atom i summed over the orbitals of d character only. The spin polarizations of s and p electrons have been neglected since they are very small (Eq. ( 2) of [START_REF] Autès | Magnetism of iron: from the bulk to the monatomic wire[END_REF]). Let us also mention that for reasons that will be explained in Section 3.5 presenting the Hartree-Fock interaction in the TB model, the local magnetic moment is the "net" magnetic moment (i.e., without including the overlap terms). In practice, since the overlaps of d orbitals are rather small, this does not significantly modify the results compared to a case where the Mulliken charges are used.

At this point, it is important to note that the Stoner potential is diagonal in spin space, i.e. there are no coupling terms involving both spins. In practice, one has to diagonalize separately two Hamiltonians, which leads to two sets of eigenvalues corresponding to the "up" ε α,↑ (k) and "down" ε α,↓ (k) spins. This procedure should obviously be performed self-consistently since, as in the case of the local charge neutrality, the potential depends on the local magnetizations. Therefore similar mixing techniques to those explained above are used to reach convergence.

Determination of the Stoner parameter

In transition metals, magnetization is very dominantly borne by d orbitals (which explains the form of our TB potential where only M i,d is taken into account) and the Stoner parameter of d orbitals I d mainly determines the magnetization of the material and we have set I s = I p = I d /10. As seen from the Stoner criterion, the density of states at the Fermi level plays a crucial role in the onset of magnetism in a material. In addition, since the width of the DOS gets narrower when the interatomic distance increases, this necessarily produces higher DOS, so that normalization is preserved. Consequently, a uniform bond stretching is almost always accompanied by an increase of the magnetization. Even non-magnetic materials will become magnetic for a large-enough lattice spacing, since most isolated atoms are magnetic. This argument can also be applied to the understanding of the general trend of the variation of magnetization with the average coordination of a system. Indeed, similarly to the argument put forward to explain the effect of bond stretching, the average DOS width is decreasing when the average coordination decreases. This basically explains why low-coordinated systems tend to have larger magnetizations.

The onset of magnetism with lattice stretching is usually rather abrupt and this explains why we have chosen to tune the Stoner parameter by trying to reproduce this variation as precisely as possible. In practice, the magnetic moment M is computed as a function of the lattice parameter a alat of a given crystallographic structure for several values of I d . These curves are compared with ab-initio data and this very simple procedure allows a precise determination of I d within a rather small energy range. A few values are presented in Table 1. For the case of Fe, we give two values 0.88 and 0.95 eV, since we found that, even though I d = 0.88 eV leads to spin moments in better agreement with ab-initio data the phase stability is better reproduced with I d = 0.95 eV as discussed in the next section.

Total energy and magnetic phase stability

Magnetism being extremely sensitive to the local environment of atoms, in particular the bond-lengths and the coordination number, it is not surprising that it plays a crucial role on the phase stability of magnetic materials. As an illustrative example let us mention the case of small magnetic clusters as well as bulk iron since they are particularly instructive.

Cluster physics is very rich and complex due to the large number of possible geometrical atomic configurations (that grows exponentially with the number of atoms in the cluster). This complexity increases when magnetism is entering into play since, besides the many metastable geometrical configurations to be considered, there also often exists several magnetic solutions [START_REF] Pastor | Size and structural dependence of the magnetic properties of small 3d-transition-metal clusters[END_REF][START_REF] Barreteau | spd tight-binding model of magnetism in transition metals: application to Rh and Pd clusters and slabs[END_REF]. In addition the bond-lengths of the relaxed structures are usually depending on the magnetic solution. As a general trend it is found that solutions with higher magnetic moments favor larger interatomic spacings, showing once again that magnetism and structure are intimately entwined.

The intrication between magnetism and atomic structure is also very important in bulk materials. One of the most emblematic example is the case of iron. Indeed, if Fe were not magnetic its crystallographic structure should have been hexagonal closed packed as found by standard non spin-polarized Density Functional calculations. Magnetism strongly stabilizes the body centered cubic structure in iron and favors it with respect to more compact structures. This is illustrated in Fig. 1 where we present the result of our TB total energy calculations for bulk Fe as a function of the atomic volume. Note that for these total energy calculations we have taken I d = 0.95 eV, since this value gives a better description of the phase stability while the magnetization is slightly overestimated. The reason for this small inconsistency is due to the fact that the energy difference between fcc and bcc is overestimated (compared to ab-initio results) in the non-magnetic phase. As a consequence, the stabilization of bcc requires a larger Stoner parameter to gain more magnetic energy.

Non-collinear magnetism

In the previous section, we have explained how magnetism can be introduced in our TB scheme and we used the terms of "up" and "down" spins, which implicitly supposed that the spin quantization axis was taken along the magnetization direction that will be denoted as z in the following. However the z axis has no reason to be related to any crystallographic direction of the crystal, since spin and orbital variables are belonging to different spaces. It will be seen how the spin-orbit coupling potential is connecting these two spaces. Before discussing this point, we will first express the magnetic TB Hamiltonian in fixed spin coordinate axes, which is more appropriate to deal with non-collinear magnetic configurations, since in such systems the spin quantization axes are changing from one site to the other.

When studying non-collinear magnetism, it is necessary to use a spin-orbital function with two-components

| ↑ (k) | ↓ (k)
where both components are usually non-zero, in contrast with the collinear case where the wave functions can be separated into "up" and "down" spin wave functions corresponding to the two "independent" spin channels. The operators now act simultaneously on both components of the spin-orbitals, such as the Pauli matrices:

σx = 0 1 1 0 ; σy = 0 -i i 0 ; σz = 1 0 0 -1 3.2.

Local and global spin coordinate axes

We will first introduce global xyz and local x y z spin axes. By definition, the local axes are chosen such that z is pointing in the direction of the local magnetization. Starting from the global xyz axes of the crystal, a new frame x y z is obtained by a rotation of an angle φ around z, finally a further rotation of an angle θ around y gives the local spin frame x y z .

The spin one-half rotation matrix [START_REF] Kübler | Theory of Itinerant Electron Magnetism[END_REF] permits to change coordinate from local to global (and vice-versa) spin frame:

U (θ, φ) = e -i φ 2 cos θ 2 -e -i φ 2 sin θ 2 e i φ 2 sin θ 2 e i φ 2 cos θ 2 = e -i φ 2 σz e -i θ 2 σy (52)
The components of the eigenvectors | ↑ loc in the global frame are given by the following matrix equation U (θ, φ) 1 0 so that:

where | ↑ loc and | ↓ loc are the spinstates that diagonalize the Pauli matrix σ z (local frame), while | ↑ and | ↓ are the spin eigenstates diagonalizing σ z (global frame).

The Stoner model in the global axes

By construction the Stoner potential matrix is diagonal in spin when expressed in the local spin axes and takes the form:

V Stoner, loc inλ, jmμ = - 1 2 I iλ (M i,d σz )δ inλ, jmμ (55)
It is straightforward to show, using the transformation V Stoner, glob = U (θ, φ) V Stoner, loc U -1 (θ, φ) that the Stoner potential expressed in the global axes is given by: V Stoner, glob inλ, jmμ = - [START_REF] Anisimov | Band theory and Mott insulators: Hubbard U instead of Stoner I[END_REF] where M i,d = M i,d (sin θ i cos φ i , sin θ i sin φ i , cos θ i ) is the magnetization vector in x, y, z coordinates, and σ = ( σx , σy , σz ) is the vector operator built from the Pauli matrices. This is a rather intuitive expression whose simple form comes from the fact that the Pauli matrix vector also transforms like a regular three-dimensional space vector under a space rotation matrix.

1 2 I iλ (M i,d . σ )δ inλ, jmμ
Note that the Stoner potential is the only one that is spin dependent, since H TB and V LCN are purely diagonal in spin space and independent of the spin, while V Stoner is diagonal in real space and non-diagonal in spin-space. The density matrix can as well no longer be split into "up" and "down" spin contributions. Its local components are now written:

ρ iλ = ρ ↑↑ iλ ρ ↑↓ iλ ρ ↓↑ iλ ρ ↓↓ iλ ( 57
)
where

ρ σ σ iλ = Re α C α iλσ C α iλσ ( 58
)
The local charges and magnetic moments are then expressed as an appropriate trace over the product of density and Pauli matrices:

N iλ = Tr( ρ iλ ) = ρ ↑↑ iλ + ρ ↓↓ iλ (59a) M x iλ = Tr( ρ iλ σx ) = ρ ↑↓ iλ + ρ ↓↑ iλ (59b) M y iλ = Tr( ρ iλ σy ) = i(ρ ↑↓ iλ -ρ ↓↑ iλ ) (59c) M z iλ = Tr( ρ iλ σz ) = ρ ↑↑ iλ -ρ ↓↓ iλ ( 59d
)
Finally an important point has to be mentioned before closing this section: in the absence of spin-orbit coupling, even for a non-collinear configuration, a global rotation of the spin magnetization does not affect the total energy of the system. The system is said to be invariant by rotation.

An example of non-collinear magnetic configuration

Non-collinear configurations usually show up in two cases: in very inhomogeneous systems, or in frustrated systems. The first situation typically occurs in small clusters. Even in strongly ferromagnetic materials such as iron, one can observe departures from collinearity which however remain small [START_REF] Hobbs | Fully unconstrained noncollinear magnetism within the projector augmented-wave method[END_REF]. The second situation generally leads to more obvious noncollinear configurations. To illustrate this situation, let us consider a very simple system where strong antagonist interactions are present: the triangular lattice of chromium.

Chromium is a material with a clear antiferromagnetic (AF) first-neighbor interactions. In bulk bcc chromium, the magnetic groundstate is antiferromagnetic [START_REF]Cr is a so-called spin-density wave (SDW) well established experimentally but calculations based on DFT[END_REF]. When Cr atoms are located on a triangular lattice, a perfect AF order cannot be established, since it is not possible that all first-neighbor pairs have magnetic moments pointing into opposite directions in a triangular network. We have carried out a series of magnetic TB calculations on the Cr triangular lattice for lattice parameters ranging from 2.1 Å to 3 Å. The magnetic groundstate is actually non collinear, forming a well-known Néel structure as the one shown in Fig. 2 resulting from our non-collinear TB calculation. A row-wise collinear anti-ferromagnetic order (see Fig. 1b of reference [START_REF] Hobbs | Fully unconstrained non-collinear magnetism in triangular Cr and Mn monolayers and overlayers on Cu(111) substrates[END_REF]) does exist in this system, but at a higher energy. In addition, above a given lattice spacing, a ferromagnetic order can also develop, but its energy is even higher than the AF one. This is in perfect agreement with ab-initio results [START_REF] Hobbs | Fully unconstrained non-collinear magnetism in triangular Cr and Mn monolayers and overlayers on Cu(111) substrates[END_REF].

Generalized Bloch theorem and spin spirals

Some systems have another type of non-collinear magnetic groundstates called spin spirals. In these magnetic structures, the magnetization is rotated by a constant angle from a unit cell to the next one, this angle being defined by some spin spiral vector q. Since there is no spin-orbit coupling, a global rotation of the spin can always bring the spin spiral structure to the case where q is the axis of rotation. The spiral is then fully characterized by its axis direction, the azimuthal angle that evolves linearly along the rotation axis and can be written at a given atom φ = φ 0 + q•R, and the constant polar angle θ = θ 0 . For convenience, the direction of the q vector is taken as the z axis and the spin spiral looks like the one presented in Fig. 3 for a monatomic linear chain with a single atom per unit cell. The magnetization vector can be written: M (R) = M(cos(φ 0 + q•R) sin θ 0 , sin(φ 0 + q•R) sin θ 0 , cos θ 0 ). The z axis has no reason to be related to any specific crystallographic axis, but to avoid lengthy notations, we will keep the same notation xyz to denote the crystallographic axes and the spin spiral axes. The following demonstration will be done in the global spin framework and we can check that the final result does not depend on the chosen axes. A brute force strategy to describe this type of systems would consist in considering a large super-cell containing an entire period of the spiral. However, this would be very time consuming and inefficient since small q vectors would necessitate extremely long super-cells. In addition, with this approach the spiral should be commensurate with the lattice. Fortunately it is possible to use the so-called generalized Bloch theorem [START_REF] Sandratskii | Noncollinear magnetism in itinerant-electron systems: theory and applications[END_REF].Considering that a lattice translation of R accompanied with a rotation q•R leaves the system invariant, a generalized Bloch wave function must verify the relation:

U (0, q•R) ↑ k,q (r + R) ↑ k,q (r + R) = e ik•R ↑ k,q (r) ↓ k,q (r) (60)
so that it can be written as

| σ (k) = n e ik•R n U -1 (0, q•R n )|nλσ (61)
Fig. 3. An example of a spin spiral in a monatomic wire of lattice spacing a. The z axis is taken along the wire. In that particular geometry the "space" and spin axes coincide. The rotation angle φ along the wire is equal to φ 0 + nq z a in that case.

According to Eq. ( 52), the rotation matrix U (0, φ) is extremely simple, since it is diagonal when it involves only a φ angle.

U (0, q•R) = e -i q•R 2 0 0 e i q•R 2 = e -i q•R 2 σz (62)
Writing the Hamiltonian in this Bloch-wave basis leads to a TB matrix:

H λσ ,μσ (k) = m e i(k+ σ 2 q)•R n H 0λσ ,mμσ (θ 0 , φ 0 ) ( 63
)
where H 0λσ ,mμσ (θ 0 , φ 0 ) is a sum of three terms H TB , V LCN and V Stoner . Since H TB and V LCN are diagonal in spinspace and V Stoner consist only of (spin-dependent) onsite elements, they are invariant under the U (0, q•R) transformation. H 0λσ ,mμσ (θ 0 , φ 0 ) are matrix elements of an Hamiltonian identical to the one of a truly periodic system with a fixed spin magnetization of orientation given by the Euler angles (θ 0 , φ 0 ). For a spin-spiral of wave vector q, the k vector is simply replaced by a k + 1 2 σ q for the spin orbital of component σ . It should be noted that the result obviously does not depend on the spin axis introduced for the demonstration. For the case of several atoms per unit cell we can introduce a set of θ i 0 angles for each atom of the unit cell and the generalization of Eq. ( 63) is straightforward. Summing up over the k vector gives the total energy (per unit cell) of a spin-spiral defined by any vector q in the irreducible Brillouin zone. [START_REF] Campo | Extended DFT+U+V method with on-site and inter-site electronic interactions[END_REF] This total energy depends on the q vector and also on the polar angle θ 0 . A minimum of the energy curve at a q vector different from zero will demonstrate the existence of a spin spiral solution. In addition, the curve often presents accidents since by varying q over the whole Brillouin zone we explore magnetic configurations ranging from the ferromagnetic case q = 0 to cases where the q vector is at the border of the Brillouin zone. It can happen that some magnetic situations are not stable and either do not "converge" or lead to non-magnetic solutions [START_REF] Soulairol | Structure and magnetism of bulk Fe and Cr: from plane waves to LCAO methods[END_REF]. Note that in the case θ 0 = 0, there is no spiral and the system is simply ferromagnetic. This can be seen from Eqs. ( 63) and ( 64), since when θ 0 = 0 the Hamiltonian matrix H 0λσ ,mμσ (θ 0 , φ 0 ) becomes diagonal relative to spin indices and one recovers the situation where the spin space can be split into up and down sub-spaces. Indeed, each independent summation over k + σ spans the whole irreducible Brillouin zone becomes independent of the q vector. Let us stress again that the generalized Bloch theorem is only applicable in the absence of SOC. In case of moderate spin-orbit coupling, it is still possible to apply this theorem and then treat SOC as a perturbation of the system as in Ref. [START_REF] Kashid | Dzyaloshinskii-Moriya interaction and chiral magnetism in 3d-5d zigzag chains: tight-binding model and ab initio calculations[END_REF].

E tot (q, θ 0 ) = α,kocc ε α (k, q) -E dc
Let us illustrate the concept of spin spiral on a very simple yet instructive system made of a monatomic Fe wire of lattice spacing a (see Fig. 3). We have calculated the total energy of spin spirals for q wave vectors along -X ( = (0, 0, 0) X = Fig. 4. Total energy of spin spirals as a function of the spiral vector q in a monatomic Fe wire at various lattice spacings and θ 0 = π /2. π a (0, 0, 1)) for θ 0 = π/2 at different lattice spacings. q = and q = X correspond to ferromagnetic and antiferromagnetic orders, respectively. The results of our calculations are shown in Fig. 4. Interestingly a clear minimum appears in the total energy curve at approximately q = X 2 for a = 2.1 Å. A shallow minimum still persists at a = 2.15 Å, but disappears at larger lattice spacings. This shows that a spin spiral corresponding to a π/2 rotation of the magnetization between adjacent atoms is the most stable structure. This structure is replaced by a simple ferromagnetic solution for lattice spacings above 2.15 Å.

Spin-orbit coupling

The spin-orbit coupling (SOC) interaction originates from relativistic effects. Starting from the Dirac equation and expanding the four-component wave functions in the small v/c limit allows one to recover the non-relativistic Hamiltonian to which several terms are added. One of these terms is the so-called spin-orbit interaction, which plays a very important role in magnetic systems since it couples the orbital and spin moments and breaks the global rotational invariance. The strength of the spin-orbit potential is the strongest close to the nucleus and it can be very well approximated by a spherical potential which takes the simple form in real space:

V SO = ξ(r) L h • Ŝ h ( 65
)
where L = r ∧ p = hl and Ŝ = h σ 2 are, respectively, the orbital and spin moment operators, and the function ξ(r) is expressed in terms of the spherical electrostatic potential V (r) ξ(r) = h2

Spin-orbit coupling interaction

The matrix elements of V SO in a crystal built from a sum of atomic-like potentials with spherical symmetry, written in the basis of atomic spin orbitals |iλσ takes the form:

iλσ | V SO | jμσ = ξ iλμ 1 2 λσ | l• σ | μσ δ i, j (67)
Only diagonal elements of the SOC potential at each site i are retained since ξ i (r) is localized near r = 0. In addition, since the angular moment operator does not couple orbitals of different natures (s, p or d) and is zero for s, the spin-orbit coefficients ξ iλ are determined by only two parameters, namely ξ p and ξ d :

ξ iλ = ∞ 0 R 2 iλ (r)r 2 ξ i (r) dr , λ= p or d ( 68
)
We recall that λ denotes the angular part of the atomic orbital. Since l acts only on the orbital and σ on the spin, the matrix elements of each component can be written as a product:

λσ | lη ση | μσ λ| | μ σ | |σ ( 69
)
where η = x, y, z are the of the operators in a (space spin) framework. The matrix elements of angular moment l η λμ = λ| lη | μ are better known in the spherical harmonics basis, but are easily derived in real spherical harmonics by simple linear algebra transformations, these formula can be found for p and d orbitals in Ref. [START_REF] Schena | Tight-binding treatment of complex magnetic structures in low-dimensional systems[END_REF]. Note that all the matrix elements of lη in the real harmonics are purely imaginary.

Determination of the spin-orbit coupling constant

The d component of the spin-orbit coupling constant ξ d can easily be determined by comparison with the band structure calculated from ab-initio codes including spin-orbit coupling effects. Since the SOC potential is very localized near r = 0, it is not sensitive to the atomic structure, therefore a band structure calculation for any simple crystallographic structure can give an excellent evaluation of the SOC parameter. We have checked on different elements and for many different crystal lattices, including low-dimensional systems, that the value extracted by this simple procedure is perfectly transferable. In Table 2, we give a few values of ξ d for various elements. In contrast, the p component cannot be obtained in the same way, since the p orbitals are much more dispersive and well above the Fermi level; therefore, the effect on the band structure is difficult to quantify. However from ab-initio calculations, it is possible to evaluate ξ p by a direct integration of Eq. ( 68). In practice, it is found that ξ p is about three to five times larger than ξ d . For example, in Fe it is found that ξ p = 0.18 eV [START_REF] Schena | Tight-binding treatment of complex magnetic structures in low-dimensional systems[END_REF]. However, despite this large value, the effect of ξ p is almost negligible on most physical phenomena.

An effect of SOC in non-magnetic systems: Rashba splitting of gold

The clearest evidence of spin-orbit coupling in atoms, molecules or extended solids is the removal of degeneracies of the one-electron levels. Let us illustrate the effect of SOC on two cases: a bulk non-magnetic cubic material and a surface.

The effect of SOC is to remove degeneracies of degenerate levels when matrix elements of V SO exist between the eigen- states of the corresponding levels. For example, in a non-magnetic cubic material the t 2g orbitals are sixfold degenerate at the point (this is the 25 irreducible group representation). When SOC is switched on the sixfold degenerate levels (including spin variable) split into fourfold and doubly degenerate levels [START_REF] Friedel | Etude du couplage spin-orbite dans les métaux de transition. Application au platine[END_REF]. In the limit of small ξ d , perturbation theory can be applied and it is easy to show that the splitting is equal to 3 2 ξ d . This result provides a simple way to determine the ξ d parameter at least in systems with sufficiently small SOC. This typically applies to 3d and 4d but not to 5d transitions metals.

In a non-magnetic system the time reversal symmetry imposes that ε α,↑ (k) = ε α,↓ (-k). In addition, if the crystal has an inversion symmetry, it comes out that ε α,↑ (k) = ε α,↑ (-k) and ε α,↓ (k) = ε α,↓ (-k). Therefore, in a system with both time reversal and inversion symmetries, the band structure must satisfy ε α,↑ (k) = ε α,↓ (k) and there is no possible spin splitting.

At the surface of a crystal, the inversion symmetry is broken and then spin splitting is allowed; however, when the crystal is non-magnetic, this splitting should not produce any spin unbalance. Such a splitting occurs in particular at the (111) surface of face centered cubic crystals where the inversion symmetry in a direction perpendicular to the surface is broken and therefore ε α,σ (k ) = ε α,σ (-k ) (where k is the component of the wave vector parallel to the surface). This leads to the so-called Rashba effect. A very obvious manifestation of this effect can be evidenced on the splitting of the Shockley surface state that appears at the (111) surface of noble metals. This is particularly significant at the (111) surface of gold for which the surface state exhibits a doublet as evidenced in the pioneering Angular-Resolved Photoelectron Spectroscopy (ARPES) measurements by Lashell et al. [START_REF] Lashell | Spin splitting of an Au(111) surface state band observed with angle resolved photoelectron spectroscopy[END_REF]. In Fig. 5, we present the dispersion of the Shockley state obtained from a TB calculation including SOC with ξ d = 0.65 eV and ξ p = 1.5 eV. This curve is in surprisingly good quantitative agreement with the one obtained from ab-initio calculations for example in Ref. [START_REF] Mazzarello | Spin-orbit modifications and splittings of deep surface states on clean Au(111)[END_REF]. Interestingly, even though the Shockley state has a strong p z character (10 times larger than the one of d z character), the Rashba splitting is essentially governed by the value of ξ d . More precisely, as pointed out in Ref. [START_REF] Lee | Role of d orbitals in the Rashba-type spin splitting for noble-metal surfaces[END_REF], ξ d and ξ p play in opposite directions. If ξ p is set to zero, the Rashba splitting obtained is slightly larger than when both ξ p and ξ d are taken into account. However, the correction due to the p component of the SOC is much smaller than the one originating from the d orbitals. The Shockley state can be described by an effective two-dimensional Rashba Hamiltonian with the following dispersion curve [START_REF] Petersen | A simple tight-binding model of spin-orbit splitting of sp-derived surface states[END_REF][START_REF] Mazzarello | Spin-orbit modifications and splittings of deep surface states on clean Au(111)[END_REF]:

E ± = E 0 + h2 2m k 2 ± γ SO k ( 70
)
where k is the modulus of the component of the wave vector in the surface plane, m is an effective mass, γ SO an effective SO parameter that quantifies the Rashba splitting and E 0 is the energy at which the two curves cross at the ¯ point. Our results are well fitted by m /m e = 0.24, γ SO = 4.2•10 -9 eV•cm and E 0 = -0.33 eV. The dispersion relation and the Rashba splitting agree extremely well with ab-initio calculations [START_REF] Mazzarello | Spin-orbit modifications and splittings of deep surface states on clean Au(111)[END_REF]. The position of the Shockley state is found slightly above the experimental value (E 0 = -0.42 eV [START_REF] Lashell | Spin splitting of an Au(111) surface state band observed with angle resolved photoelectron spectroscopy[END_REF] or E 0 = -0.49 eV [START_REF] Reinert | Direct measurements of the L-gap surface states on the (111) face of noble metals by photoelectron spectroscopy[END_REF]), but in the energy range of the ab-initio results.

SOC effect in magnetic systems: orbital moment

In a magnetic system, the magnetization has two origins: the spin and the orbital moments. The total magnetization being given by the expectation value of -μ B ( l + σ ). In a non-orthogonal atomic basis, we can define a Mulliken-like component of the orbital magnetization on a given site i:

L i = Re αk f (ε α (k)) λμσ C α iλσ (k)l λμ C α iμσ (k) ( 71
)
When SOC is neglected, the expansion coefficients can be taken as real, and since the matrix elements l λμ are imaginary, the orbital moment L i is necessarily zero. Therefore in magnetic systems the SOC is at the origin of the non-vanishing expectation-values of the orbital moment. In bulk the value of the orbital moment is usually small, since it is quenched by crystal-field effects. In low-dimensional systems or when the symmetry is broken, the orbital moment is less quenched and can reach larger values. Finally, let us point out that, except along high-symmetry directions, the orbital and spin moments have no reason to be aligned. This is indeed the case for a magnetization along an arbitrary direction of a crystal for which one can observe a small misalignment between spin and orbital magnetizations.

In magnetic systems, another important consequence of SOC is also that not only degeneracies are removed in the band structure, but they also depend on the orientation of the magnetization. Therefore, the band structure itself will be different depending on the magnetization vector. This is at the origin of the magnetocrystalline anisotropy discussed in the following section.

At this point let us note that this dependence on the orientation (θ, φ) of the magnetic moment does not show up in the expression of the SOC matrix elements itself when it is expressed in the global basis. The θ and φ dependence appears only in the Stoner potential. If we had adopted the local spin axes by construction, the magnetic potential would have been diagonal in the spin-space, but then the spin-orbit-coupling Hamiltonian ξ 2 ( lx σx + ly σy + lz σz ) takes a more complicated form (expressed in terms of the Euler angles θ and φ) since the operators ση should now be written in the local spin basis. Such an expression can be found in Refs. [START_REF] Friedel | Etude du couplage spin-orbite dans les métaux de transition. Application au platine[END_REF][START_REF] Autès | Magnetism of iron: from the bulk to the monatomic wire[END_REF]. This alternative point of view is more convenient if one wants to make use of a perturbative formula with respect to the SOC.

Magnetic anisotropy

Magnetic anisotropy denotes the dependence of the total energy of a magnetic system on the orientation of its average magnetization. The orientation that corresponds to the minimum of energy defines the so-called easy axis, the hard axis being the one of maximum energy. The magnetic stability of a system crucially depends on its magnetic energy distribution.

Magnetic anisotropy can formally be written as a sum of two terms with very different physical origins: the so-called shape anisotropy and the magneto-crystalline anisotropy. The shape anisotropy originates from classical dipole-dipole interactions and essentially depends on the global shape of the sample and not very much on the local atomic environment. It usually favors magnetization along elongated directions of the sample. For example, in films, it favors in-plane anisotropy. The magneto-crystalline anisotropy energy (MAE) is a purely quantum effect intrinsically related to the spin-orbit coupling and extremely dependent on the dimensionality and symmetry of the system, since it is sensitive to tiny details of the electronic structure. As a general trend, low dimensionality and anisotropic geometrical configurations favor larger MAE, but it can change sign depending on small variations of the local density of states.

The magnetocrystalline anisotropy

The MAE is obtained as the energy difference between two magnetization directions. The numerical values of MAE are usually extremely small, typically ranging from a few μeV per atom in the bulk to fractions of meV at a surface. A brute force method that consists in performing two self-consistent calculations to extract a very small energy difference is often not the best strategy, since it is very computationally demanding and furthermore often not more precise than the Force Theorem. In contrast, the Force Theorem is particularly well suited to this kind of calculation. Within this approximation, the MAE is obtained by the difference of band energy between two magnetic orientations, but ignoring self-consistent effects. The calculation proceeds in three steps: i) a first self-consistent collinear TB calculation (without SOC), ii) a spin rotation of the density matrix to align the magnetization along a given direction, and iii) a single-step (i.e. non-self-consistent) non-collinear TB calculation including SOC. Note that if the SOC were not included, the energy would not depend on the magnetization direction and the MAE would be zero. The MAE corresponding to the difference of energy between two spin orientations of unit vectors u 1 and u 2 can be written:

MAE ≈ E 1 F En 1 (E) dE - E 2 F En 2 (E) dE (72)
where n 1 (E) and n 2 (E) are the densities of states and E 1 F , E 2 F the Fermi levels of the spin orientations u 1 and u 2 , respectively.

The local picture

Starting from Eq. ( 72), it is tempting to decompose the MAE into a sum over the atomic sites by introducing the projected density of states on each site n i (E) such that n(E) = i n i (E). However, the atom-resolved MAE obtained from such a procedure is inappropriate, since it suffers from long-range (Friedel-like) charge oscillations (due to the neglect of selfconsistency in the Force Theorem approach) that affects the value of the local energy. However, due to the conservation of the total charge, these oscillating terms cancel out when summed up over the whole system. A more appropriate way to define the local MAE is to work within the grand-canonical ensemble, which allows us to get rid of these oscillations that produce spurious contributions to the local MAE [START_REF] Li | Magnetocrystalline anisotropy energy of Fe(001) and Fe(110) slabs and nanoclusters: a detailed local analysis within a tight-binding model[END_REF]. The atom resolved MAE is then written:

MAE i ≈ E F (E -E F ) n i (E) dE (73)
where

n i (E) = n 1 i (E) -n 2 i (E)
is the difference of DOS at site i between the two magnetic orientations, and E F is the Fermi level obtained from the self-consistent collinear TB calculation without SOC. A straightforward generalization of this formula can be applied to define an atom and orbital-resolved MAE. This approach based on the grand-canonical ensemble has been used in many different contexts, since it is also particularly well suited to a perturbative development [START_REF] Ducastelle | Order and Phase Stability in Alloys[END_REF].

We have applied our method to analyze the MAE of iron and cobalt slabs and nanocrystals [START_REF] Li | Magnetocrystalline anisotropy energy of Fe(001) and Fe(110) slabs and nanoclusters: a detailed local analysis within a tight-binding model[END_REF][START_REF] Li | Out-versus in-plane magnetic anisotropy of free Fe and Co nanocrystals: tight-binding and first-principles studies[END_REF] or the magnetic properties of FePt clusters [START_REF] Barreteau | Magnetic and electronic properties of bulk and clusters of FePt L1(0)[END_REF]. However, to illustrate the strength of this local picture analysis, let us consider the interesting case of an ultrathin film of cobalt deposited on a gold (111) surface. It is known that above five layers, cobalt grows in a hcp stacking with its bulk lattice parameter. Although the first neighbor distance in Au (d Au ) is about 15% larger than in Co (d Co ), the system adopts a configuration that allows an almost perfect strain relaxation [START_REF] Belkhou | Structure and growth mode of epitaxial Co/Au(111) magnetic thin films[END_REF].

If we note that d Au /d Co ≈ 8/7, an almost perfectly commensurate structure can been obtained since every eight Co atom falls in perfect registry with every seven Au atom along a dense (first-neighbor) atomic row. Consequently, we have built up a plausible interface structure from an 8 × 8 super-cell of Co hcp(0001) in contact with an 7 × 7 super-cell of Au fcc(111).

We have adopted a finite slab geometry, i.e., with a limited number of atomic layers. In practice, we took five layers of Au in fcc stacking in contact with five layers of Co in hcp stacking, the separation between the two slabs is equal to the average inter-layer distance between Co and Au. The total slab therefore possesses two free surfaces (one of Au and one of Co) and an interface (see inset of Fig. 6). Our aim is to compare this system with a slab made of 10 layers of Co to evaluate the influence of the gold interface.

The MAE is defined as the difference between energies of in-plane and out-of-plane magnetic configurations MAE = E -E ⊥ and, for the sake of simplicity, we have chosen the most symmetric in plane orientation. By definition, a positive MAE means that out-of-plane configurations are energetically favorable. We have calculated the atom-resolved MAE of the Au/Co system and compared it with the 10 layers of pure cobalt. The result of our calculations is presented in Fig. 6. It shows a very erratic behavior for the Au/Co system in the vicinity of the interface. The MAE of atoms belonging to the same atomic plane can change by more than 0.3 meV. Moreover, from the comparison with the pure cobalt system, it appears that the contact with gold strongly favors out-of-plane anisotropy of the cobalt film and furthermore, even though gold presents a very small magnetization at its interface layer, it still contributes significantly to the total MAE (favoring out-of-plane anisotropy). Let us add that if the SOC parameter of gold is set to zero (not shown on the graph), then the MAE of the gold interface layer vanishes.

Perturbation theory and Bruno Formula

The application of the Force Theorem to calculate the MAE is based on the assumption that the variation of the SOC potential is small and that any self-consistent correction can be ignored in the energy difference. We can go a step further by carrying out a perturbation expansion with respect to the SOC potential itself. Second-order perturbation theory predicts that the MAE is a quadratic function of the direction cosines (sin θ cos φ, sin θ sin φ, cos θ ) of the spin quantization axis with respect to the crystal axis. The coefficient of the expansion are written as complicated integrals over products of density matrices whose physical interpretation is not easy [START_REF] Autès | Magnetism of iron: from the bulk to the monatomic wire[END_REF]. At this point, let us note that it is more convenient to work in the local spin axes to derive this formula.

Starting from this second-order perturbation expansion, Bruno [START_REF] Bruno | Tight-binding approach to the orbital magnetic moment and magnetocrystalline anisotropy of transition-metal monolayers[END_REF] went a step further. By neglecting spin-flip terms (provided that the exchange splitting is large enough), he derived the following linear relation between the MAE and the variation of the orbital magnetization:

MAE i = - ξ i 4 L 1 i •u 1 -L 2 i •u 2 (74)
where

L 1 i •u 1 (L 2 i •u 2 )
is the orbital moment of site i projected along the spin magnetization direction u 1 (u 2 ). This formula first derived for fcc monolayers [START_REF] Bruno | Tight-binding approach to the orbital magnetic moment and magnetocrystalline anisotropy of transition-metal monolayers[END_REF] has been generalized by Cinal et al. [START_REF] Cinal | Magnetocrystalline anisotropy in ferromagnetic films[END_REF] to periodic systems with several atoms per unit cell. It shows that the direction corresponding to the minimum of energy (easy axis) is also the one with the largest orbital magnetization. This is almost always true and Eq. (74) gives the right trend. However, from a quantitative point of view, the linear relation between the MAE and the orbital moment variation is not perfectly fulfilled in many cases. Therefore it is safer to use the local definition of the MAE that we have presented.

Hartree-Fock in the tight-binding method: TB+U

In the TB model presented above, it is assumed that the on-site elements depend only on the spin population at each site and not on their repartition between the orbital states. However, it is well known that this charge redistribution is determined mainly by intra-atomic Coulomb interactions. Thus we will now introduce the corresponding potential, which can be written in the second quantization formalism, where c † iλ σ (c iλσ ) denotes the creation (annihilation) operators of an electron at site i in the atomic spin orbital λσ :

V int = 1 2 iλ 1 λ 2 λ 3 λ 4 σ σ U i λ 1 λ 2 λ 3 λ 4 c † iλ 1 σ c † iλ 2 σ c iλ 4 σ c iλ 3 σ ( 75
)
where:

U i λ 1 λ 2 λ 3 λ 4 = φ at iλ 1 (r), φ at iλ 2 (r )| e 2 |r -r | |φ at iλ 3 (r), φ at iλ 4 (r ) (76)
In the standard Hartree-Fock decoupling, this interaction becomes [START_REF] Desjonquères | Orbital contribution to the magnetic properties of iron as a function of dimensionality[END_REF]:

V HF int = iλ 1 λ 2 λ 3 λ 4 σ σ U i λ 1 λ 2 λ 3 λ 4 < c † iλ 1 σ c iλ 3 σ > c † iλ 2 σ c iλ 4 σ -< c † iλ 1 σ c iλ 4 σ > c † iλ 2 σ c iλ 3 σ ( 77
)
In the following we will limit ourselves to a monatomic system and the interactions are restricted to d orbitals since they are much less extended than the s and p ones. In this case, the matrix elements (Eq. ( 76)) that involve at most four different orbitals can be expressed as linear functions of three Racah parameters A, B, C , (or else three Slater integrals F 0 , F 2 , F 4 , [START_REF] Griffith | The Theory of Transition Metal Ions[END_REF]). If cubic harmonics (real combination of spherical harmonics) are used, the terms with three and four orbitals are all proportional to B and the non-vanishing terms with two orbitals are of two types:

Coulomb integrals U λ 1 ,λ 2 = U λ 1 λ 2 λ 1 λ 2 and exchange integrals J λ 1 λ 2 = U λ 1 λ 2 λ 2 λ 1 = U λ 1 λ 1 λ 2 λ 2 ,

and the following relationship is verified

J λ 1 λ 2 = (U λ 1 λ 1 -U λ 1 λ 2 )/2.
It is customary to use an other set of parameters that simplifies the expression of the Coulomb and exchange integrals, i.e., the average values of these integrals, which are independent of λ indices:

U = 1 4 λ 1 λ 2 λ 1 =λ2 U λ 1 λ 2 = A -B + C ; J = 1 4 λ 1 λ 2 λ1 =λ2 J λ 1 λ2 = 5 2 B + C (78)
Moreover, the parameters U , J , B are more physically transparent than A, B and C , since U is linked to the Coulomb integral, J to the exchange integral and B to the orbital dependence of the electronic interactions [START_REF] Autès | Transport électronique polarisé en spin dans les contacts atomiques de fer[END_REF]. Let us now comment on the choice of the basis set. Obviously, if all the matrix elements of U are included (U , J , B model), the results do not depend on this choice. This is also the case if we set B = 0 (U , J model), which in the real basis set amounts to neglect the terms involving three and four orbital since they are relatively small due to the angular dependence of the atomic orbitals that point in different directions of space. However, they should be kept to study orbital magnetism. When using a spherical harmonics basis set denoted by the value of the quantum number m, these terms are functions of B and C . Thus if they are neglected without changing the terms involving one and two atomic orbitals, the rotational invariance is destroyed unless we set B = C = 0, in which case the Coulomb integral U = A is completely isotropic and J vanishes.

Thus the corresponding Hamiltonian is probably oversimplified to study orbital magnetism, but could be sufficient for a description of spin magnetism.

If in the expression of the full Hartree-Fock Hamiltonian (U , J , B model), we make now the following approximations:

(i) the intra-atomic density matrix is assumed to be diagonal with respect to both orbital and spin indices, (ii) for each spin and at each site, the exact population of the λσ spin orbital is replaced by its average value over all orbitals of spin σ , i.e., n iλσ → (1/5) λ n iλσ = N idσ /5 then knowing that U λ 1 λ 2 λ 3 λ 4 vanishes when three indices are equal and that the elements involving three different orbitals obey the relations λ 2 U λ 2 λ 1 λ 2 λ 3 = 0 and λ 2 U λ 2 λ 1 λ 3 λ 2 = 0, it is easy to show that:

V iλ 1 σ ,iλ 1 σ = λ 2 U λ 2 λ 1 - J λ 2 λ 1 2 N i,d 5 -σ λ 2 U λ 2 λ 1 M i,d 10 (79)
where

N i,d = N id ↑ +N id ↓ (M i,d = N id ↑ -N id ↓
) are the net d population (spin moment) at site i and σ = +1(-1) for majority (minority) spin. Carrying out the summation over λ 2 leads to:

V iλσ ,λσ = 9U -2 J 10 N i,d - σ 2 U + 6 J 5 M i,d (80)
We recognize the Stoner-like potential, where I = (U + 6 J )/5 is the Stoner parameter and (9U -2 J )/10 is an effective Coulomb integral. To summarize starting from the full Hartree-Fock (U , J , B model) interaction terms that accounts for orbital and spin magnetisms, we have derived the (U , J ) model Hamiltonian that correctly describes spin magnetism. Then using further approximations, we recovered the Stoner model [START_REF] Desjonquères | Orbital contribution to the magnetic properties of iron as a function of dimensionality[END_REF]. Actually at the origin of DFT+U (see below), Anisimov and collaborators [START_REF] Anisimov | Band theory and Mott insulators: Hubbard U instead of Stoner I[END_REF] claim that the standard local spin density approximation (LSDA) is controlled by a kind of effective Stoner parameter and adding a Hubbard-type interaction allows one to solve some of the known weaknesses of DFT.

In all these models, a double counting term E dc =< V int > /2 must be subtracted from the sum of occupied one-electron energies in order not to count twice electron-electron interactions in the total energy and the determination of the average intraorbital and interorbital elements of the density matrix should be done self-consistently. Finally, we must not forget that some electronic interactions are already included in the tight-binding Hamiltonian since this Hamiltonian has been parameterized by fitting the results of LDA or GGA calculations. Following the treatment done in the "atomic-limit" theory, we must subtract from V int (ρ σ σ iλμ ) the quantity V int ( 1 5 N i,d δ λ 1 λ 2 δ σ σ ), ρ being the matrix density and N i,d the average charge on site i. Let us point out the similarities that hold between our TB+U theory in its formulation and the LSDA+U one, i.e. the density functional theory in the local spin density approximation [START_REF] Liechtenstein | Density-functional theory and strong interactions: orbital ordering in Mott-Hubbard insulators[END_REF] to which an intra-atomic Hubbard U has been added. In this approach, a set of localized orbitals can be identified. The electronic interaction between these orbitals are treated in a self-consistent Hartree-Fock (mean-field) manner. The modified one-electron Hamiltonian can therefore be written as:

H LSDA+U = H LSDA + H HF int -H dc (81)
where H LSDA is the usual local spin-density functional H HF int is the electronic interaction between the localized orbitals and H dc is a double counting term that accounts as far as possible for the electron-electron interaction already included in H LSDA . The expression of H HF int is the same as in TB+U (see Eq. ( 77)) and is treated as in the UJB or UJ model [START_REF] Shick | Implementation of the LDA+U method using the full-potential linearized augmented plane-wave basis[END_REF][START_REF]Note that most DFT+U implementations are expressed in the spherical[END_REF] in order to have a rotationally invariant formulation [60]. Several expressions of H dc can be found in the literature, but the closest to our TB+U approach is the so-called atomic limit expression [START_REF] Shick | Implementation of the LDA+U method using the full-potential linearized augmented plane-wave basis[END_REF] since H dc can be obtained by considering the energy of an isolated atomic shell with N = N ↑ + N ↓ electrons. If the exchange interaction is included, it is then written as:

E dc = U 2 N(N -1) - J 2 [N ↑ (N ↑ -1) + N ↓ (N ↓ -1)] (82)
Let us finally point out that the reference system in our TB theory is the non-magnetic system, since the determination of the parameters is done by a fit on the band structures (and total energies) on the non-magnetic system, whereas in LSDA+U the reference system, i.e. H LSDA , is the magnetic one.

Interatomic Coulomb interaction: TB+U+V

Let us now show how the tight-binding model can be generalized to take into account the influence of the interatomic Coulomb interactions. For the sake of simplicity, we first consider the case of a narrow s band all atoms being of the same chemical species. Moreover we assume that the set of atomic orbitals centered at each site i is orthogonal. Following Hirsch [START_REF] Hirsch | Metallic ferromagnetism in a single-band model[END_REF], we consider the following Hamiltonian in which the intersite interactions are limited to the first nearest neighbor: where n iσ = c † iσ c iσ , and -t is the hopping integral between nearest neighbors. The Coulomb interactions are described by the on-site term ∝ U , and the two-site terms: charge-charge interactions ∝ V , exchange interactions ∝ J inter , and the "pair hopping" term ∝ J inter .

Applying the Hartree-Fock decoupling to the two body term leads to:

t σ = t + (V -J inter )B σ -( J inter + J inter )B -σ (85) ε σ = Z (V -J inter )n + (U + Z J inter )n -σ (86)
E dc stands for the double counting energy terms, Z is the number of nearest neighbors, n σ the occupation number for electrons of spin σ , n = σ n σ is the total band filling and B σ = ρ iσ , jσ = c † jσ c iσ is the intersite density matrix between two first neighbors. It can be seen on Eq. (85) that the hopping integrals are renormalized, this renormalization depending not only on the spin polarization but also on the atomic lattice, as pointed out by Hirsch [START_REF] Hirsch | Metallic ferromagnetism in a single-band model[END_REF]. Consequently, the Stoner criterion is changed and becomes:

I eff n 0 (E 0 F ) > 1 (87)
with:

I eff = U + Z J inter + Z (V + J inter ) E 0 F Zt 2 ( 88
)
where n 0 (E 0 F ) is the density of states per spin in the PM state at the Fermi level E 0 F . This generalizes the criterion derived by Hirsch [START_REF] Barreteau | Effects of intersite Coulomb interactions on ferromagnetism: application to Fe, Co, and Ni[END_REF] for the particular case of a constant density of states. Thus the influence of the inter-site exchange integral J inter is to act in favor of the FM state for any band filling, since it increases I eff and decreases the bandwidth of the PM state (see Eqs. ( 85) and ( 88)). Let us now examine the effect of V and J inter . At low and high band fillings, the renormalization of the hopping integral in the PM state tends to zero since B 0 = B(n/2) vanishes. As a consequence, due to the term proportional to V + J inter in I eff , the PM state is more easily destabilized for low values of n since, in this case, the ratio (E 0 F /Zt) 2 is close to unity, the bottom of the band being at E = -Zt. This is also true when n approaches n = 2 for simple and body-centered cubic lattices and this tendency is weakened for the face-centered cubic lattice since (E 0 F /Zt) 2 = 1/9.

We have generalized this analysis to an s, p, d tight-binding model [START_REF] Barreteau | Effects of intersite Coulomb interactions on ferromagnetism: application to Fe, Co, and Ni[END_REF], but keeping only the most important Coulomb interactions since the intersite exchange integrals are much smaller. The calculations that we have done on the electronic structure (band structure, densities of states, magnetic moment) of ferromagnetic Fe, Co, and Ni are in excellent agreement with the calculations performed by the ab-initio method. In particular, the modifications of the hopping integrals lead to a width of the majority spin band smaller than that of the minority spin one (see for example Fig. 6 of reference [START_REF] Barreteau | Effects of intersite Coulomb interactions on ferromagnetism: application to Fe, Co, and Ni[END_REF]). Note that recently an extension of DFT+U (so-called DFT+U+V) has also been proposed to include inter-site electronic interactions [START_REF] Campo | Extended DFT+U+V method with on-site and inter-site electronic interactions[END_REF].

Conclusion

We have presented a parameterized magnetic s, p and d tight-binding model. Its initial form was introduced by Mehl and Papaconstantopoulos [START_REF] Mehl | Applications of a tight-binding total-energy method for transition and noble metals: elastic constants, vacancies, and surfaces of monatomic metals[END_REF] to describe non-magnetic and monatomic systems. We have successfully extended it to multi-component systems and introduced spin polarization by adding a Stoner-like term. Spin-orbit coupling as well as non-collinear magnetism have also been included into our model. In practice, our method only requires a careful fit on bulk monatomic band structures and total energies obtained from non-magnetic ab-initio calculations to describe the distance dependence of the usual hopping integrals, overlap integrals, and onsite terms. A single Stoner parameter and a spin-orbit coupling constant per chemical element are sufficient to entirely determine the magnetic properties of a material. It can also be further refined into a TB+U-like method by taking into account an intra-atomic Coulomb interaction with orbital-dependent Coulomb and exchange integrals, which can be important when dealing with low-dimensional or strongly asymetric systems. We finally briefly discuss inter-site Coulomb interactions leading to a so-called TB+U+V model.

Our model has been applied to a wide range of magnetic materials and atomic structures. It is able to reproduce accurately the magnetic phase stability of materials as well as the magneto-crystalline anisotropy of complex systems containing thousands of atoms, or non-collinear magnetic configurations in quantitative agreement with ab-initio results. Let us finally mention that our TB Hamiltonian has also been incorporated into an electronic transport code based on a Landauer approach within a standard Green function formalism [START_REF] Autès | Giant orbital moments are responsible for the anisotropic magnetoresistance of atomic contacts[END_REF][START_REF] Autès | Electronic transport in iron atomic contacts: from the infinite wire to realistic geometries[END_REF].

Fig. 1 .

 1 Fig.1. Total energy per atom as a function of the atomic volume for bulk Fe in the body centered cubic (bcc), face centered cubic (fcc), and hexagonal closed pack (hcp) structure. For each structure, we have considered the ferromagnetic (FM) and the non-magnetic (NM) configurations.

Fig. 2 .

 2 Fig. 2. Non collinear Néel magnetic structure of a triangular chromium lattice obtained from a TB calculation. Right: Total energy per atom of the non-magnetic (NM), ferromagnetic (FM), row-wise antiferromagnetic order (AF) and non-collinear Néel phase as a function of the first-neighbor distance of the free standing Cr triangular lattice.

Fig. 5 .

 5 Fig. 5. Rashba splitting of the Shockley state at the (111) surface of gold. The origin of energy is set to the Fermi level.

Fig. 6 .

 6 Fig. 6. Atomically resolved MAE (E i, -E i,⊥ ) of the system made of five atomic layers of Co(0001) in hcp stacking in contact with five layers of fcc stacking Au(111) (full red line), compared with a 10-layer hcp(0001) slab of pure cobalt (full black line). Each atomic layer of the Au/Co unit cell is made of 64 atoms in the case of Co and 49 in the case of Au, since we have considered a super cell built from a 8 × 8 cobalt unit cell in contact with a 7 × 7 Au unit cell which gives an almost perfect lattice coincidence. (The super cell is presented in the inset, periodic boundary conditions are applied in the plane of the surface.) The total number of atoms per unit cell is 565 for the Au/Co system and 640 for the pure Co system. Note that the pure Co system could have been modelled with a (1 × 1) surface unit cell but for the sake of comparison of the two systems we have preferred to use a much larger 8 × 8 slab with ten layers. We have only shown the gold interface layer and the 5 outermost layers of Cobalt. Layer 5LCo corresponding to the surface.

  =i,σ t σ c † iσ c jσ + iσ ε σ n iσ -E dc (84)with the following spin-dependent hopping integrals and orbital energies,

  x , p y , p z , d xy , d yz , d xz , d x 2 -y 2 , d z 2 , etc.,

  We have been able to include magnetism and local charge neutrality corrections within a non-orthogonal spd TB model and our Hamiltonian matrix takes the following expression:

		Cr	Fe	Co	Ni	Pt
	I d (eV)	0.82	0.88-0.95	1.10	1.05	0.6
	3.1.3. The magnetic spd TB model

Table 2

 2 SOC parameters of d orbitals in eV for various elements obtained from a comparison with ab-initio band structures.

	element	Cr	Fe	Co	Ni	Pt	Au
	ξ d (eV)	0.035	0.06	0.08	0.1	0.45	0.65

q vectors (for σ = ±1), where k

Acknowledgements

We wish to thank Jacques Villain and Hélène Bouchiat for having given us the opportunity to present our work in this dossier of Comptes rendus Physique dedicated to the memory of Jacques Friedel. We would also like to associate several colleagues and students who have participated at some point in the elaboration of the tight-binding model presented in this article: Ricardo-Guirado López for his contribution to the magnetism of transition metal clusters, Andrzej Oleś for his expertise in the Hartree-Fock Hamiltonian and Coulomb interaction, Gabriel Autès for the implementation of Green function formalism in the context of electronic transport using our magnetic TB model, Chu Chun Fu and Romain Soulairol for introducing us to the complex world of magnetism and defects in Fe-Cr alloys, Timo Schena for his contribution to the implementation of Dzyaloshinskii-Moriya interaction in a TB scheme, Dongzhe Li and Alexander Smogunov for their important contribution to the modelling of magneto-crystalline anisotropy in transition-metal nanostructures. Finally CB is particularly indebted to François Ducastelle for his numerous sound advices based on his deep understanding of physical phenomena.