Prediction of the Curie temperature considering the dependence of the phonon free energy on magnetic states

Prediction of the Curie temperature is of significant importance for the design of ferromagnetic materials. Even though the Curie temperature has been estimated using the Heisenberg model, magnetic exchange coupling parameters widely used is thus far based on first-principles calculations at zero temperature. In the explicit consideration of temperature effects, it is important to minimise the total free energy, because the magnetic and phonon free energies correlate with each other. Here, we propose a first-principles thermodynamic approach to minimise the total free energy considering both the influences of magnetism on phonons and the feedback effect from phonons to magnetism. By applying our scheme to bcc Fe, we find a significant reduction of the Curie temperature due to the feedback effect. This result inevitably enforces us to change our convention as follows: we should use exchange coupling constants for the disordered local moment state, not for the ferromagnetic state, in the prediction of the Curie temperature. Our results not only change the fundamental understanding of finite-temperature magnetism but also provide a general framework to predict the Curie temperature more accurately.

The Curie temperature (T C ) is one of the essential properties of ferromagnetic materials because it characterises their applicability and performance 1,2 . The method of predicting T C is, therefore, important not only for a fundamental understanding of ferromagnetic materials but also for the material design for applications. A typical technique for predicting T C is a downfolding method from first-principles calculations to an effective lattice model as below: Firstly, one derives exchange coupling constants (J ij ) by applying Green's function-based methods 3,4 or a frozen magnon approach 5,6 . Secondly, one builds an effective lattice model such as the Heisenberg model and assign J ij to the model. Finally, one solves the model analytically or numerically and estimates T C . This technique is applied to a broad range of materials, such as 3d transition metals 4,7-11 and rareearth magnets [12][13][14][15][16][17] . Such many studies demonstrated that the prediction technique has some predictive accuracy.
Such a technique usually does not include temperature effects on magnetic interactions. Moreover, temperatureinduced interactions between magnetism and other excitations such as phonons sometimes make the accurate prediction of T C difficult. At a high-temperature range around T C , there are two kinds of interaction between magnetism and phonons. One is the effect of thermal atomic displacements on J ij [18][19][20] . The change in J ij obviously modifies T C . The other interaction is the effect of magnetic disordering on phonon frequencies. Some ferromagnetic materials such as bcc Fe and Pd 3 Fe shows phonon softening at elevated temperatures near T C 21-23 . Some research groups approached this phenomenon by different theoretical methods [23][24][25][26][27][28][29][30] and achieved the same conclusion: The phonon softening is due to magnetic disordering near T C . Regarding the predictive accuracy of T C , the importance of the former interaction is easily understandable, whereas the latter interaction does not apparently seem to be related to T C . However, we will recognise the phonon softening due to magnetic disordering is closely related to T C by standing a thermodynamic viewpoint.
Thermal equilibrium states at finite temperature correspond to the minimum of the total free energy at given conditions. This is usually called as the minimum principle for the free energy. Usual procedures to study finite-temperature magnetism is constructing a magnetic Hamiltonian and deriving thermodynamic quantities such as the magnetic energy and the magnetisation. This series of procedures is equal to interpret that equilibrium magnetic quantities are determined through the magnetic free energy only. This interpretation, however, collapses in the systems that exhibit the phonon softening due to magnetic disordering. The phonon frequencies are directly related to the phonon free energy. Thus the phonon softening due to magnetic disordering means that magnetic states affect the phonon free energy as well as the magnetic free energy. As a result, equilibrium magnetic states should be determined through not only the magnetic free energy but also the phonon free energy, according to the minimum principle for the free energy. We call this effect of phonons on equilibrium magnetic states through the change of the phonon free energy as a thermodynamic feedback effect. This feedback effect surely affects T C as a consequence of the change of equilibrium magnetic states. However, the significance of the feedback effect on T C is unclear because the existence of the effect has been overlooked.
In this article, we propose a thermodynamic formulation to treat the feedback effect from phonons to magnetism. The formulation results in a simple optimisation problem for the total free energy. The ingredients to solve the problem are evaluated by first-principles phonon calculations and Monte Carlo simulations based on the Heisenberg model. By applying the formulation to bcc Fe, we demonstrate that T C of bcc Fe significantly decreases by nearly 580 K. This result proves the feedback effect is crucial for accurate prediction of T C . We also discuss the relationship between the predictive accuracy of T C and reference magnetic states in the derivation of J ij . Remarkably, we find a significant overestimation of T C in a paramagnetic disordered local moment (DLM) state is rather a correct tendency. Quantitative description of finite-temperature magnetism plays an important role in both basic and applied materials science. Therefore, our results have an impact on the fundamental understanding of magnetism and materials design for ferromagnetic materials.
We organise the following part of this paper as below. Firstly, we introduce a new thermodynamic formulation to treat the thermodynamic feedback effect. Our formulation based on the minimum principle for the free energy is justified through the Legendre transformation and results in a simple optimisation problem. Next, we evaluate the magnetic entropy and the phonon free energy of bcc Fe as functions of the magnetic energy. These functions are needed to solve the optimisation problem. Finally, we evaluate the equilibrium magnetic energy around T C by solving the optimisation problem. The shift of T C of bcc Fe is estimated from the results of the equilibrium magnetic energy.

THERMODYNAMIC FORMULATION FOR MAGNETIC MATERIALS
In conventional thermodynamic approaches for magnetic materials, the phonon and magnetic contributions are assessed independently. We start from this typical case for comparison with our formulation. The fundamental relation is written as where E is the energy and S is the entropy. The subscripts tot, ph and mag represent total, phonon and magnetic, respectively. Here, we consider the Gibbs free energy, where T represents the temperature, p the pressure, V the volume, M the magnetisation, H the external magnetic field and µ 0 the vacuum permeability. In the following, we derive the formalism for p = 0 and H = 0, but the discussion remains unchanged for the case with finite external fields p and H. The Gibbs free energy G T = T 0 equilibrium E mag by the Heisenberg model equilibrium E mag by the minimisation of the total free energy G ph(T 0 , Emag) +G mag(T 0 , Emag) Free energy Schematic image of the free energy minimisation at a temperature T0. Within a common framework using the Heisenberg model, the equilibrium magnetic energy (Emag) is corresponding to the minimum of the magnetic free energy, Gmag (blue line). On the other hand, the equilibrium magnetic energy in our scheme is corresponding to the minimum of the total free energy, G ph + Gmag (orange line).
is derived by applying the Legendre transformation.
Here we used the one-to-one correspondence between energy and entropy for fixed other thermodynamic parameters such as V and M . The independent assessment in conventional approaches is based on this trivial formulation. Next, we incorporate the dependence of the phonon free energy on magnetic states. We assume that the magnitude of the interaction between magnetic disordering and phonon frequencies can be written as thermodynamic quantities of the magnetic part. Körmann et al.
proposed a solid treatment with this assumption 27 . They treated the forces on each atom as a function of the magnetic energy. As a result, the phonon frequencies, consequently the phonon free energy, have the dependence on the magnetic energy (see Methods). Thermodynamically speaking, their treatment means the phonon energy depends not only on the phonon entropy but also on the magnetic entropy. The fundamental relation thus can be written as In principle, E mag also has a dependence on S ph . This dependence can be regarded as influences of thermal atomic displacements on J ij 20 . If we want to incorporate this effect into the thermodynamic formulation, we have to express the magnitude of the effect as a thermodynamic quantity such as S ph . However, the correspondence between the thermal displacements and S ph is not obvious. We thus focus only the dependence of E ph on S mag .
We apply the Legendre transformation as before.
Note that the entropy (or energy) and the temperature can be treated as independent variables during the minimisation procedure. The thermodynamic relationship between the entropy and the temperature, such as ∂G/∂T = −S, holds after the minimisation, i.e. after the Legendre transformation. The last expression is very intuitive from a thermodynamic viewpoint: The equilibrium magnetic energy at a temperature T 0 is determined to minimise the total free energy ( Fig. 1) as = argmin

EVALUATIONS OF THE MAGNETIC ENTROPY AND THE PHONON FREE ENERGY
We demonstrate the significance of the dependence of the phonon free energy on magnetic states for bcc Fe as an example. As a starting point, we evaluate the magnetic entropy and the phonon free energy depending on the magnetic energy (S mag (E mag ) and G ph (T, E mag )), in order to solve the minimisation problem in equation (12).
To obtain S mag (E mag ), we carried out the rescaled Monte Carlo method 31 based on the Heisenberg model. This method brings thermodynamic quantities derived from classical Monte Carlo simulations closer to those from quantum Monte Carlo simulations. The exchange coupling constants (J ij ) in the Heisenberg model are derived from the paramagnetic disordered local moment (DLM) state 3,32 (see Methods). The magnetic energy and entropy as functions of lattice-model temperature T are shown in Fig. 2 (a). The theoretical T C (1522 K) is higher than the experimental value (1043 K). Such overestimation has also been reported in previous studies 3,20,33 using the DLM state. The overestimation has been recognised as a disadvantage of the DLM state, and it will be discussed later associated with our results. Since this magnetic system does not show the first-order phase transition, the one-to-one correspondence holds between not only E mag and S mag but also T , We constructed the function S mag (E mag ) ( Fig. 2 (b)) by using this relationship. Phonon frequencies depending on the magnetic energy (G ph (T, E mag )) can be calculated by using first-principles phonon calculations and Monte Carlo simulations following the previous research 27 (see Methods). The phonon dispersions and the phonon density of states of bcc Fe depending on the magnetic energy are shown in Fig. 3. The dependence of the frequencies on the magnetic energy is represented through the parameter α (see Methods). The calculated phonon dispersions in the ferromagnetic (FM, α = 1) and paramagnetic DLM (PM, α = 0) limits are consistent with the previous research 27 . Once the phonon frequencies at various magnetic energies (i.e. at various α) are calculated, the phonon free energy can be evaluated from the analytical form, where k B represents the Boltzmann constant, ω qj (E mag ) the phonon frequency of the j-th branch at the wave number q as a function of E mag and N q the total number of q points. Note that the more disordered the magnetic state is, the lower the phonon frequencies are. This tendency means the phonon free energies of paramagnetic states are smaller than that of the ferromagnetic state because of the monotonicity of the phonon free energy for the phonon frequency. Consequently, paramagnetic states are thermodynamically stabilised by the phonon softening effect.

TOTAL FREE ENERGY MINIMISATION
We are now able to proceed to the total free energy minimisation in equation (12) by using the functions G ph (T, E mag ) and S mag (E mag ). The minimisation procedures are simple. Firstly, we fix the temperature at T 0 . Secondly, we calculate the total free energy (G ph (T 0 , E mag )+E mag −T 0 S mag (E mag )) for various E mag values. The variable range of E mag is from the ferromagnetic limit to the paramagnetic limit. Thirdly, we find E mag corresponding to the minimum total free energy. The orange line in Fig. 1 is a visualisation of these steps. Finally, repeat these steps for a temperature range around T C .
The equilibrium magnetic energies of bcc Fe obtained by two difference methods are shown in Fig. 4: One is the minimisation of the total free energy G mag + G ph ; the other is the Monte Carlo simulations based on the Heisenberg model (the result is the same as the blue line in Fig. 2 (a)). Note that the result from the latter method is corresponding to that of considering only G mag in the minimisation of the free energy. The equilibrium magnetic energies obtained by the minimisation of the total free energy are larger than those of considering only G mag . This is, as mentioned before, due to the stabilisation of paramagnetic states by the phonon softening effect, and the magnitude of the stabilisation indicates that the phonon contribution is not negligible at all in the determination of equilibrium magnetic states around T C The stabilisation of paramagnetic states leads to a decrease in T C . As shown in Fig. 4, T C in the results of the minimisation of the total free energy (946 K) is lower than that of considering only G mag (1522 K), and the magnitude of the decrease reached nearly 580 K. Notably, T C of considering both G mag and G ph is dramatically closer to experimental value (1043 K) than that of considering only G mag , i.e. T C in the Heisenberg model. Although T C = 946 K is still underestimated the experimental value to some extent, the anharmonicity of phonons probably compensates for the deviation. Heine, Hellman and Broido 29 investigated the phonon softening phenomenon in bcc Fe with including anharmonic effects. They show that at 1043 K, where the anharmonicity is effective, the differences between the frequencies of the ferromagnetic and paramagnetic states are reduced compared with those at 300 K. Thus the difference of the phonon free energies between the ferromagnetic and paramagnetic states is also reduced. This consequently makes the degree of the decrease in T C smaller than our result. The underestimation in our result is, therefore,  Fig. 2 (a)). The Curie temperature TC in the minimisation of the total free energy is defined as the temperature with the same magnetic energy as the Heisenberg model. This definition is reasonable because the same magnetic energy gives the same magnetic ordering as long as Jij values do not vary. qualitatively correct.
The substantial decrease in T C gives a doubt on the usual recognition of the accuracy in prediction techniques for T C . Roughly speaking, there are three reference magnetic states in the derivation of J ij : ferromagnetic state, paramagnetic DLM state and conical spin-spiral states. The former two are used within Green's function-based methods 3,4 , whereas conical spin-spiral states are used within the frozen magnon approach 5,6 . We can summarise the relationship between the reference states and predictive accuracy of T C in bcc Fe: The ferromagnetic state and spin-spiral states give T C near the experimental value [4][5][6]8,11,34,35 , whereas the paramagnetic DLM state overestimates T C significantly 3,20,33 . Therefore, the ferromagnetic state and spin-spiral states have been recognised to have an enough predictive accuracy of T C regarding bcc Fe. However, our study clearly shows that this recognition is questionable because the contribution of the phonon free energy decreases T C of bcc Fe significantly. The substantial decrease in T C suggests the DLM state shows correct tendency regarding T C prediction, rather than the ferromagnetic state and spin-spiral states. Note that this suggestion is of great importance for theory of finite-temperature magnetism as follows. In the development of the theory, T C of bcc Fe has been recognised as a touchstone: Whether the predicted T C of bcc Fe agrees with the experimental value or not has been an element to examine the validity of a new theory. Our result, however, indicates such an examination way is inappropriate. Instead, an appropriate judgment criterion is as follows: Without considering the phonon softening, a theory that accurately describes finite-temperature magnetism must overestimate T C of bcc Fe.
Our thermodynamic formulation becomes complete if we incorporate the dependence of E mag on S ph . This dependence may be related to the effect of thermal atomic displacements on J ij . Ruban and Peil 20 studied this effect by combining J ij calculations and molecular dynamics. They clearly showed J ij values of bcc Fe were reduced by atomic displacements, and consequently, T C was also largely decreased compared with the case of excluding thermal atomic displacements. The dependence of E mag on S ph is thus important and intriguing from a thermodynamic viewpoint. However, a concrete expression of this dependence is yet to be obtained.

CONCLUSIONS
We have quantitatively evaluated the thermodynamic feedback effect from phonons to magnetism on T C regarding bcc Fe. The phonon softening due to magnetic disordering lead to the stabilisation of paramagnetic states. As a result, T C of bcc Fe was decreased by nearly 580 K from the value in the case of ignoring the feedback effect, i.e. the value for the Heisenberg model. This deviation in bcc Fe is of great importance because bcc Fe is recognised as a touchstone for the study of finite-temperature magnetism. We stress two important knowledge regarding the prediction of T C : (i) An appropriate theory of magnetism without considering the contribution of the phonon free energy must overestimate T C in bcc Fe, contrary to conventional understanding. (ii) We should use J ij for the DLM state rather than for the ferromagnetic state in the accurate description of T C .
Finally, we mention the applicability of our thermodynamic formulation. We focused on bcc Fe in this study, but our formulation is not restricted to it. It is intriguing to apply the formulation to other magnetic materials such as permanent magnets in which T C is critically important. In addition, the core concept of the formulation can be applied to other interacting excitation phenomena, not only the interaction between phonons and magnetism: If a contribution (X) affects other contribution (Y) and changes the free energy of Y, the thermal equilibrium state of X is also affected through the minimum principle for the free energy. In our study, X is magnetic states, and Y is phonons. Therefore, the concept of our thermodynamic formulation can be applied to other interacting excitations if one can express the magnitude of the interaction as thermodynamic quantities (E mag in our study). The formulation will be helpful for a quantitative description of the finite-temperature properties of materials.

First-principles phonon calculations.
All of the phonon calculations were carried out within the harmonic approximation. To evaluate the phonon frequencies at an intermediate magnetic ordering, we employed a force-averaging method 27 . In this method, the atomic forces at an intermediate magnetic ordering are determined by mixing the forces at the ferromagnetic (FM) and paramagnetic (PM) DLM states. Following the reference 27 , the atomic forces at an intermediate magnetic ordering can be written as where F i is the atomic force vector on i-th atom and α is a mixing parameter. They also proposed a solid expression of α by using the magnetic energy (E mag ) as below: where E PM mag (E FM mag ) is the magnetic energy at high (low) temperature limit in the Heisenberg model. In the original paper 27 , they assumed the temperature dependence of E mag is determined by the Monte Carlo results only. Therefore, α was treated as a function of temperature (α = α( T )). This is equivalent to that the equilibrium magnetic energy at a temperature is determined to minimise the magnetic free energy, not total free energy. On the other hand, in our study, α is not regarded as a function of temperature but is interpreted as a function of energy (α = α(E)). This interpretation allows that the phonon free energy G ph can be regarded as a function of the magnetic energy (G ph = G ph (T, E mag )). The temperature dependence of E mag is determined after the minimisation of the total free energy in equation (12). This interpretation is the most important key for solving the minimisation problem in the minimum principle for the free energy.
The paramagnetic DLM state in the phonon calculations was mimicked by a special quasirandom structure 36 on the spin configuration (up and down) as obtained from the ATAT package 37 . The atomic forces were calculated by the direct method 24,38 . We used the 3×3×3 cubic supercell (54 atoms) for the force calculations in both ferromagnetic and paramagnetic conditions. The employed lattice constant a = 2.86Å was derived by combining the relaxed lattice constant and experimental lattice expansion ratio at T = 1043 K 39 . Although such determination procedure of lattice constant probably gives some pressure even in the framework of the quasiharmonic approximation, we assume its effect is minor and fixed the volume. First-principles calculations were based on density functional theory within the projector augmented wave method 40 , as implemented in the VASP code 41,42 . For the exchange-correlation functional, the generalised gradient approximation parametrised by Perdew, Burke and Ernzerhof 43 was used. The cutoff energy 400 eV and 9×9×9 k-point grid for the supercell were used for the force calculations. The derivation of force constants and the calculations of the phonon free energy were performed by using the ALAMODE code 44 .
Calculations of exchange coupling constants.
Exchange coupling constants J ij in the Monte Carlo simulations were derived with magnetic force theorem 4 and the Korringa-Kohn-Rostoker (KKR) Green's function method along with the coherent potential approximation (CPA) 45,46 , implemented in the AkaiKKR code 46,47 .
The exchange-correlation functional was treated within the local density approximation 48 . The lattice constant was set to be the same one as in the phonon calculations. Paramagnetic DLM state 3,32 was employed as a reference magnetic state in the derivation of J ij . Calculated J ij values are listed in Table I.

Monte Carlo simulations.
To evaluate the magnetic entropy as a function of the magnetic energy, we carried out the classical Monte Carlo simulations based on the Heisenberg model where J ij denotes the exchange coupling constant and e i is the unit vector on site i. We included up to the third nearest neighbour pairs as interacting shells. The classical Monte Carlo simulations were performed by using the ALPS code 49 . To obtain more accurate results, we employed the rescaled Monte Carlo method 31 which reproduces the quantum specific heat from the classical specific heat. The magnetic energy and entropy were derived by integrating the specific heat. The spin quantum number S = 1.07 for the DLM condition as calculated by KKR-CPA was used in the rescaled Monte Carlo method. The Monte Carlo simulations were carried out using a 16×16×16 sites and involve 300,000 steps for equilibration and 2,700,000 steps for averaging. Temperature grids of 0.1 and 0.2 mRy were used in the range of near T C and other ranges, respectively. Note that the entropy in the rescaled Monte Carlo method does not go to zero at T → 0. Thus this method is not suitable to describe thermodynamic quantities at a low-temperature range. Our thermodynamic formulation, however, needs only the result at a temperature range around T C . Thus the shortcoming does not matter in this study.