Lattice dynamics effects on finite-temperature stability of R_{1-x}Fe_{x} (R = Y, Ce, Nd, Sm, and Dy) alloys from first principles
Guangzong Xing, Takahiro Ishikawa, Yoshio Miura, Takashi Miyake, Terumasa Tadano
LLattice dynamics effects on finite-temperature stability of R − x Fe x ( R = Y, Ce, Nd,Sm, and Dy) alloys from first principles Guangzong Xing,
1, 2, ∗ Takahiro Ishikawa, Yoshio Miura, Takashi Miyake,
3, 2 and Terumasa Tadano
1, 2, † Research center for Magnetic and Spintronic Materials,National Institute for Materials Science, Tsukuba, Ibaraki, 305-0047, Japan Elements Strategy Initiative Center for Magnetic Materials,National Institute for Materials Science, Tsukuba, Ibaraki, 305-0047, Japan Research Center for Computational Design of Advanced Functional Materials,National Institute of Advanced Industrial Science and Technology, Tsukuba, Ibaraki, 305-8568, Japan (Dated: February 5, 2021)We report the effects of lattice dynamics on thermodynamic stability of binary R − x Fe x (0 < x <
1) compounds ( R : rare-earth elements, Y, Ce, Nd, Sm, and Dy) at finite temperature predictedby first-principles calculation based on density functional theory (DFT). We first demonstrate thatthe thermodynamic stability of R − x Fe x (0 < x <
1) alloys cannot be predicted accurately bythe conventional approach, where only the static DFT energy at T = 0 K is used. This issuecan be overcome by considering the entropy contribution, including electronic and vibrational freeenergies, and we obtained convex hull plots at finite temperatures that successfully explain thethermodynamic stability of various known compounds. Our systematic calculation indicates thatvibrational entropy helps stabilize various R − x Fe x compounds with increasing temperature. Inparticular, experimentally reported R Fe compounds are predicted to become thermodynamicallystable above ∼
800 K. We also show that thermodynamic stability is rare-earth dependent and discussits origin. Besides the experimentally reported structures, the stability of two new monoclinic R Fe structures found by Ishikawa et al. [Phys. Rev. Mater. , 104408 (2020)] based on a geneticalgorithm are investigated. These monoclinic phases are found to be dynamically stable and havelarger magnetization than the ThMn -type R Fe . Although they are thermodynamically unstable,the formation energies decrease significantly with increasing temperature, indicating the possibilityof synthesizing these compounds at high temperatures. I. INTRODUCTION
As one of the essential functional materials in the in-dustry, permanent magnets (PM) have widespread appli-cations, such as the electrical motors in vehicles, genera-tors in the wind turbines, and hard disk drives in infor-mation storage [1, 2]. The performance of the PMs canbe measured by the maximum energy product ( BH ) max ,which is mainly governed by two quantities, the satu-ration magnetization ( M s ) and coercivity ( H c ) that isrelated to the magnetocrystalline anisotropy. The chal-lenge here is to find materials with large magnetizationand strong uniaxial magnetocrystalline anisotropy. Al-loys composed of rare-earth (RE) and transition metals(TM) are a class of PMs satisfying these requirementsand possessing the highest performance. The strong mag-netocrystalline anisotropy mainly arises from the strongspin-orbit coupling of RE-4 f orbital shells [3] and the3 d electrons of the TM elements can provide large M s .However, RE elements are not only resource critical andexpensive but also detrimental to M s . As a consequenceof this, TM-rich phases of RE alloys are the best candi-dates for high-performance PMs.Currently, the high-performance PM in commercialuse is the Nd-Fe-B-based ternary alloys, the main phase ∗ [email protected] † [email protected] of which is Nd Fe B [4–7] with the magnetization andmagnetocrystalline anisotropy of 1273 emu/cm and 67kOe [7] at room temperature, respectively. However, theCurie temperature ( T c ) of 585 K as measured by Sagawa et al [8] is marginally above the practical use of the high-temperature edge of 450–475 K. T c and magnetocrys-talline anisotropy field can be improved through the sub-stitution of Nd with Dy [9]. However, Dy is consideredas a less abundant element with a high price on earth.Thus, considerable efforts have been made to design orsearch for RE-lean or RE-free magnets [10–14].In this context, two promising phases, ThMn -type R Fe and R Fe ( R : rare earth) with hexagonal andrhombohedral symmetries, shown in Fig. S1 of the SI [15],have attracted renewed interest, owing to the high con-tent of Fe with respect to the best performance mag-net Nd Fe B. Unfortunately, although the large mag-netization could be realized in the R Fe class of com-pounds, the low T c and plane magnetic crystal anisotropylimit their application as high-performance PMs. Coeyand co-workers reported [16, 17] that small interstitialatoms such as N or C introduced into binary Sm Fe and Y Fe could elevate their T c by factor two and alterthe magnetic anisotropy direction. The major problemhere is R Fe X δ ( X =C, N, B) is thermodynamicallyunstable, making them easy to decompose to more stablebinary phases and bcc Fe at high temperatures [18, 19].Besides the interstitial atoms, substitution of RE or Fe isanother way to improve the performance of the R Fe a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b compounds. Both results from experiments and theoret-ical study prove that the substitutes from the RE site orFe site can stabilize R Fe X δ phase and enhance theirperformance. Chen and coworkers [20] reported thatCr substitute helps in forming the Sm (Fe,Cr) alloyswith a strong uniaxial magnetocrystalline anisotropy.Sm Fe (C,N) could be stabilized through substitutingLa or Ce for Sm without affecting the saturation magne-tization too much [21]. These strategies mentioned abovesuggest that R Fe compounds are attractive candidatesfor permanent magnets.When it comes to ThMn -type R Fe compounds, themajor obstacle with technological application is that bulkphase is difficult to be synthesized because of the ther-modynamic instability, making it easy to decompose intothe stable R Fe phase with lower energy. Suzuki and co-workers [22] reported YFe as a metastable phase ob-tained through the rapid quenching method. This phaseis not favorable for high-performance PM as it is challeng-ing to convert this material into highly textured magnetsfully utilizing the M s . Although bulk binary ThMn -type R Fe compounds do not exist, a thin film phaseof SmFe [23–26] and NdFe [27, 28] has been success-fully synthesized and widely investigated. However, aPM needs to have a certain thickness along its magneticflux line to induce the magnetic flux outside of the mate-rial, which means thin films are not suitable as PMs [29].Bulk phase of R Fe can be stabilized via the substitu-tion of a small amount of Fe with stabilizer elements,such as M = Ti, Si, V, Cr and Mo to form the ternar-ies R Fe − x M x [30–40]. The main drawback of stabilizeris that magnetization decreases fast with the increasedcontent of M , which has a strong site preference and ex-clusively replace only one of three nonequivalent Fe sites.Up to now, the problem of phase stability for ThMn -type R Fe is not yet eliminated, especially the mecha-nism of forming a thin-film phase and metastable bulkphase of YFe [22] are not clear. Formation of a thinfilm or metastable phase means that the driving forcefor phase separation is exceedingly week, and diffusionis kinetically inhibited in the nonequilibrium phase, asdiscussed in the context of other metastable materi-als [41, 42]. This can not be elucidated when usingstatic formation energy obtained by conventional first-principles calculation approach based on density func-tional theory (DFT) at T = 0 K, as discussed inRefs. [43–46], of which the formation energy of R Fe lies highly above the convex hull tie line. Motivatedby this, we report in this paper the lattice dynamics ef-fects on thermodynamic stability of the binary R x Fe − x (0 < x <
1) ( R = Y, Ce, Nd, Sm, and Dy) compoundsevaluated by considering vibrational entropy contribu-tion at finite temperature. By going beyond the conven-tional prediction conducted at T = 0 K, we show that theformation energy in convex hull plot decreases dramat-ically with increasing temperature, which may explainthe existence of known thin-film and metastable phasesas mentioned above. This entropy-driven stabilization of various R − x Fe x compounds can be explained by the x -dependence of phonon frequencies. We also show that theformation energy is RE dependent, and the RE elementshaving smaller ionic radius are more effective in stabiliz-ing the ThMn -type R Fe at finite temperatures.The structure of this paper is organized as follows.In the next section, we describe our theoretical meth-ods in detail. The approximation methods of calculat-ing Helmholtz free energy and computational details aregiven in Secs. II A and II C. Our strategy of system-atically generating the studied alloy compounds is ex-plained in Sec. II B. In Sec. III, we show our main re-sults. Ground state properties of promising ThMn -type R Fe and R Fe compounds are discussed in Sec. III A.The dynamically stable R − x Fe x (0 < x <
1) com-pounds are identified by performing phonon calculationsin Sec. III B. In Sec. III C, the thermodynamic stabil-ity is evaluated by calculating formation energy at finitetemperatures. Two new types of monoclinic R Fe arestudied in Sec. III D. We summarize this study in Sec. IV. II. METHODSA. Helmholtz free energy
To evaluate the finite-temperature phase stability ofthe R − x Fe x alloy compounds in this study, we employthe Helmholtz free energy defined as F ( T ) = E ( c ) + F el ( c , T ) + F vib ( c , T ) , (1)where E ( c ) is the static internal energy of electrons ob-tained by a DFT calculation, F el ( c , T ) is the electronicfree energy, and F vib ( c , T ) is the vibrational free energyat temperature T and cell parameter c , which is a vec-tor comprising six parameters defining the Bravais lat-tice. In general, the thermal expansion of material is fi-nite and thus makes c temperature-dependent, i.e., c ( T ),which is determined as c ( T ) = arg min c F ( c , T ). How-ever, predicting c ( T ) for all of the studied R − x Fe x al-loy compounds is extremely challenging because most ofthe studied alloys display a non-cubic lattice, for whichthe quasiharmonic (QH) phonon calculation needs to beconducted on a multidimensional grid of the lattice pa-rameters. Hence, for all temperatures, we use the latticeparameters at 0 K estimated as c = arg min c E ( c ).The electronic free energy F el ( T ) is obtained based onthe fixed density-of-states (DOS) approximation [47] as F el ( T ) = U el ( T ) − T S el ( T ) , (2) U el ( T ) = (cid:88) k nσ w k (cid:15) k nσ [ f k nσ − θ ( (cid:15) F − (cid:15) k nσ )] , (3) S el ( T ) = − k B (cid:88) k nσ w k [ f k nσ ln f k nσ + (1 − f k nσ ) ln (1 − f k nσ )] . (4)Here, U el and S el are the internal energy and the elec-tronic entropy, respectively, and (cid:15) k nσ represents theKohn–Sham eigenenergy of the n th band with spin σ atmomentum k . f k nσ is the Fermi–Dirac distribution func-tion defined as f k nσ = f ( (cid:15) k nσ , T ) = [ e ( (cid:15) k nσ − µ ) /k B T +1] − with k B and µ being the Boltzmann constant and tem-perature dependent chemical potential, respectively, and θ ( x ) is the step function. The summation over k is per-formed with the irreducible set of k points in the firstBrillouin zone (FBZ) with its weight w k . The secondterm on the right hand side of Eq. (3) subtract E , whichis already included in Eq. (1).The vibrational free energy F vib ( T ) is evaluated basedon the harmonic approximation as F vib ( T ) = 1 N q (cid:88) q ν (cid:26) (cid:126) ω q ν k B T ln (cid:2) − exp ( − (cid:126) ω q ν k B T ) (cid:3)(cid:27) , (5)where (cid:126) is the reduced Planck constant and N q is numberof q points in the FBZ. The phonon frequency ω q ν isobtained by diagonalizing the dynamical matrix D αβκκ (cid:48) ( q ) = 1 √M κ M κ (cid:48) (cid:88) (cid:96) Φ αβ (0 κ ; (cid:96)κ (cid:48) ) e i q · r ( (cid:96) ) , (6)where M κ is the atomic mass of the κ th atom andΦ αβ (0 κ ; (cid:96)κ (cid:48) ) is the second-order interatomic force con-stant between the κ th atom in the unit cell at the centerand the κ (cid:48) th atom in the (cid:96) th unit cell shifted from theorigin by r ( (cid:96) ). B. Studied alloy compounds
The key point of evaluating thermodynamic stabilityis to include as many competing phases as possible in theconvex hull plot. Since the composition space, includingternary and quaternary alloy compounds, is too vast toexplore by ab initio calculations, we restrict our attentionto the binary compounds in this study.To make a comparison between our prediction andexperiment, we systematically generated the competingphase structures as follows; first, experimentally-reportedstructure prototypes of binary R m Fe n ( R = Y, Ce, Nd,Sm, Dy) compounds were obtained from the InorganicCrystal Structure Database. The obtained compositionsare R Fe , R Fe , R Fe , R Fe , R Fe , R Fe , R Fe , R Fe , and R Fe . Two different structures reportedfor R Fe and R Fe were also considered. By contrast,we excluded R Fe and R Fe because their structuresare highly disordered or complicated. Next, for the se-lected prototypes, we substituted R with Y, Ce, Nd, Sm,or Dy. In this way, we generated 45 different materialsin total.In addition to these experimentally-reported struc-tures, we also consider two new monoclinic ( C /m )phases of R Fe that have recently been reported byIshikawa et al. based on a structure-search genetic al-gorithm (GA) in Ref. [48]. While these structures havebeen shown to be less stable than the ThMn -type R Fe at 0 K, they are predicted to have higher Curie temper-ature and larger magnetization than the ThMn -type R Fe , which are favorable as a high-performance PM.Since dynamical stability and finite-temperature thermo-dynamic stability are of vital importance for realizingthese phases in the laboratory, two monoclinic ( C /m )structures, which are labeled as type-I and type-II fol-lowing Ref. [48], are included in our convex hull plot.The complete list of the investigated binary R m Fe n com-pounds is shown in Table S1 of the SI [15]. C. Computational details
All of the DFT calculations in this study were per-formed by using the projector augmented wave (PAW)method [49], as implemented in the Vienna ab initio simulation package (VASP) [50], within the Perdew–Burke–Ernzerhof (PBE) generalized-gradient approxima-tion (GGA) [51]. The 4 f electrons for all the RE elementsare frozen in the core assuming trivalent configuration.Lattice constants and atomic positions of all of the stud-ied R − x Fe x compounds were carefully optimized with akinetic-energy cutoff of 400 eV for the plane-wave expan-sion, and the k -point mesh was generated automaticallyin such a way that the mesh density in the reciprocalspace becomes larger than 450 ˚A − . To test the con-vergence with respect to the k -mesh density, we used amuch denser mesh density of 800 ˚A − ; the difference ofthe total energies was less than 0.1 meV/atom, indicat-ing convergence with the mesh density of 450 ˚A − . Forthe structural optimization and phonon calculations, weused the Methfessel–Paxton smearing method [52] withthe width of 0.2 eV. On the other hand, the tetrahedronmethod with the Bl¨ochl correction [53] was used for cal-culating E .The collinear spin-polarization was applied for all theDFT calculations, including the phonon calculations.The initial local moment of 3 µ B was set for all Fe atomsalong the [001] direction. We tried both the spin momentof RE elements parallel and anti-parallel with that of Featoms with the corresponding initial moment of 0.3 and-0.3 µ B and calculated the total energy of R − x Fe x com-pounds. We obtained lower energy for the anti-parallelcase and used this magnetic configuration for all of thecalculations. The optimized maximum, minimum andaverage moments of different site of Fe atoms in R − x Fe x compounds are listed in Table S1 of the SI [15].The harmonic phonon calculations were conducted bycomputing second-order interatomic force constants us-ing the finite-displacement method, as implemented in phonopy [54]. The size of the supercell for each struc-ture was chosen so that the number of atoms became ∼
100 and larger, which was sufficient to reach the conver-gence of F vib within the error of 1 meV/atom at 1200K. Here we choose the maximum temperature of 1200K as the high-temperature limit to evaluate the entropycontribution to thermodynamic stability. This tempera-ture is slightly above the annealing temperature used byexperiments to obtain the R Fe phases. The adaptedsupercell sizes for all compounds are summarized in Ta-ble S1 of the SI [15]. III. RESULTS AND DISCUSSIONA. Ground state properties
Table I. Calculated lattice constants of R Fe and ThMn -type R Fe phases compared with experimental data at roomtemperature obtained from Refs. [22, 24, 55–59]. m ave is theaverage of the calculated spin moment of Fe atoms in eachalloy.Composition Calc. Expt. m Fe a (˚A) c (˚A) a (˚A) c (˚A) m ave ( µ B ) R Fe - h ( P / mmc )Dy 8.413 8.247 8.453 8.287 2.30Y 8.425 8.245 8.465 8.298 2.31Sm 8.490 8.283 2.36Nd 8.529 8.315 2.38Ce 8.578 8.362 2.41 R Fe - r ( R m )Dy 8.416 12.357 2.30Y 8.428 12.353 2.31Sm 8.517 12.425 8.558 12.447 2.38Nd 8.549 12.477 8.579 12.461 2.41Ce 8.584 12.536 8.496 12.415 2.42 R Fe ( I mmm )Dy 8.415 4.672 2.20Y 8.418 4.674 8.440 a a b b a metastable, Ref. [22] b thin-film, Ref. [24] According to the high content of Fe, R Fe andThMn -type R Fe shown in Fig. S1 of the SI are twopromising compounds for high-performance PMs amongthe compounds studied here. We compare the calcu-lated and experimental lattice constants of these promis-ing compounds in Table I. Here and in what follows, thecalculated results are shown in the ascending order ofthe RE ionic radius, namely, in Dy, Y, Sm, Nd, and Ce,to highlight the RE dependence better. It is reportedby experiments that R Fe phase possesses two differ-ent types of crystal structures: (i) Th Zn -type struc-ture with rhombohedral symmetry ( R m space group)and (ii) Th Ni -type structure with hexagonal symme-try ( P /mmc space group). The structure of R Fe phase is mainly determined by the rare earth elements.Experimentally-reported Dy Fe and Y Fe mainlydisplay the hexagonal structure while the rhombohedralone is found for Sm Fe , Nd Fe , and Ce Fe . Wenote that rhombohedral Dy Fe , Y Fe and hexagonal Ce Fe structures were also reported in Refs. [60–62]based on different temperature treatments.As shown in Table I, the calculated values show rea-sonable agreement with the experimental data. For mostof the studied materials, the calculation slightly underes-timates the experimental lattice constants at room tem-perature, which can be attributed to thermal expansionas all the experimental lattice constants shown in Ta-ble I were measured at room temperature. The latticeconstants of ThMn -type R Fe and R Fe reported inRefs. [43, 44, 63–65] obtained via different DFT codes arealso underestimated to some extent, which are consistentwith our results. The exception for cerium occurs pre-sumably because Ce in the rhombohedral Ce Fe adoptsa valency of 4+ instead of 3+ [44], which is quite distinctfrom the other RE elements. In our study, the valency ofCe in rhombohedral Ce Fe is 3+. Harashima et al. [44]reported the lattice constants of rhombohedral Ce Fe with both 3+ and 4+ valencies. Their calculated latticeconstants of Ce Fe with 4+ valency ( a = 8 .
459 ˚A and c = 12 .
513 ˚A) agree well with the experimental valuesshown in Table I. Moreover, our results with 3+ valencyare consistent with theirs, where the lattice constants arereported as a = 8 .
605 ˚A and c = 12 .
557 ˚A, indicating ourresults are reasonable. A slightly larger discrepancy ob-served for the ThMn -type SmFe can be attributedto an interface effect or presence of TbCu and Th Zn sub-phases in the thin-film sample [24].The highest and lowest magnetic moments ( m Fe ) basedon different site of Fe and volume dependent magnetiza-tion ( M ) for R Fe and R Fe are shown in Table S3 ofthe SI. To quantify m Fe of different compounds, we alsosummarize the weighted average moments m ave basedon all different sites of Fe shown in Table I. It is clearthat m ave of R Fe is larger that that of correspondingThMn -type R Fe compounds. However, according toa much higher concentration of Fe, M of ThMn -type R Fe is superior to that of R Fe compounds. More-over, the different magnetic moment ( m [i]) of nonequiv-alent Fe sites in ThMn -type R Fe , with a trend of m [Fe(8 i )] > m [Fe(8 j )] > m [Fe(8 f )], was observed in ac-cord with the previous calculation [66]. B. Dynamical stability
When it comes to phase stability, the dynamical stabil-ity should be considered first, which can be evaluated bythe phonon dispersion curves. At equilibrium, where thefirst-order force constant equals 0, a dynamically stablecrystal structure means that its potential energy alwaysincreases with any small displacements of atoms. In otherwords, it is equivalent to the condition that no imagi-nary phonon modes ( ω q ν ≥
0) are present in the phonondispersion curves. The phonon dispersion curves for allbinary compounds investigated here are listed in Fig. 1and Figs. S2 and S3 of the SI [15]. As shown in thesefigures, no imaginary phonon modes were detected, indi-
Γ X M Γ Z Z MX PN Γ Γ X M Γ Z Z MX PN Γ Γ X M Γ Z Z MX PN Γ Γ X M Γ Z Z MX PN Γ Γ X M Γ Z Z MX PN Γ
Γ XU K Γ L W X Γ XU K Γ L W X Γ XU K Γ L W X Γ XU K Γ L W X Γ XU K Γ L W X (a) R Fe Dy Y Sm Nd Ce F r e q u e n c y ( m e V ) F r e q u e n c y ( m e V ) (b) R Fe Figure 1. Calculated phonon dispersion curves of ThMn -type R Fe (a) and R Fe (b) compounds. The curves for R Fe and R Fe compounds with different rare earth elements of Dy, Y, Sm, Nd, and Ce are labeled, respectively. Note no unstablebranches are present, which means all the compounds are dynamically stable. cating that all of the studied structures are dynamicallystable at 0 K. As evidently seen in Fig. 1(b), the high-est frequency of the optical phonons gradually decreaseswith going from Dy to Ce, which can be attributed tothe increase in the lattice constants and the associateddecrease in the interatomic force constants. Moreover,YFe has the largest acoustic-phonon frequencies becausethe atomic mass of Y is the lightest among the studiedRE elements. By contrast, the RE-dependence is lessobvious in R Fe compounds (Fig. 1(a)). These signifi-cant or weak RE dependencies of the phonon frequenciescharacterize the RE-dependence of the vibrational freeenergy via Eq. (5) and thereby influence the formationenergy at finite temperature, as will be discussed in thefollowing subsection. C. Thermodynamic stability
The thermodynamic stability of materials can be eval-uated conveniently by creating a convex hull plot. Inmany cases, unary compounds are selected as terminalcompositions of a convex hull plot. However, for thepure RE containing 4 f electrons, an accurate evaluationof E has been challenging due to the strong correla-tion of localized 4 f orbitals, for which the accuracy andconvergence of DFT-GGA calculation are limited, espe-cially when the 4 f bands appear near the Fermi level.The open-core treatment adopted in this study helps im-prove the convergence, but we observed that all binarycompounds R − x Fe x (0 < x <
1) were predicted to beunstable against decomposition into pure RE and bcc Fefor R = Sm and Nd, which does not explain the reality. This issue may be resolved by using the DFT+U method,but an optimal U value would be different for differentphysical properties. To mitigate these technical difficul-ties, we use bcc Fe and R Fe as terminal compositionsin this study. The open-core treatment is appropriate tobinary R − x Fe x compounds since 4 f electrons do not lo-calize near the Fermi level due to the strong hybridizationbetween RE-4 f and Fe-3 d orbitals.Figure 2 shows the convex hull plots calculated for R =Sm and Nd with and without the vibrational free en-ergy. As can be seen in the upper panel, the calculationwithout F vib incorrectly predicts R Fe phases to be un-stable; the hull distance reaches over ∼
14 meV/atom.In addition, several other reported compositions, includ-ing Sm Fe and Nd Fe , are also predicted to be ther-modynamically unstable. These results clearly highlightthe limitation of the conventional approach based on E .By considering the vibrational free energy, the hull dis-tance of these reported phases reduces gradually withincreasing temperature, as shown in the middle and bot-tom panels of Fig. 2. For all studied RE elements, the R Fe phase becomes thermodynamically stable at ∼
800 K, which is reasonably close to the annealing tem-perature used to synthesize the main phase in the labo-ratory. Moreover, the hull distance of the ThMn -type R Fe also decreases as the temperature increases andshows an RE dependence, which will be discussed later.In the convex hull plots shown in Fig. 2 and Fig. S4of the SI, we did not include the electronic free energy, F el , because its composition dependence is expected to besmaller than that of the vibrational free energy. To con-firm this point quantitatively, we compare the calculatedcomposition dependence of F el and F vib in Fig. 3 and Figure 2. Change of relative formation energies per atom with respect to R Fe ( R = Nd and Sm) and bcc Fe for different R − x Fe x (0 < x <
1) compounds obtained without vibrational free energy at 0 K (a) and with vibrational free energy at 500K (b) and 1200 K (c), respectively. The solid line shows the convex hull and the star and circle denote the synthesized phaseby experiments and hypothesized phase, respectively. The stable and unstable phases are shown in blue and red.
Fig. S5 of the SI. Due to the entropic term − T S , the freeenergy decreases monotonically with increasing temper-ature, at least, within the constant DOS approximationor the harmonic approximation. Hence, the compositiondependence of the free energy becomes more significant inthe high-temperature region. For example, we observedthat the maximum difference in F el at 1200 K, whichwas found between bcc Fe and R Fe , was no more than 3meV/atom. By contrast, the maximum difference in thevibrational free energy reaches 160 meV/atom, as shownin Fig. 3(b). These results evidently indicate the crucialimportance of the vibrational free energy and the rela-tively minor effect of F el . The physics behind this starkcontrast can be understood as follows; within the fixedDOS approximation, the temperature-dependence of F el is manifested by the Fermi–Dirac distribution function, f ( (cid:15), T ). The deviation of f ( (cid:15), T ) from f ( (cid:15), f k nσ − θ ( (cid:15) F − (cid:15) k nσ ) in Eq. (3), as well as f k nσ ln f k nσ in Eq. (4) are nonzero only in the range of ± (cid:15) F even at 1500 K. In this narrow energywindow, the DOS is dominated by the Fe 3 d orbitals and DOS/atom does not depend much on the fraction of Fe, x , leading to the weaker composition dependence of F el per atom (Fig. 3(a)). By contrast, since the energy scaleof phonons is less than ∼
50 meV in the studied alloys, allphonon modes contribute to F vib significantly even below1000 K. It can be shown easily that F vib is an increasingfunction of ω q ν ; thus, the negative value of F vib /atom be-comes larger in its magnitude with decreasing the phononfrequency. By considering that the phonon frequency isroughly proportional to M − / κ (Eq. (6)) and assumingthat the atomic mass can be replaced with the averagemass M avg ( x ) = (1 − x ) M R + x M Fe , it is easy to showthat phonon frequency tends to increase (thereby F vib tends to decrease) with increasing x because M R > M Fe is satisfied for all the studied RE elements. This x depen-dence of F vib /atom explains why the largest differencewas observed between R Fe ( x = ) and Fe ( x = 1).Interestingly, if we follow the above simple argumentbased on M avg ( x ) , it is easy to show that F vib /atom is aconvex function of x . Indeed, we observed that the calcu-lated data of F vib /atom roughly follows the convex curve, -20-15-10-5 0 F e l ( m e V /a t o m ) (a) -600-500-400-300-200-100 0 100 0 200 400 600 800 1000 1200 F v i b ( m e V /a t o m ) Temperature (K)
SmFe NdFe SmFe NdFe Fe (b) Figure 3. Calculated electronic free energy F el (a) and vi-brational free energy F vib (b) with respect to temperature ofThMn -type R Fe , R Fe ( R = Sm and Nd) phases and bccFe. while some deviation from the curve was observed likelydue to the x dependence of interatomic force constants.The approximately convex shape of the F vib /atom curveexplains qualitatively why many compounds are stabi-lized in the convex hull plot (Fig. 2) at finite tempera-tures.Next, we discuss the RE dependence of the forma-tion energy, particularly focusing on the Fe-rich phases, R Fe and R Fe . As shown in Fig. 2 and Fig. S4 of theSI, the convex hull plots for the five different RE looksomewhat similar to each other. Still, after a detailedinvestigation, we found a clear RE dependence. For ex-ample, the formation energy of R Fe defined as∆ F | R Fe ← F [ R Fe ]+13 F [Fe] = F [ R Fe ] − (2 F [ R Fe ] + 13 F [Fe])= ∆ E + ∆ F vib + ∆ F el (7)is shown in Fig. 4 as a function of r calc R , which is theradius of the RE element estimated as the half of theshortest bond length in the pure hexagonal RE phases.As already discussed in Sec. III A, R Fe displays twodifferent structures, namely, hexagonal or rhombohedralphase, depending on the RE element. Hence, the for- r calc R (˚A) −400−2000200400 ∆ F | R F e ← F [ R F e ] + F [ F e ] ( m e V ) DFT 0Kw/ph+ele 700Kw/ph+ele 1200K
Dy Y Sm Nd Ce
Figure 4. Formation energy defined for R Fe ( R = Dy,Y, Sm, Nd, and Ce) in Eq. (7) as a function of atomic ra-dius at different temperatures . The filled and open symbolsdenote the formation energies of R Fe with hexagonal andrhombohedral structure, respectively. r calc R (˚A) −50050100150 ∆ F | R F e ← F [ R F e ] + F [ F e ] ( m e V ) DFT 0Kw/ph+ele 700Kw/ph+ele 1200K
Dy Y Sm Nd Ce
Figure 5. Formation energy defined for ThMn -type R Fe ( R = Dy, Y, Sm, Nd, and Ce) in Eq. (8) as a function ofatomic radius at different temperatures. The filled and opensymbols denote the formation energies related to R Fe withhexagonal and rhombohedral structure, respectively. mation energy [Eq. (7)] was calculated for the hexag-onal and rhombohedral phases and shown in Fig. 4by filled and open symbols, respectively. Since thefree energy difference between the two structures wassmall, we considered F el in addition to F vib . Whenwe considered E alone, the formation energy ∆ E be-came positive for all the studied RE elements, and ittends to increase with increasing r calc R , as shown by thesquare symbols in Fig. 4. With increasing temperature,∆ F | R Fe ← F [ R Fe ]+13 F [Fe] changes the sign, and the REdependence becomes weaker, which is most notable at1200 K. The weaker RE dependence at a high tempera-ture can be understood as follows; first, the temperature-dependent part of Eq. (7), i.e., ∆ F vib ( T ) + ∆ F el ( T ), de-creases with increasing r calc R particularly for R = Sm, Nd,and Ce. This tendency was commonly observed in bothhexagonal and rhombohedral phases. Since the r calc R de-pendence of ∆ F vib ( T ) + ∆ F el ( T ) is opposite to that of∆ E , the sum of these two terms shows a weaker r calc R dependence. Second, we observed for all the studied REthat the free energy difference between the hexagonaland rhombohedral structures, F [ R Fe - h ] − F [ R Fe - r ],keeps increasing with increasing the temperature, whichcan be barely inferred from Fig. 4 and is more clearlyshown in Fig. S6 of the SI. For Dy Fe and Y Fe ,the hexagonal phase is more stable than the rhombohe-dral one in the low-temperature range, while the rhom-bohedral phase is more stable for the other R Fe com-pounds. With increasing the temperature, the rhombo-hedral phase acquires more energy gain than the hexag-onal phase, thus contributing to the overall weak RE de-pendence of ∆ F | R Fe ← F [ R Fe ]+13 F [Fe] at 1200 K.Finally, we discuss the thermodynamic stability of theThMn -type R Fe at finite temperature. A much morecompetitive reference phase comparing to R Fe for eval-uating formation energy of ThMn -type R Fe com-pounds is R Fe , since all compounds of this phase stud-ied here are realized by experiments and have a relativelylower formation energy than R Fe phase. Hence, we eval-uated the formation energy (meV/f.u.) by∆ F | R Fe ← F [ R Fe ]+ F [Fe] = F [ R Fe ] − (cid:18) F [ R Fe ] + 72 F [Fe] (cid:19) (8)for different ThMn -type R Fe compounds with respectto the hexagonal and rhombohedral R Fe phase andbcc Fe. The clear RE dependent trend is shown in Fig. 5;the smaller the atomic radius is, the more stable theThMn -type R Fe tends to be. Our calculation pre-dicts DyFe and YFe to be thermodynamically stableabove ∼
700 K and SmFe to be nearly stable when thetemperature reaches 1200 K. These results are consistentwith the fact that only YFe has been synthesized asbulk by the high-temperature annealing followed by rapidquenching [22]. Recently, Harashima et al. [43] reportedbased on DFT calculations that the hydrostatic pressureof ∼ andSmFe by ∼
50 meV/f.u. and ∼
25 meV/f.u., respec-tively. They also reported that the pressure did notreduce the formation energy of DyFe . On the otherhand, the entropic stabilization at 1200 K shown in Fig. 5amounts to 70 meV/f.u., 60 meV/f.u., and 35 meV/f.u.for DyFe , SmFe , and NdFe , respectively. These re-sults indicate that the temperature is more effective instabilizing the ThMn -type R Fe phase than pressure,particularly when the atomic radius is small. In addition,our result indicates that the finite-temperature stabilityof NdFe and SmFe may be enhanced by partially sub-stituting (Nd, Sm) with (Y, Dy). D. Monoclinic R Fe Recently, Ishikawa et al. reported two new metastablephases of YFe [48] named Type-I and Type-II withmonoclinic C m structures using the structure searchmethod based on a scheme of GA [67]. These newstructures were reported to possess larger magnetiza-tion M and higher Curie temperature T c than corre-sponding ThMn -type YFe , which is attractive as ahigh-performance PM material if the structure can beformed as a (meta-)stable phase. Knowing the details ofthe structure information, we generated four analogs of R Fe ( R = Dy, Sm, Nd, and Ce) and investigated theirphase stability and magnetic properties. Figure 6(a)and (b) shows the schematic crystal structures of theType-I and Type-II monoclinic R Fe compounds, whosestructure parameters are listed in Table S2 of the SI.Compared with YFe , the optimized lattice parametersand angles of R Fe analogs change slightly. We alsoperformed phonon calculations of these monoclinic com-pounds. As shown in Figs. 6(c),(d) and Fig. S7 of the SI,all of these compounds are dynamically stable withoutexhibiting any unstable modes. The phonon frequency isless RE dependent than that of R Fe phase.Thermodynamic stability of these monoclinic struc-tures at finite temperatures are evaluated and shownin Fig. 2 and Fig. S4 of the SI. We observe that allthese monoclinic R Fe compounds are thermodynamicunstable. Without the vibrational free energy contri-bution, the formation energies exceed more than 30meV/atom above the R Fe -Fe tie line. Fortunately,they decrease dramatically as increasing temperature andbecome less than 15 meV/atom above tie line for alltype-I R Fe compounds at 1200 K. Although the type-I monoclinic R Fe compounds are thermodynamicallyunstable, these dynamically stable phases are, at least,metastable. Thus, it is still possible to synthesize themonoclinic phase by experiments through applying a con-jugated field (pressure, temperature, or surface area) asdiscussed in Ref. [68]. The magnetization M of new type R Fe compounds is shown in Table S3. As seen, thepredicted Type-I phases of R Fe possess the largest M ,which is consistent with that of Type-I YFe compoundpredicted by Ishikawa et.al. [48]. Given the lack of f elec-trons in YFe and high cost of DyFe , the new struc-tures of SmFe and NdFe would be good candidates asa PM whose performance is superior to that of the corre-sponding ThMn -type compound if the phase instabilityproblems can be solved, for example, by chemical substi-tution. IV. SUMMARY
To summarize, we investigated the dynamical and ther-modynamic stability of R − x Fe x ( R = Y, Ce, Nd, Sm,and Dy) compounds by using first-principles calculationsbased on DFT with the effect of the vibrational and elec- Figure 6. The low-energy R Fe compounds: the lowest energy structure, type-I (a) and the second lowest energy structure,type-II (b) as obtained from the GA search. (c),(d) Calculated phonon dispersion curves for type-I and Type-II types of Sm-and NdFe , respectively. Note the predicted structures are dynamically stable. tronic entropies. By performing phonon calculations sys-tematically, we showed that all compounds studied hereare dynamically stable. We also demonstrated that theinclusion of the vibrational entropy significantly improvesthe prediction accuracy of the thermodynamic stabilityat a finite temperature compared to the conventional ap-proach based on the static DFT energy.The ThMn -type R Fe compounds, which arepromising for PM applications, were thermodynamicallyunstable because of the presence of competing R Fe and R Fe phases. Nevertheless, the formation energies of R Fe decreased significantly with increasing tempera-ture, and the ThMn -type R Fe phases were predictedto become stable in the high-temperature range, partic-ularly for DyFe and YFe that have relatively smallerlattice constants than the other three RE cases. Weshowed that the observed stabilization and its RE de-pendence could be explained by the difference in thephonon frequencies and thereby in the vibrational freeenergy. Moreover, two new metastable phases of mon-oclinic R Fe predicted using a scheme of genetic algo-rithms were included in this study. These monoclinicphases showed larger magnetization, which is superiorto that of corresponding ThMn -type R Fe . Althoughthey were found to be thermodynamically unstable, the formation energies decreased dramatically with heatingand became less than 15 meV/atom for the Type-I struc-ture at 1200 K.While our prediction that incorporates the vibrationaland electronic entropies should be, in principle, moreaccurate than the conventional approach based on thestatic energy, it was still difficult to reach perfect agree-ments between theory and experiment. Since some otherfactors, including magnetic and mixing entropies, andlattice anharmonicity, are still missing in the present cal-culation, we expect the inclusion of these factors will im-prove the prediction accuracy even further, which is leftfor a future study. ACKNOWLEDGMENTS
This work was partially supported by the Ministryof Education, Culture, Sports, Science and Technol-ogy (MEXT) as the Elements Strategy Initiative Cen-ter for Magnetic Materials (ESICMM), Grant NumberJPMXP0112101004 and as “Program for Promoting Re-searches on the Supercomputer Fugaku” (DPMSD).0 [1] J. M. Coey,
Magnetism and magnetic materials (Cam-bridge university press, 2010).[2] L. H. Lewis and F. Jim´enez-Villacorta, Perspectives onpermanent magnetic materials for energy conversion andpower generation, Metall. Mater. Trans. A , 2 (2013).[3] T. Miyake and H. Akai, Quantum theory of rare-earthmagnets, J. Phys. Soc. Jpn. , 041009 (2018).[4] M. Sagawa, S. Fujimura, N. Togawa, H. Yamamoto, andY. Matsuura, New material for permanent magnets on abase of Nd and Fe, J. Appl. Phys. , 2083 (1984).[5] J. Herbst, R Fe B materials: Intrinsic properties andtechnological aspects, Rev. Mod. Phys. , 819 (1991).[6] J. J. Croat, J. F. Herbst, R. W. Lee, and F. E. Pinkerton,Pr-Fe and Nd-Fe-based materials: A new class of high-performance permanent magnets, J. Appl. Phys. , 2078(1984).[7] S. Hirosawa, Y. Matsuura, H. Yamamoto, S. Fujimura,M. Sagawa, and H. Yamauchi, Magnetization and mag-netic anisotropy of R Fe B was measured on single crys-tals, J. Appl. Phys. , 873 (1986).[8] M. Sagawa, S. Hirosawa, H. Yamamoto, S. Fujimura,and Y. Matsuura, Nd–Fe–B permanent magnet materi-als, Jpn. J. Appl. Phys. , 785 (1987).[9] L. Li, J. Yi, Y. Peng, and B. Huang, The effect of com-pound addition Dy O and Sn on the structure and prop-erties of NdFeNbB magnets, J. Magn. Magn. Mater. ,80 (2007).[10] K. Skokov and O. Gutfleisch, Heavy rare earth free, freerare earth and rare earth free magnets-vision and reality,Scripta Mater. , 289 (2018).[11] I. Poenaru, A. Lixandru, S. Riegg, B. Fayyazi, A. Taubel,K. G¨uth, R. Gauß, and O. Gutfleisch, Ce and La as sub-stitutes for Nd in Nd Fe B-based melt-spun alloys andhot-deformed magnets: a comparison of structural andmagnetic properties, J. Magn. Magn. Mater. , 198(2019).[12] R. McCallum, L. H. Lewis, R. Skomski, M. Kramer, andI. Anderson, Practical aspects of modern and future per-manent magnets, Annu. Rev. Mater. Res. , 451 (2014).[13] M. KuzMin, K. Skokov, H. Jian, I. Radulov, and O. Gut-fleisch, Towards high-performance permanent magnetswithout rare earths, J. Phys.: Condens. Matter ,064205 (2014).[14] G. Trumpy, E. Both, C. Djega-Mariadassou, andP. Lecocq, Mssbauer-effect studies of iron-tin alloys,Phys. Rev. B , 3477 (1970).[15] See Supplemental Material at [URL] for crystal structureinformation of R Fe , ThMn -type, and new mono-clinic R Fe phases; supercell sizes chosen for phonon cal-culations and magnetic information of R − x Fe x phases;phonon dispersion curves and convex hull plots of therest R − x Fe x phases besides the main text.[16] J. Coey and H. Sun, Improved magnetic properties bytreatment of iron-based rare earth intermetallic com-pounds in anmonia, J. Magn. Magn. Mater. , L251(1990).[17] H. Sun, J. Coey, Y. Otani, and D. Hurley, Magnetic prop-erties of a new series of rare-earth iron nitrides: R Fe N y (y approximately 2.6), J. Phys.: Condens. Matter , 6465(1990).[18] X. Kou, R. Gr¨ossinger, M. Katter, J. Wecker, L. Schultz, T. Jacobs, and K. Buschow, Intrinsic magnetic propertiesof R Fe C y N x compounds:(R= Y, Sm, Er, and Tm),J. Appl. Phys. , 2272 (1991).[19] G. Sadullaho˘glu and B. Altuncevahir, Influence of boronaddition on magnetic properties of Sm Fe (2019).[20] Z. Chen and G. Hadjipanayis, Effects of Cr substitu-tion on the formation, structure and magnetic propertiesof Sm (Fe,Cr) C x alloys, IEEE Trans. Magn. , 3856(1997).[21] T. Pandey, M.-H. Du, and D. S. Parker, Tuning the mag-netic properties and structural stabilities of the 2-17-3magnets Sm Fe X (X= C, N) by substituting La or Cefor Sm, Phys. Rev. Appl. , 034002 (2018).[22] H. Suzuki, Metastable phase YFe fabricated by rapidquenching method, AIP Adv. , 056208 (2017).[23] F. Cadieu, H. Hegde, A. Navarathna, R. Rani, andK. Chen, High-energy product ThMn Sm-Fe-T and Sm-Fe permanent magnets synthesized as oriented sputteredfilms, Appl. Phys. Lett. , 875 (1991).[24] D. Wang, S.-H. Liou, P. He, D. J. Sellmyer, G. Had-jipanayis, and Y. Zhang, SmFe and SmFe N x filmsfabricated by sputtering, J. Magn. Magn. Mater. , 62(1993).[25] H. Sepehri-Amin, Y. Tamazawa, M. Kambayashi,G. Saito, Y. Takahashi, D. Ogawa, T. Ohkubo, S. Hi-rosawa, M. Doi, T. Shima, et al. , Achievement of highcoercivity in Sm(Fe . Co . ) anisotropic magnetic thinfilm by boron doping, Acta Mater. , 337 (2020).[26] Y. Hirayama, Y. Takahashi, S. Hirosawa, and K. Hono,Intrinsic hard magnetic properties of Sm(Fe − x Co x ) compound with the ThMn structure, Scripta Mater. , 62 (2017).[27] Y. Hirayama, Y. Takahashi, S. Hirosawa, and K. Hono,NdFe N x hard-magnetic compound with high magneti-zation and anisotropy field, Scripta Mater. , 70 (2015).[28] T. Sato, T. Ohsuna, M. Yano, A. Kato, and Y. Kaneko,Permanent magnetic properties of NdFe N x sputteredfilms epitaxially grown on V buffer layer, J. Appl. Phys. , 053903 (2017).[29] S. Hirosawa, M. Nishino, and S. Miyashita, Perspectivesfor high-performance permanent magnets: applications,coercivity, and new materials, Adv. Nat. Sci.: Nanosci.Nanotechnol. , 013002 (2017).[30] B.-P. Hu, H.-S. Li, and J. Coey, Relationship betweenThMn and Th Ni structure types in the YFe − x Ti x alloy series, J. Appl. Phys. , 4838 (1990).[31] F. De Boer, Y.-K. Huang, D. De Mooij, and K. Buschow,Magnetic properties of a series of novel ternary inter-metallics (RFe V ), J. Less-Common Met. , 199(1987).[32] I. Tereshina, P. Gaczy´nski, V. Rusakov, H. Drulis,S. Nikitin, W. Suski, N. Tristan, and T. Palewski, Mag-netic anisotropy and m¨ossbauer effect studies of YFe Tiand YFe TiH, J. Phys.: Condens. Matter , 8161(2001).[33] O. Isnard, S. Miraglia, M. Guillot, and D. Fruchart, Hy-drogen effects on the magnetic properties of RFe Ticompounds, J. Alloy. Compd. , 637 (1998).[34] L. Ke and D. D. Johnson, Intrinsic magnetic propertiesin R(Fe − x Co x ) TiZ (R= Y and Ce; Z= H, C, and N), Phys. Rev. B , 024423 (2016).[35] C. Skelland, T. Ostler, S. Westmoreland, R. Evans,R. Chantrell, M. Yano, T. Shoji, A. Manabe, A. Kato,M. Ito, et al. , Probability distribution of substituted ti-tanium in RT (R= Nd and Sm; T= Fe and Co) struc-tures, Ieee T. Magn. , 1 (2018).[36] Y. Harashima, K. Terakura, H. Kino, S. Ishibashi, andT. Miyake, First-principles study of structural and mag-netic properties of R (Fe,Ti) and R(Fe,Ti) N (R= Nd,Sm, Y), in
Proceedings of Computational Science Work-shop 2014 (CSW2014) (2015) p. 011021.[37] Y.-c. Yang, X.-d. Zhang, L.-s. Kong, Q. Pan, and S.-l. Ge,New potential hard magnetic material—NdTiFe N x ,Solid State Commun. , 317 (1991).[38] Y.-c. Yang, X.-d. Zhang, S.-l. Ge, Q. Pan, L.-s. Kong,H. Li, J.-l. Yang, B.-s. Zhang, Y.-f. Ding, and C.-t. Ye,Magnetic and crystallographic properties of novel Fe-richrare-earth nitrides of the type RTiFe N − δ , J. Appl.Phys. , 6001 (1991).[39] K. Ohashi, Y. Tawara, R. Osugi, and M. Shimao, Mag-netic properties of Fe-rich rare-earth intermetallic com-pounds with a ThMn structure, J. Appl. Phys. , 5714(1988).[40] A. Sch¨onh¨obel, R. Madugundo, O. Y. Vekilova, O. Eriks-son, H. C. Herper, J. Barandiar´an, and G. Hadjipanayis,Intrinsic magnetic properties of SmFe − x V x alloys withreduced V-concentration, J. Alloy. Compd. , 969(2019).[41] H. Sellinschegg, S. L. Stuckmeyer, M. D. Hornbostel, andD. C. Johnson, Synthesis of metastable post-transition-metal iron antimony skutterudites using the multilayerprecursor method, Chem. Mater. , 1096 (1998).[42] G. Xing, Y. Li, X. Fan, L. Zhang, W. Zheng, and D. J.Singh, Sn Se : A conducting crystalline mixed valentphase change memory compound, J. Appl. Phys. ,225106 (2017).[43] Y. Harashima, T. Fukazawa, H. Kino, and T. Miyake,Effect of R-site substitution and the pressure on stabilityof RFe : A first-principles study, J. Appl. Phys. ,163902 (2018).[44] Y. Harashima, T. Fukazawa, and T. Miyake, Ceriumas a possible stabilizer of ThMn -type iron-based com-pounds: A first-principles study, Scripta Mater. , 12(2020).[45] H. ˙I. S¨ozen, S. Ener, F. Maccari, K. P. Skokov, O. Gut-fleisch, F. K¨ormann, J. Neugebauer, and T. Hickel, Abinitio phase stabilities of Ce-based hard magnetic mate-rials and comparison with experimental phase diagrams,Phys. Rev. Mater. , 084407 (2019).[46] M. Matsumoto, T. Hawai, and K. Ono, (Sm,Zr)Fe − x M x (M= Zr, Ti, Co) for permanent-magnet applications: Abinitio material design integrated with experimental char-acterization, Phys. Rev. Appl. , 064028 (2020).[47] X. Zhang, B. Grabowski, F. K¨ormann, C. Freysoldt, andJ. Neugebauer, Accurate electronic free energies of the3d, 4d, and 5d transition metals at high temperatures,Phys. Rev. B , 165126 (2017).[48] T. Ishikawa, T. Fukazawa, and T. Miyake, MonoclinicYFe phases predicted from first principles, Phys. Rev.Mater. , 104408 (2020).[49] G. Kresse and D. Joubert, From ultrasoft pseudopoten-tials to the projector augmented-wave method, Phys.Rev. B , 1758 (1999). [50] G. Kresse and J. Furthm¨uller, Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set, Phys. Rev. B , 11169 (1996).[51] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalizedgradient approximation made simple, Phys. Rev. Lett. , 3865 (1996).[52] M. Methfessel and A. T. Paxton, High-precision samplingfor brillouin-zone integration in metals, Phys. Rev. B ,3616 (1989).[53] P. E. Bl¨ochl, O. Jepsen, and O. K. Andersen, Improvedtetrahedron method for brillouin-zone integrations, Phys.Rev. B , 16223 (1994).[54] A. Togo and I. Tanaka, First principles phonon calcula-tions in materials science, Scripta Mater. , 1 (2015).[55] A. Van der Goot and K. Buschow, The dysprosium-iron system: Structural and magnetic properties ofdysprosium-iron compounds, J. Less-Common Met. ,151 (1970).[56] T. Jacobs, K. Buschow, G. Zhou, X. Li, and F. De Boer,Magnetic interactions in R Fe − x Al x compounds (R=Ho, Y), J. Magn. Magn. Mater. , 220 (1992).[57] A. Teresiak, O. Gutfleisch, N. Mattern, B. Gebel, andK.-H. M¨uller, Phase formation and crystal structure ofSm Fe − y Ga y compounds during hydrogen dispropor-tionation and desorption recombination (hddr-process),J. Alloy. Compd. , 235 (2002).[58] N. Bouchaala, M. Jemmali, T. Bartoli, K. Nouri, I. Hen-tech, S. Walha, L. Bessais, and A. B. Salah, Influenceof Fe-substitution on structural, magnetic and magne-tocaloric properties of Nd Fe − x Co x solid solutions, J.Solid State Chem. , 501 (2018).[59] P. Teplykh, A. Pirogov, A. Kuchin, and A. Teplykh, Realcrystal structure and magnetic state of Ce Fe com-pounds, Physica B Condens. Matter , E99 (2004).[60] P. ´Alvarez-Alonso, P. Gorria, J. L. S. Llamazares, G. J.Cuello, I. P. Orench, J. S. Marcos, G. Garbarino, M. Reif-fers, and J. A. Blanco, Exploring the magneto-volumeanomalies in Dy Fe with unconventional rhombohedralcrystal structure, Acta Mater. , 7931 (2013).[61] K. Buschow, The crystal structures of the rare-earth com-pounds of the form R Ni , R Co and R Fe , J. Less-Common Met. , 204 (1966).[62] K. Buschow and J. Van Wieringen, Crystal structure andmagnetic properties of cerium-iron compounds, Phys.Stat. Sol. , 231 (1970).[63] D. Odkhuu, T. Ochirkhuyag, and S. C. Hong, Enhanc-ing energy product and thermal stability of SmFe byinterstitial doping, Phys. Rev. Appl. , 054076 (2020).[64] T. Fukazawa, H. Akai, Y. Harashima, and T. Miyake,First-principles study of intersite magnetic couplings inNdFe and NdFe X (X= B, C, N, O, F), J. Appl. Phys. , 053901 (2017).[65] P. Delange, S. Biermann, T. Miyake, and L. Pourovskii,Crystal-field splittings in rare-earth-based hard magnets:An ab initio approach, Phys. Rev. B , 155132 (2017).[66] T. Miyake, K. Terakura, Y. Harashima, H. Kino, andS. Ishibashi, First-principles study of magnetocrystallineanisotropy and magnetization in NdFe , NdFe Ti, andNdFe TiN, J. Phys. Soc. Jpn. , 043702 (2014).[67] T. Ishikawa, T. Miyake, and K. Shimizu, Materials infor-matics based on evolutionary algorithms: Application tosearch for superconducting hydrogen compounds, Phys.Rev. B , 174506 (2019). [68] W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain,W. D. Richards, A. C. Gamst, K. A. Persson, and G. Ceder, The thermodynamic scale of inorganic crys-talline metastability, Sci. Adv.2