First principles and atomistic calculation of the magnetic anisotropy of Y₂Fe₁₄B

Rare-earth (RE) transitional metal permanent magnetic materials play a critical role in hybrid and electric vehicles and electric power generation.¹ Recent concern about the impact of climate change has renewed interest in understanding and optimizing these materials to improve the energy efficiency of these key technologies. The most technologically important Nd–Fe–B magnet consists of 2:14:1 phase Nd$2$Fe$14$B. This Nd–Fe–B magnet has the highest energy product among all known permanent magnet materials.¹ The magnetocrystalline anisotropy (MAE) is a key factor for understanding the high coercivity of RE$2$TM$14$M (RE = rare earth, TM = transition metals, and M = B, C, N) permanent magnets.² These elements form stoichiometric compounds in the 2:14:1 phase of rare earths, transition metals, and metalloids, respectively, which allows the study of different magnetic couplings RE–RE and RE–TM between the different sites.^3–5 Related with the MAE is the effect of strain on these materials since the manufacturing process could promote variations in their lattice parameters and therefore some residual strain.^6–8 

In addition, rare earths have localized 4$f$ electrons, adding further complexity to their study. These $f$ electrons are relatively insensitive to their environment in contrast with 3$d$ electrons that are quite sensitive to lattice changes due to their itinerancy. As proposed by Torbatian et al.,⁸ the compound Y$2$Fe$14$B is suitable for the study of Fe–3$d$ electrons since Y does not have $f$ electrons, being a prototypical $f0$ rare earth element keeping the same geometrical structure as the other RE-TM bulk materials.

Due to the structural symmetry of RE$2$Fe$14$B series alloys, studies of Y$2$Fe$14$B play an essential role in understanding the contribution of the Fe sublattice to the overall MAE for the case of more technologically relevant alloys including Nd or Dy. In the majority of RE$2$Fe$14$B series alloys, the rare earth ions dominate the MAE at low temperatures.⁹ However, the combination of the strong temperature dependence of the RE sublattice magnetization and higher order contributions to the MAE lead to a strong temperature dependence of the MAE, so that at room temperature, the Fe sublattice contributes a significant fraction of the total MAE.⁹ A common feature of the so-called “non-magnetic” RE$2$Fe$14$B alloys, where RE = Y, Ce, Th is an increase of the anisotropy field with increasing temperature, in direct contrast to the usual expectation of a reduction of the MAE due to spin fluctuations.^10–12 Bolzoni et al.¹³ suggested that the origin of this anomalous increase of the MAE may be due to asymmetric temperature dependent lattice expansion,¹⁴ causing an increase of the effective MAE at elevated temperatures. Torbatian et al. calculated the MAE in Y$2$Fe$14$B for two cases: an equilibrium and a compressed lattice, finding that the MAE is enhanced upon lattice compression. However, understanding the influence of lattice changes on the MAE in Y$2$Fe$14$B for a wider range of $a$ and $c$ values is essential to ascertain whether temperature-induced lattice expansion can account for the observed increase in the anisotropy field at elevated temperatures.

Here, we systematically investigate the effects of the lattice parameters $a$ and $c$ on the magnetic properties of Y$2$Fe$14$B, including changes in the local magnetic moments, electronic structure, and MAE. The paper is structured as follows. In Sec. II is presented in brief the theoretical tools employed to perform the calculations. In Sec. III A, a geometric analysis of the optimized bulk structures is presented. The projected density of states on each atomic species and site is presented in Sec. III B together with the charge transfer between different atomic species. The analysis of the magnetic moments is presented in Sec. III C, and the magnetic anisotropy calculated values are presented in Sec. III D. The dependence of the anisotropy with the temperature is shown in Sec. IV. Finally, in Sec. V we summarize the main conclusions of the work.

We have undertaken geometrical, electronic, and magnetic structure calculations of the Y$2$Fe$14$B bulk alloy by means of DFT using the SIESTA code.^15,16 To describe the core electrons, we have used fully separable Kleinmann–Bylander¹⁷ and norm-conserving pseudopotentials (PPs) of the Troulliers–Martins¹⁸ type. As exchange correlation (XC) potential we have employed the generalized gradient approximation (GGA) following the Perdew, Burke, and Ernzerhof (PBE) version.¹⁹ To have a better description of magnetic systems, pseudocore (pc) corrections were used to include in the XC terms not only the valence charge density but also the core charge.²⁰ In general, the correction will only be significant in the range where valence and core charges overlap. As a basis set, we have employed double-$ζ$ polarized (DZP) strictly localized numerical atomic orbitals (AO). The electronic temperature—kT in the Fermi-Dirac distribution—was set to 50 meV. Real space integrals are computed over a three-dimensional grid with a resolution of 1600 Ry, a mesh fine enough to ensure convergence of the electronic/magnetic properties. To obtain the charge distribution, we have used the Mulliken partitioning scheme.²¹ 

The MAE is defined as the difference in the total self-consistent energy between hard and easy magnetization directions. MAE values are commonly of the order of meV and hence it is necessary to perform an accurate calculation through the convergence of relevant DFT parameters such as the number of $k$-points. For Y$2$Fe$14$B, we find that the MAE values are particularly sensitive to the $k$-point sampling. We therefore carried out an exhaustive analysis of the MAE convergence in order to achieve a total energy tolerance below 10$−5$ eV. We employed more than 500 $k$ points in the calculations for each geometric configuration. To obtain the MAE, we have used a recent implementation of the off-site spin–orbit coupling (SOC)^22–25 in the SIESTA code. This approximation takes into account not only the local SOC contributions to the total energy but also the neighboring interactions between atoms to obtain the total self-consistent energy.

Figure 1 shows a schematic view of the crystal structure of Y$2$Fe$14$B. It is composed of a tetragonal unit cell (UC) with a space group of $P$4$2$/$mnm$ and the UC has two different Y sites (⁠$f$⁠, $g$⁠), six distinct Fe sites (⁠$j1$⁠, $j2$⁠, $k1$⁠, $k2$⁠, $e$ and $c$⁠), and one B site, with a total of 68 atoms in the unit cell. In this work, we have performed the study of seven in-plane lattice values: $a=8.48Å+k⋅δa$⁠, with $δa=+0.10$ Å and $k=0,…,6$⁠. Variations of the $a$ lattice parameter will also naturally lead to a variation of the $c/a$ ratio, and so it was necessary to optimize the tetragonal out-of-plane distortions for each value of $a$⁠. Accordingly, $c/a$ values are found to be: 1.45, 1.42, 1.40, 1.37, 1.34, 1.31, and 1.29, respectively.

All of the geometric configurations studied in the present work were subject to an optimization process. For the range of study, we chose in-plane lattice constants around the experimental room temperature value of 8.75 Å (see the end of Sec. II).¹⁴ The optimization process was carried out for fixed in-plane lattice constant $a$ by performing self-consistent energy calculations for a series of out-of-plane $c$ values. Plotting the total energy as a function of $c$ we are able to obtain local minima, geometrically characterizing each unit cell in terms of the $a$ and $c$ parameters. In this respect, the ground state (GS) energy value was achieved for the pair $(a,c)=(8.78Å,12.00Å)$⁠. As noted in Sec. II, we are using the generalized gradient approximation as the XC functional and so the bonds between atoms will be slightly larger than the experimental ones.^5,14 However, our optimized lattice value is within 0.03 Å of the values obtained from previous studies of Y$2$Fe$14$B. As an example of this process, we show in Fig. 2 three different optimization curves that show the energy dependence with respect to the out-of-plane $c$ for fixed $a$⁠. All the values have been shifted with respect to the lowest energy value pair (filled black dot curve), which shows clearly the energy difference of $∼0.25$ eV between the GS and the other values.

Figure 3 shows the radial distribution function (RDF) of the bond distances between the Fe atoms located in the upper part of the unit cell and those (denoted Fe') located in the lower side, A and B, respectively. A visualization of the atom groupings is shown in the inset in Fig. 3. The plotted range represents small bond values, below 4.5 Å, ensuring that only nearest and next nearest neighbors are taken into account. In order to encompass the geometrical evolution of the UC for different in-plane lattice values, we have chosen a representative $a$ value for each zone: 8.48, 8.78, and 9.08 Å, small, intermediate, and large, respectively. On each graph two different peaks appear around $∼2.5$ Å and $∼4.25$ Å. Clearly, these peaks represent bonds between the first and second neighbors of a specific Fe atom, respectively. It is beyond the scope of this work to discriminate inequivalent Fe distances between different groups located in the upper and lower UC. However, after careful inspection of each configuration, we are able to give a qualitative understanding of the geometric expansion/contraction of the bulk Y$2$Fe$14$B alloy. We find different structural properties for Fe atoms at different points in the unit cell, corresponding to Fe and Fe” sites in Fig. 3. In particular, the small bond-length range represents bonds mainly between Fe atoms at different out-of-plane level and the longer to those located almost in the same horizontal plane. The overall trend for d$Fe−Fe$ and d$Fe′−Fe′$ is quite similar, having a small positive deviation around each peak as $a$ increases. For example, focusing on Fig. 3(a) we observe that the shift is more pronounced around 4.25 Å than for the smaller bond-length, with the Fe–Fe and Fe’–Fe” bond distances increasing by 0.2 Å with $a$ expansion. Then, $a$ variations cause the Fe atoms located at the same plane level (⁠$∼4.2$ Å) to lie further away from each other while those Fe species at different planes are positioned almost at the same distance. In summary, we conclude that locally the variation of $a$ promotes that the bonds between different groups of Fe atoms behave differently on whether they are located at the lower or upper part of the UC. However, the overall out-of-plane contraction when $a$ increases is mainly due to the upper and lower Fe blocks approaching. This makes that the UC keeps constant its volume.

In Fig. 4 (left), the spin-resolved projected density of states (PDOS) is shown for three different in-plane lattice configurations. In each column, the Fe DOS were projected to show the contributions from different sites, namely, $k$⁠, $j$ and $c,e$ (see Fig. 1 for inequivalent Fe sites within the unit cell). In general, after inspection of these different Fe sites, we observe that only the $k1$ and $k2$ DOS share a similar form. The $j$⁠, $c$⁠, and $e$ Fe sites have their $d$-band peaks at different energy positions, for example, the Fe(⁠$j1$⁠) up-states (solid line) have two peaks at 0.8 and 2 eV and in (⁠$j2$⁠) (dashed line) these peaks have shifted to the lower energy values of 3.5 and 2 eV, promoting a population of the Fe(⁠$j2$⁠)-$d$ band. We can see this behavior in detail in Fig. 5 where the Fe(⁠$j1$⁠) (empty circles) have lost more charge than Fe(⁠$j2$⁠), with respect to the isolated atoms. In the same fashion, we are able to explain the charge transfer of $d$-bands in Fe(⁠$e$⁠) and Fe(⁠$c$⁠) (last column). In a column, the shape of the DOS curves shares a common trend: as $a$ increases, upper to lower panels, the $d$-band width becomes slightly smaller and higher, having the states concentrated in a smaller energy range. For instance, for Fe(⁠$k$⁠), the energy extends from zero to almost $−6$ eV for the smaller $a$ and this width is reduced for bigger lattice up to $−5$ eV. In other cases, even though the broadness is similar for small and large lattice, the center of the $d$-band tends to increase, concentrating more states in the range of $−3$ to $−2$ eV. It is also interesting to note in Fig. 4 (right) the two peaks around $−8$ and $−9$ eV in Fe(⁠$e$⁠) and Fe(⁠$k1$⁠), which permit interaction with the Y and B states.

The PDOS for the Y(⁠$g$⁠), Y(⁠$f$⁠), and B(⁠$g$⁠) species are depicted in Fig. 4 (right). As we will see in Sec. III C, Y and B are antiferromagnetically coupled to Fe, having an excess of down-states in the population. In general, the in-plane lattice constant expansion for these species does not significantly affect the shape of the Y PDOS. However, around the Fermi level, both Y sites differ in height, implying a slightly different amount of up-charge contributing to different imbalance in the up-/down-populations. As previously mentioned, there are two peaks; at $−8$ and $−9$ eV belonging to the B species, mainly of $s$ character, that permit the hybridization of Fe(⁠$k1$⁠) and Fe(⁠$e$⁠) with the B states. We note that these Fe atoms are located close to B species allowing this interaction.

To obtain further insight into the charge distribution on each atomic species in Fig. 5, we analyze the increase or reduction of the total charge (not spin-resolved) on each atom compared to its isolated counterpart. A negative value means that the species has lost charge and positive means that the atom has been charged. The upper panel shows how the charge on the different Fe sites changes with the lattice spacing. Only Fe(⁠$k1$⁠) is charged positively, however its net charge tends to be, for bigger $a$ values, close to that of an isolated atom. The remaining Fe sites have transferred charge to the neighbors by different amounts and, depending on the site that they occupy in the unit cell, their net charge will be almost constant or will decrease with the lattice expansion. This is exactly what happens at Fe(⁠$c$⁠) that reduces its charge by $∼0.20$e/at when $a$ changes from 8.48 to 9.08 Å. Fe(⁠$k2$⁠), Fe(⁠$j1$⁠) and Fe(⁠$j2$⁠) remain nearly constant in all this $a$ range. Analysis of the Y(⁠$f$⁠) and B(⁠$g$⁠) charge (lower panel in Fig. 5) shows an enhancement of their net charges mainly due to the Fe neighbors and from the B species as well.

We now present the total and local magnetic moment (MM) values of the Y$2$Fe$14$B bulk alloy for different (⁠$a,c$⁠) pairs. Our total MM predictions assign bigger MM values to $a=9.08$ Å and smaller to $a=8.48$ Å, depicting an enhancement of the MM as the volume of the UC expands. The complete MM/f.u. values sequence for each (⁠$a,c$⁠) pair is: 28.70, 29.77, 30.60, 30.99, 31.24, 31.43, and 31.79, the first one corresponding to the smaller in-plane $a$⁠. Therefore, our GS configuration (8.78 Å) presents a total MM of 30.99 $μB$/f.u., in good agreement with other experimental and theoretical works.^4,5,8 In addition, it is of huge importance to understand how different Fe, Y, and B sites contribute to the total MM.

In Fig. 6(a), we show the average magnetic moment (MM) for different Fe sites for a common lattice value. The MM has been calculated as MM$M/NM$/N$M/NM$⁠, where N$M/NM$ are the total number of magnetic (M) or non-magnetic (NM) species. All the Fe sites have a clear tendency of their MM values to increase as the in-plane lattice constant increases. The largest growth (of 0.55 $μB$/at) takes place for the Fe(⁠$j1$⁠) and the smallest variation is for Fe(⁠$c$⁠). The remaining sites exhibit an increase to a greater or lesser extent. We are able to observe this trend qualitatively by inspecting the spin-resolved projected DOS shown in Fig. 4 (left). For example, in Fe(⁠$j1$⁠) (black solid curve), the total up-states increase with $a$⁠; this reflects that the initial small shoulder at the Fermi level displaces inwards toward the valence band allowing an increase in the spin up charge. Conversely, the down-states act in a different manner, reducing the down charge. As a result, the imbalance in the up/down charge will be more pronounced for the expanded lattices and will promote an enhancement of the on-site Fe(⁠$j1$⁠) MM values when $a=9.08$ Å. Through the $j2$ PDOS curves, it is also easy to understand why they have the higher MM values among all of the Fe sites. They present bigger imbalance in their up-/down-states in comparison with the other Fe sites.

The Y(⁠$f$⁠), Y(⁠$g$⁠), and B(⁠$g$⁠) MM values are presented in Fig. 6(b). The net MM per atom of these two species are negative for all the lattice constant values [labeled in Fig. 6(b) by filled and empty squares, triangles, circles and empty diamonds]. The Y’s MM values exhibit different behavior depending on whether the lattice constant changes seen by its vertical dispersion between $−1.05$ and $−0.8μB$/at. Both, Y(f) and Y(g) location sites, tend to behave similarly. Regarding the variation with $a$⁠, the tendency is to increase the net MM with lattice expansion, by around $−0.3μB$/at for both Y(⁠$f$⁠) and Y(⁠$g$⁠). On the other hand, for the Y atoms located at different sites, the MM are slightly bigger for Y(⁠$f$⁠) than for Y(⁠$g$⁠), the difference being at most of 0.04 $μB$/at. B MM values behave in a different manner since they do not change with the lattice expansion having a common value of $−0.20μB$/at. In general, the down-states population exceeds that of the up-states, resulting in a negative MM for the Y and B atoms. We also note the presence of a significant magnetic moment on the Y sites, in contradiction to the usual assumption of a non-magnetic RE.⁹ 

We now present the calculations of the total magnetic anisotropy energy.

As was pointed out in Sec. II, the MAE values are of the order of a few meV. Consequently, the calculation of the self-consistent total energies involved have to be sufficiently accurate. In the same fashion as Kitagawa and Asari,⁵ we performed convergence tests for the $k$-points sampling and other relevant DFT parameters, resulting in sufficient accuracy for the quantities under study. In Fig. 7 is shown the MAE convergence for the Y$2$Fe$14$B bulk alloy with increasing $k$-points. The dashed black lines indicate our required accuracy in the MAE values, of $5×10−5$ eV. We observe that in the range from 567 $k$-points (indicated in red in Fig. 7) to 1400 all calculated values lie within the required tolerance. Consequently, we decided to use 567 $k$-points to perform the calculation of the MAE for all the configurations. It is worth noting that, although Isao et al. used a code based on a linear combination of pseudo-atomic-orbital, similar to the SIESTA code, the convergence of each code will depend not only on the scheme/formalism to solve the Kohn–Sham equations but also on the pseudopotential and many other parameters involved.

In Fig. 8, we show the MAE as a function of the optimized $c/a$ ratio. As was pointed out at the end of Sec. II, an increase (reduction) of $a$ means a reduction (increase) of the out-of-plane constant, thus from right to left along the X axis will mean an in-plane lattice expansion. The MAE’s dispersion over the studied $c/a$ range is of 2.8 meV. Specifically, the MAE tends to zero when $a$ increases and hence $c/a$ decreases while it presents closer values to 2.8 meV when the in-plane $a$ constant decreases promoting larger $c/a$ values. The ground state geometry configuration depicts a MAE value of 1.74 meV/u.c., that is, in reasonable agreement with previously reported works of Torbatian et al ⁸ of 2.6 meV/u.c., Miura et al.²⁶ of 0.622 meV/u.c. and Isao et al. ⁵ 3.2 meV/u.c. The present trend of the bulk MAE of Y$2$Fe$14$B shown here is similar to the FePt–L1$0$ on MgO, where the MAE reduces its value as the in-plane expansion takes place.²⁷ In FePt, the tetragonality of the system contributes to the MAE, and we argue that the similar qualitative behavior seen for Y$2$Fe$14$B is due to the same physical effect due to the anisotropy of the local electronic environment.

Finally, we turn to the motivation of this work, which is to understand the unusual (non-monotonic) temperature dependence of the MAE in Y$2$Fe$14$B. For clarity, we refer to the MAE as the intrinsic anisotropy value from DFT calculations and the anisotropy as the observed value at non-zero temperature. The latter is the free energy of the system, which can have two contributions: first, the thermodynamic decrease driven by spin fluctuations,¹¹ and second any change in the MAE arising from lattice expansion. Clearly, this can potentially lead to non-monotonic temperature dependence of the anisotropy and as such is the first option considered here.

For this calculation, we employ an atomistic model, which introduces the thermally driven spin fluctuations which, in the absence of a temperature dependent intrinsic anisotropy, are responsible for the reduction in the observed anisotropy value at non-zero temperatures.

Using the atomistic model, we calculate the temperature variation of the magnetisation using conventional Metropolis Monte Carlo methods and the temperature dependence of the anisotropy using a constrained Monte Carlo algorithm¹² with adaptive updates,²⁸ both implemented in the vampire software package.^29,30

First, we consider whether the non-monotonic behavior of the measured anisotropy can be explained through our predicted change in the MAE due to changes in the $c/a$ ratio coupled with the thermally induced spin fluctuations. However, our calculations (not shown) give a conventional monotonic decrease of anisotropy with temperature.

We can interpret this as follows. We can estimate the reduction in anisotropy over a given temperature range using the change in magnetization and the Callen–Callen law scaling of anisotropy with $Mn$ with exponent $n=3$ for uniaxial anisotropy. From the temperature dependence of the magnetization, shown in Fig. 9 it can be seen that there is approximately a 20% reduction in $M(T)$ at $T=300$ K. Consequently, to compensate for the corresponding decrease in anisotropy due to spin fluctuations we require at least a factor 2 increase in the MAE. This is the minimum required to give rise to an increase of the anisotropy with temperature. Qualitatively from the DFT calculations, we find an increase in the intrinsic MAE with increasing $c/a$⁠. However, the experimental range¹⁴ for the $c/a$ ratio in the bulk is typically 1.37–1.38 up to $TC$⁠, which suggests a relatively modest increase in the MAE from our calculations. Therefore, while we obtain the correct qualitative behavior, the calculations are so far unable to explain the cause of the large increase in MAE. In particular, we note that the temperature variation of the MAE is governed by a combination of spin fluctuations, leading to a decrease in MAE, coupled with any change in the intrinsic MAE due to the lattice expansion as investigated here. While our calculations predict an increase of MAE with lattice expansion, the consequent increase with T is not sufficient to overcome the effects of spin fluctuations, leading to a monotonic decrease in the MAE, in contrast to experiments. Thus, our calculations show that the previous interpretation where the increase is attributed to lattice expansion effects^13,14 seems to be unlikely.

We finally consider the possibility of a combination of opposing single-ion and two-ion anisotropy as the origin of the large increase in anisotropy. We note that the two-ion anisotropy is essentially an exchange anisotropy which is not accessible to our DFT calculations. Recent investigations³¹ suggest that the origin of an increase in anisotropy with temperature may be due to a competition between single-ion and exchange anisotropy. To extract the two-site anisotropy, one uses a generalized tensorial form of the exchange

from which one obtains a two-site (exchange) anisotropy as

Our ab initio calculations cannot access the tensorial exchange. Consequently, for this investigation, we introduce a Hamiltonian including a phenomenological two-site anisotropy while considering the 2-14-1 crystallographic positions of the Fe sites,³² 

Equation (3) has three terms; the first being the Heisenberg exchange term and the second and third terms representing two-ion and single-ion anisotropies, respectively. The exchange interactions are parameterized with an exponential distance dependence with a cutoff of 5 Å and normalized to achieve a Curie temperature of 560 K in agreement with experimental measurements. We note that the single- and two-ion terms have different temperature dependences. Specifically, in terms of the temperature scaling $K(T)/K(T=0)=(M(T)/M(T=0))n$⁠, the exponent $n=2.28,3$ for the two- and single-ion terms, respectively.^11,31 Thus, for $kij$ and $ki$ of opposite sign, a non-monotonic variation $K(T)$ of the net anisotropy becomes possible.

Figure 10 shows the simulated temperature dependence of the effective magnetic anisotropy of Y$2$Fe$14$B using the constrained Monte Carlo algorithm¹² with adaptive updates²⁸ implemented in the vampire software package.^29,30 We apply temperature rescaling to achieve a correct description of spin fluctuations and magnetization with temperature.³³ The optimized anisotropy constants to achieve agreement with the experimental data are $kij=+1.83546×10−22$ J/atom and $ki=−1.743687×10−22$ J/atom, noting the large values and competing signs of the anisotropy constants. While the agreement between the simulation and experimental data is good, there are some small systematic differences that may be accounted for with strain effects determined from our ab initio simulations, or indeed due to the temperature dependence of the Hamiltonian parameters themselves.

We wish to note the existence of previous explanations of the increasing anisotropy with temperature from crystal-field effects suggested by Bolzoni et al.¹³ This explanation was not based on a microscopic description but assumed arbitrary higher order anisotropy terms that are fictitious in nature and the model could not perfectly explain the experimental behavior. The crystal-field analysis by Miura et al.³⁵ concludes that the underlying origin of anisotropy in the Fe sublattice remains an unresolved problem. In comparison, our analysis of mixed anisotropy is based on a microscopic model and provides a close quantitative agreement with the observed increase in the magnetic anisotropy energy. Further ab initio calculations may be able to determine the relative contributions from two-ion and single-ion anisotropy,³⁶ though these calculations are still excessively expensive in terms of computational complexity due to the large cell size and number of electrons.

We have found that tuning the unit cell geometry of the Y$2$Fe$14$B bulk alloy, specifically changing the in-plane/out-of-plane constants $a$ and $c$⁠, leads to a significant change in the magnetocrystalline anisotropy energy. In the present work, we have undertaken a geometric, electronic, and magnetic study of seven different configurations of the Y$2$Fe$14$B alloy where the in-plane $a$ range was varied from 8.48 to 9.08 Å and the corresponding out-of-plane $c$ parameter was optimized for each case. The ground state configuration has an in-plane value of 8.78 Å with $c=12.00$ Å and its total MM is 30.99 $μB$/f.u. in good agreement with previous experimental and theoretical results. Breaking up the total MM for individual Fe, Y, and B inequivalent sites, we have demonstrated that the higher MM value is for the Fe(⁠$j2$⁠) sites and the higher dispersion acts on Fe(⁠$j1$⁠), having a lower value for smaller $a$ increasing for the larger $a$ value. The Y and B species are antiferromagnetically coupled and only the Y changes similarly to the pair (⁠$a$⁠,$c$⁠), decreasing their net MM value as the lattice is compressed. It is clear that the strain has an important impact on the magnetic anisotropy of these alloys so that a volume reduction of the UC promotes higher values of the MAE. The MAE dispersion value for the present work is around 2.8 meV, having almost zero for the bigger unit cell and 2.8 meV for the (8.48 Å,1.46) pair. Our calculations of the variation of the MAE with lattice expansion is not sufficient to explain the non-monotonic variation of MAE with temperature observed experimentally. Atomistic simulations including contributions from competing two-ion and single-ion anisotropies are able to reproduce the observed increase in anisotropy with increasing temperature and offer an alternative explanation for the phenomenon in RE-TM intermetallic alloys. The electronic origins of anisotropy require further investigation due to the essential role played by Y$2$Fe$14$B in understanding the contribution of the Fe sublattice to the overall MAE for the case of more technologically relevant alloys including Nd or Dy.

This work is based on results obtained from the future pioneering program “Development of magnetic material technology for high-efficiency motors” commissioned by the New Energy and Industrial Technology Development Organization (NEDO).

The data that support the findings of this study are available from the corresponding author upon reasonable request.