{'fraction_non_alphanumeric': 0.0327306643707849, 'fraction_numerical': 0.00937006744152676, 'mean_word_length': 4.605775418275418, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exchange interactions and local moment fluctuation corrections in ferromagnets at finite temperatures based on non-collinear density-functional calculations

We explore the derivation of interatomic exchange interactions in ferromagnets within density-functional theory (DFT) and the mapping of DFT results onto a spin Hamiltonian. We delve into the problem of systems comprising atoms with strong spontaneous moments together with atoms with weak induced moments. All moments are considered as degrees of freedom, with the strong moments thermally fluctuating only in angle and the weak moments thermally fluctuating in angle and magnitude. We argue that a quadratic dependence of the energy on the weak local moments magnitude, which is a good approximation in many cases, allows for an elimination of the weak-moment degrees of freedom from the thermodynamic expressions in favor of a renormalization of the Heisenberg interactions among the strong moments. We show that the renormalization is valid at all temperatures accounting for the thermal fluctuations and resulting in temperature-independent renormalized interactions. These are shown to be the ones directly derived from total-energy DFT calculations by constraining the strong-moment directions, as is done e.g. in spin-spiral methods. We furthermore prove that within this framework the thermodynamics of the weak-moment subsystem, and in particular all correlation functions, can be derived as polynomials of the correlation functions of the strong-moment subsystem with coefficients that depend on the spin susceptibility and that can be calculated within DFT. These conclusions are rigorous under certain physical assumptions on the measure in the magnetic phase space. We implement the scheme in the full-potential linearized augmented plane wave method using the concept of spin-spiral states, considering applicable symmetry relations and the use of the magnetic force theorem. Our analytical results are corroborated by numerical calculations employing DFT and a Monte Carlo method.


I. INTRODUCTION
In recent years there has been increasing activity in the prediction of high-temperature magnetic properties of solids, especially regarding critical magnetic transition temperatures. The theoretical approach is founded in many cases on two assumptions: (i) that the magnetic excitations of a system can be phenomenologically described within a classical or quantum Heisenberg model, and (ii) that the exchange parameters entering the model (i.e., the excitation energies) can be microscopically derived from total energy results of e.g. density-functional theory calculations. This "magnetic multiscale modelling" has proven succesful in many cases, including itinerant elemental ferromagnets such as Fe and Co, [1][2][3][4] localized-moment systems such as Gd, 5,6 magnetic alloys as NiMnSb or Co 2 MnSi, 7 or diluted magnetic semiconductors as (Ga,Mn)As. 8,9 In these systems, the estimated Curie temperatures T C are within 10-15% of the experimental values, showing that the approach is reliable for practical purposes.
The derivation of magnetic excitation energies from density-functional calculations requires additional assumptions, since density-functional theory is, in principle, valid only for the description of the ground state of many-electron systems. Mainly, one relies on an adiabatic hypothesis, which conjectures that the slow motion of low-energy magnetic excitations can be decoupled from the fast motion of intersite electron hopping, so that the local electronic structure has time to relax under the constraint that a magnon traverses the system. Then the magnons are regarded as practically static or "frozen" objects, and constrained density functional theory 10 is employed. This adiabatic approximation, together with the realization that the local moments persist until above the Curie temperature, constitutes a paradigm which has proven fruitful. Quite a few theories and methods have been based on it, aimed at the calculation of thermodynamic quantities, including T C , by harvesting the model parameters from first principles calculations and working out the thermodynamics within Monte Carlo methods or other suitable approaches.
Among the first to discuss the adiabatic approximation in connection to DFT calculations were Gyorffy and coworkers in the development of a mean-field theory of magnetic fluctuations, 11 even though the concept was applied earlier on the level of solutions to many-body model Hamiltonians (see e.g. Small and Heine 12 ). Later Antropov et al. 13 and Halilov et al. 14 derived equations for adiabatic abinitio spin dynamics. Further elaboration came from Niu and collaborators 15 and Gebauer and Baroni, 16 who showed that the Berry curvature enters the adiabatic dynamics equations; they also demonstrated mathematically how the Born-Oppenheimer method can be generalized for adiabatic spin dynamics without the requirement that there should be a large mass difference defining the time-scale difference.
Based on the adiabatic approximation, a proper parametrization of the excitation spectrum is called for. In many, but not all, ferromagnetic systems the magnitude of the local atomic moments is relatively robust under modest rotations, so that mainly the pair exchange constants entering a Heisenberg model are required. From the point of view of methodology, two approaches are commonly used for the determination of the pair exchange constants. The first is based on multiple-scattering theory and Green function methods, frequently in the approximation of infinitesimal rotations. 17 It is widely used in methods where the Green function is available. 1,2,8,[18][19][20][21][22][23][24] A variant of this approach is the use of the disorder-local-moment (DLM) reference state 11 as an approximation to the magnetic disorder at the Curie temperature, from which one obtains exchange parameters 25,26 in a manner analogous to the method of infinitesimal rotations. 17 Higher-order interactions (e.g. the biquadratic, 4-th order, or anisotropic exchange) can be also obtained in a systematic way either from the ferromagnetic ground state 27,28 or from the DLM state. 23 The second approach, which we follow in this work, is based on reciprocal-space calculation of spin-wave excitation energies by constraining the system to spin spirals of given wavevectors q. This requires non-collinear calculations, restricted to the primitive chemical unit cell (in the absence of spin-orbit coupling) by virtue of a generalized Bloch theorem. 33 A subsequent integration in q-space (in form of a back-Fourier transformation) yields the pair exchange constants. Although the force theorem 34 is in principle not necessary here, practically it is often used in order to avoid the numerical load of a self-consistent calculation for each vector q. This second approach is well-suited for electronic structure methods which are based on Hamiltonian diagonalization rather than Green functions, and has been developed, for example, for the augmented spherical wave (ASW) method 36,37 or the LMTO method. 14 Concerning the second approach, we should also mention that a proper energy fit to a large number of constrained-moment-directions has been used to extract model parameters, inspired by alloy-theory methods, for example by applying the Connolly-Williams theory 29 to fit a number of antiferromagnetic states 30 or in a more refined way by applying a systematic spin-cluster expansion including higher-order interactions 31,32 for a fit to non-collinear states.
The results of the two approaches (infinitesimal rotations and spin-wave spectra) on the values of exchange interactions agree well, at least as long as the density-functional calculations are done within the same electronic structure method. 22 There is also good agreement in the spin-wave spectra of the two approaches if these are calculated in the first approach from a Fourier-transformation of the real space coupling parameters, as long as a sufficiently large number of atomic shells is taken in the Fourier series. 19 The present paper contains two major parts. In the first part, Sections II and III, we start by presenting the formalism for the calculation of pair exchange constants which we have implemented in the full-potential augmented linearized plane wave (FLAPW) method. 38,39 The approach of spin spirals and inverse Fourier transformations 14 is used to arrive at formulae for the pair exchange constants in the case of one or more magnetic atoms per unit cell. Symmetry relations are derived that reduce the selection of spin-spirals to the irreducible part of the Brillouin zone. Then, we focus on accuracy tests concerning the use of the force theorem at finite rotations; we also address the problem of subtracting the contribution of the magnetization in the interstitial. As a test, we calculate the Curie temperature of certain compounds by a Monte Carlo method. In the second part, Sec-tion IV, we discuss a way to parametrize the energetic contribution of longitudinal changes in the atomic spin moment so as to include them in a Monte Carlo method in the case that we are faced with a compound containing strong-moment and induced-moment sublattices. We obtain a scheme for the study of temperature-dependent magnetic properties and derive equations that allow a simplification of the computational method and a reduction of the computational cost under certain frequently-met physical conditions; importantly, we show the hitherto unnoticed result that the fluctuating weakmoment degrees of freedom can be analytically eliminated in favour of renormalized, temperature-independent strongmoment interactions, and that the thermodynamics of the full system (strong plus weak moments with bare interactions) can be derived from the thermodynamics of the strong-moment system only with renormalized interactions (a general derivation is given in the Appendix). A numerical demonstration of this rigorous result is presented in the case of the half-Heusler alloy NiMnSb. Finally, in Sec. V we place our work in a wider context comparing with the treatment of the weak moments or the concept of renormalized exchange interactions presented by other authors and we state the basic physical assumptions underlying our result. We then conclude with a summary.
Adopting the adiabatic approximation, the magnetic interactions are modelled by a classical Heisenberg Hamiltonian. The part of the total energy due to these interactions is then obtained from the expression where R mα(nβ) ≡ R m(n) + τ α(β) . Here, R m(n) are the lattice vectors and τ α(β) are the basis vectors specifying the positions of the atoms within the unit cell. M mα(nβ) are the atomic magnetic moments at the sites R mα(nβ) , while J(R mα , R nβ ) is the exchange coupling constant for the pair of atoms situated at these sites, and is the quantity to be calculated. The summations over the indices n, m are carried out over all lattice vectors, and the ones using indices α, β, over all the atoms in the unit cell. The factor 1/2 takes care of the double counting and the on-site term (R mα = R nβ ) is left out. The constants J(R mα , R nβ ) contain the information about the inter-site interaction due to the exchange coupling. The knowledge of these exchange interactions is essential for the description of thermal excitations in magnetic solids and their deriving from ab-initio calculations is the core problem in the attempt to describe the system with the Heisenberg Hamiltonian. The correspondence between the ab-initio theory and the Heisenberg model is established by using the ansatz which follows from Eq. (1), as a defining relation of J(R mα , R nβ ) within density-functional theory. Here, δM nα and δM nβ are to be understood as small differences with respect to the direction only, not the magnitude. I.e., an appropriate energy functional (usually within the local density or generalized gradient approximation (LDA or GGA)) E DFT [ρ, m] of charge density ρ(r) and magnetization density m(r) is used in the place of E in Eq. (2). When evaluating (2), it is assumed that the exchange-correlation field is collinear within each atom, so that the derivative with respect to atom-cell integrated moments M mα = cell mα m(r) d 3 r is meaningful. The intra-atomic collinearity is an approximation, justified by the energetic dominance of the moment formation (usually of the order of eV) compared to the formation of ferromagnetic order (usually of the order of less than 0.1 eV). From these comments it is also evident that we do not require that the local moments are quantized either in the form M = S(S + 1) µ B or M z = S µ B with S integer or half-integer, as, e.g., would be the case in ferromagnetic semiconductors. 40 Rather, Eq. (1) represents the lowest term in an expansion of the total energy in terms of the magnetization direction, neglecting longitudinal enhancement or suppression of the moments, and Eq. (1) represents a classical Heisenberg model, valid after local quantum spin fluctua-tions have been averaged out (see, e.g., the discussion in [14]). There are known cases when these approximations are insufficient, and we discuss such a case in Section IV. Next, the collective transverse magnetic excitations are approximated by static spin spirals, the energy of which is calculated within the non-collinear FLAPW method. 39,41 For the Fourier and back-Fourier transformations that are needed we follow the formalism of Halilov et al.. 14 In the case of a spin spiral with wave vector q, the azimuthal angle of the magnetic moment of an atom at the position R nα is given by ϕ nα = q · R nα . The magnetic moment of an atom at the position R nα is where θ α is the so-called cone angle, a relative angle between the final and initial direction of the local magnetic moment (here chosen along the z axis; this choice does not limit the generality in absence of spin-orbit coupling), and φ α an eventual phase factor, also called phase angle. Taking advantage of the translational invariance we define R ≡ R n − R m and Here, the energy E is a function of the spin-spiral vector q, as well as of the cone and phase angles of the magnetic moments on all the atoms of the unit cell. The dependence on these angles is collectively expressed by Θ for the set of all cone angles {θ α } and by Φ for the set of all phase angles {φ α } in the argument of E. To account for the condition τ α = τ β − R under which the sum in Eq. (4) is conducted (and is indicated by a prime), from now on we set J(τ α , τ α ) ≡ 0, for all the atoms α in the unit cell.
With the aim to obtain the exchange interaction constants J(τ α , τ β − R) at the minimum of computational expense, we define in the following a set of expressions which are evaluated computationally. We first define the Fourier transform It is straightforward to show that with the use of this Fourier transform, Eq. (4) becomes Finally, we note that from the definition (10) it is clear that J µν (q) satisfies the same symmetry relations as the coefficients J µν (q). The described calculations can be time consuming, since they involve the determination of small energy differences (typically of the order of a few mRyd). Due to the oscillatory phase in Eq. (12), the appropriate size of the q point set increases with the distance between the two atoms for which the interaction constant is being calculated, the grid fineness being basically determined by the inverse of the quantity |τ µν − R| that enters the exponential in Eq. (12). In the spirit of the Nyquist-Shannon sampling theorem, 42 and for |R| ≫ |τ µν |, the q-grid fine spacing in the direction of R should be at most half the value 2π/|R| . Additionally, sufficient accuracy requires larger plane-wave basis sets and a finer k-point grid compared to a simple ground-state calculation. A rule of thumb for increased accuracy is that, given a q-grid the k-grid should be twice as fine per spatial dimension of the lattice in order to avoid spurious oscillatory behaviour of period of the grid-spacing δk inJ µν (q). A self-consistent calculation of all energies needed here is computationally very demanding. Fortunately, in many cases the spin spiral can be considered a small enough perturbation that the force theorem 34,35 can be used to calculate the energy differences. We discuss this in the following subsection.

D. Test of the Applicability of the Force Theorem
The magnon energy is calculated as the difference between the total energy of the system with a spin spiral (this is an excited state), and the ground state, which is ferromagnetic in the systems under study here. A self-consistent calculation of the spin spiral requires use of a constraint, in the form of an external spiraling magnetic field, which forces the magnetization to take the form given in Eq. (3). On the other hand, an application of the force theorem requires a position dependent rotation of the exchange-correlation field B xc so that its direction acquires the form (3), i.e., is the self-consistent exchangecorrelation field of the collinear calculation at atom type α (V xc σα is the exchange-correlation potential dependent on spin σ =↑, ↓). In the FLAPW method, this rotation is applied in the muffin-tin spheres, i.e., non-overlapping spheres around the atomic nuclei where the potential is expanded in radial functions and spherical harmonics, as well as in the interstitial space between these spheres. Using the field of Eq. (13) in the Kohn-Sham equations yields a (non-self-consistent) sum of eigenvalues of the occupied levels; according to the force theorem, the magnon energy is approximated by the difference of this sum to the sum of eigenvalues of occupied levels in the self-consistent ferromagnetic ground state.  The approximation is expected to be better for smaller perturbations, i.e., for smaller cone angles and/or smaller magnon wavevectors |q|. As a test, in the left panel of Figure 1 the dispersion curve is shown for a spin spiral in NiMnSb, defined by a cone angle θ = 30 o , and a spin-wave vector along [001] direction. The force-theorem and self-consistent calculations agree rather well. The right panel of Fig. 1 shows the dependence of the magnon energy on the squared sinus of the cone angle (sin 2 θ) for a fixed spin-spiral vector q = (0, 0, 1)π/a lat . We see a a maximal deviation of the order of 7-8% for the unfavourable case of θ = 90 o and q = π/a lat , while the deviation starts becoming visible at cone angles larger than θ max ∼ 50 o .
A general conclusion is that if one wants to use the force theorem to obtain the spin-spiral dispersion within the whole first Brillouin zone of the crystal, a cone angle θ = 30 o seems to be a reasonable choice, since the energy differences are not too large for the magnon to stop being a small perturbation, but are also not too small so that one would have to employ a very large basis, or k-point set. This cone angle was used in the calculation of the Heisenberg interaction constants J presented in Section III.
We close this section with a comment on a spurious effect of the exchange-correlation field in the interstitial region when the force theorem is used. (The interstitial region is covered by empty spheres or empty cells in some ab-initio methods, e.g. in the ASW 36,37 or LMTO method, 14 if the structures are relatively open, as are for instance the half-Heusler or the zinc-blende structure.) In a force-theorem application, the trial exchange-correlation potential in the interstitial is a smooth periodic function. However, in a self-consistent spinspiral calculation, the resulting exchange-correlation potential is in many cases a much less smooth, although still periodic, function. 43 Then the interstitial magnetization in the forcetheorem calculation can cause a serious overestimation of the spin-spiral energy. This spurious energy contribution can be circumvented by setting the magnetic part of the exchangecorrelation potential in the interstitial to zero. Depending on the volume filling of the touching "muffin-tin" atomic spheres of the system, the spurious energy can be considerable, becoming larger for open systems. In Fig. 2a we show a calculated example for NiMnSb. The self-consistent spiral energy and the force-theorem energy (here normalized to sin 2 θ) are practically identical, if we set B xc = 0 in the interstitial in the force theorem calculation (in the self-consistent calculation, B xc is non-zero also in the interstitial); if this correction is not used, then the spiral energy is strongly overestimated. Figs. 2b-d show the same for Fe 2 MnSi in the full-Heusler structure, FeMnSi in the half-Heusler structure, and FeSi in the zincblende structure. All three structures are based on an fcc lattice and were calculated here with the same lattice parameter a lat = 5.663Å, but differ in the number of atoms in the unit cell and therefore in the volume of the interstitial region. Fe 2 MnSi has the smallest interstitial region, and evidently the interstitial magnetization has almost no effect; both force theorem calculations practically coincide with the selfconsistent result. For FeMnSi, however, the interstitial volume is larger and the spiral energy is strongly overestimated in the force theorem calculation, if the interstital magnetization is not set to zero. The effect is even stronger for FeSi, where the interstitial volume is largest.

III. EXCHANGE INTERACTION PARAMETERS AND CURIE TEMPERATURE
Following the prescription of Sec. II we calculated Heisenberg exchange interaction parameters (J ij ; i, j are abbrevia-tions of mα, nβ) of Co 2 MnSi and NiMnSb, shown in Fig. 3 and used a Monte Carlo method to obtain an estimate of the Curie temperatures for these compounds. The calculations for both compounds were performed within the generalized gradient approximation 44 on a 4096 k-point mesh with 2744 qpoints in the full Brillouin zone. The planewave cutoff was k max = 3.8 au −1 . The convergence was checked with respect to the above parameters.
In Co 2 MnSi, both Co and Mn atoms have strong and stable magnetic moments, whose interaction is described with J parameters depicted in Fig. 3a. In NiMnSb, though the magnetic moment of Ni is small and actually induced by the Mn surrounding atoms, both Ni and Mn were treated as magnetic atoms and the parameters of Mn-Mn, Ni-Ni, and Mn-Ni interaction were calculated (Fig. 3b). As will be discussed in Sec. IV, this treatment gives an insight into the thermal behaviour of the Ni sublattice and with the use of an extended Heisenberg model more useful information can be obtained.
Co 2 MnSi and NiMnSb are half-metallic ferromagnets, i.e., the density of states (DOS) in one spin direction (here majority spin) is metallic, while in the other spin direction the DOS shows a band gap around the Fermi level.
For a non-half-metallic ferromagnet, the interaction constants J ij as a function of distance follow a decaying oscillating behavior that in the simplest case is of the Ruderman-Kittel-Kasuya-Yoshida (RKKY) type, decaying with |R i − R j | −3 , but in general can have many superimposed periods and a different decay power law depending on the details of the Fermi surface. 2 On the other hand, for half-metallic ferromagnets the gap at one spin direction leads to an imaginary wave vector and an exponential decay with distance. 2 In both cases shown in Fig. 3, we notice this very fast decay of the interaction parameters with the distance between the atoms. For a comparison, in Fig. 4, the exchange interaction parameters of fcc Co (which is not half-metallic) are shown. The decay here is much slower.
In a short digression, we compare the results of the spinspiral approach within the FLAPW method to the approach of infinitesimal rotations 17 calculated with the full-potential Korringa-Kohn-Rostoker Green function method (KKR). 45 In KKR we use both the ferromagnetic state and the disordered local moment (DLM) state as reference points (see also the discussion in Sec. V B). We see in Fig. 4 that the agreement is good for fcc Co, with a deviation of 3% nearest-neighbor coupling, if the ferromagnetic state is used as a reference, while the deviation is large if the DLM state is used as reference (also the atomic moment decreases from 1.65 µ B in the ground state to 1.08 µ B in the DLM state). For Co 2 MnSi (not shown here) we find that the difference between KKR and FLAPW is larger, with the Mn-Co nearest-neighbor coupling being approximately 10% weaker in the KKR calculation, while if the DLM state is used as reference, the Co moment vanishes altogether.
Once the exchange interactions are known, the effective Heisenberg model can be solved for the Curie temperature. While the mean field approximation is computationally the simplest method to this task, it is known (and has been shown in practice) that the resulting T C is overestimated. The ran- Spiral energy E(q, θ) as function of sin 2 θ in NiMnSb for the cases of (i) a self-consistent calculation (including the self-consistent, spiraling polarization of the interstitial) (ii) a force-theorem calculation including the polarization of the interstitial, and (iii) a force-theorem calculation excluding the polarization of the interstitial. Evidently, case (ii) stongly overestimates the spiral energy, since the interstitial polarization acts as an effective magnetic field. Case (iii) on the other hand is a good approximation to the exact result (i). The spiral wavevector is q = (0, 0, 0.15) (π/a lat ). (b,c,d): Same as in (a) but for Fe2MnSi, FeMnSi, and FeSi, in full-Heusler, half-Heusler, and zinc-blende structures, respectively, to demonstrate the effect of increasing the volume of the interstitial. The structures correspond to an fcc lattice with four, three, and two atoms per unit cell, respectively. All three were calculated in the same lattice parameter (a lat = 5.663Å) and at a spiral vector of q = (0, 0, 0.25) (π/a lat ). dom phase approximation (also known as Tyablikov approximation), on the other hand, is rather accurate. 46 The most accurate, but numerically more expensive, way to calculate T C of a classical Heisenberg model is the Monte Carlo method, in particular when taking advantage of the cumulant expansion to account for the finite size of the simulation supercells. In this work we applied the Monte Carlo method, locating T C by the peak in the temperature-dependent, static susceptibility. As the simulation supercells are rather large, finite-size corrections to T C are small. We also show Monte Carlo magnetization curves for Co 2 MnSi in Fig. 5. 47 The classical Heisenberg Hamiltonian, Eq. (1) was used to model the systems, assuming that the magnetic moments can change their orientation, but not their length. NiMnSb is discussed in more detail in Sec. IV. In the calculations the Metropolis method was employed. A supercell of 2048 unit cells (4096 magnetic atoms) was used for NiMnSb and one of 2197 unit cells (6591 magnetic atoms) for Co 2 MnSi; interactions to neighbors up to a distance of four lattice constants were included. For each temperature the number of Monte Carlo sampling events was 5000, after allowing an initial relaxation time of 1000 steps and taking one sampling event every 10 sweeps of the lattice.
The magnetization curves (Fig. 5) do not go to a sharp zero at T C , but rather have a tail, as a result of the finite supercell. The peak of the susceptibility is, on the other hand, rather sharp and its position can be used to determine the Curie temperature. From the positions of these peaks, we estimated for Co 2 MnSi (Fig. 5, left) T C =980 K (experimental 48 985 K), and for NiMnSb (see Sec. IV) T C =870 K (experimental 48 value being 730 K). Note that here, for NiMnSb, the Ni moment was also taken into account within the Heisenberg model (but see also the discussion in Sec. IV, in particular IV C). For fcc Co we find a Curie temperature of 1200 K, while the experimental result is 1403 K (the high-temperature stable phase of Co is fcc).

IV. CONSIDERATION OF LONGITUDINAL MOMENT FLUCTUATIONS
So far we have discussed a formalism and examples of calculation of exchange parameters for a Heisenberg model, assuming that the magnitude of the local magnetic moments is rigid. However, as is long known, certain systems are weakly magnetic and the moment formation takes place at a rela-tively low energy scale, comparable to the magnon energies. Moriya's quantum mechanical spin fluctuation theory 49 already factors in these effects. An introduction of a classical Hamiltonian including longitudinal and transverse degrees of freedom has been done by Uhl and Kübler 50 by parametrizing ab-initio total energy results; thermodynamic quantities were then obtained within a mean-field approach which coupled longitudinal and transverse degrees of freedom. 37,50 More elaborate approaches by a classical fit of the energy to densityfunctional results, including transversal and longitudinal degrees of freedom, together with Monte-Carlo or classical spin-dynamics calculations, were presented e.g. by Rosengaard and Johansson, 51 Ruban et al., 20 Ma and Dudarev, 52 or Derlet. 53 Particularly interesting are compounds of a strongly magnetic and a weakly magnetic subsystem, where the weak moments cannot be treated as having rigid-magnitude, but on the other hand are large enough that they cannot be neglected. Examples are FePt, FePd, or NiMnSb, where the strongly magnetic atoms are Fe and Mn, while the weakly magnetic atoms are Pt, Pd and Ni. Then the assumption of rigid Heisenberg spins does no longer seem plausible. One way to circumvent this is to use only the strong magnetic moments as independent variables, but with their pair exchange parameters renormalized by the exchange among the weakly magnetic atoms and by the enhanced susceptibility. Such an approach was, for example, used by Mryasov et al. 54,55 to model the temperature dependence of the magnetic anisotropy in FePt. Sandratskii et al. 56 also find improvement in the theoretical results if the weak moments are not treated as independent, rigid Heisenberg spins, but rather fully adjust to the magnetization of the strong-moment sublattice. There is, however, the point of view of treating the weak moments not as fully enslaved to their strong-moment environment, but as independent degrees of freedom whose fluctuation can affect the thermodynamics. This point of view has been pushed for example in order to understand the antiferromagnetic-to-ferromagnetic transition of FeRh. 57,58 In this section we develop a method of treating such systems, emphasizing the choise of the energetically relevant magnetic degrees of freedom at temperatures up to the critical temperature, as well as the practical implementation of these degrees of freedom to models that can be solved via Monte Carlo simulations. The additional challenge, compared to the theory developed in the previous sections, stems from two facts. First, a treatment of the weak moments as independent degrees of freedom must include the energy scale of the moments magnitude as well as their direction. Second, if the interactions of the strong moments are calculated as described earlier, by considering non-collinear magnetic configurations, then the calculated energies include a spurious contribution because the magnitude of the weak moments changes during the application of the non-collinear constraint. This spurious contribution must be accounted for. However, at the end we see that, in certain commonly occuring situations, the degrees of freedom of the weak-moment atoms can be eliminated in the statistical-mechanical calculations by using renormalized exchange parameters, while the statistical averages of the weak moments follow directly from the averages of the strong moments. We show that this is an exact result if the weak degrees of freedom give only quadratic contributions to the energy. Here we make a step-by-step approach for the special case where the weak moments do not interact with each other, while we present the more general case when they interact in the Appendix.

A. General approach
The magnetically strong atoms are the main carriers of spin moment that also acts as an effective field H eff polarizing the magnetically weak atoms. The polarization of the latter is strengthened by an enhanced local susceptibility, which (neglecting complications from the band-structure) is χ = χ Pauli /(1 − Iχ Pauli ) (with χ Pauli the Pauli paramagnetic susceptibility and I the exchange integral). Let us denote by M the rigid-magnitude local moment of the magnetically strong atoms and by µ the supple local moment of the magnetically weak atoms. A simple expression that parametrizes the energy in terms of the local moment µ and the polarizing field H eff is It is implied that H eff = κ n∈neighb M n + H ext , where M n is the summed moment of the neighboring atoms inducing a polarization, κ is a phenomenological parameter which encapsulates all microscopic processes (in particular electron hoppings) that couple µ to M n , while H ext is an external magnetic field that we henceforth set to zero. 59 In the absence of a field H eff Eq. (14) contains only even powers of µ, because symmetry requires E loc [0, µ] = E loc [0, −µ]. In Eq. (14) we always must have b ≥ 0. The case of spontaneous polarization is described by a < 0 and b > 0, while the case of induced magnetic moments is described by a > 0. In this case, the on-site susceptibility is Thus the enhancement of the susceptibility is contained in the parameter a.
Eq. (14), if the fourth-order term is negligible, can be written as To elaborate on these ideas we use NiMnSb as a concrete example. In this case the Mn subsystem is magnetically strong and the Ni subsystem magnetically weak. We have also found that in this case it is a good approximation to set b = 0 (as we show below), which we adopt. Combining Eq. (14) at b = 0 with the Heisenberg-model energy expression for the Mn subsystem results in the extended model Hamiltonian where n(l) M n is the sum of the moments of the Mn nearest neighbors of the l-th Ni atom, n(l), and M i and µ l refer to the Mn and Ni moment respectively. The superscript "b" in the Mn-Mn exchange interaction is placed in anticipation that it does not coincide with the quantity calculated within the method presented in Sec. II, but with a "bare" quantity, as will be discussed below. It has been also assumed that the Ni-Ni interaction can be neglected, as there are no Ni-Ni nearest neighbor pairs, while from the calculations it becomes evident that more distant Ni-Ni interactions are negligible (the case where the interactions among the weak-sublattice moments are non-negligible is discussed in the Appendix). The above equation can be rewritten in a somewhat more convenient form, if we define (with an obvious notation) Then Eq. (17) takes the form This form is appealing, as it contains a Heisenberg-like expression among all magnetic atoms plus a local correction term that accounts for the longitudinal degree of freedom of the Ni moments. However, it is important to remember that the Ni moments included in the scalar products of the RHS can change in size as well as angle, deviating from the traditional Heisenberg model. We can derive the parameters κ and a from DFT calculations as follows. From Eq. (16) it follows that the equilibrium value of the Ni moment at T = 0 and for zero external field is linearly dependent to the polarizing neigboring Mn moments: where N c = 4 is the coordination number of a Ni atom with respect to Mn neighbors. From this expression one can deduce the ratio κ/a once M and µ eq have been calculated in the ferromagnetic state by an ab-initio method. On the other hand, a can be determined by introducing into the DFT calculation a constraint on the Ni moment, either in the form of a longitudinal magnetic field constraining the magnitude or in the form of a transverse magnetic field constraining the angle θ of µ with respect to the moment directions at the Mn neighbors, but allowing the magnitude to relax to µ eq (θ). The former method results in a parabolic (for b = 0) moment-dependence of the total energy, as suggested by Eqs. (14,16). The latter method results in the following dependence, as can be easily found by an energy minimization of Eq. (16) at a given θ: If the extended Heisenberg model provides a good approximation to the energetics of the system, then Eqs. (16) and (23) should both yield a good approximation to the calculated DFT energy difference, when a constraint is applied on the Ni moment, and the extracted parameter a for the two cases should be approximately the same. We find that this is the case in NiMnSb (see Fig 6). 60 This allows us to extract the value of a from the DFT calculation and from this the value of κ (equivalent to J b Mn−Ni ) via Eq. (21). Thus the ingredients of formula (20) are accessible.
There remains, however, a correction to be made, connected to a "renormalization" of the Mn-Mn exchange constants which have been calculated within the spin-spiral method described in Section II. To clarify the problem we remind the reader that, when calculating spin spirals, one normally constrains the direction of the moments but not their magnitude. For the Ni atoms, where the rigid-moment approximation is not valid, the magnitude µ is reduced by the spin-spiral formation due to the canting of the neighboring Mn moments in different directions. This effect provides an additional, indirect energy contribution to the spin spiral and to the Mn-Mn interaction, compared to the case that the Ni-moments magnitude would have been kept constant. It is the spurious contribution that we mentioned in the introduction to Sec. IV. Let us call this contribution J Mn−Ni−Mn . The calculated Mn-Mn interaction consists thus of two parts: with J b Mn−Mn the sought-after bare interaction entering Eqs. (18) and (20) Here, n and n ′ run through the Mn moments. The second step follows under the condition that, in the DFT calculation, the Ni moment relaxes to the particular equilibrium value that is dictated by the neighboring Mn-moment directions, i.e., µ = κ 2a n M n . From Eq. (25) we obtain the indirect interaction: At this point it is important to note that, whether one chooses to work with J Mn−Mn or J b Mn−Mn , depends on the choice of the degrees of freedom. If only the Mn moments are chosen as degrees of freedom, then J Mn−Mn has to be used. If, however, the Ni moments are also chosen as independent degrees of freedom, then J b Mn−Mn has to be used together with J b Mn−Ni = κ. In the latter case, according to the prescription leading to Eqs. (18) and (26), we arrive at the extended Heisenberg Hamiltonian (20) with the appropriate bare parameters in the nearest-neighbor coupling.

B. Analytical elimination of weak degrees of freedom at T > 0
In the previous subsection we discussed the bare and renormalized parameters of the model arising from total-energy calculations of the constrained ground state. In this section we examine the case of thermodynamic quantities at T > 0, where it is not a priori obvious that the same renormalization is still valid if the weak moments are allowed to fluctuate. We find that, in the absence of fourth order terms in E loc [b = 0 in Eq. (14)], also at T > 0 the weak moments can be eliminated in favour of the same renormalized parameters as the ones appearing at T = 0. Our conclusion is based on an exact analytical integration of the weak-moment part of the partition function in the case b = 0.
First we observe that the energy functional (16) under the action of H eff has the same quadratic form as the one with H eff = 0, only with the minimum shifted to µ eq = H eff /(2a). This simplifies the integration of the partition function. To see this, we first transform the Hamiltonian (17) in such a way that the renormalized interactions J ij appear explicitly. We start by rewriting Eq. (17) as: The indices i, j run over the Mn atoms, the index l over the Ni atoms, and n(l) over the Mn neighbours of the l-th Ni atom. The last term is now split in an interatomic contribution and an on-site contribution: In the last step we have converted the sum l n(l),n ′ (l)n =n ′ to a sum over i, j, by introducing a combinatorial factor c(i, j) that counts how many common Ni neighbors the i-th and j-th where the Heisenberg Hamiltonian involving only the Mn sublattice with the renormalized interactions, H r = − 1 2 ij J ij M i · M j , has been introduced for convenience as part of the full Hamiltonian. Note that there is no renormalization for the Ni-Mn interactions: either they should be completely omitted within H r that includes only the Mn moments, or they should be given by Eq. (18) in the generalized model. The spin-spiral result for J Mn−Ni does not enter in our theory.
The partition function is: where dΩ 1 · · · dΩ N denotes an integration over the Mn-moment solid angles. The integration over the Ni-moments d 3 µ 1 · · · d 3 µ N contains only the exponential of a complete square and has been analytically integrated to (πk B T /a) 3N/2 , i.e., it corresponds to the partition function of 3N uncoupled harmonic degrees of freedom and is independent of the value of M n . Z r is the partition function corresponding to the Hamiltonian H r . The magnetization M of the Mn sublattice also turns out to depend only on the renormalized partition function (and Hamiltonian). The same is true for any moment-moment correlation function of the Mn sublattice, (M 1x ) m1x (M 1y ) m1y · · · (M N z ) mNz , α ∈ (x, y, z), where the m nα ≥ 0 are integer exponents defining the order of the correlation function. To see this we consider explicitly where the Hamiltonian in the exponential has again been split in two terms according to Eq. (29) and the integration over d 3 µ 1 · · · d 3 µ N has again been carried out analytically, cancelling out the factor (πk B T /a) 3N/2 in the full partition function. But expression (35) is just the expression for the correlation function of only the Mn sublattice with the renormalized Heisenberg interactions. Thus we see that, as far as the Mn-sublattice magnetization is concerned, one can use the exchange interactions calculated by the spin-spiral method that are renormalized by construction, while neglecting the Ni moments and the Mn-Ni interaction. If the Ni moments are to be included, then the bare interactions have to be used in an extended model, however, the Mn-sublattice properties in the two cases will be the same. Note, finally, that this exact result is based on the fact that the exponent in the integration over d 3 µ 1 · · · d 3 µ N contains a complete square, i.e., that the "harmonic" approximation, b = 0, is valid; if b = 0 and one proceeds to an elimination of the weak degrees of freedom, then the resulting renormalized strong-moment Hamiltonian can have higher-order terms (biquadratic, four-spin, etc.) with temperature-dependent parameters. A semi-analytical result follows in an analogous way also for the average Ni moment and fluctuation amplitude per atom: I.e., there is a "harmonic" part and a part induced by the fluctuation of the Mn sublattice. The former is independent of the number of atoms, while the latter, proportionally to the number of magnetic atoms in the system for T < T C , as in normal Heisenberg systems. 62 The longitudinal susceptibility χ Ni of the Ni sublattice can be found in a completely analogous way: i.e., there is again a harmonic part and a part proportional to the Mn sublattice susceptibility. (In Eq. (39), z is implied to be the direction of magnetization.) According to these results, as long as the approximation b = 0 holds, one can deduce the thermodynamic properties of the weak sublattice by merely a calculation on the strong sublattice, avoiding the extra numetical cost that a full Monte Carlo simulation entails. Actually this procedure can be carried out for higher-order correlation functions of the Ni moments, (µ 1x ) m1x (µ 1y ) m1y · · · (µ N z ) mNz (α ∈ (x, y, z)), where once again the m lα ≥ 0 are integer exponents defining the order of the correlation function. In the resulting formula the integration over d 3 µ l can be carried out analytically, yielding a sum of correlation functions of the M i of the form: where a change of variables µ lα → µ lα + κ 2a n(l) M nα has been performed, the binomial expansion of (µ lα + κ 2a n(l) M nα ) m lα has been used and we have taken into account that integrals of the type  (38) and (40) are special cases of application of this formula. Obviously, mixed correlation functions between strong and weak moments can also be reduced to strong-moment correlation functions of the renormalized Hamiltonian by first eliminating the weak moments following the same prescription.
We should stress that choosing to eliminate the weak moments in favour of renormalized interactions does not mean that the weak moments are physically less valid as degrees of freedom; to do so is merely a matter of mathematical or computational convenience, especially since efficient methods exist for the calculation of thermodynamic quantities within the classical Heisenberg model.

C. Calculations in NiMnSb
We exemplify the above results with calculations in NiMnSb. First we establish which interactions can be neglected. To this end we performed calculations of the exchange coupling parameters by the spin-spiral method. As it turns out, the Ni-Ni interactions are negligible. About the Mn-Mn interactions, Fig. 3b suggests that it should be enough to include up to 2nd neigbors, but this is misleading; we find by Monte Carlo calculations a 2nd neighbor approximation leads to an overestimation of approximately 100 K in T C compared to the value including also more distant neighbors, therefore we include Mn-Mn interactions up to a distance of three lattice parameters. Among the Mn-Ni interactions (Fig. 3b) only the nearest-neighbor coupling has some small influence, changing T C by a modest 20 K. Given the conclusions of the previous subsection, the Mn-Ni interaction should be excluded for an estimation of T C , since the Mn-Mn interactions are already renormalized (in the corresponding spin-spiral calculations the magnitude of the Ni moment was allowed to relax). The Mn-Ni interaction must be included if one wishes to extract information on the behavior of the Ni magnetiza-tion; then, however, the interaction type and strength has to be corrected, since the Ni atoms belong to a soft-magnetic sublattice.
After deriving the exchange parameters J ij by the spinspiral method, utilizing the force theorem which according to Fig. 2a is accurate enough for this purpose, we performed constrained-angle calculations for the Ni moments with full self-consistency. From the total-energy results shown in Mn−Ni is significantly different than the value of 0.92 mRyd/µ 2 B for the Mn-Ni interaction, which was derived from the spin-spiral calculation (Fig. 3B). Finally, from Eqs. (24) and (26)  Mn−Ni M µ eq = 1.87 mRyd; i.e., the Mn-Mn bare interaction is weaker than the Mn-Ni, which is not surprising, since the Ni polarization stems from a direct, nearest-neighbor exchange interaction with Mn.
Next we present a series of Monte-Carlo-simulation results, examining the effect and importance of the bare interactions. We performed simulations in the framework of the traditional model (only transverse fluctuations allowed) as well as the extended model (longitudinal fluctuations also allowed), with the renormalized and bare parameters. Putting it all together, we must substitute the above-found parameters Mn−Ni , and a, into the extended Hamiltonian (20) which is to be treated with a Monte Carlo method where the (weak) Ni moments µ l should be allowed to vary in length and direction while the (strong) Mn moments should be allowed to vary in direction only. In this particular example, a distinction between bare and renormalized interactions was made only for the nearest-neighbour Mn-Mn and Mn-Ni coupling. Hamiltonian containing only the Mn moments that vary in direction, where now J ij are the Mn-Mn renormalized parameters, i.e., obtained from the spin-spiral method presented in Sec. II and shown (multiplied by the moments) in Fig. 3b.
In both cases the Mn-Mn interaction was included up to a distance of three lattice parameters. The first set of results is contained in Fig. 7, where magnetization curves are shown, calculated within different assumptions. Here we omit showing the susceptibiliy, but the Curie temperature can be recognized by the characteristic inflection point of the magnetization curve. The central results are included in curves (a) and (e) of the top panel, as well as in the bottom panel, while curves (b), (c) and (d) merely show that neglecting some of the degrees of freedom or some of the bare constants leads to almost arbitrary results. In the simulations including longitudinal moment fluctuations, more Monte Carlo steps were necessary, typically by an order of magnitude, in order to arrive at the same quality as in the simulations including only transverse fluctuations.
Curve (a) shows the result of the traditional model with spin-spiral-derived (i.e., renormalized) Mn-Mn exchange parameters, including spin-spiral-derived Mn-Ni interactions. This results in T C = 870 K. Curve (b) includes the bare Mn-Ni parameters, but stays within the traditional model. T C increases by a modest amount of 60 K, since the bare Mn-Ni parameters are stronger than the spin-spiral-calculated ones. Curve (c) also stays within the traditional model, but includes the bare Mn-Mn interaction, which is weaker by a factor 1/3 compared to the renormalized value; T C drops significantly to 730 K. This coincides with the experimental value, but the coincidence is probably fortuitous: the longitudinal fluctuations at Ni, that are essential to the derivation of the bare parameters, are yet unaccounted for in the simulation (see, however, the discussion on the phase space measure in subsection V). Next comes curve (d), where the longitudinal fluctuations are allowed for in the simulation, but taking the bare Mn-Ni and the renormalized Mn-Mn exchange. The increase in T C (1090 K) with respect to curve (b) (same parameters but rigid Ni moment) is striking. The difference stems from the fact that the Ni local moment becomes larger at high temperatures (see discussion on Fig. 8 below), thus the ferromagnetic contribution J Mn−Ni M · µ to the Hamiltonian is effectively strengthened. Finally, if we account also for the (weaker) bare Mn-Mn coupling in the extended model, we obtain curve (e), which is very close to the original curve (a), also with a very similar T C = 850 K. This value, however, is the same that one obtains if the traditional model is used, but with the Ni moments completely neglected. In fact, the corresponding magnetization curve falls exactly on top of the Mn-contribution to curve (e), and the same is true for the Mn-sublattice susceptibility. This striking agreement is also demonstrated in the bottom panel of Fig. 7, and is expected on the basis of Eq. (35) and the discussion thereafter.
Let us consider now the contrast between the traditional model with spin-spiral calculated interactions and the extended model with bare interactions. Monte-Carlo results on these are shown in Fig. 8a and b, representing the traditional and extended model respectively. The main difference lies in the behavior of the Ni-sublattice magnetization µ z and susceptibility χ Ni . For better comparison we have scaled up these quantities. In Fig. 8a we see that the traditional model results in a Ni magnetization that drops rather fast with temperature, much faster than the Mn magnetization. This is due to the relatively weak coupling between the Mn and Ni sublattices. Although the difference in energy scale is not directly obvi-ous from Fig. 3, recall that a Mn atom is surrounded by 4 only Ni atoms but by 12 Mn atoms at distance a lat / √ 2 and 6 Mn atoms at distance a lat , where the exchange coupling is still appreciable. The Ni susceptibility, χ Ni shows a large plateau below T C and follows the critical peak of Mn at T C rather weakly. The behavior of the Mn magnetization and susceptibility, on the other hand, is characteristic of a ferromagnetic phase transition.
In Fig. 8b we see results within the extended model. Here, the Ni magnetization follows the behavior of Mn. The susceptibility χ Ni starts off at a finite value at T = 0, which coincides with the value given by Eq. 15, in contrast to the vanishing susceptibility at T = 0 for rigid-moment systems. At T C the Ni susceptibility shows a peak following the critical behavior. The Ni magnetization and susceptibility have been scaled up by appropriate factors to show that Eqs. (36) and (40) are reproduced by the Monte Carlo simulation. A technical point to be mentioned for accuracy is that, in calculating µ z and µ 2 z in the simulation one should take the projection of l µ l in the direction of the total moment at each Monte Carlo step, instead of making the approximation to calculate the average and variance of | l µ l |.
It is also interesting to see that the average magnitude of the local Ni moment (dotted line) increases monotonically with temperature. This effect, pointed out by Sandratskii 61 in an analytical low-temperature approximation for NiMnSb, is connected to the fact that even above T C there is some shortrange order in the system, so that the equilibrium value of µ = κ 2a N c M is non-vanishing. On top of this, the fluctuations of µ shift the value of |µ| = (µ 2 x + µ 2 y + µ 2 z ) 1/2 to higher values. The increase is expected to cease when µ reaches such high values that the approximation b = 0 is no more valid (so that the fluctuations are moderated by the fourth-order term); this correction does not apply for NiMnSb, however, at least at temperatures as high as T C , since the constrained-DFT calculation shown in Fig. 6 (top panel) yields a quadratic energy dependence at values of µ comparable to the ones close to T C in Fig. 8b.
Experimentally, the temperature dependent magnetization of NiMnSb follows Bloch's T 3/2 law up to about 70-100 K, 63,64 which cannot be reproduced within classical Heisenberg models. Element (or sublattice) specific experiments were done using neutron scattering 65 and X-ray magnetic circular dichroism (XMCD). 66 Neutron scattering 65 at 15 K and 260 K shows a thermally stable Ni moment (0.18 µ B ), but a decreasing Mn moment (dropping from 3.79 µ B at 15 K to 3.55 µ B at 260 K). Contrary to this, according to the XMCD results of Ref. [66] both the Ni and Mn moments drop rapidly at 80 K to half their ground-state values, and then level off up to at least 250 K. Such behavior would be counter-intuitive; the authors in [66] write that surface effects possibly complicate the interpretation of MCD data.

A. Remarks on the treatment of the weak moments and on the concept of renormalization
Recently, Wysocki, Glasbrenner and Belashchenko 67 (WGB) presented a study of a classical spin-fluctuation model. Their model Hamiltonian is similar to the one that we use here, with the difference that in the WGB paper all atoms can change the magnitude of their moments, the 4thorder term is not neglected, and in practice only one atomic species is considered in their calculations. WGB point out that the magnetic moments are not canonical variables, therefore there is no obvious way how to choose the phase space measure, which should therefore be given as part of the model together with the Hamiltonian. In the present work we have chosen what they call uniform phase space measure, which basically amounts to dividing the µ-space in equal-volume infinitesimal cells with equal integration weight, and which is the most common choice in the literature. 68 It also amounts to taking a simple d 3 µ integration with no further weight in Eqs. (30,34,41). Different choices of measure can lead to qualitatively different results, e.g. a fast drop of magnetization in the weakly-magnetic sublattice, as also shown by Sandratskii,61 or even a first-order transition, as WGB find. 67 The correct measure can only be certified by the best classical approximation to the full quantum-mechanical solution, which, however, remains an open problem. However, if it is assumed that the correct measure is not uniform but e.g. proportional to µ −2 (this was one alternative choise by WGB), then the weak-sublattice magnetization will drop fast at low temperatures, so that close to T C only the bare parameters of the strong sublattice would be relevant.
Mryasov et al. 54 have discussed the idea of renormalized interactions in the case of FePt alloys. Related is also the work by Polesya et al. 21 who adopt a model for FePd and CoPt alloys. In these works, the weak moments of Pd or Pt are determined from the strong moments of their neighbors via the susceptibility. However, even at higher temperatures the weak moments are not treated as independent variables that can fluctuate (either in direction or in length) but rather as enslaved quantities to their immediate neighbourhood i.e., their role in the thermodynamics is only to mediate an additional interaction between the strong moments. Our idea of renormalization has an analogous starting point but goes along a different path, since we consider the weak moments in the spin Hamiltonian as independent fluctuating variables.
In addition, the main novelty of the present study in this respect is the formal and numerical proof of the equivalence of two approaches: all thermodynamic quantities can be derived by using either the extended model with bare interactions or the traditional model with renormalized interactions. Furthermore it is shown in the present study that the renormalized intercations are actually the ones that are harvested by the spinspiral approach within density-functional calculations with no further manipulations; however, manipulations are necessary if one wishes to extract the bare parameters. The present result holds under the assumption of a uniform phase space measure and of a quadratic on-site energy. The latter seems to hold true e.g. for NiMnSb and for FePt, 54,55,61 but not, for example, in the case of FeRh 55,58 where higher-order corrections are necessary. In case that either of these two requirements is not met, most probably a hypothesis of temperature-independent renormalized interactions cannot be justified on the grounds of fluctuating weak moments, but should instead be conjectured as an ad-hoc hypothesis within the model.
Bruno 70 also introduces a concept of renormalized exchange parameters. However, his approach encapsulates different physics than our present approach. Bruno's renormalization corrects for a systematic error, mainly due to the difference of the constraining-field direction to the resulting moment direction when the force theorem is applied. Our renormalization on the other hand concerns the error due to the reduction of the weak-moment magnitude when the strong moments are tilted.
Yet another concept of renormalization is described by Lounis and Dederichs. 71 Using a multiple-scattering approach, they consider the energy expansion as a function of the angles between moments. As they find, at high angles corrections are necessary to the phenomenological Hamiltonian (e.g., biquadratic or four-spin terms). However, at low angles these corrections can be partly included in the Heisenberg model via a renormalization of the exchange parameters, recovering correct energy scales and Curie temperatures.

B. Remarks on the prediction of the Curie temperature
The Curie temperatures calculated within the mainstream approach to the adiabatic spin-dynamics are in many cases in agreement with experiment to within 10-15%, but with no obvious systematics toward over-or underestimation. The main source of error is not clear. Considering the most serious approximations made, error can stem from • (i) the use of local density functional theory (LDA or GGA) for total energy calculations without further corrections for electron exchange and correlation.
• (ii) the use of the adiabatic approximation, • (iii) the assumption of a classical, rather than quantum, Heisenberg model, • (iv) the assumption that the exchange constants do not change as a function of temperature, • (v) the assumption of rigid spins of the strong moments in the Heisenberg model.
In general, these factors have possibly different weight for different materials. Concerning point (i), theories that provide a better treatment of correlations exist, e.g., the LSDA+U or LSDA combined with dynamical mean field theory at zero or finite temperatures. Also within such theories exchange parameters can be derived (see, e.g., Ref. 18). However, parameters are required (as is the Coulomb repulsion U ), and it is usually not obvious how to determine these uniquely.
As for point (ii), there are promising developments in the calculation of magnon spectra (including magnon lifetime effects) within time-dependent density functional theory [72][73][74] or many-body perturbation theory 75 . These can prove very useful in the future (they can also be extended to finite temperatures), however at this point they are computationally too demanding for systematic calculations of the Curie temperature.
Point (iii) (the classical assumption), can be improved upon by using the random phase approximation to solve the quantum Heisenberg model. However, for itinerant electron systems the local moment M does not corespond to some integer or half-integer value of the spin S either in the form M = S(S + 1) µ B or in the form M z = S µ B . In fact, calculations of Heusler compounds in Ref. [7] have shown that for a reasonable choice of S the Curie temperature is strongly overestimated, while the classical limit of the random phase approximation, with a choice of large S (with appropriate normalization of the exchange parameters so that the product J ij S i S j remains constant), results in reasonable values of T C . Therefore, a quantum Heisenberg model is perhaps a better for a correct description of the shape of the magnetization curve M (T ), but a poor choice for a correct T C , at least in itinerant electron systems, if the exchange constants are calculated within the adiabatic approximation. This conclusion is in accordance to the spirit of adiabatic spin dynamcics, 13 where the effective interactions J nn ′ correspond to the equation of motion of the expectation value of the local moments, i.e., to classical quantities, not operators.
Corrections to point (iv) can be treated within local density functional theory if the exchange constants are calculated starting from a disordered local moment state. This requires use of the coherent-potential approximation (CPA), and has been proposed for example in Ref. [3]. The use of the CPA for the description of the disordered local moment state at T C underestimates the existence of magnetic short range order (which is known to be present); however, it constitutes a promissing approach, since it can be systematically improved e.g. by the use of a non-local CPA. 76 Finally, point (v) becomes a serious approximation in systems of weak magnetic moments, such as ferromagnetic Ni, and has been widely discussed in the literature as we noted in Sec. IV. Corresponding corrections for multiple-scattering based methods have been recently proposed, e.g., by Bruno 70 and Shallcross and co-workers. 3 In a more recent work by Ruban et al., 20 based on an expansion of the energy within the disordered local moment state, promising results were obtained showing the fundamental importance of longitudinal corrections to the local moment for T C in ferromagnetic Ni.

VI. SUMMARY AND CONCLUSIONS
In the first part of this work we have investigated the calculation of interatomic exchange constants that we implemented in the FLAPW method based on the concept of adiabatic spin dynamics.
The exchange constants are harvested by an inverse Fourier transformation involving static spin-spiral energies. 14 Symmetry relations obeyed by the spin spiral energies have been found to greatly reduce the numerical effort, in particular regarding confinement of the inverse Fourier transformation in the irreducible wedge of the Brillouin zone. Furthermore, the force-theorem approximation has been tested and found to be adequate for small cone angles of the spin spirals. However, we have shown that application of the force theorem requires special treatment of the intersitial region, namely setting there the magnetic part of the exchange-correlation field to zero.
In the second part of the present work we have proposed a way to explore multicomponent systems where a magnetically strong sublattice coexists with a magnetically weak sublattice, necessiating a consideration of longitudinal and transverse changes of the weak local magnetic moments while they are still treated as independent variables. We find the rigorous result that, under the frequently met condition of a parabolic dependence of the energy on the weak-moment magnitude, the weak moments and their interactions can be eliminated via an analytical integration of the partition function in favour of the strong moments with renormalized, temperature-independent exchange constants, with the renormalization accounting for the weak-moment fluctuations at non-zero temperatures. We also show that the renormalized constants are actually the ones probed by constrained spin-spiral calculations of the strongmoment subsystem, thus simplifying calculations. Finally we show that the thermodynamic correlation functions of the full system including the strong and weak moments can be derived as polynomials of the correlation functions of the system of strong moments only but with renormalized interactions.
This renormalization will affect various quantities such as temperature-dependent magnetization, susceptibility, or spin-stiffness constant. The method can prove useful for systems comprising 3d atoms with strong moments together with 4d or 5d atoms with weak moments, such as transitionmetal alloys, Heusler alloys or 3d overlayers deposited on 4d or 5d metal surfaces. ered, assuming that Eqs. (44)(45)(46) are a good approximation to the DFT energy landscape. The answer is straightforward if we calculate the total energy of the constrained Heisenberg Hamiltonian, i.e., without an integration over the M i . The constrained partition function is just the last term of Eq. (48), Z 1 = d 3 µ 1 · · · d 3 µ Nw e −(Hw+Hint)/kB T . The total energy at T → 0 is q.e.d. This means that constrained (spin spiral) DFT calculations on the strong sublattice are already corresponding to the renormalized Hamiltonian and are therefore by virtue of Eqs. (48,50) sufficient for the calculation of the strongsublattice thermal averages, without the need to calculate the bare parameters or A ll ′ . However, if the weak-sublattice thermal averages are to be calculated, one must additionally gain knowledge on the matrix A ll ′ of Eq. (45) as well as of the bare parameters J b of Eq. (45) and of κ in Eq. (46). The scheme presented in Sec. IV A shows how this can be done in the rather simple case of NiMnSb (where A ll ′ = aδ ll ′ is diagonal), but in general this problem must be solved according to the geometry and other factors in each case. In particular for the calcula-tion of the off-diagonal A ll ′ , probably it is easiest to calculate directly the susceptibility matrix (A −1 ) ll ′ by applying in the DFT calculation a longitudinal external field on one atom and probing the response in the moment of the neighboring atoms, or to make a transformation to the normal modes in Fourier space and calculate A(q).
If all the ingredients are available, then one can calculate thermal averages and correlation functions of the weak moments either by a direct Monte-Carlo calculation, or by reducing these correlation functions to correlations of the strong moments in the presence of only the renormalized Hamiltonian, in an analogous way to the case of Eq. (41), where the parameters A ll ′ and κ must be inevitably contained in the expansion coefficients. We present an outline of how this is achieved in practice. The complication here compared to Eq. (41) is that the matrix A ll ′ is not diagonal. It is, however, real and symmetric, therefore one proceeds by bringing it to a diagonal form D with elements D l . If R is the diagonalization matrix, then D = RDR T , and we define the transformation of the moment-coordinatesμ qα = l R ql µ lα andM qα = n R qn M nα . Also since A ll ′ is real and symmetric the phase-space element is unchanged: d 3μ = d 3 µ. (This transformation is analogous to one bringing a system of coupled harmonic oscillators in a normal-mode representation.) Then one has, for the arbitrary correlation function among the xyz-components of the N w weak moments, (µ 1x ) m1x · · · (µ Nwz ) mN w z , the following integral: The last expression can be handled by expanding the term ( k R T lkμ kα ) m lα and proceeding to an analytic integration of powers ofμ times the exponential just as in Eq. (41). Mixed correlation functions between strong and weak moments can also be reduced to strong-moment correlation functions in the renormalized model via this prescribed route.