Monte Carlo calculations of Curie temperatures of Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$ pseudobinary system

The close-packed AB$_2$ structures called Laves phases constitute the largest group of intermetallic compounds. In this paper we computationally investigated the pseudo-binary Laves phase system Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$ spanning between the YFe$_2$, YCo$_2$, GdFe$_2$, and GdCo$_2$ vertices. While the vast majority of the Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$ phase diagram is the ferrimagnetic phase, YCo$_2$ along with a narrow range of concentrations around it is the paramagnetic phase. We presented results obtained by Monte Carlo simulations of the Heisenberg model with parameters derived from first-principles calculations. For calculations, we used the Uppsala atomistic spin dynamics (UppASD) code together with the spin-polarized relativistic Korringa-Kohn-Rostoker (SPR-KKR) code. From first principles we calculated the magnetic moments and exchange integrals for the considered pseudo-binary system, together with spin-polarized densities of states for boundary compositions. Furthermore, we showed how the compensation point with the effective zero total moment depends on the concentration in the considered ferrimagnetic phases. However, the main result of our study was the determination of the Curie temperature dependence for the system Y$_{1-x}$Gd$_x$(Fe$_{1-y}$Co$_y$)$_2$. Except for the paramagnetic region around YCo$_2$, the predicted temperatures were in good qualitative and quantitative agreement with experimental results, which confirmed the ability of the method to predict magnetic transition temperatures for systems containing up to three different magnetic elements (Fe, Co, and Gd) simultaneously. For the Y(Fe$_{1-y}$Co$_y$)$_2$ and Gd(Fe$_{1-y}$Co$_y$)$_2$ systems our calculations matched the experimentally-confirmed Slater-Pauling-like behavior of T$_C$ dependence on the Co concentration.


Introduction
The largest group of intermetallic compounds are the Laves phases [1]. They are binary close-packed structures with a formula AB 2 , which are found in three different types -hexagonal MgZn 2 -type (C14), cubic MgCu 2 -type (C15), and hexagonal MgNi 2 -type (C36), see Fig. 1. In this work, we study theoretically the pseudo-binary cubic Laves phase Y 1−x Gd x (Fe 1−y Co y ) 2 . Its boundary cases are quite well known. GdFe 2 and GdCo 2 are ferrimagnets with Curie temperature (T C ) equal to 805 K [2] and 390 K [3], respectively. YFe 2 is a ferromagnet with Curie temperature equal to 545 K [2] and YCo 2 is an exchangeenhanced Pauli paramagnet [4] that undergoes a metamagnetic transition in a field of 70 T at 10 K [5,6]. In Fig. 2 we present a summary of the experimental Curie temperature -concentration relationships for all boundary concentrations [2,3,[7][8][9]. We expect the Y 1−x Gd x (Fe 1−y Co y ) 2 YFe 2 has been identified as a hydrogen storage material [6,10,11]. YCo 2 and its alloys have been considered for permanent magnet applications [12]. For the system Y 1−x Gd x Co 2 , the magnetocaloric [13], thermopower [14], and electronic transport [15] properties have been examined. For the system Y(Fe 1−y Co y ) 2 , electrical resistivity and Mössbauer effect have been investigated [7]. Moreover, rare-earth compounds with large localized magnetic moment, such as Gd, Tb or Er, are some of the most promising candidates for magnetic cooling among Laves phases [16].
Our previous experimental and theoretical studies on the Y 1−x Gd x (Fe 1−y Co y ) 2 system have included such topics as: effect of YCo 2 doping [17], magnetocaloric effects in Y 1−x Gd x Co 2 [13], Curie temperature of Y(Fe 1−y Co y ) 2 [18], structural disorder in YCo 2 [19], and electronic specific heat coefficient of Y(Fe 1−y Co y ) 2 [20]. In this work we show the Curie temperatures determined from first principles for pseudo-binary Laves phases Y 1−x Gd x (Fe 1−y Co y ) 2 . The two-site chemical disorder is modeled using Monte Carlo (MC) simulations. The T C 's are obtained by fashioning the Heisenberg model Hamiltonian with MC simulations using parameters from firstprinciples calculations. We study the dependence of the exchange integrals on the interatomic distance and analyze the behavior of the total and partial magnetic moments calculated from first principles. We also investigate the valence bands of the stoichiometric boundary compositions. Furthermore, the MC method also allows us to predict the temperature dependence of magnetization (M ) and magnetic susceptibility (χ). T C calculations were performed  for systems containing up to four different elements simultaneously, including three magnetic elements. Although we have previously presented T C dependence on concentration, e.g. for Y(Fe 1−y Co y ) 2 [18], in this work for the first time our calculations are based fully on MC simulations and not on the disordered local moment method as before.

Calculations' details
One method for determining the Curie temperature is Monte Carlo simulations of the Heisenberg model. The simulations allow to determine the temperature dependence of magnetization and magnetic susceptibility. The location of the peak in magnetic susceptibility allows us to determine T C with an accuracy of ±10 K. We obtained the values of the magnetic moments and the exchange integrals necessary to perform the MC simulations from first-principles calculations utilizing the spin polarized relativistic Korringa-Kohn-Rostoker (SPR-KKR) code, version 7.7 [21,22]. We used the generalized gradient approximation (GGA) in the Perdew, Burke, and Ernzerhof parametrization [23] and the atomic-sphere approximation (ASA) in the fully-relativistic approach. We used the coherent potential approximation (CPA) [24] to simulate chemical disorder. We used the basis functions up to l = 4, a 45 × 45 × 45 k-mesh, and 40 energy points. Exchange integrals were obtained using the method of Liechtenstein et al. [25] with respect to the ferrimagnetic ground state. Although in previous work we have used GGA + U corrections to describe both d [26] and f [27] valence electrons, due to the exploratory nature of this work we decided not to extend the current model beyond standard GGA. The effect of on-site Hubbard-type interactions on the values of the exchange integrals has recently been discussed elsewhere [28].  For Monte Carlo simulations of the Heisenberg Hamiltonian we used the Uppsala atomistic spin dynamics (Up-pASD) code [29,30]. The simulated system consisted of ∼13000 atoms with periodic boundary conditions. The radius of the exchange integrals cutoff sphere in the Heisenberg Hamiltonian was set to 1.5 lattice parameters. The exchange integrals were assumed to be temperature independent in the MC simulations. The simulations were performed using classical (Boltzmann) statistics, but it is worth noting that it has recently been possible to use quantum (Bose-Einstein) statistics using the UppASD code [31]. The magnetic moment on Y, since it is induced and expected to vanish with temperature, was set equal to zero in the MC simulations (first-principles calculations gave a moment of ∼0.5 µ B ). In the MC simulations, we used exchange integrals calculated with respect to the ferrimagnetic ground state. We justify this on the grounds that MC simulations of GdCo 2 with exchange integrals obtained with respect to the paramagnetic ground state have yielded that material to be a paramagnet. This can be explained by the fact that the Co sublattice is metamagnetic in GdCo 2 [32].
Experimental lattice parameters were used to model the boundary compositions YFe 2 (7.36), YCo 2 (7.22), GdFe 2 (7.38), and GdCo 2 (7.24 Å) [7,33]. For intermediate Y 1−x Gd x (Fe 1−y Co y ) 2 concentrations, we assumed a linear behavior of the lattice parameters. The full range of Co and Gd concentrations were prepared with a step of 0.1, leading to a total of 121 cases considered (11 × 11). The space group and atomic positions of the C15 Laves phase are shown in Table 1. The picture of the unit cell, generated using the VESTA code [34], is shown in Fig. 1.

Results and Discussion
In Fig. 3 we presents a Curie temperature dependence of Gd and Co concentrations for Y 1−x Gd x (Fe 1−y Co y ) 2 determined from Monte Carlo simulations with parameters from first-principles calculations. For comparison, the four curves (red, green, blue, and yellow) are plotted based on experimental results [2,3,[7][8][9]. Along the Co concentrations (y), our MC simulations preserve the characteristic Slater-Pauling behavior for T C (y). For Y 1−x Gd x Fe 2 , we reproduce the linear behavior of the dependence of T C on Gd concentration, please compare with the bottom panel of Fig. 2. Our T C estimates agree relatively well with experimental results in all boundary regions except for the YCo 2 neighborhood. The reason for the observed discrepancy is a combination of two factors: (a) the paramagnetic nature of the YCo 2 remaining on the verge of satisfying the Stoner criterion (N (E F )I ≃ 0.9) [20,32], and (b) the previously raised problem of GGA's failure to correctly describe the magnetic ground state of Co [18]. Leaving aside the overestimation of the results near YCo 2 , the map obtained for a system containing three magnetic elements shows that a realistic first-principles T C analysis for complex magnetic materials is possible allowing, for example, to determine the concentrations of the individual components leading to a specific T C value, e.g. 300 K, as sought for magnetocaloric materials for use in magnetic coolers.
T C values were estimated from magnetic moments and exchange integrals calculated from first principles. Figure 4 shows the calculated exchange integrals as a function of normalized interatomic distance for considered ferrimagnetic boundary compositions. It is easy to see that in each case the dominant contribution to the Heisenberg Hamiltonian comes from first-neighbor interactions at the Fe/Co sites. We verified that restricting the MC simulation to this dominant exchange interaction would lead to an underestimation of the Curie temperature on the order of hundreds of kelvins compared to simulations that include other exchange integrals. In addition to the values of exchange interactions just discussed, the second factor that significantly affects the T C values obtained are the magnetic moments. Figure 5 presents the calculated magnetic moments in the ground state (T = 0 K) for the Y 1−x Gd x (Fe 1−y Co y ) 2 system. While the atom-resolved values of magnetic moments do not differ by more than ∼ 0.2 µ B with concentrations, the total magnetic moment varies from about −4 µ B to +4 µ B . This range is due to the antiparallel alignment of moments at the Fe/Co and Gd sites. The calculated values of the magnetic moments are about 2.2 and 1.3 µ B for Fe and Co and about -6.6 µ B for Gd. These results are in good agreement with our previous theoretical results for the systems Y 1−x Gd x Co 2 [13] and Y(Fe 1−y Co y ) 2 [20]. The calculated total magnetic moments are also in good agreement with the experimental results for the systems Y 1−x Gd x Co 2 [13] and Y(Fe 1−y Co y ) 2 [35]. Interestingly, for certain concentrations we observe a zero total moment, implying a complete compensation of the opposite moments present on the different elements. We would like to recall that the results in the upper left corner of the magnetic moment map, describing the nearest neighborhood of YCo 2 , are incorrect because YCo 2 is in fact a magnetically disordered phase whose ground state is not correctly described in GGA.
For the positive (red) region of the total magnetic moment, due to the presence of uncompensated opposite magnetic moments, the relation M (T ) shows a positive slope at low temperatures, see Fig. 4(d). For the example Y 0.8 Gd 0.2 Co 2 concentration, our MC predicted increase in magnetization at low temperatures is confirmed by experimental (zero-field cooled) results [13]. For considered compounds with large positive total magnetic moment, we expect a large inverse magnetocaloric effect to occur in a region with positive M (T ) slope. The physics responsi-ble for the shape of M (T ) curves can be understood by considering exchange integrals for different pairs of atoms. In Fig. 4 we see that the 4f -4f and 4f -3d interactions are much weaker than the 3d-3d interactions therefore the thermal disorder first affects the 4f sublattice, while the 3d sublattice remains (at least relatively) ordered, which means that the magnetic moments of the 4f sublattice no longer compensate, or compensate to a much lesser extent, the 3d moments, hence the initial increase in magnetization. For temperatures higher than the temperature for which the magnetization maximum occurs in Fig. 4(d), the simulated material begins to behave as a ferromagnet. Figure 6 shows the valence band DOS for four stoichiometric compositions of the Y 1−x Gd x (Fe 1−y Co y ) 2 system. All panels of Fig. 6 show results with spin polarization that is proportional to the values of calculated magnetic moments. As we said, the magnetic moments are about 2.2, 1.3, and -6.6 µ B for Fe, Co, and Gd. In each case, the dominant contribution to the valence band comes from the 3d orbitals of the 3d elements (Fe and Co). For compounds with Gd, we observe an almost completely occupied one spin channel of the Gd 4f orbital and an almost empty second spin channel. The 4f bands are also much more localized than the 3d bands. The location of the occupied 4f bands is about -3.5 eV below the Fermi level. The difference in DOS between systems containing Co and Fe is mainly due to filling of the valence band with an extra electron for each Co atom relative to the Fe atoms. Similarly, as we mentioned earlier when discussing the results of Curie temperature calculations, the obtained ferromagnetic solution for YCo 2 should be treated with caution in view of the experimentally found paramagnetic ground state for this phase.

Summary and Conclusions
In this paper, we theoretically described the Curie temperature dependence of the Y 1−x Gd x (Fe 1−y Co y ) 2 pseudobinary system using the combined Monte Carlo and first- principles methods. The calculation results agree well with the previous measurements except for the vicinity of the YCo 2 system. For the latter, instead of paramagnetic ground state, we obtained a ferromagnetic one. For the Y 1−x Gd x Fe 2 subsystem we reproduced the linear dependence of T C on Gd concentration, and for the Y(Fe 1−y Co y ) 2 and Gd(Fe 1−y Co y ) 2 subsystems we reproduced the characteristic Slater-Pauling-like dependence of T C on concentration. The results presented here confirm the ability to efficiently predict Curie temperatures for magnetic systems containing up to three different magnetic elements simultaneously. Furthermore, the first-principles results show that the largest contribution to the Heisenberg Hamiltonian comes from Fe/Co nearest-neighbor exchange interactions. For the Y 1−x Gd x (Fe 1−y Co y ) 2 system, we have also shown how the occurrence of a compensation point with an effective zero magnetic moment depends on the concentration of Co and Gd. In addition, using the example of the Y 0.8 Gd 0.2 Co 2 ferrimagnetic phase, we have shown that the model used allows to predict non-trivial magnetization-temperature dependence.