Collective dynamics in atomistic models with coupled translational and spin degrees of freedom

TL;DR

This study uses combined molecular and spin dynamics simulations to explore how phonons and magnons interact in ferromagnetic iron, revealing shifts in frequencies, reduced lifetimes, and new magnetic excitations due to their coupling.

Contribution

It introduces a coupled atomistic model for simultaneous simulation of translational and spin dynamics, providing new insights into their mutual influence in ferromagnetic materials.

Findings

Phonons and magnons influence each other's frequency spectra.

Coupling reduces the lifetimes of vibrational and magnetic modes.

Additional longitudinal magnetic excitations are observed with phonon frequencies.

Abstract

Using an atomistic model that simultaneously treats the dynamics of translational and spin degrees of freedom, we perform combined molecular and spin dynamics simulations to investigate the mutual influence of the phonons and magnons on their respective frequency spectra and lifetimes in ferromagnetic bcc iron. By calculating the Fourier transforms of the space- and time-displaced correlation functions, the characteristic frequencies and the linewidths of the vibrational and magnetic excitation modes were determined. Comparison of the results with that of the stand-alone molecular dynamics and spin dynamics simulations reveals that the dynamic interplay between the phonons and magnons leads to a shift in the respective frequency spectra and a decrease in the lifetimes. Moreover, in the presence of lattice vibrations, additional longitudinal magnetic excitations were observed with the same frequencies as the longitudinal phonons.

Equations

$$\mathcal{H}=\sum _{i=1}^N\frac{m{{\mathrm{v}}_i}^2}{2}+U(\{{\mathbf{r}}_i\})-\sum _{i<j}{J}_{ij}(\left\{{\mathbf{r}}_k\right\}){\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$\frac{d{\mathbf{r}}_i}{dt}$$

$$\frac{d{\mathbf{v}}_i}{dt}$$

$$\frac{d{\mathbf{S}}_i}{dt}$$

$$U(\{{\mathbf{r}}_i\})=U_{\text{DD}}-E_{\text{spin}}^{\text{ground}},$$

$${\rho }_n(\mathbf{r},t)=\sum _i\delta \left[\mathbf{r}-{\mathbf{r}}_i(t)\right].$$

$$F_{nn}(\mathbf{q},t)=\frac{1}{N}⟨{\rho }_n(\mathbf{q},t){\rho }_n(-\mathbf{q},0)⟩,$$

$$S_{nn}(\mathbf{q},\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{F}_{nn}(\mathbf{q},t)e^{-i\omega t}dt,$$

$$F_{vv}^{L,T}(\mathbf{q},t)=\frac{1}{N}⟨{\mathbf{\rho }}_v^{L,T}(\mathbf{q},t)\cdot {\mathbf{\rho }}_v^{L,T}(-\mathbf{q},0)⟩,$$

$${\mathbf{\rho }}_s(\mathbf{r},t)=\sum _i{\mathbf{S}}_i(t)\delta \left(\mathbf{r}-{\mathbf{r}}_i(t)\right).$$

$$F_{ss}^k(\mathbf{q},t)=\frac{1}{N}⟨{\rho }_s^k(\mathbf{q},t){\rho }_s^k(-\mathbf{q},0)⟩,$$

$$F_{ss}^L(\mathbf{q},t)=F_{ss}^z(\mathbf{q},t),$$

$$F_{ss}^T(\mathbf{q},t)=\frac{1}{2}\left(F_{ss}^x(\mathbf{q},t)+F_{ss}^y(\mathbf{q},t)\right).$$

$$S(\mathbf{q},\omega )=\frac{I_0{\Gamma }^2}{(\omega -{\omega }_0{)}^2+{\Gamma }^2},$$

$$S(\mathbf{q},\omega )=I_c\exp \left(-{\omega }^2/{\omega }_c^2\right)+\frac{I_0{\Gamma }^2}{(\omega -{\omega }_0{)}^2+{\Gamma }^2},$$

$${\omega }_{ij}^{-}({\mathbf{q}}_i\pm {\mathbf{q}}_j)=\omega ({\mathbf{q}}_i)-\omega ({\mathbf{q}}_j),$$

Full text

I Introduction

For decades, dynamical simulations of atomistic models have played a pivotal role in the study of collective phenomena in materials at finite temperatures. Molecular dynamics (MD) Rapaport (2004); Frenkel and Smit (2002) utilizing empirical potentials has been extensively used in the analysis of vibrational properties in a variety of systems such as metals and alloys Petry et al. (1991); Heiming et al. (1991); Sun and Murthy (2006); Meyer and Entel (1998); Branicio et al. (2003), polymers Henry and Chen (2008); Lopez Navarrete and Zerbi (1991), carbon nanotubes Shiomi and Maruyama (2006), graphene Koukaras et al. (2015) etc. With regard to magnetic excitations, the lesser-known spin dynamics (SD) method Landau and Krech (1999); Tsai et al. (2000); Chen and Landau (1994); Gerling and Landau (1990); Watson et al. (1969) has proven to be an indispensable tool for investigating classical lattice-based spin models for which the analytical solutions are intractable. Over the years, SD simulations have expanded our understanding of spin waves and solitons in magnetic materials, leading to a number of groundbreaking discoveries, including the existence of propagating spin waves in paramagnetic bcc iron Tao et al. (2005), presence of longitudinal two-spin-wave modes Bunker and Landau (2000) that subsequently lead to experimental verification Schweika et al. (2002), and an unexpected form of transverse spin wave excitations in antiferromagnetic nanofilms Hou et al. (2015).

Study of collective dynamics in magnetic materials faces an enormous challenge due to the coupling of lattice vibrations and spin waves which is inherently neglected in the aforementioned atomistic models. In magnetic metals and alloys, the atomic magnetic moments and exchange interactions strongly depend on the local atomic environment Shimizu (1981); Sabiryanov et al. (1995); Sabiryanov and Jaswal (1999) and therefore change dynamically as the local crystal structure is distorted by lattice vibrations Yin et al. (2012). On the other hand, magnetic interactions themselves are integral for maintaining the structural stability of such systems Abrikosov et al. (1996); Ekman et al. (1998). For instance, the stabilization of the bcc crystal structure in iron is long conceived to be of magnetic origin Hasegawa and Pettifor (1983); Herper et al. (1999). Furthermore, a number of recent studies emphasize the significance of phonon-magnon coupling on various dynamical processes such as self diffusion Ding et al. (2014), thermal transport Boona and Heremans (2014), dislocation dynamics Dudarev et al. (2008), and spin-Seebeck effect Jaworski et al. (2011). The dynamics of atomic and magnetic degrees of freedom are, hence, inseparable and should be treated in a self-consistent manner.

The idea of integrating spin dynamics with molecular dynamics was pioneered by Omelyan et al. Omelyan et al. (2001) in the context of a simple model for ferrofluids. The foundation of this combined molecular and spin dynamics (MD-SD) approach lies in the unification of an atomistic potential and a Heisenberg spin Hamiltonian, with the coupling between the atomic and spin subsystems established via a coordinate-dependent exchange interaction. With the use of an empirical many-body potential and a parameterized exchange interaction, Ma et al. Ma et al. (2008) further extended MD-SD into a framework for realistic modeling of bcc iron. The parameterization developed by Ma et al. Ma et al. (2008) has since been successfully adopted to investigate various phenomena in bcc iron such as magneto-volume effects Chui and Zhou (2014a), vacancy formation and migration Wen et al. (2013); Wen and Woo (2014), and external magnetic field effects Chui and Zhou (2014b). Moreover, the method has been recently extended by incorporating spin-orbit interactions to facilitate the dynamic exchange of angular momentum between the lattice and spin subsystems Perera et al. (2016a). This, in particular, extends the applicability of MD-SD to accurate modeling of non-equilibrium processes.

The aim of this paper is to improve our understanding of phonon-magnon interactions in the ferromagnetic phase of bcc iron within the context of MD-SD. This study is an extension of our earlier preliminary work Perera et al. (2014a, b) which primarily focused on the effect of lattice vibrations on the spin-spin dynamic structure factor in the [100] lattice direction. In this paper, we provide a more in-depth analysis of the mutual influence of phonons and magnons on their respective frequency spectra and lifetimes for all three high-symmetry lattice directions: [100], [110] and [111]. This is achieved by comparing the results obtained for MD-SD simulations with those of standalone MD and SD simulations in which spin-lattice coupling is completely neglected. In Sec. II, we present the MD-SD formalism and the parameterization for bcc iron, followed by a comprehensive description of the methods we adopt for characterizing collective excitations. Sec. III.1 and  III.2, respectively, report our results on vibrational and magnetic excitations, followed by conclusions in Sec. IV.

II Methods

II.1 Combined molecular and spin dynamics

MD-SD is essentially a reformulation of the MD approach, in which the effective spin angular momenta of the atoms $\{\mathbf{S}_{i}\}$ are incorporated into the Hamiltonian and treated as explicit phase variables. For a classical system of $N$ magnetic atoms of mass $m$ described by their positions $\{\mathbf{r}_{i}\}$, velocities $\{\mathbf{v}_{i}\}$, and the atomic spins $\{\mathbf{S}_{i}\}$, the MD-SD Hamiltonian takes the form

[TABLE]

where the first term represents the kinetic energy of the atoms, and $U(\{\mathbf{r}_{i}\})$ is the spin-independent (non-magnetic) scalar interaction between the atoms. The Heisenberg-like exchange interaction with the coordinate-dependent exchange parameter and $J_{ij}(\left\{\mathbf{r}_{k}\right\})$ specifies the exchange coupling between the $i$th and $j$th spins. The aforementioned Hamiltonian has true dynamics as described by the classical equations of motion

[TABLE]

where $\mathbf{f}_{i}=-\nabla{}_{\mathbf{r}_{i}}\mathcal{H}$ and $\mathbf{H}_{i}^{\text{eff}}=\nabla{}_{\mathbf{S}_{i}}\mathcal{H}$ are the interatomic force and the effective field acting on the $i$th atom/spin. The goal of the MD-SD approach is to numerically solve the above equations of motion starting from a given initial configuration, and obtain the trajectories of both the atomic and spin degrees of freedom.

MD-SD is a generic framework that with proper parameterization, can be readily adopted for any magnetic material in which the spin interactions can be modeled classically. In this study, we adopt the parameterization introduced by Ma et al. Ma et al. (2008) for bcc iron, in which $U(\{\mathbf{r}_{i}\})$ is constructed as

[TABLE]

where $U_{\text{DD}}$ is the “magnetic” embedded atom potential developed by Dudarev and Derlet Dudarev and Derlet (2005); Derlet and Dudarev (2007), and ${E_{\text{spin}}^{\text{ground}}=-\sum_{i<j}J_{ij}(\left\{\mathbf{r}_{k}\right\})|\mathbf{S}_{i}||\mathbf{S}_{j}|}$ is the energy contribution from a collinear spin state, subtracted out to eliminate the magnetic interaction energy that is implicitly contained in $U_{\text{DD}}$. With the particular form of $U(\{\mathbf{r}_{i}\})$ given in Eq. (3), Hamiltonian (1) provides the same ground state energy as $U_{\text{DD}}$. The exchange interaction is modeled via a simple pairwise function $J(r_{ij})$ parameterized by first-principles calculations Ma et al. (2008), with spin lengths absorbed into its definition, i.e. ${J(r_{ij})=J_{ij}(\left\{\mathbf{r}_{k}\right\})|\mathbf{S}_{i}||\mathbf{S}_{j}|}$. We assume constant spin lengths $|\mathbf{S}|=2.2/g$, with $g$ being the electron $g$ factor.

We would like to point out that the fluctuation of the magnitudes of magnetic moments and spin-orbit interactions are not considered in this work. In transition metals and alloys, fluctuation of spin magnitudes may have a notable effect on the material properties, particularly at high temperatures. Ma et al. Ma and Dudarev (2012) proposed a way of incorporating longitudinal spin fluctuations into SD and MD-SD simulations via a Langevin-type equation of motion within the context of fluctuation-dissipation theorem. Numerical coefficients of the corresponding Landau Hamiltonian can be determined from ab initio calculations Ma and Dudarev (2012, 2014). An accurate depiction of spin-orbit interactions can be potentially achieved with the use of Hubbard-like Hamiltonians as the foundation for deriving the equations of motion Coury et al. (2016). A phenomenological approach for modeling spin-orbit interactions in MD-SD has also been recently proposed Perera et al. (2016a), but was not adopted in this study due to its computationally demanding nature.

II.2 Characterizing collective excitations

In MD and SD simulations, space-displaced, time-displaced correlation functions of the microscopic dynamical variables are integral to the study of the collective phenomena in the system Rapaport (2004); Lovesey (1984); Landau and Krech (1999). Fourier transforms of these quantities directly yield information regarding the frequency spectra and the lifetimes of the respective collective modes.

Let us define microscopic atom density as

[TABLE]

The spatial Fourier transform of the space-displaced, time-displaced density-density correlation function, namely, the intermediate scattering function Hansen and McDonald (2006) then takes the form

[TABLE]

where $\rho_{n}(\mathbf{q},t)=\int\rho_{n}(\mathbf{r},t)e^{-i\mathbf{q}\cdot\mathbf{r}}d\mathbf{r}=\sum_{i}e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}$. The power spectrum of the intermediate scattering function

[TABLE]

is called the “density-density dynamic structure factor” for the momentum transfer $\mathbf{q}$ and frequency (energy) transfer $\omega$. $S_{nn}(\mathbf{q},\omega)$ is directly related to the differential cross section measured in inelastic neutron scattering experiments Hansen and McDonald (2006). Local density fluctuations in a system are caused by the thermal diffusion of atoms as well as vibrational modes related to the propagating lattice waves Anento and Padró (2004). For liquid systems, the thermal diffusive mode can be identified as a peak in $S_{nn}(\mathbf{q},\omega)$ centered at $\omega=0$, whereas for solids this peak will disappear due to the absence of thermal diffusion Anento and Padró (2004). In crystalline solids, peaks in $S_{nn}(\mathbf{q},\omega)$ at non-zero frequencies can be uniquely associated with longitudinal vibrational modes with the corresponding frequencies and wave vectors. As the transverse lattice vibrations do not cause local density fluctuations towards the direction of wave propagation, $S_{nn}(\mathbf{q},\omega)$ is incapable of revealing information about these modes. Therefore, to identify transverse lattice vibrations, one needs to consider the time-dependent correlations of transverse velocity components.

With the microscopic “velocity density” defined as $\boldsymbol{\rho}_{v}(\mathbf{r},t)=\sum_{i}\mathbf{v}_{i}(t)\delta\left[\mathbf{r}-\mathbf{r}_{i}(t)\right]$, the spatial Fourier transform of the velocity-velocity correlation function takes the form

[TABLE]

where $\boldsymbol{\rho}_{v}^{L,T}(\mathbf{q},t)=\sum_{i}\mathbf{v}_{i}^{L,T}(t)e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}$, with the superscripts $L$ and $T$ respectively denoting the longitudinal and transverse components with reference to the direction of the wave propagation. Peaks in the corresponding power spectra $S_{vv}^{L}(\mathbf{q},\omega)$ and $S_{vv}^{T}(\mathbf{q},\omega)$ respectively reveal longitudinal and transverse vibrational modes of the system. It can be shown that $S_{vv}^{L}(\mathbf{q},\omega)$ is directly related to the density-density dynamic structure factor $S_{nn}(\mathbf{q},\omega)$ via the relationship $S_{nn}(\mathbf{q},\omega)=\omega^{2}/q^{2}S_{vv}^{L}(\mathbf{q},\omega)$ Hansen and McDonald (2006); Anento and Padró (2004).

Just as the time-dependent density-density and velocity-velocity correlations reveal vibrational excitations associated with the lattice subsystem, spin density autocorrelations can elucidate the magnetic excitations associated with the spin subsystem.

The microscopic “spin density” is given by

[TABLE]

Treating the spin-spin correlations along $x$, $y$, and $z$ directions separately, we define the intermediate scattering function as

[TABLE]

where $k=x,y,$ or $z$, and $\boldsymbol{\rho}_{s}(\mathbf{q},t)=\sum_{i}\mathbf{S}_{i}(t)e^{-i\mathbf{q}\cdot\mathbf{r}_{i}(t)}$. For a ferromagnetic system in the microcanonical ensemble, the magnetization vector is a constant of motion and serves as a fixed symmetry axis throughout the time evolution of the system. To differentiate between the magnetic excitations that propagate parallel and perpendicular to this symmetry axis, we redefine the coordinate system in spin space such that the $z$ axis is parallel to the magnetization vector. The components $\{F_{ss}^{k}(\mathbf{q},t)\}$ can then be simply regrouped to yield the longitudinal component

[TABLE]

and the transverse component

[TABLE]

Note that the separation of magnetic excitations into longitudinal and transverse modes is only meaningful for temperatures below the Curie temperature $T_{\text{C}}$, since above $T_{\text{C}}$, the net magnetization vanishes and all directions in spin space become equivalent.

Fourier transforms of $F_{ss}^{L,T}(\mathbf{q},t)$ yield the spin-spin dynamic structure factors $S_{ss}^{L,T}(\mathbf{q},\omega)$. Just like the density-density dynamic structure factor, the spin-spin dynamic structure factor is a measurable quantity in inelastic neutron scattering experiments Lovesey (1984); Landau and Krech (1999).

In this study, we are primarily interested in investigating wave propagation along the three principle lattice directions: [100], [110] and [111]. Let us denote the wave vectors in these directions as $\mathbf{q}=(q,0,0)$, $(q,q,0)$, and $(q,q,q)$, respectively. Due to the finite size of the simulation box, the accessible values of $q$ in each direction is constrained to a discrete set given by $q=2\pi n_{q}/La$, with $n_{q}=\pm 1,\pm 2,\dotsc,\pm,L$ for the [100] and [111] directions, and $n_{q}=\pm 1,\pm 2,\dotsc,\pm,L/2$ for the [110] direction, where $L$ is the linear lattice dimension and $a=2.8665$ Å is the lattice constant of bcc iron.

II.3 Simulation details

For integrating the coupled equations of motion presented in Eq. (2), we adopted an algorithm based on the second order Suzuki-Trotter (ST) decomposition of the non-commuting operators Omelyan et al. (2001); Tsai et al. (2005); Krech et al. (1998). To obtain a reasonable level of accuracy as reflected by the energy and magnetization conservation, an integration time step of $\Delta t=1$ fs was used.

For computing canonical averages of time-dependent correlation functions, we used time series obtained from microcanonial dynamical simulations, that are, in turn, initiated from equilibrium states drawn from the canonical ensemble at the desired temperature $T$. Averaging over the results of multiple simulations started from different initial states yields good estimates of the respective canonical ensemble averages Landau and Krech (1999).

For generating the initial states for our microcanonical MD-SD simulations, we adhere to the following procedure. First, we equilibrate the subspace consisting of positions and spins using the Metropolis Monte Carlo (MC) method Metropolis et al. (1953). As the second step, we assign initial velocities to the atoms based on the Maxwell-Boltzmann distribution at the desired temperature $T$. Finally, we perform a short microcanonical MD-SD equilibration run (typically $\sim 1000$ time steps with $\Delta t=1$ fs), which would ultimately bring the whole system to the equilibrium by resolving any inconsistencies between the position-spin subspace and the velocity distribution. Fig. 1 shows the time evolution of the instantaneous lattice and spin temperatures as observed in a microcanonical MD-SD simulation initiated from an equilibrium state generated from the aforementioned technique for $T=800$ K. Both lattice and spin temperatures fluctuate about a mean value of $T=800$ K, indicating that the lattice and the spin subsystems are in mutual equilibrium.

To characterize phonon and magnon modes, we performed simulations for the system size $L=16$ ($8192$ atoms) at temperatures $T=300$ K, $800$ K, and $1000$ K. $T=1000$ K was particularly chosen due to its vicinity to the Curie temperature of bcc iron, $T_{C}\approx 1043$ K. (A recent high resolution Monte Carlo study has revealed that the transition temperature of the particular spin-lattice model used in our study to be $T\approx 1078$ K Perera et al. (2016b)). Equations of motion were integrated up to a total time of $t_{\text{max}}=1$ ns, and the space-displaced, time-displaced correlation functions were computed for the three principle lattice directions: [100], [110] and [111]. To increase the accuracy, we have averaged these quantities over different starting points in the time series. Canonical ensemble averages were estimated using the results of $200$ independent simulations, each initiated from a different initial state. The time Fourier transform in Eq. (6) was carried out to a cutoff time of $t_{\text{cutoff}}=0.5$ ns.

As our primary goal is to understand the mutual impact of the phonons and magnons on their respective frequency spectra and lifetimes, we have also performed standalone MD and SD simulations for comparison. For MD simulations, we used the Dudarev-Derlet potential to model the interatomic interactions while completely neglecting the spin-spin interactions. SD simulations were conducted with the atoms frozen at perfect bcc lattice positions, and the exchange parameters determined from the same pairwise function used for MD-SD simulations.

III Results

III.1 Vibrational excitations

For all the temperatures considered, we observe well defined excitation peaks at non-zero frequencies in the density-density dynamic structure factor $S_{nn}(\mathbf{q},\omega)$, as well as in the longitudinal and the transverse components of the velocity-velocity dynamic structure factor: $S_{vv}^{L}(\mathbf{q},\omega)$ and $S_{vv}^{T}(\mathbf{q},\omega)$. For each $\mathbf{q}$ along [100] and [111] lattice directions, all three quantities show single peaks (See Fig. 2 for an example). The peak positions in $S_{nn}(\mathbf{q},\omega)$ and $S_{vv}^{L}(\mathbf{q},\omega)$ for the same wave vector coincide with each other as they are both associated with the longitudinal vibrational modes, and hence convey the same information. The peaks in $S_{vv}^{T}(\mathbf{q},\omega)$ are associated with the transverse lattice vibrations. Since there are two orthogonal directions perpendicular to a given wave vector $\mathbf{q}$, there are, in fact, two transverse vibrational modes for each $\mathbf{q}$. Due to the four-fold and three-fold rotational symmetry about the axes [100] and [111], respectively, the two transverse modes for the wave vectors along these directions become degenerate Dove (2005). As a result, we only observe a single peak in $S_{vv}^{T}(\mathbf{q},\omega)$ for the wave vectors along these directions.

We also observe single peak structures in $S_{nn}(\mathbf{q},\omega)$ and $S_{vv}^{L}(\mathbf{q},\omega)$ for the wave vectors along the [110] direction. However, for the case of $S_{vv}^{T}(\mathbf{q},\omega)$, one can clearly identify two distinct peaks. This is a consequence of the two transverse modes being non-degenerate due to the reduced rotational symmetry (two-fold) about the [110] axis in comparison to [100] and [111] directions Dove (2005).

To extract the positions and the half-widths of the phonon peaks, we fit the simulation results for the dynamic structure factor to a Lorentzian function of the form Tsai et al. (2000); Tao et al. (2005)

[TABLE]

where $\omega_{0}$ is the characteristic frequency of the vibrational mode, $I_{0}$ is the intensity or the amplitude of the peak, and $\Gamma$ is the half-width at half maximum (HWHM) which is inversely proportional to the lifetime of the excitation. The errors of the fitting parameters were estimated using the following procedure. The complete set of correlation function estimates obtained from $200$ independent simulations was divided into $10$ groups, and the data within each group were averaged over to yield $10$ results sets. Dynamic structure factors were independently computed for these $10$ correlation function sets. To estimate the errors in the fitting parameters, we separately performed curve fits to these $10$ independent dynamic structure factor estimates, and calculated the standard deviations of the fitting parameters. Statistical errors bars obtained in this manner were found to be an order of magnitude larger than the error bars estimated by the curve-fitting tool.

For all the temperatures and wave vectors considered, the Lorentzian lineshape given in Eq. (12) fitted well with the peaks observed in $S_{nn}(\mathbf{q},\omega)$ and $S_{vv}^{L,T}(\mathbf{q},\omega)$. Fig. 2 shows an example curve fit for the MD-SD results of $S_{nn}(\mathbf{q},\omega)$ for $\mathbf{q}=(1.1$ Å-1,$0,0)$ at $T=300$ K. To fit the two peak structure observed in $S_{vv}^{T}(\mathbf{q},\omega)$ for the [110] direction, we use the sum of two Lorentzians.

Using the peak positions obtained from the Lorentzian fits, one can construct phonon dispersion relations for the three principle lattice directions. Fig. 3 shows the the dispersion curves determined from our MD-SD simulations for $T=300$ K, along with the experimental results Minkiewicz et al. (1967); Brockhouse et al. (1967) obtained from inelastic neutron scattering. For comparison, we have also shown the results of standalone MD simulations for the same temperature. In general, for small to moderate $q$ values, both MD-SD and MD dispersion curves agree well with the experimental results, but deviations can be observed for larger $q$ values, particularly near the zone boundaries in [100] and [111] directions. Although the MD-SD and MD dispersion curves are indistinguishable within the resolution of Fig. 3, we will show later on that there are, in fact, deviations larger than the error bars.

At temperatures in the vicinity of absolute zero, due to the low occupation of vibrational modes, phonons behave as weakly interacting quasiparticles that can be treated within the harmonic approximation Maradudin and Fein (1962). In this limit, characteristic frequencies of the phonons are well defined and the lifetimes are practically infinite. As the temperature is increased, phonon occupation numbers also increase, which in turn increases the probability of mutual interactions. As a result of such phonon-phonon scattering at elevated temperatures, characteristic frequencies of the phonons may shift, and the lifetimes may shorten Maradudin and Fein (1962); Fultz (2010). In magnetic crystals, the co-existence of phonons and magnons gives rise to another class of scattering processes, namely, phonon-magnon scattering. Just as phonon-phonon scattering, phonon-magnon scattering may also lead to a shift in the characteristic phonon frequencies, as well as shortening of the phonon lifetimes. As the occupancy of both phonon and magnon modes increases with temperature, these effects will be more pronounced as the temperature is increased.

To carefully examine the changes in the phonon frequency spectrum due to magnons, we compare the characteristic frequencies determined from MD-SD simulations ($\omega_{{}_{\text{MD-SD}}})$ with the ones obtained from MD simulations ($\omega_{{}_{\text{MD}}}$) by calculating the fractional frequency shift, ${\left(\omega_{{}_{\text{MD-SD}}}-\omega_{{}_{\text{MD}}}\right)/\omega_{{}_{\text{MD}}}}$. The results for the three principle directions are shown in Figs. 4 and 5, for the longitudinal and the transverse modes, respectively. With the exception of the high frequency transverse branch along [110] direction (TA2), phonon frequencies shift to higher values in the presence of magnons. In general, the shift in frequencies becomes more pronounced as the temperature is increased. A particularly interesting behavior occurs in the longitudinal branch for the [111] direction where we observe dips in the curves for all three temperatures at the same $q$ value. For all three temperatures, the frequency shift of the vibrational mode that corresponds to the bottom of the dip is close to zero. Therefore, the frequency of this phonon mode appears to be unaffected by the presence of magnons.

Lifetimes of the phonon excitations are inversely proportional to the half-widths at half maximum of the corresponding vibrational peaks observed in $S_{nn}(\mathbf{q},\omega)$ and $S_{vv}^{L,T}(\mathbf{q},\omega)$. To study the impact of the magnons on the phonon lifetimes, we compare the half-widths obtained from MD-SD simulations with that of the MD simulations. Fig. 6 shows the results for the longitudinal phonons. For the longitudinal phonons at $T=300$ K, a marginal increase in the half-widths can be observed due to the magnons, which becomes more pronounced as the temperature is increased. For the case of transverse phonons, we did not observe any notable difference between the MD-SD and MD half-widths outside the error bars, for all the temperatures considered.

III.2 Magnetic excitations

III.2.1 Transverse magnon modes

For the temperatures $T=300$ K and $T=800$ K, our results for the transverse component of the spin-spin dynamic structure factor $S_{ss}^{T}(\mathbf{q},\omega)$ show a single spin wave peak, that can be fitted to a Lorentzian lineshape of the form Eq. (12) (See Fig. 7 (a) for an example.). For $T=1000$ K, we also observe a diffusive central peak at $\omega=0$, as observed in neutron scattering experiments Mook and Lynn (1985) and previous SD studies Tao et al. (2005). This two-peak structure can be best captured by a function of the form Tao et al. (2005)

[TABLE]

where the first term (Gaussian) corresponds to the central peak, and the second term (Lorentzian) describes the spin wave peak (See Fig. 7 (b) for an example.). For large $\mathbf{q}$ values at $T=800$ K and $T=1000$ K, spin wave peaks in $S_{ss}^{T}(\mathbf{q},\omega)$ were found to be asymmetric, and hence did not yield good fits to Lorentzian lineshapes. Therefore, one cannot obtain reliable estimates of the magnon half-widths. However, spin wave peak positions can still be determined relatively precisely, thus the magnon dispersion relations can be constructed.

Fig. 8 shows the transverse magnon dispersion relations for small $|\mathbf{q}|$ values along the three principle directions as determined from MD-SD simulations at $T=300$ K. In agreement with the experimental findings Lynn (1975); Collins et al. (1969), the three dispersion relations are isotropic when plotted as functions of the magnitude of the wave vector $|\mathbf{q}|$. Moreover, for small $|\mathbf{q}|$ values, our results agree quantitatively with the experimental results for the [110] direction Lynn (1975); Collins et al. (1969). Fig. 9 shows the complete dispersion curves determined from MD-SD and SD simulations for $T=300$ K, $T=800$ K, and $T=1000$ K. For both MD-SD and SD, the characteristic frequencies shift to lower values as the temperature is increased. This indicates increased magnon-magnon scattering at elevated temperatures. For all three temperatures, particularly near the zone boundaries, we can observe a marginal difference between the MD-SD and SD dispersion curves. This, in fact, is a result of phonon-magnon scattering.

To further investigate the magnon softening due to phonons, we calculate the fractional frequency shift of the magnons, ${\left(\omega_{{}_{\text{MD-SD}}}-\omega_{{}_{\text{SD}}}\right)/\omega_{{}_{\text{SD}}}}$. The results are shown in Fig. 10 for the three principle directions. For small $q$ values, magnon modes shift to lower frequencies in the presence of phonons. As $q$ increases, the direction of the shift is reversed. Moreover, the shift in frequencies becomes more pronounced as the temperature is increased.

Fig. 11 compares the transverse magnon half-widths obtained from MD-SD and SD simulations for $T=300$ K. Although the difference between the half-widths is negligible for small $q$ values, for moderate to large $q$ values, half-widths for the MD-SD results are significantly larger than that for the SD results. This indicates significant shortening of the magnon lifetimes due to phonon-magnon scattering.

III.2.2 Longitudinal magnon modes

Our results for the longitudinal spin-spin dynamic structure factor $S_{ss}^{L}(\mathbf{q},\omega)$ obtained from both MD-SD and SD simulations show many very low-intensity excitations peaks, for all wave vectors considered. Fig. 12 shows $S_{ss}^{L}(\mathbf{q},\omega)$ for a small system size $L=8$ at $T=300$ K, where we compare the SD results [panel (a)] with the MD-SD results [panel (b)] for $\mathbf{q}=\frac{2\pi}{La}(1,0,0)$.

In the context of classical Heisenberg models, Bunker et al. Bunker and Landau (2000) showed that the excitation peaks observed in $S_{ss}^{L}(\mathbf{q},\omega)$ are two-spin-wave creation and/or annihilation peaks which result from the pairwise interactions between transverse magnon modes. For ferromagnetic systems, only spin wave annihilation peaks are present, and their frequencies are given by

[TABLE]

where $\mathbf{q}_{i}$ and $\mathbf{q}_{j}$ are the wave vectors of the two transverse magnon modes which comprise the two-spin-wave excitation. Since the set of allowable wave vectors $\{\mathbf{q}_{i}\}$ depends on the system size $L$, the resultant two-spin-wave spectrum also varies with $L$. For a real magnetic crystal where $L$ is practically infinite, the two-spin-wave spectrum would become continuous.

To verify whether the peaks we observe in $S_{ss}^{L}(\mathbf{q},\omega)$ are two-spin-wave peaks, we chose a relatively small system size ($L=8$) so that the set of allowable wave vectors is reduced to a manageable size. Then, using MD-SD and SD simulations, we separately determined the transverse magnon frequencies that correspond to the first few $n_{q}$ values along all possible lattice directions. With this information at hand, we can predict the expected positions of all the two-spin-wave annihilation peaks using Eq. (14) for both SD and MD-SD case. As an example, let us consider the wave vector pair $\mathbf{q}_{i}=(1,1,1)$ and $\mathbf{q}_{j}=(1,1,0)$. Since $\mathbf{q}_{i}-\mathbf{q}_{j}=(0,0,1)$, they produce a spin wave annihilation peak in $S_{ss}^{L}(\mathbf{q},\omega)$ for $\mathbf{q}=(0,0,1)$ at the frequency $\omega^{-}=\omega(\mathbf{q}_{i})-\omega(\mathbf{q}_{j})$. (Note that we have ignored the common pre-factor $2\pi/La$ from the wave vectors.)

In Fig. 12 (a) and (b), we have superimposed the predicted spin wave annihilation peak positions corresponding to each case. We see an excellent match between the observed peaks and the predicted two-spin-wave peak positions, with the exception of the particular sharp peak at $\omega\approx 10$ meV which only appears in panel (b). Surprisingly, the position of this peak coincides with the frequency of the longitudinal phonon mode for the same $\mathbf{q}$ as determined from the peak position of $S_{nn}(\mathbf{q},\omega)$ or $S_{vv}^{L}(\mathbf{q},\omega)$. Similar excitation peaks were observed for all wave vectors, for all system sizes and temperatures considered. The origin of these coupled phonon-magnon modes can be explained as follows. Unlike transverse phonons, when a longitudinal phonon propagates along a certain lattice direction, it generates fluctuations in the local atom density along that direction with the corresponding phonon frequency. This, in turn, leads to fluctuations in the local density of the longitudinal components of the spins (i.e. components of the spin vectors parallel to the net magnetization). These longitudinal spin fluctuations propagate along with the phonon, yielding a sharp, coupled mode in the the longitudinal magnon spectrum.

Fig. 13 and Fig. 14 show $S_{ss}^{L}(\mathbf{q},\omega)$ for $\mathbf{q}=\frac{2\pi}{La}(1,0,0)$ at $T=800$ K and $T=1000$ K, respectively, where we compare the SD results [panel (a)] with the MD-SD results [panel (b)]. In each figure, the inset of panel (b) shows the longitudinal density-density dynamic structure factor for the same wave vector. In comparison to the results for $T=300$ K, we observe that the diffusive central peak becomes more pronounced as the temperature rises, and many of the low-intensity two-spin-wave peaks broaden and disappear into its tail. These observations are in qualitative agreement with previous SD studies of the ferromagnetic Heisenberg model Bunker and Landau (2000). The coupled phonon-magnon peak also becomes less pronounced with increasing temperature, as the diffusive central peak becomes more pronounced. At $T=1000$ K, the intensity of the peak is very low and is barely recognizable. Above the Curie temperature, spins are randomly oriented and the vector sum of spins per unit volume will be zero on average. Hence, the coupled phonon-magnon mode should entirely disappear; however, it is already so faint at $T=1000$ K that we clearly would not have sufficient resolution to test this behavior above the Curie temperature.

We would like to point out that the existence of these coupled phonon-magnon modes is a phenomenon that so far hasn’t been discovered experimentally. In fact, this is not surprising since the experimental detection of these peaks would be extremely challenging due to their very low intensities.

IV Conclusions

To investigate collective phenomena in ferromagnetic bcc iron, we performed combined molecular and spin dynamics (MD-SD) simulations at temperatures $T=300$ K, $T=800$ K, and $T=1000$ K. From the trajectories of these simulations, space- and time-displaced correlation functions associated with the atomic and spin variables were calculated. Fourier transforms of these quantities, namely, dynamic structure factors, directly reveal information regarding the frequencies and the lifetimes of the vibrational and magnetic excitation modes. For small $q$ values, the dispersion relations obtained from our simulations at $T=300$ K agree well with the experimental results, but deviations can be observed for large $q$ values, especially for the transverse magnon dispersion curves. These discrepancies can be attributed to the anharmonic effects not being faithfully captured in the embedded atom potential and the pairwise functional representation of the exchange interaction. In fact, Yin et al. Yin et al. (2012) recently pointed out that the exchange parameters in bcc iron depend on the local atomic environment in a complicated manner that may not be properly characterized through a pairwise distance-dependent function. Thus, a more accurate depiction of magnetic interactions necessitates the development of sophisticated models of exchange interactions that effectively capture the contribution of the local environment.

To understand the mutual influence of the phonons and magnons on each other, we compared our results with that of the standalone molecular dynamics and spin dynamics simulations. Due to phonon-magnon coupling, we observe a shift in the characteristic frequencies, as well as a decrease in the lifetimes. These effects become more pronounced as the temperature is increased. Moreover, the frequency shifts and the lifetime reductions that occur in magnons due to phonons are found to be far more pronounced than the corresponding effects experienced by phonons due to magnons. This is not surprising considering the fact that the energy scale associated with the spin-spin interactions is about an order of magnitude smaller than that of the atomic (non-magnetic) interactions.

A comparison of our results at different temperatures shows that the effects of spin-lattice coupling becomes more pronounced as the temperature rises. However, due to critical fluctuations, the size of the error bars for magnon properties increases rapidly as the temperature approaches the Curie temperature (See Fig. 7 (b) and Fig. 10 for example.). Therefore, obtaining reliable estimates of magnon properties becomes increasingly difficult as the Curie temperature is approached.

The unprecedented resolution provided by our simulations has allowed us to clearly identify two-spin-wave peaks in the longitudinal spin-spin dynamic structure factor with amplitudes down to six orders of magnitude smaller than that of the highest single spin wave peak observed. In addition, in the presence of lattice vibrations, we also observe additional longitudinal magnetic excitations with frequencies which coincide with those of the longitudinal phonons. This is an unexpected form of longitudinal spin wave excitations that so far has not been detected in inelastic neutron scattering experiments, presumably due to their very low intensities.

Acknowledgements.

This work was sponsored by the “Center for Defect Physics”, an Energy Frontier Research Center of the Office of Basic Energy Sciences (BES), U.S. Department of Energy (DOE); the later stages of the work of G.M.S. and M.E. was supported by the Materials Sciences and Engineering Division of BES, US-DOE. We also acknowledge the computational resources provided by the Georgia Advanced Computing Resource Center.