Investigating the magnetovolume effect in isotropic body-centered-cubic iron using spin-lattice dynamics simulations

The understanding of the magnetovolume effect lacks explicit consideration of spin-lattice coupling at the atomic level, despite abundant theoretical and experimental studies throughout the years. This research gap is ﬁlled by the recently developed spin-lattice dynamics technique implemented in this study, which investigates the magnetovolume effect of isotropic body-centered-cubic (BCC) iron, a topic that has previously been subject to macroscopic analysis only. This approach demonstrates the magnetic anomaly followed by the volumetric changes associated with the effect, each characterized by the corresponding ﬁeld-induced inﬂection temperature. The temperature of the heat capacity peaks is useful in determining the temperature for retarding the atomic volume increase. Moreover, this work shows the correlation between the effects of temperature and ﬁeld strength in determining the equilibrium atomic volume of a ferromagnetic material under a magnetic ﬁeld. C (cid:2) 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Com-mons Attribution 3.0 Unported License.


I. INTRODUCTION
The magnetovolume effect is the structural deformation induced by the change in magnetization that results from the magnetic phase transition induced by an applied magnetic field. The change in magnetization leads to numerous anomalous phenomena in the physical properties of the lattice structure, especially in ferromagnetic materials, such as body-centered-cubic (BCC) iron. 1-4 The magnetovolume effect caused by magnetic fluctuation is strongest across the ferro/paramagnetic (FM/PM) phase boundary, at which the abrupt change in magnetic order brings about drastic changes in the thermal and elastic properties of Fe-Cr alloys. 5 An example of this effect can be found in quenched face-centered-cubic (FCC) Fe 65 Ni 35 , which has a nearly zero thermal expansion coefficient in a certain temperature region. 6 Materials with this property have dimensional stability that is beneficial for applications such as magnetic shielding, precision instruments, and thermostats in electric water heaters.
A primary concern when investigating the effect is the choice of an approach that models the interaction between the spin and lattice subsystems. Macroscopic study of the magnetovolume effect has been conducted for dozens of years and has mainly suggested the importance of the distancedependent exchange integral used to describe the coupling in material size determination. To list a few results, Lee 3 has stated the interrelation among the average magnetic moment, the volume, the pressure, and the external magnetic field in a thermodynamic context, while Callen et al. 7 have deduced an anisotropic coupling relation between the spins and the external strain in response to the direction of the external magnetic field. Shiga 8 has proposed the distance dependence of the norm of a spin, which indicates the spin-lattice coupling in magnetic materials. Experimental work has also been performed to observe the effect. Tanji 9 has experimentally studied the derivative of the exchange integral with respect to interatomic distance as the key to the anomalous volume change caused by the magnetovolume effect. Oomi and Mōri 10 have observed the effect of stress on the magnetic properties, known as the reverse magnetovolume effect, yet discussion of the resulting magnetic phase transition is lacking.
From a computational perspective, atomic simulation of the magnetovolume effect evaluates macroscopic quantities based on the microscopic condition of the atoms. A notable example of such a case was performed by Grossmann and Rancourt 11 in two steps, with the lattice motion treated in terms of molecular dynamics, and the presumably less-frequent spin motion was modelled using a Monte Carlo technique. However, the direct coupling between the spin and lattice subsystems cannot be considered in this formulation. Skubic et al. 12 employed the model parameters from density functional theory computations and introduced the damping parameters of electron and lattice subsystems to achieve coupling to the spin subsystem. Wang et al. 13 formulated an ab initio scheme of incorporating spin and lattice coupling with low-energy modes in simulating multiferroic materials. Ma, Woo and Dudarev (MWD) 14,15 developed the spin-lattice dynamics (SLD), which mimics atomic behavior by treating the spin and lattice subsystems with equal weight by means of the distance-dependent exchange integral. In this regard, a more realistic and intuitive account of the magnetovolume effect is expected, together with a more reliable determination of the inflection temperature.
The aim of this study is the adoption of the SLD technique for studying the magnetovolume effect of isotropic BCC iron at elevated temperatures. We show that SLD can demonstrate both the volume change associated with the magnetic properties.
To investigate the magnetovolume effect under stress-free conditions, spin-lattice dynamics is adopted. The Hamiltonian of this approach is established, and the equations of motion relating to this Hamiltonian are shown. The NPT ensemble is employed, in which the pressure is set to zero. The field strength is varied between 0 T and 100 T, and the thermalization temperature is varied between 300 K and 1,400 K. After thermal equilibrium is reached, the time averages of observables such as the atomic magnetic moment, the spin-spin correlation functions and the atomic volume are obtained to represent the ensemble averages in accordance with the ergodicity principle.

A. Spin-lattice dynamics
Spin-lattice dynamics (SLD) 14,15 is an improvement of the conventional molecular dynamics technique in solving the equations of motion because of its incorporation of spin degrees of freedom, which are interrelated with the position and momentum degrees of freedom. Therefore, SLD is suitable for simulations of ferromagnetic materials, in which the magnetic properties affect the mechanical properties. The Hamiltonian of ferromagnetic iron used in SLD is shown as In equation (1), m i is the mass of atom i, {p i } is the momentum space, {R i } is the position space, {e i } is the classical unit spin space, and U ({R i }) is the total lattice potential. Additionally, j i j (R i j ) is the exchange integral attributed to the spin separation R i j , 14 whose expression is where j 0 = 904.90177 meV, R c = 3.75 Å, and is the Heaviside unit-step function. It can be observed that j i j is a decreasing function of R i j ; the value of the exchange integral drops with interatomic distance during thermal expansion. Besides, S is the norm of a classical spin. The first term on the right-hand side of equation (1) represents the kinetic Hamiltonian of the lattices. The second term on the right-hand side refers to the lattice potential. The third term on the right-hand side of equation (1) represents the Heisenberg exchange Hamiltonian. The energy attributable to an external field is expressed as the last term on the right-hand side of equation (1), where H ext is the where U CDD ({R i }) is a magnetic potential known as the CDD potential with version number MP-CS3-30, 16,17 and the ground-state magnetic energy (−1/2) is derived from the embedded-atom method (EAM) that parametrizes the potential in terms of interatomic distance, in the form Here, ρ i is the local electronic density of atom i in the form The embedding energy F(ρ i ) is given by A and B in equation (6) are fitting parameters listed in Ref. 16. The first term is the non-magnetic dband energy and the second term is the ferromagnetic band energy. Then, the repulsive pair potential V (r ) and pair density function f (r ) are both fifth-order knot functions in the form Here, N V and N f are number of knot points used for the potential and density functions. The coefficients V n and f n , the positions of the knot functions r V n and r f n , λ V and λ f are the fitting parameters of the MP-CS3-30 parametrization, all of which are listed in Ref. 16.
The exchange integral used in this study is isotropic, such that the effect that arises from it has no preference for any of the lattice axes. In short, the spin subsystem is connected to the lattice subsystem by introducing the exchange integral, which varies with interatomic distance.
The phase-space trajectory of the SLD technique can be evaluated from equation (1) to be: The numerical solutions to the equations of motion rely on the recently developed integration algorithm 18 based on the Suzuki-Trotter decomposition. [19][20][21][22][23] Both the lattice and spin thermostats in SLD are implemented using the Langevin thermostat 24 and the fluctuation-dissipation theorem. 25,26 The implementation of the spin thermostat can be The arrows indicate the inspected inflection points for respective field strengths in both SLD and SD cases. By comparing them, it is evident that both the thermal expansion and spin-lattice coupling contribute to the agreement with the experimental zero-field Curie point.
obtained from Ref. 14. The stress is handled in SLD by adopting the Berendsen barostat, 27 such that the fluctuation of the simulation box volume with temperature can be modeled.

B. Simulation steps
First, the CDD potential 16 was chosen for the simulation. BCC iron bulks of 54,000 atoms with spins were created at 0 K in dimensions of 30 × 30 × 30 unit cells. The lattice constant of BCC iron in the ground state was set to 2.8665 Å. The bulks were exposed to an external magnetic field of a given strength along the +z direction, followed by thermalization. One simulation box was subjected to a fixed temperature and a fixed field strength. Both the lattice and spin subsystems rely on the Langevin thermostat. 24 The time step was 1 femtosecond (fs). The atomic motion was evaluated using the Suzuki-Trotter decomposition. [19][20][21][22][23] The cut-off distance R C for the spin interaction was set to 3.75 Å so that the second-nearest-neighbor atoms could be included in the spin computations. To mimic an infinite bulk, a periodic boundary condition 28 was adopted on each face of the simulation box.
The heat capacity at constant pressure P (zero in this case) is given by 29 where Q eq is the heat absorbed at equilibrium volume at absolute temperature T . The thermal expansion coefficient α at constant pressure P is given by 30 where the temperature-dependent atomic volume V (T ) for BCC metals can be written as 31 where a(T ) is the temperature-dependent lattice constant at a certain field strength.

III. RESULTS
The change in the magnetization, which plays a key role in varying the atomic volume, is investigated first. Figure 1(a) shows the SLD-generated temperature dependence of the average magnetic moment per atom, M, along the applied field direction to express the overall spin collinearity that characterizes the long-range magnetic order. In this graph, the atomic magnetic moment that characterizes the magnetization of BCC iron decreases with increasing temperature for all field strengths. However, a stronger field inhibits the decrease of the magnetic moment at higher temperatures (see the trend beyond the Curie temperature of BCC iron, T C ), leading to a finite field-induced magnetization. This simulation result is consistent with the experimental results that indicate long-range ferromagnetic order beyond the Curie temperature in the presence of an external field. 32 It should be noted that zero magnetization should occur at some temperature far beyond the range shown in Fig. 1(a). The inflection point for the zero-field curves in Fig. 1(a) may be used to distinguish the magnetic phase transition, 33,34 beyond which point the magnetization, or the long-range magnetic order, vanishes. The inspected inflection points for the field-induced magnetization are located as follows: 1,000 K at 0 T, 1,025 K at 10 T, 1,050 K at 50 T, 1,100 K at 100 T, indicating the increase in the inflection temperature with increasing external field.
As a comparison to the SLD results, the temperature dependence of the magnetization produced by the pure spin-dynamics (SD) simulation is plotted in Fig. 1(b). The pure spin-dynamics simulation was undertaken by fixing the ground state lattice constant at 2.8665 Å and by neglecting the contributions of the lattice potential term U ({R i }) in Eqs. (1) and (10), so that the effect of lattice expansion was ignored. In Fig. 1(b), the inflection points inspected from the pure-spin system under 0 T, 10 T, 50 T, and 100 T are about 1,100 K, 1,125 K, 1,155 K, and 1,180 K respectively, which are higher than the spin-lattice coupled counterparts in Fig. 1(a). In particular, the zero-field inflection point in the pure-spin system (1,100 K) is higher than the well-known experimental result of 1,043 K. Without increasing the lattice separation during the increase in temperature, the exchange integral j i j (R i j ) which is a decreasing function would maintain at a larger value than that with thermal expansion, promoting spin alignment. In this case, a higher temperature is required to disrupt the spin alignment more severely, thereby returning a higher inflection point. The comparison thus demonstrates that the lattice component should contribute to the closeness of the zero-field inflection point with the experimental value of 1,043 K. In other words, volumetric variation is vital for modelling the inflection points of magnetization, at least for the zero-field condition.
An external field not only varies the average magnetization but also changes the temperature dependence of the spin-spin correlation functions e i · e j that characterize the short-range magnetic order. 35 The temperature dependence of the spin-spin correlation functions due to the first, second, and third nearest neighbors are depicted in Fig. 2(a) to 2(c), respectively. The Curie point of zerofield BCC iron (1,043 K) is marked, and the figure also shows the existence of the short-range order (SRO) in the presence of an external field beyond the Curie point. All three plots exhibit similar inflection points by inspection: 1,000 K for 0 T, 1,030 K for 10 T, 1,050 K for 50 T, and 1,100 K for 100 T. Accordingly, the spin-spin correlation functions exhibit similar temperature dependence to that of the average magnetization in Fig. 1(a), indicating the consistency between the long-range and short-range magnetic order in the expression of spin collinearity. The short-range order that exists under strong thermal agitation accounts for the convergence of e i · e j to a non-zero value beyond the Curie point for all fields that is seen in this work, which has been confirmed experimentally in ferromagnetic materials. 36 The similarity between the second-nearest neighbor and third-nearest neighbor correlation functions results from the computational model that just considers a cut-off distance of the exchange integral up to the second-nearest neighbors. 14 In addition to macroscopic magnetization, the heat capacity is another indicator of the magnetic order; it reaches its maximum at the onset of disorder of spins. The heat capacity is obtained by equation (12) after the atomic energy has been calculated. Figure 3 shows the heat capacity of isotropic BCC iron under stress-free conditions and under various field strengths ranging from 0 T to 100 T. The temperature of the heat capacity peaks increases as the applied field increases, yet the peak values of the heat capacity appear to be suppressed by a stronger field. From Fig. 3, the anomaly temperatures determined from the abrupt changes in heat capacity near 1,000 K at 0 T, 1,010 K at 10 T, 1,040 K at 50 T, and 1,065 K at 100 T. When these values are compared with those obtained from the inflection temperatures in Figs. 1 and 2, it is clear that the anomaly temperature in the heat capacity curve (see Fig. 3) can also represent the temperature of magnetic order stabilization. The introduction of a distance-dependent exchange integral to the Hamiltonian of ferromagnetic iron implies a change in atomic volume caused by a variation in magnetic properties. Fig. 4 demonstrates this relation by plotting the atomic volume V atom against the absolute temperature for a number of external field strengths. It is shown that the atomic volume exhibits a generally linear relation with temperature, with an extrapolated offset of approximately 11.77 Å 3 at zero temperature, corresponding to the same ground-state atomic volume for all fields simulated in this study. The volumetric anomaly for the zero-field curve occurs in a similar temperature region as in the atomic volume curve derived from the zero-field lattice constants of BCC iron measured by Ridley and Stuart. 37 In addition, the offset agrees with the ground-state atomic volume obtained from ab initio calculations by Friák et al. 38 and Ekman et al. 39 Two features can be observed in Fig. 4. First, the simulated magnetovolume effect seems to be limited, in the sense that the simulated volumetric anomaly is less pronounced than in the experimental results, which might be caused by the stiff nature of the adopted interatomic potential. Magnetic potentials of this kind focuses on reproducing the anomaly temperature, rather than the trend after the anomaly. Second, the inhibited increase in the simulated atomic volume with increasing field strength can be more clearly observed in Fig. 4 near the critical region of BCC iron (near 1,043 K), which can be attributed to the SLD technique, which relates magnetic properties to volumetric properties.
Alternatively, the temperature dependence of the thermal expansion coefficient α, which can be derived from the atomic volume data used in Fig. 4, is plotted in Fig. 5. This quantity is evaluated according to equation (13) after the equilibrium atomic volume has been computed by equation (14). Since the experimental atomic volume remains roughly unchanged just after reaching the anomaly (see the flat portion in Fig. 4), the resulting thermal expansion coefficient returns a rather small value compared to the simulated SLD results. On the other hand, the experimental thermal expansion coefficients were derived from the zero-field lattice constants recorded in Ref. 37 for comparison. It can be seen that the simulated zero-field α at 300 K returns a value of 35 × 10 −6 K −1 , similar to that derived according to Ref. 37. In this graph, a stronger field suppresses the rising trend of the thermal expansion coefficient and increases the values of the anomaly temperature located at approximately 950 K for 0 T, 1,000 K for 10 T, 1,025 K for 50 T, and 1,040 K for 100 T. These values generally coincide with those obtained from Figs. 1(a) and 2. In fact, the peak of each curve in Fig. 5 indicates the magnetic anomaly temperature under a certain applied field, and these values are consistent with the inflection temperatures in Figs. 1(a) and 2 and the anomaly temperatures in Fig. 3.

IV. DISCUSSION
This study has demonstrated the atomic volume anomaly associated with the magnetic properties of BCC iron using the spin-lattice dynamics technique. The implications of the results are discussed in detail below.
The success of the lattice potential is attributed to its inclusion of the directional spin degrees of freedom, which is able to address the thermal disorder of magnetization and the interaction between spin waves and lattice vibrations. 14 Also, the adopted potential has undergone checking of the third-order elastic constants, 16 adding its credibility of describing thermal expansion. The closeness of the present calculations to the experimental values can be regarded as a demonstration of the applicability of this simulation technique in handling thermodynamics of BCC iron. Besides, the success of describing the magnetization is attributed to the employed exchange integral, which is parametrized by fitting two sets of independent ab initio data, 40,41 considering up to 4-nearestneighbor (4nn) interactions. Even though the adopted exchange integral is a computationally pairwise one, it has absorbed the effect due to more distant atoms.
In the same way, the agreement between the calculation and the experimental data should result from the chosen interatomic potential and exchange interaction used in the Hamiltonian in equation (1). First, the adopted lattice potential is a practical treatment of the possible anharmonicity of BCC iron during thermal expansion, an issue worried about by the radiation damage community. The major aim of this potential is the reproduction the experimental bulk properties of BCC iron, including the lattice constant and elastic constants. 17 Second, according to the authors of the spinlattice dynamics technique, 14 the data used for fitting the exchange integral are obtained from ab initio density functional theory (DFT) in the local spin density approximation (LSDA), together with the tight-binding linear-muffin-tin-orbital (TB-LMTO) method in the atomic sphere approximation (ASA). [40][41][42] The difference between the method by Morán et al. and that by Sabiryanov et al. is that the former uses the frozen-magnon method, and that the latter uses scattering matrices and Green's function. In view of the methods of determining the lattice potential and the exchange integral, both empirical and theoretical aspects have been considered in the foregoing simulations.
The inflection points in Figs. 1 and 2 can be treated as the stabilizing points of decreasing magnetization. Because the interatomic distance is related to the magnetic properties of ferromagnetic materials, the inflection temperature of the magnetization should correspond to the retardation point of the increase in atomic volume with further temperature increase.
The temperature dependence of the heat capacity shown in Fig. 3 can be explained as follows. It is known from the thermodynamic relation that C P = T ( ∂ S ∂ T ) P , where S is the entropy related to the spin-lattice-coupled phase space. Accordingly, the peak of the heat capacity occurs at the point of the largest change in entropy with respect to temperature. The conclusion that an external field delays the appearance of the largest entropy, i.e., the lowest magnetic order, to a higher temperature can be deduced. Additionally, the reduction of the peak C P value with an applied field results from the enhancement in the maintenance of magnetic order caused by the field and hence from the reduction in the amount of heat absorbed, which disrupts the spin system.
From Figs. 1-3, it can be seen that the inflection temperature that characterizes the magnetic phase transition is an important factor that reflects the response of ferromagnetic iron to an external field during the course of thermal expansion. The findings concerning both the atomic volume (Fig. 4) and the thermal expansion coefficient (Fig. 5) suggest that an external magnetic field suppresses thermal expansion because of the phononphonon interaction resulting from the anharmonicity of the interatomic potential. 43 The discrepancy between the simulated and experimental α in Fig. 5, especially beyond the critical region, most likely arises because the adopted form of the magnetic potential fails to fully describe the features of the order parameter beyond the inflection point.
It can be inferred that a magnetic effect should be responsible for the trends of V atom and α because Figs. 1 to 5 display anomalies at approximately the same temperatures. The success of modelling these phenomena is attributed to the SLD simulation approach, which links both lattice and spin degrees of freedom.
Figs. 4 and 5 illustrate the compromise between temperature and the external field in the magnetovolume effect at the microscopic level, which can be explained by using Figs. 1 and 2 together with the Hamiltonian in equation (1). First, the atoms gain kinetic energy during thermal agitation, leading to the increase in the kinetic Hamiltonian. Simultaneously, the short-range spin alignment is reinforced by increasing the magnetic field strength such that the spin-spin correlation function e i · e j approaches its maximum value of 1. The increase in spin-spin correlation enforces the short-range magnetic order, thereby lowering the Heisenberg exchange Hamiltonian in the third term of equation (1). However, the presence of an external field preserves the long-range magnetic order and increases the average magnetization, leading to a decrease in the Hamiltonian component that is contributed by the external field. In this manner, the total Hamiltonian of the spin subsystem in equation (1) is reduced. As a result, the restoring force acting on the atoms decreases in magnitude, followed by a shortening of the interatomic distance R i j , which leads to bulk shrinkage in an external field. This shortened interatomic distance enhances the exchange coupling j i j due to its decreasing nature, which further lowers the Heisenberg exchange Hamiltonian in the third term of equation (1). The bulk continues to shrink because of the decreasing total Hamiltonian until the kinetic Hamiltonian that is responsible for thermal expansion can no longer compensate for the magnetic shrinkage.
An experiment regarding the magnetovolume effect in ferromagnetic materials has been explained in terms of the disappearance of the magnetic domains beyond the Curie temperature, 44 accompanied by a prominent change in volume. Because a typical magnetic domain of ferromagnetic iron has a dimension on the order of 1 μm to 100 μm, much larger than the simulation box dimension used in this work (approximately 8.5 nm), the magnetovolume effect at the atomic level that is demonstrated in this study might be the underlying mechanism at larger scales.

V. CONCLUSION
Magnetovolume effect can be demonstrated by the recent simulation approach of spin-lattice dynamics (SLD), even though it appears limited due to the stiff interatomic potential used in this work. Firstly, the magnetic anomaly, reflected by the inflection temperature, can be identified through the Hamiltonian used in SLD, which has incorporated spin degrees of freedom in addition to those for position and momentum. Secondly, it is realized from the various results that the magnetic anomaly corresponds to the volumetric anomaly at the corresponding anomaly temperature for a given field, implying that SLD is able to demonstrate the phenomenon of the magnetovolume effect arising from variation of magnetic properties. Besides, it is the compromise between thermal expansion and magnetic shrinkage that is responsible for the magneto-volume effect in thermal equilibrium.

ACKNOWLEDGMENTS
C. P. C. wishes to thank the Hong Kong Polytechnic University for providing CUDA computing servers.