Realistic first-principles calculations of the magnetocaloric effect: applications to hcp Gd

Abstract

We present an efficient computational approach to evaluate field-dependent entropy of magnetocaloric materials from ab-initio methods. The temperature dependence is reported for the entropy change, specific heat and magnetization for hcp Gd. To obtain optimal accuracy in the calculations, a mixed-scheme for magnetic Monte Carlo simulations is proposed and found to be superior to using pure quantum or classic statistics. It is demonstrated that lattice and magnetic contributions play a role in the entropy change and that the dominating contribution comes from the magnetic contribution. The total calculated entropy change agrees with measurements at room temperature.

GRAPHICAL ABSTRACT

IMPACT STATEMENT

Demonstration of the accuracy of ab-initio theory, coupled to statistical methods, for accurate calculations of the total entropy variation associated with the magnetic transition of Gd. Reproduction of experimental data of entropy change.

Keywords:

• Magnetocaloric
• Gd
• entropy
• Monte Carlo
• mixed statistics

1. Introduction

Historically, Gd-based compounds were in the center of the development of magnetic refrigeration devices operating at room temperature. In fact, metallic Gd was the magnetic refrigerant of the first magnetic refrigerator device, a choice motivated by its Curie temperature ($T_C\approx 290$ K) close to room temperature and a great magnetocaloric response [1]. Another important mark of this field is the observation in 1997 of a giant magnetocaloric effect (MCE) in Gd${}_5$(Si${}_2$Ge${}_2$) by Pecharsky and Gschneidner [2], which captured the attention of the scientific community for the magnetic cooling technology and pushed the field forward. Due to these landmark results, and the availability of multiple studies in the literature, the magnetocaloric community takes bulk Gd as a reference for comparing the performance of the MCE in materials [3–6]. The entropy variation ($\mathrm{\Delta }S$) is the most important parameter to characterize the MCE and the performance of magnetic materials in magnetic cooling devices.

The MCE is a response of an applied field, and the magnetic contribution to the entropy change ($S_{\mathrm{m}\mathrm{a}\mathrm{g}}$) is expected to be the dominant term, although it is established that the lattice ($S_{\mathrm{l}\mathrm{a}\mathrm{t}}$) and electronic ($S_{\mathrm{e}\mathrm{l}\mathrm{e}}$) contributions are also relevant [7–9]. For this reason, the entropy can be accurately described as the sum of three different contributions viewed as independent terms: $\mathrm{\Delta }S=\mathrm{\Delta }{S}_{\mathrm{e}\mathrm{l}\mathrm{e}}+\mathrm{\Delta }{S}_{\mathrm{m}\mathrm{a}\mathrm{g}}+\mathrm{\Delta }{S}_{\mathrm{l}\mathrm{a}\mathrm{t}}$. However, note that for simplicity, couplings between the system's degrees of freedom are not included in this approach, despite their relevance for magnetocalorics.

In this work, we study the MCE in the hexagonal closed packed form (hcp) of Gd, focusing on the calculation of $\mathrm{\Delta }S$ from first-principles calculations and Monte Carlo (MC) simulations. Along with these calculations, we introduce and test the performance of a new scheme for MC simulations, which combines quantum and classic statistics to obtain an improved description over a wide temperature range.

2. Materials and methods

2.1. Mixed-statistics Monte Carlo method

Traditionally, magnetic MC simulations use the Boltzmann distribution to emulate and regulate the thermal fluctuations at finite temperatures. While this approach is suitable to simulate order-disorder transitions, it neglects quantum effects on the fluctuations, thus becoming less accurate for lower temperatures. A well-known example of such failure is the resulting heat capacity of 1 $k_B$ in the limit $T\to 0$ (see Figure 1(a)), which is in accordance with the classical equipartition theorem¹, but in disagreement with experimental observations. At low temperatures, quantum effects become relevant and must be included in fluctuations of any appropriate approach.

Figure 1. Temperature-dependent magnetization derived from MC simulations of hcp Gd using different statistical approaches. Temperature is normalized according to $T_C$ of the experimental (293 K) and calculated (315 K) data. Experimental data extracted from Ref. [25].

With that intent, previous works have developed methods to capture the correct behavior at low temperatures. For example, in Ref. [10], a quantum scaling factor of the probability density at temperature T is introduced so the transition probability between MC states is driven by quantum statistics. Effectively [11,12], this approach corresponds to an effective rescaling of the simulation temperature (T) defined from the average energy of the magnons $〈E_{\mathrm{m}\mathrm{a}\mathrm{g}}〉$, calculated using the density of states of the magnons (MDOS) and Bose–Einstein statistics at a certain T, as: $T_{\mathrm{q}\mathrm{t}}=〈E_{\mathrm{m}\mathrm{a}\mathrm{g}}〉\left(T\right)/k_B$.

Despite a significant improvement at low temperatures (shown in Figure 1), this quantum treatment sharpens the systems behavior around $T_C$, as is evident in the resulting magnetization (Figure 1(b)) and heat capacity peak (Figure 1(a)). Such effect is not exclusive of Gd as similar results were obtained in analogous calculations [10,11,13]. An explanation for these results might come from the use of a quasi-harmonic approximation (see more details in Refs. [10,11]) to calculate a temperature-dependent MDOS from the adiabatic magnon spectra. This approximation should only be valid at low temperatures, while not being accurate close to $T_{\mathrm{C}}$ [11].

Comparing the methods based on classic and quantum statistics, we note that they are suitable in the $T\to T_{\mathrm{C}}$ and $T\to 0$ limits, respectively. However, both have reduced accuracy when used outside these regions. These limitations constitute a problem for simulations over a wide temperature range. In particular, for quantities derived from integration over the whole temperature range, such as the entropy: (1) $S_{\mathrm{m}\mathrm{a}\mathrm{g}}\left(T\right)={\int }_0^T\frac{C_{\mathrm{m}\mathrm{a}\mathrm{g}}}{T}\mathrm{d}T.$(1) To tackle this issue, we propose a simulation scheme that combines quantum and classic statistics to improve the practical issues observed.

Generally, in MC simulations, the sampling of magnetic configurations is performed using the transition probability between MC states (W). This probability depends on the sampling method. Commonly in MC simulations, the Metropolis algorithm is considered, where the state transition probability between configurations is defined as: (2) $W=\{\begin{array}{ll}\exp (-{\displaystyle \frac{\mathrm{\Delta }E}{k_{B}T}}),& \mathrm{i}\mathrm{f}\mathrm{\Delta }E>0\\ 1,& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$(2) with $\mathrm{\Delta }E$ as the energy difference between states.

In the mixed method used here, we combine the transition probability calculated in the quantum scheme, $W(\mathrm{\Delta }E,T_{\mathrm{q}\mathrm{t}})$, with the transition probability calculated in the classical approach $W(\mathrm{\Delta }E,T)$ as a weighted average: (3) $\begin{array}{r}{W}_{mix}(\mathrm{\Delta }E,T)=\alpha W(\mathrm{\Delta }E,T_{\mathrm{q}\mathrm{t}})+(1-\alpha )W(\mathrm{\Delta }E,T)\end{array}$(3) where α is the weight factor. This weight factor should be temperature-dependent so that at low temperatures, the $W(\mathrm{\Delta }E,T_{\mathrm{q}\mathrm{t}})$ quantum contribution is dominant, and near $T_{\mathrm{C}}$ the $W(\mathrm{\Delta }E,T)$ classic contribution becomes dominant. Furthermore, α should be defined by a smooth function to avoid a sharp transition between the two statistics schemes. With this in mind, we propose an α with the following temperature behavior: (4) $\alpha =\{\begin{array}{ll}1-{\displaystyle \frac{T}{T_{\mathrm{C}}}},& T\le T_{\mathrm{C}}\\ 0,& T\ge T_{\mathrm{C}}\end{array}$(4) This choice is the most straightforward function that satisfies the necessary conditions with the advantage of not needing additional free parameters than $T_{\mathrm{C}}$, which is also required for the quantum simulations and can be calculated using classic MC simulations. The factor α also enters in the rescaling of the MDOS used for the pure quantum statistics method [11].

2.2. Computational details

The half-filled 4f shell of Gd is responsible for its large magnetic moment and strongly correlated electronic structure. However, since these orbitals are highly localized, the overlap between neighboring atomic sites is negligible, being the magnetic interactions ruled by the 6s, 6p and 5d states [14]. Despite the high localization of the 4f states, the correct treatment of these orbitals, in first-principles calculations, is a topic of discussion since standard calculations for the ferromagnetic configuration (FM) predict a non-physical hybridization of these states with the remaining valence states near the Fermi energy. Some of the most efficient approaches in the literature are treating the f-states as spin-polarized core states or using the LDA+U correction [14,15].

In the paramagnetic configuration (PM), the magnetization averages to zero while retaining the local magnetic moment. Non-spin-polarized calculations fail to reproduce this picture since they constrain the degree of freedom of the spin. Usually, more sophisticated approaches as the disordered local moments (DLM) approximation are necessary to describe this magnetic state. However, for Gd, it is reasonable to assume that a calculation with the f-states treated as non-spin-polarized core states is a good approximation since the dispersing valence electron states are treated without a polarizing field generated by the highly localized f-states, which is consistent with the zero averaged magnetization of the PM state. Comparison of the density of states (DOS) calculated in this setup against a DLM calculation supports the validity of this approach (see Supplement).

Additionally, we know from the results of similar calculations [15] that in Gd the electronic and magnetic properties are better described by LDA exchange-correlation potentials. In contrast, the structural properties benefit from the use of GGA.

The structural and magnetic properties of hcp Gd were determined from first-principles calculations of the electronic structure using the VASP package [16–19] and RSPt software [20]. The PHONOPY [21] package was utilized to compute the vibrational density of the states. For calculating the electronic structure of disordered magnetic phases, we considered the DLM method by utilizing the SPR-KKR code [22]. From the DFT calculations, the magnetic moments and the $J_{ij}$ exchange couplings (within the Lichtenstein–Katsnelson–Antropov–Gubanov formalism [23]) were determined and used as input for the MC simulations. We performed the magnetic MC simulations using an atomistic Heisenberg Hamiltonian: (5) $\mathcal{H}=-\sum _{i\ne j}{J}_{ij}{\mathbf{m}}_{\mathbf{i}}\cdot {\mathbf{m}}_{\mathbf{j}}-\sum _i\mathbf{H}\cdot {\mathbf{m}}_{\mathbf{i}}$(5) in the UppASD code  [24]. The Hamiltonian in Equation (5(5) $\mathcal{H}=-\sum _{i\ne j}{J}_{ij}{\mathbf{m}}_{\mathbf{i}}\cdot {\mathbf{m}}_{\mathbf{j}}-\sum _i\mathbf{H}\cdot {\mathbf{m}}_{\mathbf{i}}$(5) ) describes the pair exchange interactions between atomic magnetic moments (m) in different sites and the interaction of the magnetic moments with the magnetic field (H), respectively. In the Supplement, we provide more complete details on the setup of the calculations.

3. Results and discussion

3.1. Monte Carlo simulations

The proposed mixed scheme captures the behavior of the quantum scheme at low temperatures and it correctly describes the classical limit showing the same behavior as the classical approach at $T\to T_C$ as can be seen from the MC simulation of the temperature-dependent magnetism, see Figure 1. Such convergence of results from the mixed to the quantum scheme at $T\to 0$ and convergence to the classic scheme at $T\to T_C$ is guaranteed by the mixing weight choice (see Equation (4(4) $\alpha =\{\begin{array}{ll}1-{\displaystyle \frac{T}{T_{\mathrm{C}}}},& T\le T_{\mathrm{C}}\\ 0,& T\ge T_{\mathrm{C}}\end{array}$(4) )).

The success of the mixed method becomes clear by comparing the magnetization curves (Figure 1(b)) with experimental results. In particular, it should be noted that at intermediate temperatures the mixed-scheme calculations are very close to the experimental measurements [25].

Comparing the magnetization obtained from simulations with the experimental data allows to determine the temperature ranges in which the quantum and classic descriptions are more adequate, see Figure 1. From this figure, the existence of an intermediate regime in which the actual magnetization lies between the classic and quantum regimes can be seen. This implies that in this region the origin of the magnetic fluctuations is different, and possibly that both quantum and classic excitations co-exist. The success of the mixed-statistics scheme suggests that it is the latter case, or at least that the experimental observations can be more satisfactorily be reproduced by this description.

Since the primary goal of this study is to calculate the entropy change of hcp Gd when an applied magnetic field is present, it is worth discussing how the proposed scheme compares with the conventional classic statistics (see Figure 2). We calculated the variation of the entropy of the MCE from the heat capacity difference between a simulation without an applied magnetic field and a simulation with an applied field of 2T: (6) $\mathrm{\Delta }{S}_{\mathrm{m}\mathrm{a}\mathrm{g}}(H,T)={\int }_0^T\frac{C_{\mathrm{m}\mathrm{a}\mathrm{g}}(H,T)-C_{\mathrm{m}\mathrm{a}\mathrm{g}}(0,T)}{T}\mathrm{d}T$(6) As Figure 2 shows there is in all three methods a peak of $\mathrm{\Delta }{S}_{\mathrm{m}\mathrm{a}\mathrm{g}}$ at the ordering temperature. This property is of primary interest to this investigation. Although the different approaches used here give different results of the M(T) curves and the specific heat (see Figure 1) we note that the observed variations of the entropy change (Figure 2) are small. This suggests no significant reasons to adopt the mixed scheme described above, over the conventional MC simulations. Nevertheless, it is important to note the existence of substantial changes in the entropy at $T<T_C$, in the quantum and mixed regimes, which can be crucial to consider in the study of phase stability in more complex materials.

Figure 2. Temperature-dependent magnetic entropy variation from MC simulations using different statistics schemes. For a detailed description of the mixed scheme, see text.

3.2. Electronic and lattice entropy contributions

Based on the conclusions from a previous study [26], the electronic entropy, $S_{\mathrm{e}\mathrm{l}\mathrm{e}}$, was calculated from the Sommerfeld approximation while the lattice entropy, $S_{\mathrm{l}\mathrm{a}\mathrm{t}}$, was obtained from full phonon calculations using the ground state properties of each magnetic phase.

It is necessary to consider all intermediate states between the fully ordered and fully disordered magnetic configurations to describe the entropy variation accurately. To achieve this, we determined the magnetization dependence of $S_{\mathrm{e}\mathrm{l}\mathrm{e}}$ and $S_{\mathrm{l}\mathrm{a}\mathrm{t}}$ such that these quantities could be related to the magnetization-temperature curves from the MC simulations (Figure 1). For this purpose, we calculated the Sommerfeld coefficient for intermediate magnetic configurations, by adjusting the net spin moments of DLM calculations to coincide with the data of Figure 1. Doing this, we obtained a continuous value of the Sommerfeld coefficient, as a function of the magnetization, allowing us to calculate $S_{\mathrm{e}\mathrm{l}\mathrm{e}}\left(M\right)$, and from there $S_{\mathrm{e}\mathrm{l}\mathrm{e}}\left(T\right)$. To avoid the appearance of numerical artifacts in $S_{\mathrm{e}\mathrm{l}\mathrm{e}}$, we applied a trend-conserving smoothing filter on the original data, see more details in the Supplement.

We performed a similar treatment for $S_{\mathrm{l}\mathrm{a}\mathrm{t}}$, but in this case, for practical reasons of treating magnetic configurations of intermediate disorder, we made use of the virtual crystal approximation. This approximation is often used to calculate phonon properties in alloys since it simplifies the calculation of the properties of the alloy to compositionally weighted averages of the properties of the pure systems [27]. Similarly to chemically disordered alloys, the intermediate magnetic configurations between the PM and FM phases can be treated as mixing of these phases (FM${}_x$ + PM${}_{1-x}$). Within this picture, we determine $S_{\mathrm{l}\mathrm{a}\mathrm{t}}$ as: (7) $\begin{array}{r}{S}_{\mathrm{l}\mathrm{a}\mathrm{t}}(M,T)=\frac{M\left(T\right)}{M_0}{S}_{\mathrm{F}\mathrm{M}}\left(T\right)+(1-\frac{M\left(T\right)}{M_0})S_{\mathrm{P}\mathrm{M}}\left(T\right)\end{array}$(7) where $S_{\mathrm{F}\mathrm{M}}$ and $S_{\mathrm{P}\mathrm{M}}$ were calculated from phonon density of states of the ferromagnetic spin-polarized and the unpolarized DFT calculations, respectively. Having derived a formalism for the $S_{\mathrm{e}\mathrm{l}\mathrm{e}}\left(M\right)$ and $S_{\mathrm{l}\mathrm{a}\mathrm{t}}\left(M\right)$, we computed the entropy variation using the simulated magnetization curves for the H=0 and H=2 T cases: $\mathrm{\Delta }S=S(H=2)-S(H=0)$. The results are shown in Figure 3. While the electronic contribution now becomes negligible in agreement with the expectations, the lattice contribution is still surprisingly large for a 2nd order transition without structural changes.

Figure 3. Temperature-dependent electronic and lattice contributions for the entropy variation.

3.3. Total entropy variation

Comparison between the distinct entropy contributions, see Figure 4, reveals that the magnetic term is dominant, followed by a lattice contribution that is almost as large. Together with the small electronic contribution, they add up to the total entropy, that has a peak at the transition temperature. When comparing to experimental data, we note a decent agreement: theory gives a value of 7.56 J/kg/K at the ordering temperature (which is $\sim 320$ K in the present calculations) while the measured value is 5.2 J/kg/K, at the experimental Curie temperature (292 K) [6].

Figure 4. Variation of the total entropy change between the PM and FM configuration at a magnetic field change of 2 T. Temperature is normalized according to the entropy peak temperature of the experimental (291 K) and calculated (317 K) data. Experimental data from Ref. [6].

A surprising result of these calculations is the large lattice contribution, and it is tempting to associate the overestimate of the total entropy change, when comparing to observations, to a theory that results in a too large value of $\mathrm{\Delta }{S}_{\mathrm{l}\mathrm{a}\mathrm{t}}$. One can relate this overestimation of $\mathrm{\Delta }{S}_{\mathrm{l}\mathrm{a}\mathrm{t}}$ to the consequence of calculating the entropy from lattice properties at T = 0 K. In Ref. [28], measurements of the elastic constants for Gd from 0 to 360 K for the cases $H=0\text{T}$ and $H=2.5\text{T}$, allow us to examine the importance of finite temperature effects on the elastic properties of the FM phase, and therefore also to estimate $\mathrm{\Delta }{S}_{\mathrm{l}\mathrm{a}\mathrm{t}}$.

Applying the Debye model with the experimental data of elastic constants as input, with and without an applied field and at 300 K [28], we estimate a value of $\mathrm{\Delta }{S}_{\mathrm{l}\mathrm{a}\mathrm{t}}=0.986$ J/kg/K, at the Curie temperature (for a more detailed comparison between our values and experimental estimates, see the Supplement). When added to the calculated electronic and magnetic contributions to the entropy change (gray diamond in Figure 4), we obtain a $S_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ of 5.41 J/kg/K nearly perfect agreement with the observed experimental value. This close agreement highlights the importance of including all the entropy contributions for calculations of the magnetocaloric effect under the risk of underestimating $\mathrm{\Delta }S$ [29] the performance of the material.

Despite some controversy in its calculation and use, the relative cooling power (RCP) is a widely used magnetocaloric materials performance parameter [6]. It can be calculated from the refrigerant capacity (q) [30]: (8) $q={\int }_{T_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}}^{T_{\mathrm{h}\mathrm{o}\mathrm{t}}}\mathrm{\Delta }S\left(T\right)\mathrm{d}T$(8) with $T_{\mathrm{h}\mathrm{o}\mathrm{t}}$ and $T_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}$ approximated by the temperatures for full width at half maximum of $\mathrm{\Delta }S\left(T\right)$ [4,31]. Using this definition on the calculated entropy with experimental elastic constants, we estimate an RCP of approximately 219 J/kg in close agreement with experimental calculations (226.9 J/kg [6]), attesting to the accuracy of the calculated $\mathrm{\Delta }S$ curve.

4. Conclusions

We propose a new scheme for magnetic MC simulations that reproduce the experimental measurements of magnetism of hcp Gd, in a wide temperature range. The success obtained stresses the importance of including quantum and classic fluctuations to describe the magnetic state correctly in temperatures between limits $T\to 0$ and $T\to T_{\mathrm{C}}$, and from there an appropriate value of the magnetic entropy change. Furthermore, based on the simulated temperature-dependence of the magnetization we calculate the entropy of the electronic and lattice subsystems. Calculations of all entropy contributions, in the absence and presence of an external magnetic field, allow to estimate the magnetocaloric effect of Gd, with encouraging agreement with experiments.

Note, that first-principles calculations of the $\mathrm{\Delta }S$ contributions are not a trivial task since lifting approximations imply a drastic increase in required computational resources. In addition, it is not always obvious which finite temperature mechanisms are relevant. For example, finite temperature effects could be included in the lattice entropy calculation through the quasi-harmonic approximation. However, other phenomena like phonon-phonon interactions might be needed for an accurate description. Moreover, there are hysteresis losses in experimental measurements, which one cannot account for in a first-principles approach dealing with ideal systems.

Despite methodological and numerical challenges and the need to introduce approximations, the present work establishes a realistic procedure to calculate from first-principles theory coupled to methods of statistical physics, the entropy variation associated with the magnetocaloric effect. The method is applied to hcp Gd and is demonstrated to accurately reproduce observations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Swedish Foundation for Strategic Research, within the project Magnetic Materials for Green Energy Technology (ref. EM16-0039), STandUPP, eSSENCE, the Swedish Energy Agency (Energimyndigheten) as well as Vetenskapsrådet. The computations performed in this work were enabled by resources provided by the CSC – IT Center for Science, Finland, for computational resources and by the Swedish National Infrastructure for Computing (SNIC) at NSC and PDC centers, partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

Notes

1 At low temperatures, the magnetic fluctuations can only happen in the directions transverse to the magnetization, i.e. the system has 2 degrees of freedom.