Spin-lattice dynamics simulation of external field effect on magnetic order of ferromagnetic iron

Modeling of field-induced magnetization in ferromagnetic materials has been an active topic in the last dozen years, yet a dynamic treatment of distance-dependent exchange integral has been lacking. In view of that, we employ spin-lattice dynamics (SLD) simulations to study the external field effect on magnetic order of ferromagnetic iron. Our results show that an external field can increase the inflection point of the temperature. Also the model provides a better description of the effect of spin correlation in response to an external field than the mean-field theory. An external field has a more prominent effect on the long range magnetic order than on the short range counterpart. Furthermore, an external field allows the magnon dispersion curves and the uniform precession modes to exhibit magnetic order variation from their temperature dependence.


I. INTRODUCTION
Magnetic materials have found diverse applications, such as the magnetic data storage in hard drives and power controlling components in electric machinery. For high temperature applications, materials in the magnetic phase, such as ferritic steels, are suitable for building nuclear fission and fusion reactors which require withstanding neutron bombardment. However, the physical properties of ferritic steels are vulnerable to phase changes, which would severely affect the applicability of their designed functions.
An important characteristics of magnetic materials is the ferro/paramagnetic (FM/PM) phase transition, characterized by the transition temperature which is determined by the quantum exchange interaction between the lattice and spin subsystems. Magnetic order can be classified as long range and short range, determined by the magnetic order at varying separation of classical spins. The change in magnetic order is important to determine the magnetic phase transition; a material is ferromagnetic if the overall magnetic order exists, but is paramagnetic if the order vanishes. Magnetization is a physical quantity describing the long range order as it represents the average atomic spin orientation in a given bulk material. Likewise, spin-spin correlation and effective magnetic field reflect the strength of the short range magnetic order.
The external field effect on magnetization and structural deformation has long been investigated. For example, Holstein and Primakoff 1 have treated ferromagnetic materials in domain scale and expressed the Hamiltonian using the exchange interaction plus the magnetic energy due to an external field. Wojtowicz and Rayl 2 have analysed the external field dependence of magnetization and heat capacity using mean field theory, suggesting that an external field can shift the magnetic phase transition point. Choi et al. 3 have employed the mean field theory to find the field-induced magnetization of iron under fields of tens of Tesla, and reported that the transition temperature between austenite and ferrite state increased with the applied field. Koch 4 has also mentioned the higher shift of phase transition temperature under an external field, and the structure of materials can be controlled by an external field. These theoretical and experimental studies indicated that the variation of magnetic order would result in changes of material phase.
Furthermore, studies of the zero-field magnon dispersion curves of iron, which can also reflect the magnetic order, could be identified as a measure of spin stiffness due to collective thermal excitation of magnons. For example, neutron diffraction scattering has been performed by Shirane et al. 5 on Fe (4% of Si), suggesting the necessity of higher order correction terms in spin stiffness. Experimentally, Lynn 6 has studied the temperature dependence of zero-field magnetic excitation using the neutron coherent inelastic scattering, and provided magnon dispersion relations beyond the Curie temperature. Oguchi et al. 7 have presented the zero-field magnon dispersion relations beyond the Curie temperature based on the coherent potential approximation (CPA) and local spin density functional (LSDF) formalism, with short range magnetic order neglected. Katsnelson and Lichtenstein 8 have performed the LSDF with magnetic correlations in electronic materials considered, showing that such a correlation is essential for modeling ferromagnetic iron. Surprisingly, few of the related studies focus on the field-induced magnon dispersion curves of ferromagnetic iron. This research gap has to be filled for a better description of the external field effect.
In addition to the macroscopic and ab-initio study of magnetic properties, molecular simulation helps determine the macroscopic conditions from microscopic level in the course of atomic vibrations. Another simulation tackling both spin and lattice subsystems has been carried out by Grossmann and Rancourt, 9 who have performed molecular dynamics (MD) on lattice subsystem and Monte Carlo (MC) simulation on spin subsystem on the basis that lattice vibration is more frequent than that of spin. However, it lacks an explicit treatment of direct coupling between the two subsystems. Antropov et al. 10 have attempted an approach that combines finite-temperature spin dynamics to first-principle molecular dynamics. However, only the simulation results at zero temperature has been presented. 11 Ma, Woo and Dudarev (MWD) 11,12 have developed spin-lattice dynamics (SLD) by incorporating these two subsystems on an equal footing, using the exchange integral varying with lattice separations. The atomic spins vary with instantaneous lattice dimensions, and the corresponding magnetic order is more pertinent to real ferromagnetic substances.
The external field effect on ferromagnetic materials is yet to be studied. In view of this problem, this paper studies the field-induced thermo-mechanical properties of body-centred-cubic (BCC) ferromagnetic iron. We focus on the temperature at which the trend of magnetic order starts to stabilize, based on the orientation of atomic spins.

A. Mean field theory
According to the mean field theory (MFT), the temperature-dependent magnetization in ferromagnetic materials can be derived from the exchange field. [13][14][15] The reduced mean field magnetization about the z-component, e z = M(T )/M(0), is equal to where ζ = gμ B H eff /k B T . In addition, H eff is the effective magnetic field strength found by In BCC structure such as ferromagnetic iron in ground state, the total exchange integral can be expressed by the lattice separation R i j in terms of the lattice constant a as with √ 3a/2 representing the first nearest neighbor (1nn) distance and a the second nearest neighbor (2nn) distance in a BCC structure.

B. Magnon dispersion relations and spin stiffness
The magnon dispersion relation is found from the Fourier transform of the autocorrelation function C e of the unit spin vector e, given by ensemble averaging as since e i (q, t) is a real-valued function. Then after the time Fourier transform, one can get where q is the wavevector and ω q is the frequency at q. P can be expressed in terms of q by summing the contributions from all frequencies: The smallest interval of a wavevector for each Cartesian direction is, according to the reciprocal lattice vectors in the first Brillouin zone, where N is the unit cell dimension. In this case, the wavevectors along the three major directions of a reciprocal lattice, i.e. , , and , are : q k = 2π Na (0, k, 0) for k = 0, 1, 2, · · · , N ; : q k = 2π Na (k, k, k) for k = 0, 1, 2, · · · , N 2 ; : The magnon dispersion relations were obtained using 16,384-point Fast Fourier Transform, and the wavevectors have a general form where N 1 , N 2 , and N 3 are unit cell dimensions in x-, y-, and z-directions, respectively, and l 1 , l 2 , and l 3 are positive integers smaller than N 1 , N 2 , and N 3 , respectively. For the simulation boxes used in this paper, The spin stiffness, characterized by the curvature, can be evaluated by fitting the points with the dispersion relation for small wavevector, known as where ω is the magnon energy, is the offset of the magnon dispersion from zero energy, D is the fitted spin stiffness, and q is the norm of the wavevector q. Specifically when q = 0, the spin waves are all in phase. The corresponding mode is known as the uniform precession mode, 16,17 being totally responsible for the magnon energy for q = 0. Also by |q| = 2π/λ, where λ is the wavelength of a spin wave, the magnon energy at zero wavevector corresponds to an infinite wavelength obtainable in an infinite bulk.

C. Spin-lattice dynamics
This technique 11,12 is based on a typical Hamiltonian of ferromagnetic materials: where m i is the mass of atom i, {p i } is the momentum space, {R i } is the position space, and {e i } is the space of classical unit spins, and S is the norm of a spin vector. U ({R i }) is the total lattice potential, and j i j (R i j ) is the exchange integral subsuming the norm S of spins i and j, i.e.
The exchange integral used in SLD is isotropic, implying that the effect of spin correlation is non-directional. In this formalism, the last term represents the energy contributed by an external field. Note in this definition that the atomic magnetic moment is opposite to the classical spin, i.e. the atomic magnetic moment M i = −gμ B S i . The lattice potential adopted in the simulation was the Chiesa-Derlet-Dudarev (CDD) potential U CDD ({R i }), 18 minus the ground state After removing the ground state energy, the lattice potential can reflect the condition of finite temperature.
The phase space trajectory of the SLD technique is evaluated as For equation (17), H eff i is the effective magnetic field experienced by atom i (see equation (21)), and is the reduced Planck constant. The solutions to each of these equations can be obtained from the integration algorithm for SLD, 19 based on Suzuki-Trotter decomposition. [20][21][22][23][24] Note that the spin equation of motion in equation (17) is derived from the Poisson bracket, treating spins quantum mechanically: where

D. Simulation steps
The BCC iron bulks contained 54,000 atoms with spins, arranged in 30 × 30 × 30 unit cells. The lattice constant at ground state was set as 2.8665 Å The bulks were subject to a constant magnetic field along the +z direction, while undertaking thermalization under NPT condition, with the pressure set to be zero. The field ranged from 0 T to 100 T, whereas the thermalization temperature ranged from 300 K to 1,400 K. Langevin thermostats, 25  used for both spin and lattice subsystems. Readers may refer to Ref. 11 for the mechanism of the spin thermostats. Berendsen barostat 28 is used to control the stress of the simulation box. For each atom, its stress is evaluated by the virial theorem. 29 The simulation was carried out with time step of 1 femtosecond (fs). Suzuki-Trotter decomposition was used to solve the equations of motion. The cut-off distance of spin interaction was set as 3.75 Å in order to include the second nearest neighbor atoms during simulation. Periodic boundary condition 29 was adopted to avoid the surface effect. After the equilibrium volume had been reached, the equilibrium lattice constants were determined, which were used to implement simulation in NVT condition for finding the equilibrium magnetization. The atomic magnetization was given by where the angle brackets · · · T eq ,V eq represent the ensemble average obtained at equilibrium temperature T eq and equilibrium volume V eq , n is the number of atoms in a system, S is the magnitude of a classical spin, and e z i is the z-component of the unit spin vector of atom i.

A. Atomic magnetization
Magnetization, characterising long range order, is often expressed within the mean field theory (MFT) by assuming that the exchange field of spins originates from one single external magnetic field -the mean field. In the presence of just one mean field, the MFT has neglected spin correlations along the plane of spin precession, resulting in a less accurate representation of spin correlation. Recently, the spin-lattice dynamics 11, 12 method has addressed this limitation since the precession of each atomic spin can be computed at each simulation time step. Therefore, the effect of spin correlation on the magnetic order can be demonstrated. Figure 1 shows the temperature dependence of reduced magnetization by SLD. The corresponding MFT results are also plotted for comparison. When an external field is applied, the reduced magnetization does not totally disappear. Rather, some finite values of magnetization can be recorded, known as the field-induced magnetization. According to the literature, 30, 31 the approximate Curie temperature under the zero-field condition, T C , can be represented by the inflection points. Then one may consider the inflection points for field-induced magnetization as well, known as the inflection temperature T inflection . Hence T inflection = T C if the external field is absent. For SLD, Fig. 1 returns the inspected inflections at about 1,000 K for 0 T, 1,025 K for 10 T, 1,050 K for 50 T, 1,100 K for 100 T, respectively. Similarly for MFT results, the inspected values are about 1,250 K for 0 T, 1,275 K for 10 T, 1,300 K for 50 T, and 1,350 K for 100 T, respectively. These values show that the inflection temperatures from both SLD and MFT increase with the field. Their deviation from the Curie temperature is more pronounced when the external field increases. One can quantify the discrepancy between SLD and MFT by scaling the MFT effective field H eff with a factor A and then by finding the scale with the best fit to the SLD data, such that Here, e z is the unit spin vector along the z-component, g is the gyromagnetic ratio, μ B is the Bohr magneton, J i j is the exchange integral experienced in atom j due to neighbor atom i, and H ext z is the effective field along the z-component.
Three convergence criteria are chosen to find A: (i) most approximate area under the reduced magnetization curve, (ii) least algebraic sum of residuals between the two curves, and (iii) least squared sum of residuals between the two curves. The fitting results of A against the external field strength are depicted in Fig. 2, which shows that A increases with the external field regardless of the convergence criteria. The result A ≤ 1 implies that the mean field lowers its contribution along the z-component (i.e. the applied field direction) and contributes to the spin precession components instead. Furthermore, we can realize that A → 1 if the external field gets stronger because the uncorrelated portion of the effective field (i.e. due to the external field) dominates the molecular field, as MFT assumes. In other words, MFT overestimates the temperature for the onset of the magnetization stability. The cause of this overestimation is because MFT lacks consideration of the correlation in spin precession planes, but SLD simulation explicitly takes into account the three dimensional Heisenberg model. In addition, both the MFT and SLD results would converge at about zero when the field is well beyond the Curie temperature. This phenomenon is expected because the internal field is disrupted by such a high temperature, making the external field a dominant component of the effective field.
The implication of the scaling factor A between MFT and SLD can be further refined. The factor A can be regarded as the proportion of the mean field magnetization which should lie along the external field direction. Indeed, 1 − A can be regarded as the proportion of the mean field which contributes to the precession direction, in order for MFT to be reduced to SLD. For the proportion 1 − A to occur, the orientations perpendicular to the external field need possess some atomic spin components. The field-induced inflection points deserve more concern because it represents the temperature at which the average magnetization begins to become stable. Figure 3 shows the overall effect of the external field by plotting the reduced magnetization against the reduced temperature, in which are derived by inspecting T inflection in Fig. 1. The experimental data under zero-field condition by Crangle and Goodman 32 are also presented for comparison. The simulated and experimental values are departed because SLD treats spins classically instead of quantum mechanically with just spin-up and spin down states. 11 Apart from of the departure, the simulated reduced magnetization curves have nearly the same temperature dependence trend below about 0.90T inflection . This means that the external field has a weaker effect to change the reduced magnetization below the inflection temperatures. However, as the temperature increases to about 0.90T inflection , of the corresponding field, the reduced magnetization starts to diverge and to maintain above zero in the presence of an external field. The graph also shows that increasing field strength retards the drop of magnetization in the critical region more effectively. In another perspective, we can realize that the field-induced inflection points are analogous to the Curie temperature in the zero-field case that represents the retarding change of magnetization near T inflection .
The total effective magnetic field is an indicator of the short range order, since it is calculated from the spin correlation among the nearest neighboring spins. Fig. 4(a) shows the temperature dependence of the effective magnetic field along the applied field direction. Given the simulation results, the curves tend to converge at about 625 T at 1,400 K, probably as evidence of the remaining short range order that exists under strong magnetic fields (up to 100 T according to this study) in spite of thermal excitation. This result coincides with the computations by Tao et al. 33 that demonstrates the spin wave excitations beyond the Curie temperature. Fig. 4(b) shows the reduced temperature dependence of the effective field at various external field strengths. By adopting the inspected inflection temperatures from Fig. 4(a), the field-induced effective fields exhibit a divergence at about 0.82T inflection as compared to 0.90T inflection for macroscopic magnetization in Fig. 1. Even though the proportionality constant of T inflection is slightly different, Fig. 4(b) suggests that the inflection temperature from the effective field is the temperature beginning to exhibit a levelling trend in short range order.
The findings infer that the external field contributes insignificantly to the effective field. Consider the general expression of the effective field of atom i: The first term on the right hand side of equation (21) represents the internal magnetic field contributed by the temperature-dependent spin moment S j , and the second term is the external field vector H ext . The effective field is thus a vector sum of the internal and the external fields. In Fig. 4, the simulation calculates the total effective field to be thousands of Tesla, significantly larger than any of the external fields tried in this paper. Therefore, the external field investigated in the paper has a much weaker contribution than the internal field in the effective field, despite the relative orientation between the external field and the internal field. In other words, the internal field would dominate the effective field if the external field is not excessively strong. If the spin unit vector e i of atom i is known, the ensemble average of the normalized spin correlation function e i · e j can serve as another indicator of short range magnetic order, 34 which would change with external field and temperature. Fig. 5 shows the reduced temperature dependence of the (a) first-neighbor (1nn) and (b) second-neighbor (2nn) spin-spin correlation functions when an external field of various strengths is present. Here, the inflection temperatures inspected from the effective field curve (see Fig. 4) are employed to calculate the reduced temperatures. The spin-spin correlations exhibit a similar divergence near the inflection points (at 0.96T inflection ) and a convergence beyond the inflection points, just as the average magnetization characterising long range order indicates. Besides, the remaining short range order under strong thermal agitation accounts for the convergence of e i · e j (0.128 for 1nn and 0.0753 for 2nn) at high temperature for all simulated fields.

B. Effect on magnons
The simulation results show that thermal excitation of magnons can be inhibited by an external field, which in turns varies the magnetic order. Fig. 6 shows the magnon dispersion relations at (a) 300 K, (b) 900 K, (c) 1,000 K and (d) 1,100 K, respectively, for four applied magnetic fields. The data points are joined by the fitted lines according to Eq. (12). Since the experimental magnon dispersion curve of pure iron is lacking in the literature, the magnon dispersion relation of Fe (12% of Si) 6 is taken for comparison. An external field can hardly vary the spin stiffness at temperatures below 1,000 K, i.e. before the Curie temperature (T C ) of BCC iron at 1,043 K. It means that the long range magnetic order cannot be reinforced by the magnetic field appreciably. However, at temperatures near or above the critical region of iron (i.e. 1,000 K and 1,100 K), an external field can increase the spin stiffness significantly. It is noted that the data points near the critical region (Fig. 6(c) and 6(d)) are less aligned with the fitted trend of the solid lines. Such a misalignment demonstrates the disappearance of long range magnetic order, as the atomic magnetic moment in Fig. 1 already demonstrates. Another noticeable point is the existence of non-zero magnetic energy at large wavevector when the temperature is high and the external field is absent. Due to the short range order (see the plots for 0 T beyond |q| = 0.6 in Fig. 6(d)), the magnetic energy is still present beyond the Curie temperature of BCC iron (1,043 K). Also, the experimental result 6 shows a slightly smaller curvature (230 meVÅ 2 ) than the simulated results at room temperature (300 K), possibly due to the introduction of impurities (Si) in the experimental settings, as opposed to the adoption of pure BCC iron in the simulations. Compared to Figs. 1 and 4, the temperature dependence of the spin stiffness resembles that of the macroscopic magnetization and that of the effective field. In addition, the presence of an external field is useful to maintain the spin stiffness at elevated temperatures, in particular beyond the inflection points. Fig. 7 shows the temperature dependence of spin stiffness D, according to the magnon dispersion curves in Fig. 6. The data points beyond the zero-field Curie point (1,000 K) do not agree with the fitted curves very well regardless of the field, due to the strong scattering of magnons at this temperature range. On the other hand, the spin stiffness at room temperature (300 K) agrees fairly well with the ab-initio results of You et al. 35 (237 meVÅ 2 ) and Liechtenstein et al. 36 (294 meVÅ 2 ), respectively. The spin stiffness decreases slowly with temperature from 275 meVÅ 2 at 300 K. Then it begins to vanish a temperature that exceeds the zero-field Curie point (T C ).The temperature at which the field-induced spin stiffness value begins to stabilize can be regarded as the transition temperature. One can inspect the inflection points as indicated by the arrows in the inset of Fig. 7. Though the temperature of these points differ from those observed from the atomic magnetic moments (Fig. 1), these values can demonstrate the increase in the transition point with the applied magnetic field. Far beyond the critical region, the external field has inappreciable effect on the spin stiffness values. Instead the stiffness values tend to converge at around 25 meVÅ 2 regardless of the field strengths, probably due to the short range order that remains beyond the inflection temperature.
The curves in Fig. 6 at a fixed temperature seems to be generally shifted upward by the external field, suggesting that the vertical intercept known as the uniform precession mode should also depend  Fig. 6. The Curie temperature T C (1,043 K) is located for comparing the field-induced inflection temperatures with it. In fact, the temperature dependence can be an indicator of the stabilising trend of magnetic order, indicated by the corresponding inflection point (see arrows on the inset) at which the spin stiffness becomes stable at elevated temperatures. The zero-field spin stiffness evaluated by You et al. 35 and Lichtenstein et al. 36 have been marked as comparison to our results at room temperature (300 K).
FIG. 8. Temperature dependence of the uniform precession mode, given external field strengths. The Curie temperature T C of BCC iron is marked as a reference. Data points are from SLD simulations, and the fitted lines indicate the trend of the temperature dependence being another indicator of magnetic order, just as the temperature dependence of the spin stiffness shows us. The theoretical values of the uniform precession mode are marked in dashed lines. The inspected transition points are located by arrows, which show their increase with the external field.
on the field strength. Fig. 8 plots the uniform precession mode varying with temperature for certain field strengths. The uniform precession mode increases with the applied field, and it has a declining trend of temperature dependence similar to that of other quantities. Below the critical region of iron (800 K), the uniform precession mode drops slightly with temperature for all external fields simulated. But the falling trend starts to stabilize at the onset of the critical region. Again, this trend is consistent with that of other figures, indicating the disappearance of long range order and the prominence of the short range order beyond the inflection temperatures. Some negative values are found in Fig. 8 for the zero-field case, which might be attributed by the random error in locating the peak frequency in a scattered magnon spectrum when a temperature exceeds the inflection point.
Since results from one of the spin wave components (q = 0) in a given bulk, it shows the characteristics of spin waves, e.g. the magnetic order. Therefore, the temperature dependence of the uniform precession mode can be regarded as an indicator of magnetic order. Also, the temperature at which the trend levels off can be treated as the inflection temperature of BCC iron at the corresponding external field. They are located in Fig. 8 for exhibiting their increase with the external field.
Theoretically, the uniform precession mode occurs when each spin undergoes precession at the same frequency, known as the Larmor frequency as follows independent of the norm of the classical spin. Guirreiro and Rezende 37 attempted to explain that the magnon energy intercept (See Fig. 7) has a relation to the applied field H ext as according to quantum theory with no temperature effect. Here is the reduced Planck constant. Then one may compare the simulated uniform precession modes with either equation (22) or (23), which are plotted in Fig. 8 as horizontal dashed lines. Before entering the critical region at about 1,000 K, the uniform precession modes from SLD simulations are about 2 meV larger than the theoretical values. But once beyond the inflection point, the simulated values become consistent with the theoretical values. Such a discrepancy of 2 meV below the inflection temperature might be due to the fact that SLD approach to incorporate temperature effect before the inflection point.

IV. CONCLUSION
The paper studies the external field effect on BCC iron in terms of long and short range order using spin-lattice dynamics, a current technique of atomistic simulation of spin-carrying particles. In general, an external field can reinforce both the long and short range magnetic order at elevated temperatures. An external field can resist the drop of average magnetization, so that a higher temperature is required to disrupt magnetic order. With the technique of spin-lattice dynamics by explicitly considering the spin precessions, the external field effect on ferromagnetic iron can be better studied and understood than the conventional mean field theory. On the other hand, an external field can reinforce the spin correlation and the effective magnetic field. However it requires an extremely strong field to bring about the effect on the short range magnetic order controlled by the internal field of the atom. An external field can also exhibit its effect on magnon dispersion, showing an increase of the spin stiffness and a collective upward shift of the magnon dispersion relations. At higher temperatures, an external field can help exhibit the disappearance of long range magnetic order and the prominence of short range order. In addition, the temperature dependence of the field-induced spin stiffness and the uniform precession mode can serve as indicators of the retarding change of magnetic order. In particular, the uniform precession mode found by classical spins can reflect the temperature effect that could not be given by quantum treatment of spins.