Unveiling a Scaling and Universal Behavior for the Magnetocaloric Effect in Cubic Crystal Structures: A Monte Carlo Simulation

The magnetocaloric effect and the universal character for the magnetic entropy change regarding the cubic crystal structures (SC, BCC, FCC) were investigated, in a qualitative way, using Monte Carlo simulations. A classical Heisenberg Hamiltonian with nearest neighbors, and next nearest neighbors interactions was implemented. In order to compute the critical temperature of the system depending on the coordination number, it was calculated the dependence of the magnetization and magnetic susceptibility as a function of temperature. Magnetic field dependence on the magnetization for isothermal processes was performed considering a magnetocrystalline anisotropy term. In this way, the magnetic entropy change (ΔSm) was computed. Results show that the rescaled ΔSm as well as the exponent (n) characterizing the field dependence of the magnetic entropy change curves, collapse onto a single curve for the studied crystal structures. By this reason, it can be assured that ΔSm exhibits a universal behavior regarding the strength and contribution of the magnetic exchange energy to the total magnetic energy.

The magnetocaloric effect (MCE) consists in the thermal response of a material when subjected to a magnetic field change and it is an intrinsic property of all magnetic materials. MCE in magnetic systems with second-or first-order phase transitions, near room temperature, has attracted significant interest in recent years due to their potential applications on magnetic refrigeration 1,2 . Magnetic refrigeration techniques have a larger efficiency and a smaller footprint on the environment, compared to conventional compression/expansion of gas techniques. Therefore, it is important to properly describe, both theoretically and experimentally, the field and temperature dependence of the magnetic entropy change (ΔS m ) in magnetic materials. Materials suitable for magnetic refrigeration devices should operate at low, or moderate magnetic fields, in a range (1 − 2T) 3 . Perovskites are materials of great interest since they exhibit a large MCE. This type of materials crystallize into a distorted form of a body-centered cubic (BCC) crystal structure 4,5 . Likewise, there are other materials which crystallize into a cubic type structure, and that also exhibit MCE [6][7][8] . Both, experimental and computational studies, have been performed in order to comprehend the MCE in these materials [9][10][11][12][13] .
In recent years, it has been invested a great effort on correlating ΔS m with the critical exponents of the ferromagnetic-paramagnetic phase transition [14][15][16][17] ; then, by knowing the denominated phenomenological universal curve for ΔS m , it is possible to predict the behavior of ΔS m for ferromagnetic alloys of the same compositional series or for the same material at different applied magnetic fields.
There are several cases where standard techniques like the Kouvel-Fisher method do not work properly 18,19 . Therefore, scaling laws for ΔS m have been proposed for determining the critical exponents of a material 2,[14][15][16][17][18]20 . The theoretical background of the MCE universal behavior relies on the scaling of second-order phase transitions. However, from a theoretical point of view, obtaining the universal curve requires the knowledge of the critical exponents of the material and its equation of state, which is a situation rather unlikely, especially when first undergoing with a new material in the laboratory 16 . Franco et al. 20 proposed a phenomenological procedure to describe the universal behavior of ΔS m in compounds with second-order phase transition, which can be followed without this a priori knowledge. In this way, it is possible to employ the universal curve in the characterization of new materials, such as a screening procedure of the performance of materials, as well as a method for making extrapolations to temperatures or fields not available in the simulated range or in the laboratory.
For materials of the same compositional series, it is likely to obtain different crystal structures and different exchange energy constants by slightly changing its composition. Then, the contribution of this work is focused on the study, by means of Monte Carlo simulations, of the universal behavior and scaling laws of ΔS m for different cubic crystal structures with different coordination number, and to relate them by means of the exponent of the field dependence of ΔS m .
This work is presented as follows: In section II, we describe the model and method, the thermodynamics of the magnetocaloric effect as well as the formulas to compute the physical quantities of interest. In section III, we present the numerical results and discussion and finally, in section IV, we present the conclusions.

Model and Method
Three cubic structures were studied: simple cubic (SC), body-centered cubic (BCC), and face-centered cubic (FCC) lattices. Periodic boundary conditions were imposed to the system. The size of each lattice was chosen so that all the simulated structures had approximately the same amount of ions, because the magnetocaloric properties depend on the size and the amount of ions within the structure 21 . The total number of ions was set at N ≈ 4000, which corresponded to a system size of L × L × L, where L was fixed at 16, 13, and 10 unit cells (uc), being N = 4096, 4394, and 4000 ions, for SC, BCC, and FCC lattices, respectively.
The magnetic behavior of the system was simulated using a classical Heisenberg Hamiltonian, which included terms for the exchange interaction of the magnetic moments with their nearest neighbors (NN) and their next nearest neighbors (NNN), magnetocrystaline anisotropy, and the interaction of the magnetic moments with an external magnetic field. It is worth to mention that, regarding the crystal structure, the coordination number (z) takes the values of 6, 8, and 12 NN for the SC, BCC, and FCC lattices, respectively. Regarding the NNN interactions, each structure has a total number of 12, 6, and 6 NNN for the SC, BCC, and FCC lattices, respectively. The Hamiltonian can be expressed as where S i and S j are the directions of the magnetic moments in the sites i and j, respectively, with |S i | = |S j | = 1.
Note that the first sum runs over the NN, and the second one runs over the NNN, where J NN and J NNN are the respetive exchange constants, with J NN = 1.0. H is the magnetic field intensity, k is the unitary vector in the z-direction, k an is the anisotropy constant, and n i is the anisotropy vector, which was chosen in a perpendicular direction with respect to the applied magnetic field, x-direction. The simulations were carried out using V EGAS 22 , which is an open source software package for the atomistic simulations of magnetic materials employing the Monte Carlo method based on the Metropolis algorithm. The simulations were accomplished using 1 × 10 4 Monte Carlo steps (MCS), rejecting the first 0.5 × 10 4 MCS for relaxation. In order to compute error bars, the simulations were performed using 5 different initializations.
The magnetic susceptibility (χ) was calculated as where k B is the Boltzmann constant and T is the absolute temperature. The critical temperature (T c ) was calculated as the temperature at which the magnetic susceptibility exhibits a maximum value. The total magnetization was calculated as i N i and the total normalized magnetization was evaluated as The change in the magnetic entropy, using the Maxwell relations, can be expressed as 23 where H f is the superior limit of the applied magnetic field intensity. In the case of a discrete field, ΔS m can be approximated by 24 m H H 1 2 In order to compute ΔS m employing the Eq. 6, it is necessary to calculate the magnetization of the material as a function of the applied magnetic field at small discrete steps, for several isothermal processes.
www.nature.com/scientificreports www.nature.com/scientificreports/ The aforementioned method proposed by Franco et al. 20 consists in normalizing the ΔS m curves with respect to their maximum and rescaling the temperature axis in a different way, by means of a variable θ, below and above T c , by imposing that the position of two reference points in the curve corresponds to θ = ±1, as follows where T r 1 and T r 2 are the temperatures of the two reference points. For this study, they have been selected as those where ΔS m pk is the value of the maximum of the |ΔS m | curve.
It can be assumed that the field dependence of ΔS m follows the power law with the field as 20 being n a parameter which depends on temperature and can be locally calculated as At low temperatures, well below T c , n should have a value that tends to 1, which indicates that although the magnetization curves depend on temperature, at these temperatures, this dependence is essentially field independent. At temperatures well above T c , n should tend to 2 as a consequence of the Curie-Weiss law. At T = T c , n shows a global minimum, corresponding to the value of n(T c ). Oesterreicher and Parker predicted that n(T c ) = 2/3 for the mean field case 16,25 . However, for other cases such as Heisenberg model, the value of n(T c ) is different from the one deduced for the mean field and it is closely related to the critical exponents of the model and material 20 . It was found that the value of n(T c ) is equal to Results and Discussion Figure 1 shows the temperature dependence of the magnetization and magnetic susceptibility curves for SC, BCC, and FCC crystal lattices taking into account the interactions of NN without anisotropy, NNN, and NN with anisotropy, respectively. It can be observed that, as the coordination number increases, the critical temperature increases. This is due to the fact that the exchange energy increases as the coordination number increases and, thus, more thermal energy is required to break the magnetic ordering imposed by the exchange interactions and to reach the paramagnetic regime. It can also be seen that the peak of the susceptibility curves decreases when the interaction of the NNN (Fig. 1b) is taken into account. This behavior could be explained from the fact that the role of the NNN interactions is to increase the magnetic ordering in the system by increasing the exchange energy contribution to the total magnetic energy, thus, reducing the thermal fluctuations of the magnetization. As a consequence, even more thermal energy is required to break the magnetic ordering imposed by the exchange interactions. Moreover, in Fig. 1c, it can be observed that the peak height of the magnetic susceptibility is greater when considering the magnetocrystalline anistropy term in the hamiltonian. Also, the error bars are larger, meaning that thermal fluctuations gain significance since there is a competition between both exchange and anisotropic energies.
Isothermal simulations of the magnetization as a function of the applied magnetic field were performed in order to compute ΔS m . The insets of Fig. 2a,b show ΔS m as a function of the temperature for SC, BCC, and FCC structures, taking into account the interactions of the ions with their NN without anisotropy, NNN, and NN with anisotropy, respectively. It can be noted that the values of ΔS m pk were found to be around T c , as it is expected. This is because near T c , the magnetic moments pass from having a magnetic structure, governed by the exchange interactions between ions, to be in a paramagnetic regime, which is a disordered magnetic configuration. As a result, −ΔS m has a maximum due to the abrupt change in the magnetization of the material. On the other hand, the value of ΔS m pk decreases as the coordination number increases. Also, the values of ΔS m pk are smaller when considering NNN interactions, when compared with NN interactions. This is because, since H f was fixed at the same value for the −ΔS m curves = . J ( / 1 0) www.nature.com/scientificreports www.nature.com/scientificreports/ cant as the structure has lower coordination number, or if the interactions with the next nearest neighbors are not considered. As a result, the contribution of the exchange energy will be lower within the structure. Moreover, below T c , −ΔS m increases monotonically and, for values above T c , −ΔS m decreases monotonically, which is in agreement with the standard behavior of magnetic refrigerants.
The main plot of Fig. 2a-c show the universal behavior of the magnetocaloric effect for different cubic crystal structures, considering the interactions of NN without anisotropy, NNN, and NN with anisotropy, respectively. It is possible to observe that, for every figure, all the rescaled curves converge onto the same curve. This means that the second-order phase transition in these studied crystal structures are driven by the same phenomenological behavior. In this way, there is not only a universality and scaling laws for ΔS m regarding the magnetic field, as proposed by Franco et al. 3,[15][16][17][18]20,26 , but there is also a universality and scaling laws for materials regarding the contribution and the strength of the magnetic exchange energy to the total magnetic energy, as long as these ferromagnetic alloys exhibit a second-order phase transition. Furthermore, it can be seen that the master curve of ΔS m has a few discrepancies at low temperatures. This can be explained from the fact that, at low temperatures, the magnetic behavior of the system is mainly driven by the exchange interaction. Figure 3a-c show the scaling and universal behavior for the exponent characterizing the field dependence of ΔS m , considering the interactions of the ions with their NN without anisotropy NNN, and NN with anisotropy, respectively. It can be seen that, at temperatures well below T c , n → 1. This indicates that, although the magnetization curves depend on temperature at these temperatures, this dependence is essentially field independent. Furthermore, for values well above T c , n → 2. This indicates that the change in the magnetic entropy has a quadratic field dependence, as a consequence of the well known Curie-Weiss law 27 . For T = T c , n → 0.639 (solid horizontal line in Fig. 3). This result is in agreement with the critical behavior reported, using the values of the critical exponents β and γ for the Heisenberg model. The aforementioned values were found to be β = 0.367 and γ = 1.388 26 . Therefore, using Eq. 10, n(T c ) = 0.639. It can also be seen from Fig. 3b that the scaling and universal behavior of n holds up even when considering the interaction of the ions with their next nearest neighbors. Figures 2c and 3c show the universality and scaling behavior for the temperature dependence of ΔS m and n, respectively. The simulations were carried out rejecting the interactions of the magnetic moments with their NNN (J NN /J NN = 0.0), and considering the influence of the magnetocrystalline anisotropy term = .
The www.nature.com/scientificreports www.nature.com/scientificreports/ uniaxial anisotropy axis was chosen in the x-direction, perpendicular to the applied magnetic field, with the aim of avoid blocking of the magnetic moments at low temperatures. The mentioned figures show that both temperature dependence of ΔS m and n keep the universality and scaling behavior. Therefore, as the anisotropic term is weak in energy, for this systems, this term does not play an important role, and does not affect the universality nor the scaling behavior studied in this work.

Conclusions
The universal behavior of the magnetocaloric effect regarding the crystal structure, related with the strength and contribution of the magnetic exchange energy, was investigated by means of Monte Carlo simulations using a Heisenberg model with nearest neighbors and next nearest neighbors interactions. Simulations of the magnetization as a function of the applied magnetic field were carried out at isothermal processes to find the temperature dependence of −ΔS m for different values of the magnetic applied field, for SC, BCC, and FCC structures. It was found that the value of ΔS m pk decreases as the coordination number increases. It was showed that the rescaled −ΔS m curves for different cubic crystal structures follow the same behavior, meaning that there exists a scaling law for −ΔS m for the crystal structures investigated in this work. Finally, the temperature dependence of n was obtained for each structure and different anisotropy constants, moreover, the rescaled n(T) curves lies in the same curve, and the exponent n(T c ) is in agreement with the critical exponents for the used model,confirming the universality and scaling behavior of −ΔS m for the studied magnetic systems.

Data Availability
The datasets generated during the current study are available from the corresponding author on reasonable request.