Simultaneous tuning of the magnetic anisotropy and thermal stability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ^{''}$$\end{document}α′′-phase Fe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}$$\end{document}16N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2

Simultaneously enhancing the uniaxial magnetic anisotropy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u$$\end{document}Ku) and thermal stability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ^{''}$$\end{document}α′′-phase Fe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}$$\end{document}16N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2 without inclusion of heavy-metal or rare-earth (RE) elements has been a challenge over the years. Herein, through first-principles calculations and rigid-band analysis, significant enhancement of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u$$\end{document}Ku is proposed to be achievable through excess valence electrons in the Fe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}$$\end{document}16N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2 unit cell. We demonstrate a persistent increase in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u$$\end{document}Ku up to 1.8 MJ m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text {-}3}$$\end{document}-3, a value three times that of 0.6 MJ m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text {-}3}$$\end{document}-3 in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ^{''}$$\end{document}α′′-Fe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}$$\end{document}16N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2, by simply replacing Fe with metal elements with more valence electrons (Co to Ga in the periodic table). A similar rigid-band argument is further adopted to reveal an extremely large \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u$$\end{document}Ku up to 2.4 MJ m\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\text {-}3}$$\end{document}-3 in (Fe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{0.5}$$\end{document}0.5Co\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{0.5}$$\end{document}0.5)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{16}$$\end{document}16N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_{2}$$\end{document}2 obtained by replacing Co with Ni to Ga. Such a strong \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u$$\end{document}Ku can also be achieved with the replacement by Al, which is isoelectronic to Ga, with simultaneous improvement of the phase stability. These results provide an instructive guideline for simultaneous manipulation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_u$$\end{document}Ku and the thermal stability in 3d-only metals for RE-free permanent magnet applications.

Alpha-phase iron has been known for its extraordinary magnetic properties, including high saturation magnetization ( µ 0 M s ) and Curie temperature ( T c ), in addition to its relatively simple fabrication and low price. These intriguing features make it a potential champion ever for high-performance permanent magnet applications [1][2][3][4] . However, the main drawback of α-Fe is its negligible uniaxial magnetic anisotropy ( K u ) on the order of µ eV per atom [1][2][3][4] . Two practical approaches to enhance K u in α-Fe are (1) alloying with heavy-metal (HM) or rare-earth (RE) elements, with the most prominent examples being FePt 5 and Nd 2 Fe 14 B 6, 7 , and (2) reducing the crystal symmetry from the cubic ( c/a = 1 ) to tetragonal phase ( c/a = 1) [8][9][10] . In (1), 5d or 4f electrons possess inherently large spin-orbit coupling (SOC) and orbital angular momentum (L), but the inclusion of these HM and RE elements is not desirable in terms of price and is detrimental to µ 0 M s and T c . In (2), the energy levels of the 3d orbitals evolve in tetragonal symmetry, particularly around the Fermi level ( E F ), which in turn enhances K u 11 .
The tetragonal phase of c/a = 1 is now accessible in epitaxial Fe films with a diverse choice of lattice-mismatched substrates. Nevertheless, such tetragonal distortion is feasible only for limited film thicknesses of a few nanometers 12,13 . In contrast, a bulk-scale tetragonal structure of the α ′′ -phase (c/a = 1.1) is favored when 12.5 at.% N is embedded into the α-Fe structure with octahedral interstitial sites, forming a 16:2 (Fe 16 N 2 ) stoichiometry 14 . Since a surprisingly large magnetic moment of 2.6−3 µ B per Fe atom was reported [15][16][17] , α ′′ -phase Fe 16 N 2 has received enormous attention as a possible 3d-only permanent magnet. However, the practical implementation of α ′′ -Fe 16 N 2 in obtaining monophasic samples is quite difficult as α ′′ -Fe 16 N 2 decomposes into α-Fe and γ ′ -Fe 4 N at a low temperature of approximately 500 K 14 . Numerous efforts have been made to improve the thermal stability of α ′′ -Fe 16 N 2 ; the most successful approach is Ti addition but the magnetic properties are greatly suppressed 18,19 . In addition to the weak thermal stability, another major obstacle that hampers practical applications is the still insufficient K u , which ranges from 0.4 to 1 MJ m -3 , depending on the sample preparation and film thickness [20][21][22] . In the research community, search for enhancing K u while improving the thermal stability of α ′′ -Fe 16 N 2 in the bulk has been thus very intensive and remains unresolved.
In this article, we propose a possible mechanism of tuning the number of valence electrons to simultaneously enhance the thermal stability and K u by a few times in Fe 16 N 2 and (Fe 0.5 Co 0.5 ) 16 N 2 apart from the aforementioned approaches (1) and (2), using first-principles calculations and rigid-band model analysis. We predict a persistent increase in K u up to 2.4 MJ m -3 , which is four times that (0.6 MJ m -3 ) of Fe 16 N 2 , by replacing Fe with metal elements with more valence electrons (Co to Ga and Al in the periodic table). Such a supreme K u is discussed in  14 . In the unit cell, 16 Fe atoms occupy 3 inequivalent sites at the Wyckoff positions of 4e, 8h, and 4d, while 2 N are at the octahedral interstices with the 4 Fe(8h) coordination 14 . These 4e, 8h, and 4d sites differ in magnetic moment (Table 1): 2.14 (2.33), 2.36 (2.45), and 2.82 (3.05) µ B in the present theory (experiment 17 ), respectively. The mechanism is associated with the nonidentical Fe-N bond lengths: Fe(4e)-N: 1.83, Fe(8h)-N: 1.95, and Fe(4d)-N: 3.24 Å. Furthermore, the energy levels of their 3d orbitals evolve near E F , particularly in the spin-down channel (Fig. 1b).     The obtained H f values of Fe 16 N 2 are -3.14 kJmol -1 against ( α-Fe)+N 2 and -0.80 kJ mol -1 against ( α-Fe)+(γ ′ -Fe 4 N) decomposition. The small negative value of the latter implies that the α ′′ -Fe 16 N 2 phase is stable at a low temperature but most likely decomposes into the α-Fe and γ ′ -Fe 4 N phases at an elevated temperature, as observed experimentally 14 . From the ternary Fe−M−N phase diagram, as an example for M = Al in Fig. 1c, Fe 4 N and Fe 3 M phases are identified as the most competitive binary decomposable phases to Fe 15 M 1 N 2 . For M = Ni (Cu and Zn), Fe 4 N+FeNi(Cu/Zn)+Fe decomposition has been considered since Fe 3 Ni (Fe 3 Cu/Zn and FeCu/Zn) is unstable. Figure 1d presents the calculated H f of Fe 15 M 1 N 2 (M = Co−Ga and Al) for M(4e), M(8h), and M(4d). All the M elements prefer 4d-site replacement, which in turn splits the neighboring Fe(4d) sites into Fe(4d) c along the c axis and Fe(4d) ab on the ab plane with dissimilar magnetic properties, as addressed in the following paragraphs. The α ′′ -phase becomes unstable upon Co to Cu replacements. In contrast, the replacement of Zn, Ga, and Al improves the α ′′ -phase stability with H f enhanced by 0.04−0.8 kJ mol -1 in magnitude, as their nitrides (ZnN, GaN, and AlN) have higher standard enthalpies of formation, in the range of -100 to -320 kJ mol -123-25 , than FeN (-47 kJ mol -1 ) 25 . The completely filled d-orbitals of the Zn and Ga elements provide extra stability to the system, as these elements have a symmetrical distribution of electrons and larger exchange energies than Fe 26 .

Results and discussion
We further investigate the structural stability at an elevated temperature using ab initio molecular dynamic (AIMD) simulation. Figure 2a,b present the fluctuations of the total free energy of the selected Fe 16 N 2 and Fe 15 Al 1 N 2 phases for given temperatures 300, 500, and 600 K, respectively. The total energy of Fe 16 N 2 decreases immediately within a few fs by 1.3 (at 300 K)−2.2 eV/f.u. (at 600 K), which is associated with the thermal motion and relocation of atomic coordinates. Furthermore, the energy variation during the AIMD simulation increases with temperature and reaches 3.2 eV/f.u. at 600 K, where the α ′′ -phase structure is largely distorted, as indicated in the insets in Fig. 2a. On the other hand, the energy fluctuation of Fe 15 Al 1 N 2 phase is rather small within 1.5 eV/f.u. even at 600 K (Fig. 2b). In particular, the α ′′ phase tends to maintain up to 600 K (insets in Fig. 2b), although marginal phonon vibrations and atomic coordinate distortions occur during the AIMD simulation. In Fig. 3a, an even more notable finding is the persistent increase in K u as M changes from Ni to Ga, reaching the largest value of 1.85 MJ m -3 for Ga replacement. This value is more than 3 times the enhancement attained for α  -phase is associated with tetragonal distortion. The Jahn−Teller-like d-orbital level splitting when α → α ′′ offers more electronic energy level degrees of freedom 10 . More specifically, the tetragonal distortion splits the cubic e g and t 2g levels around E F : singlets a 1 ( d x 2 -y 2 ) and b 1 ( d 3r 2 -z 2 ), and singlet b 2 ( d xy ) and doublet e ( d yz,xz ), respectively. Evidently, the energy levels near E F of the spin-down electrons, especially the a 1 and b 2 states, differ at the 4e, 8h, and 4d sites (Fig. 1b). Further analyses indicate that the difference comes from their dissimilar hybridization with N-2p orbitals. Our analysis of atom resolved magnetocrystalline anisotropy energy (MAE) in Fig. 3c indicates that MAE in this α ′′ -phase distributes unequally over the unit cell: -0.14, 0.18, and -0.03 meV at the 4e, 8h, and 4d sites, respectively. Here, MAE is scaled down to the microscopic atomic level (meV/atom), rather than the macroscopic energy density (MJ m -3 ). From the k-resolved (minority-spin) eigenvalue analysis in Fig. 1b, featured bands with cone-like shapes occur: the minimum of a 1 (parabolic) and maximum of b 2 (reverse parabolic) dispersions touch at the Ŵ point. In particular, for Fe(8h), such a conical a 1 − b 2 pair appears right at E F . In association with their reduced eigenvalue difference across E F , the SOC matrix term in the Hamiltonian can thus increase the positive contribution to MAE, according to the perturbation theory 11 .
According to Fig. 3d, the Fe(4d) c site plays a major role in the M(Cu to Ga)-induced enhancement in MAE ( MAE ) rather than the Fe(4d) ab site. The contributions to MAE from the other 4e and 8h sites cannot be ignored although minor. Meanwhile, the conical a 1 − b 2 pair of the Fe(4d) c site moves gradually toward E F with the Cu to Ga replacement (Fig. 4a), which reflects the rigid-band model. A similar phenomenon is not present for the Fe(4d) ab site because of its longer separation (4 Å) from M than Fe(4d) c (3.1 Å). We therefore attribute MAE in Fe 15 M 1 N 2 to the joint effects of the Jahn−Teller level splitting and the supplied-electron-induced level changes of the d-orbitals.
In the rigid-band picture, the shift of the electronic states is related to the change in the energy of the Bloch state with M ( �ε k ), as ρ(ε) = ρ 0 (ε) − [∂ρ 0 (ε k )/∂ε k ]�ε k 27 , where ρ(ε) and ρ 0 (ε) are the density of states (DOS) of Fe(4d) c in Fe 15 M 1 N 2 and Fe 16 N 2 , respectively. For a small amount of M, �ε k is also small and thus independent of k, where the shape of the band structure remains the same but displaced by �ε k . Eventually, for the Ga replacement, the conical a 1 − b 2 pair of the Fe(4d) c site shifts down and appears near E F , which in turn enhances MAE. Here, the Jahn−Teller argument is not applicable, as the c/a (1.1) of α ′′ -Fe 16 N 2 remains almost the same upon M replacement (Table 1). www.nature.com/scientificreports/ In accordance with the a 1 − b 2 shift (Fig. 4a), the replacement element that can maximize MAE is Ga, in line with the obtained K u in Fig. 3a. To support this scenario, we explore the replacement element Al, which is isoelectronic to Ga. Remarkably, we find K u value of 1.80 MJ m -3 for the Al replacement (Fig. 3a), similar to that (1.85 MJ m -3 ) for the Ga replacement. Accordingly, similar electronic features of the a 1 − b 2 bands and �ε k at the Fe(4d) c site are identified in Fig. 4a for the same group elements, Al and Ga. Furthermore, as mentioned early in Fig. 1d, the inclusion of Al in Fe 16 N 2 greatly improves the α ′′ -phase stability beyond the other M-replacements. From a practical viewpoint, the Ga and Al replacements (particularly, Al) for Fe are desirable for RE-free permanent magnets because of their abundances on earth.
In line with the argument we outlined thus far, an even larger K u may be achieved if more Fe in Fe 16 N 2 are replaced with metal elements with more valence electrons. To test this scenario, we replace half of the Fe in Fe 16 N 2 with Co. From the total energy minimization, all the 4d and 4e sites are occupied by Co, forming the B2-phase, while 2 N prefer the 4 Fe(8h) coordinated octahedral interstices on the same ab plane. As Co has 1 more electron and stronger SOC than Fe, we find that the K u in (Fe 0.5 Co 0.5 ) 16 N 2 is 1.65 MJ m -3 . This value is more than double that (0.6 MJ m -3 ) of α ′′ -Fe 16 N 2 . A similar argument can be applied for other replacements such as Ni and Zn (not shown here). Furthermore, the enhanced c/a (1.17) of (Fe 0.5 Co 0.5 ) 16 N 2 , compared with 1.1 in Fe 16 N 2 , is clearly an additional cause of the large K u 8 .
Remarkably, the M-replaced (Fe 0.5 Co 0.5 ) 16 N 2 compounds exhibit a trend similar to, but with notably enhanced numerical values, that in Fe 16 N 2 : a nearly linear increase in K u from Cu to Ga (Fig. 3b). Eventually, Ga replacement leads to a K u as high as 2.44 MJ m -3 . Such supreme value of K u can also be achieved for Al (2.41 MJ m -3 ). These values are more than 4 times that of Fe 16 N 2 and more than half the value of 4.5 MJ m -3 of the typical REmagnet Nd 2 Fe 14 B 28 . Similar to in Fe 16 N 2 , Co next to M on the c axis (denoted Co c ) produces the largest MAE compared with other sites (Fig. 3f), although MAE is larger for Fe than for Co (Fig. 3e). At this Co c site, it is manifested that the main mechanism of enhancing K u is the displacement of the (unoccupied) a 1 band toward E F as M changes from Cu to Ga and Al in Fig. 4b.
We believe that the present argument is rather general and can be applied to other magnetic materials. To better justify the excess-electron-induced enhancement in K u , we forcibly increase the number of valence electrons in Fe 16 N 2 and (Fe 0.5 Co 0.5 ) 16 N 2 . This approach reflects an excess electron that is uniformly accumulated over all Fe rather than at a specific site neighboring the M replacement. For both Fe 16 N 2 and (Fe 0.5 Co 0.5 ) 16 N 2 , K u increases linearly as the number of excess electrons ( ) increases (Fig. 5a). Nearly the same values of K u of 1.3 MJ m -3 in Zn-replaced Fe 16 N 2 and 2.4 MJ m -3 in Ga-replaced (Fe 0.5 Co 0.5 ) 16 N 2 are reproduced at = 0.2 e/atom. From the simplified DOS analyses in Fig. 5b for Fe 16 N 2 , the unoccupied a 1 bands of all the Fe sites displace toward E F upon an increase in , while the occupied b 2 state is rather insensitive. This result again reveals that the d-orbital level change induced by supplied electrons is the main mechanism of the K u enhancement.
We now would like to highlight the intrinsic hard magnetic properties, including maximum theoretical energy product (BH) max , anisotropic field µ 0 H a , and hardness parameter κ , of the present compounds. The  31,32 . Furthermore, as listed in Table 1, the obtained T c values (804−855 K) of the stable compounds (M = Zn, Ga, and Al) are sufficient to fulfill the basic requirement (no less than 550 K) of permanent magnets 33 .

Conclusion
In summary, we show, using first-principles calculations and rigid-band model analysis of α ′′ -phase Fe 16 N 2 , that K u can be scaled up by a few times upon the substitution of metal elements with more valence electrons than Fe, from Co to Ga, without the inclusion of RE and HM elements. More remarkably, the replacement by simple metals (Al and Ga) has potential for simultaneous enhancements of K u and the thermal stability, which would make α ′′ -Fe 16 N 2 a possibly RE-free permanent magnet, along with its high Curie temperature and low materials price. Furthermore, we demonstrate that a similar argument, as a general rule, is applicable to suitable systems to achieve enhanced intrinsic hard magnetic properties and improved thermal stability. We hope that our results can be used as a guideline for subsequent experimental investigations of RE-free high-performance permanent magnetic materials.

Methods
The density-functional theory (DFT) calculations were performed using the projector augmented wave (PAW) method 34 , as implemented in the Vienna ab initio simulation package (VASP) 35 . The exchange-correlation interactions are treated with the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) 36 . We used an energy cutoff of 500 eV and a 11 × 11× 11 Brillouin zone (BZ) k-point mesh to relax the lattice parameters and atomic coordinates until the largest force decreased to below 10 -2 eV/Å. The total energy method is applied to obtain K u , which is expressed as K u = (E a − E c )/volume , where E a and E c are the total energies with www.nature.com/scientificreports/ magnetization along the a and c axes, respectively. To obtain well-converged K u , we impose a denser k-point mesh of 15 × 15× 15 with a smaller smearing of 0.05 in the Gaussian method, where the convergence of K u with respect to the number of k points and smearing parameter is ensured. In tetragonal symmetry, K u is expressed as K u ≈ K 1 sin 2 θ + K 2 sin 4 θ , where K 1 and K 2 are the magnetic anisotropy constants and θ is the polar angle between the magnetization vector and the easy axis (c axis in the present system). For θ = π/2 , K u = K 1 + K 2 . It is a formidable task to ensure numerical results of K u with all electron methods, if we start from scratch. To this end, we have also performed full-potential calculations using the WIEN2K package 37 , adopting the optimized lattice constants and ionic positions obtained from the VASP calculations. The two methods produce consistent results. In the AIMD simulation, we adopted the Nosé-thermostat algorithm to model a canonical ensemble 38 . A time step of 1 fs and 10000 ionic steps were used for the total simulation time of 10 ps with the Ŵ-point BZ integration, where the lattice parameters and atomic coordinates are allowed to relax at constant volume. The numerical calculations for magnetization dynamics and T c were carried out using Monte Carlo simulation based on the Heisenberg model in the VAMPIRE package 30 . Here, the Heisenberg spin Hamiltonian is defined by where J ij is the exchange interaction between two spins S i at the i site and S j at the j site. The exchange interaction parameters, from the first to the third nearest neighbor atoms, were estimated by the constrained local moment approach in the VASP calculations. More detailed methodology is provided in Ref. 9 .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. (1)