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

0.94 MB·

by  Dilina Perera

Collective dynamics in atomistic models with coupled translational and spin degrees
of freedom
Dilina Perera,1,2,∗ Don M. Nicholson,3 Markus Eisenbach,4 G. Malcolm Stocks,4 and David P. Landau1
1Center for Simulational Physics, The University of Georgia, Athens, Georgia 30602, USA
2Department of Physics and Astronomy, Mississippi State University, Mississippi State, Mississippi 39762, USA
3University of North Carolina at Asheville, Asheville, North Carolina 28804, USA
4Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
7 Using an atomistic model that simultaneously treats the dynamics of translational and spin de-
1 grees offreedom, weperform combined molecular and spin dynamicssimulations toinvestigate the
0 mutualinfluenceofthephononsandmagnonsontheirrespectivefrequencyspectraandlifetimesin
2 ferromagnetic bcciron. Bycalculating theFouriertransforms ofthespace- andtime-displaced cor-
n relationfunctions,thecharacteristicfrequenciesandthelinewidthsofthevibrationalandmagnetic
a excitationmodesweredetermined. Comparisonoftheresultswiththatofthestandalonemolecular
J dynamics and spin dynamics simulations reveal that the dynamic interplay between the phonons
and magnons leads to a shift in the respective frequency spectra and a decrease in the lifetimes.
7
2 Moreover, in the presence of lattice vibrations, additional longitudinal magnetic excitations were
observed with thesame frequencies as thelongitudinal phonons.

]
i
c
s I. INTRODUCTION of recent studies emphasize the significance of phonon-
- magnoncouplingonvariousdynamicalprocessessuchas
l
r For decades, dynamical simulations of atomistic mod- self diffusion [29], thermal transport[30], dislocation dy-
t
m els have played a pivotal role in the study of collective namics [31], and spin-Seebeck effect [32]. The dynamics
. phenomena in materials at finite temperatures. Molec- ofatomicandmagneticdegreesoffreedomare,hence,in-
at ular dynamics (MD) [1, 2] utilizing empirical poten- separableandshouldbetreatedinaself-consistentman-
m tials has been extensively used in the analysis of vibra- ner.
- tional properties in a variety of systems such as metals The idea of integrating spin dynamics with molecu-
d and alloys [3–7], polymers [8, 9], carbon nanotubes [10], lar dynamics was pioneered by Omelyan et al. [33] in
n graphene [11] etc. With regard to magnetic excitations, the context of a simple model for ferrofluids. The foun-
o thelesser-knownspindynamics(SD)method[12–16]has dation of this combined molecular and spin dynamics
c
[ provento be an indispensable toolfor investigatingclas- (MD-SD) approach lies in the unification of an atom-
sical lattice-based spin models for which the analytical istic potential and a Heisenberg spin Hamiltonian, with
1 solutions are intractable. Over the years, SD simula- thecouplingbetweentheatomicandspinsubsystemses-
v tions have expanded our understanding of spin waves tablished via a coordinate-dependent exchange interac-
6
and solitons in magnetic materials, leading to a number tion. With the use of an empirical many-body potential
0
9 of groundbreaking discoveries, including the existence of and a parameterizedexchange interaction, Ma et al. [34]
7 propagating spin waves in paramagnetic bcc iron [17], further extended MD-SD into a framework for realistic
0 presence of longitudinal two-spin-wave modes [18] that modeling of bcc iron. The parameterization developed
1. subsequently lead to experimental verification [19], and by Ma et al. [34] has since been successfully adopted
0 an unexpected form of transverse spin wave excitations to investigate various phenomena in bcc iron such as
7 in antiferromagnetic nanofilms [20]. magneto-volume effects [35], vacancy formation and mi-
1 Study of collective dynamics in magnetic materials gration [36, 37], and external magnetic field effects [38].
v: facesanenormouschallengeduetothecouplingoflattice Moreover,the method has been recently extended by in-
i vibrationsandspinwaveswhichisinherentlyneglectedin corporating spin-orbit interactions to facilitate the dy-
X the aforementioned atomistic models. In magnetic met- namic exchange of angular momentum between the lat-
r als and alloys, the atomic magnetic moments and ex- tice and spin subsystems [39]. This, in particular, ex-
a
change interactions strongly depend on the local atomic tends the applicability of MD-SD to accurate modeling
environment [21–23] and therefore change dynamically of non-equilibrium processes.
as the local crystal structure is distorted by lattice vi- The aim of this paper is to improve our understand-
brations [24]. On the other hand, magnetic interactions ing of phonon-magnon interactions in the ferromag-
themselvesareintegralformaintainingthestructuralsta- netic phase of bcc iron within the context of MD-SD.
bilityofsuchsystems[25,26]. Forinstance,thestabiliza- This study is an extension of our earlier preliminary
tion ofthe bcc crystalstructure in ironis long conceived work[40,41]whichprimarilyfocusedontheeffectoflat-
tobeofmagneticorigin[27,28]. Furthermore,anumber tice vibrations onthe spin-spindynamic structure factor
in the [100]lattice direction. In this paper, we provide a
morein-depthanalysisofthemutualinfluenceofphonons
and magnons on their respective frequency spectra and
∗ [email protected] lifetimes for all three high-symmetry lattice directions:
2
[100], [110] and [111]. This is achieved by comparing where U is the “magnetic” embedded atom poten-
DD
the results obtained for MD-SD simulations with those tial developed by Dudarev and Derlet [42, 43], and
of standalone MD and SD simulations in which spin- Eground =− J ({r })|S ||S | is the energy contri-
spin i<j ij k i j
lattice coupling is completely neglected. In Sec. II, we bution fromPa collinear spin state, subtracted out to
present the MD-SD formalism and the parameterization eliminate the magnetic interaction energy that is im-
for bcc iron, followed by a comprehensive description of plicitly contained in U . With the particular form
DD
the methods we adopt for characterizing collective exci- of U({r }) given in Eq. (3), Hamiltonian (1) pro-
i
tations. Sec. IIIA and IIIB, respectively, report our vides the same ground state energy as U . The ex-
DD
results on vibrationalandmagnetic excitations,followed change interaction is modeled via a simple pairwise
by conclusions in Sec. IV. function J(r ) parameterized by first-principles calcula-
ij
tions [34], with spin lengths absorbed into its definition,
i.e. J(r )=J ({r })|S ||S |. We assume constant spin
ij ij k i j
II. METHODS lengths |S|=2.2/g, with g being the electron g factor.
We would like to point out that the fluctuation of the
A. Combined molecular and spin dynamics magnitudesofmagneticmomentsandspin-orbitinterac-
tions are not considered in this work. In transition met-
MD-SD is essentially a reformulation of the MD ap- als and alloys, fluctuation of spin magnitudes may have
proach, in which the effective spin angular momenta of a notable effect on the material properties, particularly
the atoms {S } are incorporated into the Hamiltonian at high temperatures. Ma et al. [44] proposed a way
i
and treated as explicit phase variables. For a classical of incorporating longitudinal spin fluctuations into SD
system of N magnetic atoms of mass m described by andMD-SD simulationsvia a Langevin-typeequationof
theirpositions{r },velocities{v },andtheatomicspins motion within the context of fluctuation-dissipationthe-
i i
{S }, the MD-SD Hamiltonian takes the form orem. Numerical coefficients of the corresponding Lan-
i
dau Hamiltonian can be determined from ab initio cal-
N mv 2 culations [44, 45]. An accurate depiction of spin-orbit
H= i +U({r })− J ({r })S ·S , (1) interactions can be potentially achieved with the use of
i ij k i j
2
Xi=1 Xi<j Hubbard-like Hamiltonians as the foundation for deriv-
ing the equations of motion [46]. A phenomenological
where the first term represents the kinetic energy of approachfor modeling spin-orbit interactions in MD-SD
the atoms, and U({ri}) is the spin-independent (non- hasalsobeenrecentlyproposed[39],butwasnotadopted
magnetic) scalar interaction between the atoms. The in this study due to its computationally demanding na-
Heisenberg-likeexchangeinteractionwiththecoordinate- ture.
dependent exchange parameter and J ({r }) specifies
ij k
the exchange coupling between the ith and jth spins.
The aforementioned Hamiltonian has true dynamics as B. Characterizing collective excitations
described by the classical equations of motion
In MD and SD simulations, space-displaced, time-
dr
i =v (2a) displacedcorrelationfunctionsofthemicroscopicdynam-
i
dt ical variables are integral to the study of the collective
dv f phenomena in the system [1, 12, 47]. Fourier transforms
i i
= (2b)
dt m of these quantities directly yield information regarding
the frequency spectra and the lifetimes of the respective
dS 1
i = Heff×S (2c) collective modes.
dt ~ i i
Let us define microscopic atom density as
where fi =−∇riH andHeiff =∇SiH are the interatomic ρn(r,t)= δ[r−ri(t)]. (4)
force and the effective field acting on the ith atom/spin. X
i
The goalof the MD-SD approachis to numerically solve
The spatial Fourier transform of the space-displaced,
the above equations of motion starting from a given ini-
time-displaced density-density correlation function,
tialconfiguration,andobtainthetrajectoriesofboththe
namely, the intermediate scattering function [48] then
atomic and spin degrees of freedom.
takes the form
MD-SD is a generic framework that with proper pa-
1
rameterization, can be readily adopted for any magnetic F (q,t)= hρ (q,t)ρ (−q,0)i, (5)
nn n n
material in which the spin interactions can be modeled N
classically. In this study, we adopt the parameterization where ρ (q,t) = ρ (r,t)e−iq·rdr = e−iq·ri(t). The
introducedbyMaetal.[34]forbcciron,inwhichU({r }) n n i
i power spectrum oRf the intermediate scPattering function
is constructed as
1 +∞
S (q,ω)= F (q,t)e−iωtdt, (6)
U({ri})=UDD−Esgprionund, (3) nn 2π Z−∞ nn
3
is called the “density-density dynamic structure factor” system in spin space such that the z axis is parallel to
for the momentum transfer q and frequency (energy) the magnetization vector. The components {Fk(q,t)}
ss
transfer ω. S (q,ω) is directly related to the differ- can then be simply regrouped to yield the longitudinal
nn
ential cross section measured in inelastic neutron scat- component
tering experiments [48]. Local density fluctuations in a
system are caused by the thermal diffusion of atoms as FL(q,t)=Fz(q,t), (10)
ss ss
well as vibrationalmodes relatedto the propagatinglat-
tice waves [49]. For liquid systems, the thermal diffusive and the transverse component
mode can be identified as a peak in S (q,ω) centered
nn 1
at ω =0, whereas for solids this peak will disappear due FT(q,t)= (Fx(q,t)+Fy(q,t)). (11)
ss 2 ss ss
to the absence of thermal diffusion [49]. In crystalline
solids, peaks in S (q,ω) at non-zero frequencies can be
nn Note that the separation of magnetic excitations into
uniquely associated with longitudinal vibrational modes
longitudinal and transverse modes is only meaningful
withthecorrespondingfrequenciesandwavevectors. As
for temperatures below the Curie temperature T , since
C
the transverse lattice vibrations do not cause local den-
above T , the net magnetization vanishes and all direc-
C
sity fluctuations towards the direction of wave propa-
tions in spin space become equivalent.
gation, Snn(q,ω) is incapable of revealing information FouriertransformsofFL,T(q,t)yieldthespin-spindy-
ss
aboutthesemodes. Therefore,to identify transverselat- namicstructurefactorsSL,T(q,ω). Justlikethedensity-
ss
ticevibrations,oneneedstoconsiderthetime-dependent
density dynamic structure factor, the spin-spin dynamic
correlations of transverse velocity components.
structurefactorisameasurablequantityininelasticneu-
With the microscopic “velocity density” defined as
tron scattering experiments [12, 47].
ρ (r,t)= v (t)δ[r−r (t)], the spatial Fourier trans-
v i i i In this study, we are primarily interested in investi-
form of thPe velocity-velocity correlation function takes gating wave propagation along the three principle lat-
the form
tice directions: [100], [110] and [111]. Let us denote
1 the wave vectors in these directions as q = (q,0,0),
FL,T(q,t)= ρL,T(q,t)·ρL,T(−q,0) , (7)
vv N v v (q,q,0), and (q,q,q), respectively. Due to the finite
(cid:10) (cid:11)
size of the simulation box, the accessible values of q in
where ρLv,T(q,t) = iviL,T(t)e−iq·ri(t), with the super- each direction is constrained to a discrete set given by
scripts L and T resPpectively denoting the longitudinal q = 2πn /La, with n = ±1,±2,...,±,L for the [100]
q q
and transverse components with reference to the direc- and[111]directions,andn =±1,±2,...,±,L/2forthe
q
tionofthewavepropagation. Peaksinthecorresponding [110] direction, where L is the linear lattice dimension
powerspectraSvLv(q,ω)andSvTv(q,ω)respectivelyreveal and a=2.8665˚A is the lattice constant of bcc iron.
longitudinalandtransversevibrationalmodesofthesys-
tem. ItcanbeshownthatSL(q,ω)isdirectlyrelatedto
vv
the density-density dynamic structure factor Snn(q,ω) C. Simulation details
via the relationship S (q,ω)=ω2/q2SL(q,ω) [48, 49].
nn vv
Just as the time-dependent density-density and
For integrating the coupled equations of motion pre-
velocity-velocity correlations reveal vibrational excita-
sented in Eq. (2), we adopted an algorithmbased on the
tions associated with the lattice subsystem, spin density
second order Suzuki-Trotter (ST) decomposition of the
autocorrelations can elucidate the magnetic excitations
non-commuting operators [33, 50, 51]. To obtain a rea-
associated with the spin subsystem.
sonable level of accuracy as reflected by the energy and
The microscopic “spin density” is given by
magnetization conservation, an integration time step of
∆t=1fs was used.
ρ (r,t)= S (t)δ(r−r (t)). (8)
s i i For computing canonical averages of time-dependent
X
i
correlation functions, we used time series obtained from
Treating the spin-spin correlations along x, y, and z di- microcanonial dynamical simulations, that are, in turn,
rectionsseparately,wedefinetheintermediatescattering initiated from equilibrium states drawn from the canon-
function as ical ensemble at the desired temperature T. Averaging
over the results of multiple simulations started from dif-
1
Fsks(q,t)= N ρks(q,t)ρks(−q,0) , (9) ferentinitialstatesyieldsgoodestimatesoftherespective
(cid:10) (cid:11) canonical ensemble averages [12].
where k = x,y, or z, and ρ (q,t) = S (t)e−iq·ri(t). Forgeneratingtheinitialstatesforourmicrocanonical
s i i
ForaferromagneticsysteminthemicroPcanonicalensem- MD-SD simulations, we adhere to the following proce-
ble,themagnetizationvectorisaconstantofmotionand dure. First,weequilibratethesubspaceconsistingofpo-
servesasafixedsymmetryaxisthroughoutthetimeevo- sitionsandspinsusingtheMetropolisMonteCarlo(MC)
lution of the system. To differentiate between the mag- method[52]. As the secondstep, weassigninitialveloci-
netic excitations that propagate parallel and perpendic- tiestotheatomsbasedontheMaxwell-Boltzmanndistri-
ular to this symmetry axis, we redefine the coordinate butionatthedesiredtemperatureT. Finally,weperform
4
MDsimulations,weusedtheDudarev-Derletpotentialto
model the interatomic interactions while completely ne-
glecting the spin-spin interactions. SD simulations were
conducted with the atoms frozen at perfect bcc lattice
positions,andtheexchangeparametersdeterminedfrom
the same pairwise function used for MD-SD simulations.
III. RESULTS
A. 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 (q,ω), as
nn
well as in the longitudinal and the transverse compo-
FIG. 1. Time evolution of the instantaneous lattice and spin
nents of the velocity-velocity dynamic structure factor:
temperatures as observed in a microcanonical MD-SD simu-
lation for the system size L=16 at temperature T =800K. SvLv(q,ω) and SvTv(q,ω). For each q along [100] and
The initial state for the time integration was generated from [111] lattice directions, all three quantities show single
the procedure described in Sec. IIC. The spin tempera- peaks (See Fig. 2 for an example). The peak positions
ture was measured using the formula developed by Nurdin in S (q,ω) and SL(q,ω) for the same wave vector co-
nn vv
et al. [53]. incide with each other as they are both associated with
thelongitudinalvibrationalmodes,andhenceconveythe
same information. The peaks in ST (q,ω) are associated
vv
a short microcanonical MD-SD equilibration run (typi- with the transverse lattice vibrations. Since there are
cally ∼ 1000 time steps with ∆t = 1fs), which would two orthogonal directions perpendicular to a given wave
ultimately bring the whole system to the equilibrium by vector q, there are, in fact, two transverse vibrational
resolving any inconsistencies between the position-spin modes for each q. Due to the four-fold and three-fold
subspace and the velocity distribution. Fig. 1 shows the rotational symmetry about the axes [100] and [111], re-
time evolutionofthe instantaneouslatticeandspintem- spectively,thetwotransversemodesforthewavevectors
peratures as observedin a microcanonicalMD-SD simu- along these directions become degenerate [55]. As a re-
lationinitiatedfromanequilibriumstategeneratedfrom sult, we only observe a single peak in ST (q,ω) for the
vv
the aforementioned technique for T = 800K. Both lat- wave vectors along these directions.
tice andspintemperatures fluctuate about amean value We also observe single peak structures in S (q,ω)
nn
of T = 800K, indicating that the lattice and the spin and SL(q,ω) for the wave vectors along the [110] direc-
vv
subsystems are in mutual equilibrium. tion. However, for the case of ST (q,ω), one can clearly
vv
To characterize phonon and magnon modes, we per- identify two distinct peaks. This is a consequence of the
formed simulations for the system size L = 16 (8192 two transverse modes being non-degenerate due to the
atoms) at temperatures T = 300K, 800K, and 1000K. reduced rotational symmetry (two-fold) about the [110]
T = 1000K was particularly chosen due to its vicinity axis in comparison to [100] and [111] directions [55].
to the Curie temperature of bcc iron, TC ≈ 1043K. (A To extract the positions and the half-widths of the
recent high resolution Monte Carlo study has revealed phonon peaks, we fit the simulation results for the dy-
that the transition temperature of the particular spin- namic structure factor to a Lorentzian function of the
lattice model used in our study to be T ≈ 1078K [54]). form [13, 17]
Equations of motion were integrated up to a total time
of t = 1ns, and the space-displaced, time-displaced I Γ2
max S(q,ω)= 0 , (12)
correlationfunctions werecomputedforthe threeprinci- (ω−ω )2+Γ2
0
ple lattice directions: [100], [110] and [111]. To increase
the accuracy,we have averagedthese quantities overdif- whereω isthecharacteristicfrequencyofthevibrational
0
ferent starting points in the time series. Canonical en- mode, I is the intensity or the amplitude of the peak,
0
semble averages were estimated using the results of 200 andΓisthehalf-widthathalfmaximum(HWHM)which
independent simulations, each initiated from a different is inverselyproportionalto the lifetime ofthe excitation.
initial state. The time Fourier transform in Eq. (6) was Theerrorsofthefittingparameterswereestimatedusing
carried out to a cutoff time of t =0.5ns. the following procedure. The complete set of correlation
cutoff
As our primary goal is to understand the mutual im- function estimates obtained from 200 independent simu-
pact of the phonons and magnons on their respective lations was divided into 10 groups, and the data within
frequency spectra and lifetimes, we have also performed each group were averaged over to yield 10 results sets.
standalone MD and SD simulations for comparison. For Dynamicstructurefactorswereindependentlycomputed
5
0.0004 Γ H P Γ N
simulation 40
Lorentzian fit
TA LA
0.0003
ω = 30.543 ± 0.006 meV
0 30
ω) Γ = 0.268 ± 0.006 meV LA TA2
q, 0.0002 I = (3.24 ± 0.08) x 10-4 V)
S(nn 0 Γ I0 ω (me 20 TA LA
I/2
0
0.0001 TA1
10 Exp
MD-SD
MD
0 ω
28 29 30 31 32 33
0
0
ω (meV) 0 0.5 1 0.5 0 0.5
[q 0 0] [q q q] [q q 0]
FIG. 2. Density-density dynamic structure factor for q = q (2π/a)
(1.1˚A-1, 0,0) obtained from MD-SD simulations for L = 16
at T = 300K. The symbols represent simulation data while FIG.3. Comparisonofthephonondispersioncurvesobtained
the solid line is a fit with the Lorentzian lineshape given in from MD-SD simulations (L = 16) with the experimental
Eq. (12). results [56, 57] for T = 300K. Results obtained from pure
MDsimulationsarealsoplottedforcomparison. LAandTA,
respectively,denotethelongitudinalandtransversebranches.
for these 10 correlation function sets. To estimate the
errorsinthefitting parameters,weseparatelyperformed
curvefitstothese10independentdynamicstructurefac- and the lifetimes are practically infinite. As the tem-
tor estimates, and calculated the standard deviations of perature is increased, phonon occupation numbers also
the fitting parameters. Statistical errors bars obtained increase, which in turn increases the probability of mu-
in this manner were found to be an order of magnitude tual interactions. As a result of such phonon-phonon
larger than the error bars estimated by the curve-fitting scattering at elevated temperatures, characteristic fre-
tool. quenciesofthephononsmayshift, andthelifetimes may
For all the temperatures and wave vectors considered, shorten [58, 59]. In magnetic crystals, the co-existence
theLorentzianlineshapegiveninEq.(12)fittedwellwith of phonons and magnons gives rise to another class of
the peaks observed in S (q,ω) and SL,T(q,ω). Fig. 2 scattering processes, namely, phonon-magnon scatter-
nn vv
shows an example curve fit for the MD-SD results of ing. Just as phonon-phonon scattering, phonon-magnon
S (q,ω) for q = (1.1˚A-1,0,0) at T = 300K. To fit the scattering may also lead to a shift in the characteristic
nn
two peak structure observed in ST (q,ω) for the [110] phonon frequencies, as well as shortening of the phonon
vv
direction, we use the sum of two Lorentzians. lifetimes. Asthe occupancyofbothphononandmagnon
UsingthepeakpositionsobtainedfromtheLorentzian modes increases with temperature, these effects will be
fits, one can construct phonon dispersion relations for more pronounced as the temperature is increased.
the three principle lattice directions. Fig. 3 shows the To carefully examine the changes in the phonon fre-
the dispersion curves determined from our MD-SD sim- quency spectrum due to magnons, we compare the char-
ulations for T = 300K, along with the experimental re- acteristic frequencies determined from MD-SD simula-
sults [56, 57] obtained from inelastic neutron scattering. tions (ω ) with the ones obtained from MD simula-
MD-SD
For comparison, we have also shown the results of stan- tions (ω ) by calculating the fractionalfrequency shift,
MD
dalone MD simulations for the same temperature. In (ω −ω )/ω . Theresults forthe threeprinciple
MD-SD MD MD
general, for small to moderate q values, both MD-SD directions are shown in Figs. 4 and 5, for the longitu-
and MD dispersion curves agree well with the experi- dinal and the transverse modes, respectively. With the
mental results, but deviations can be observedfor larger exception of the high frequency transverse branch along
q values, particularly near the zone boundaries in [100] [110]direction(TA2), phononfrequencies shift to higher
and [111] directions. Although the MD-SD and MD dis- values in the presence of magnons. In general, the shift
persioncurvesareindistinguishablewithintheresolution in frequencies becomes more pronounced as the temper-
of Fig. 3, we will show later on that there are, in fact, ature is increased. A particularly interesting behavior
deviations larger than the error bars. occurs in the longitudinal branch for the [111] direction
Attemperaturesinthevicinityofabsolutezero,dueto whereweobservedipsinthecurvesforallthreetempera-
thelowoccupationofvibrationalmodes,phononsbehave turesatthesameqvalue. Forallthreetemperatures,the
as weakly interacting quasiparticles that can be treated frequency shift ofthe vibrationalmode that corresponds
within the harmonic approximation [58]. In this limit, to the bottom of the dip is close to zero. Therefore, the
characteristicfrequenciesofthephononsarewelldefined frequency of this phonon mode appears to be unaffected
6
(a) (a)
2.0 0.6
[100] [100]
%) T = 300 K %) T = 300 K
(MD 1.5 TT == 810000 0K K (MD 0.4 TT == 810000 0K K
ω ω
) /MD 1.0 ) /MD 0.2
ω ω
- D - D
D-S 0.5 D-S 0.0
M M
ω ω
( (
0.0 -0.2
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
(b) (b)
2.0 0.8
%) [110] TT == 380000 KK %) [110] (TA1)
(D 1.5 T = 1000 K (D 0.6
M M
ω ω
) /MD 1.0 ) /MD 0.4
ω ω
- D-SD 0.5 - D-SD 0.2 TTT === 381000000 0KK K
M M
ω ω
( (
0.0 0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5
(c) (c)
1.5 0.2
%) [111] T = 300 K %) [110] (TA2) T = 300 K
ω (MD 1.0 TT == 810000 0K K ω (MD 0.1 TT == 810000 0K K
) /MD ) /MD 0.0
ω ω
- MD-SD 0.5 - MD-SD-0.1
ω ω
( (
0.0 -0.2
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 0.5
q (2π/a)
(d)
1.0
FIG. 4. The fractional shift in longitudinal phonon frequen- %) T = 300 K [111]
T = 800 K
cies due to magnons for L = 16 at T = 300K, T = 800K, (D 0.8 T = 1000 K
and T =1000K in the (a) [1 0 0], (b) [1 1 0], and (c) [1 1 1] M
ω
lattice directions. ) /D 0.6
M
ω 0.4
by the presence of magnons. - D
S
Lifetimes of the phonon excitations are inversely pro- MD- 0.2
portionalto the half-widths athalf maximumofthe cor- ω
(
responding vibrational peaks observed in Snn(q,ω) and 0.0
SL,T(q,ω). To study the impact of the magnons on the 0.0 0.2 0.4 0.6 0.8 1.0
vv
phonon lifetimes, we compare the half-widths obtained q (2π/a)
from MD-SD simulations with that of the MD simu-
lations. Fig. 6 shows the results for the longitudinal FIG.5. Thefractional shift in transversephonon frequencies
phonons. For the longitudinal phonons at T = 300K, due to magnons for L = 16 at T = 300K, T = 800K, and
a marginal increase in the half-widths can be observed T = 1000K in the (a) [1 0 0], (b) [1 1 0] (TA1), (c) [1 1 0]
due tothemagnons,whichbecomesmorepronouncedas (TA2),and (d) [1 1 1] lattice directions.
the temperature is increased. For the case of transverse
phonons, we did not observe any notable difference be-
7
(a) B. Magnetic excitations
T = 300 K (MD-SD)
1.2 T = 300 K (MD)
T = 800 K (MD-SD) 1. Transverse magnon modes
T = 800 K (MD)
meV) 0.8 TT == 11000000 KK ((MMDD-)SD) For the temperatures T = 300K and T = 800K, our
M ( rneasmulitcssftorrutchtuerterafancstvoerrsSeTc(oqm,pωo)nsehnotwofatshiengslpeins-psipninwadvye-
H ss
W [100] peak, that can be fitted to a Lorentzian lineshape of
H 0.4 the form Eq. (12) (See Fig. 7 (a) for an example.). For
T = 1000K, we also observe a diffusive central peak at
ω =0,asobservedinneutronscatteringexperiments[60]
and previous SD studies [17]. This two-peak structure
0.0 can be best captured by a function of the form [17]
0.0 0.2 0.4 0.6 0.8 1.0
(b) S(q,ω)=I exp −ω2/ω2 + I0Γ2 , (13)
T = 300 K (MD-SD) c (cid:0) c(cid:1) (ω−ω0)2+Γ2
T = 300 K (MD)
1.2
T = 800 K (MD-SD) where the first term (Gaussian) corresponds to the cen-
T = 800 K (MD)
tralpeak,andthesecondterm(Lorentzian)describesthe
V) T = 1000 K (MD-SD)
e T = 1000 K (MD) spinwavepeak(SeeFig.7(b)foranexample.). Forlarge
m
M ( 0.8 qinvSaTlu(eqs,aωt)Twe=re8f0o0uKndatnodbTe a=sy1m00m0eKtr,icsp,ianndwahveencpeeadkids
H ss
W not yield good fits to Lorentzian lineshapes. Therefore,
[110]
H one cannotobtain reliable estimates of the magnonhalf-
0.4
widths. However, spin wave peak positions can still be
determined relatively precisely, thus the magnon disper-
sion relations can be constructed.
0.0 Fig. 8 shows the transverse magnon dispersion re-
0.0 0.1 0.2 0.3 0.4 0.5
lations for small |q| values along the three principle
(c) 1.5 directions as determined from MD-SD simulations at
T = 300 K (MD-SD)
T = 300K. In agreement with the experimental find-
T = 300 K (MD)
T = 800 K (MD-SD) ings [61, 62], the three dispersion relations are isotropic
T = 800 K (MD) when plotted as functions of the magnitude of the wave
T = 1000 K (MD-SD) vector |q|. Moreover, for small |q| values, our results
eV)1.0 [111] T = 1000 K (MD) agreequantitativelywiththeexperimentalresultsforthe
m
M ( [110]direction[61,62]. Fig.9showsthecompletedisper-
sioncurvesdeterminedfromMD-SDandSDsimulations
H
W for T = 300K, T = 800K, and T = 1000K. For both
H0.5
MD-SD and SD, the characteristic frequencies shift to
lower values as the temperature is increased. This in-
dicates increased magnon-magnonscattering at elevated
temperatures. For all three temperatures, particularly
0.0
nearthe zoneboundaries,wecanobserveamarginaldif-
0.0 0.2 0.4 0.6 0.8 1.0
ference between the MD-SD and SD dispersion curves.
q (2π/a)
This, in fact, is a result of phonon-magnonscattering.
To further investigate the magnon softening due to
FIG. 6. Half-width at half maximum (HWHM) of the longi-
phonons,wecalculatethefractionalfrequencyshiftofthe
tudinal phonons at T = 300K, T = 800K, and T = 1000K
obtainedfrom MD-SDandMDsimulations forL=16inthe magnons,(ωMD-SD −ωSD)/ωSD. The results areshownin
(a) [1 0 0], (b) [1 1 0], and (c) [1 1 1] lattice directions. 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 be-
comesmore pronouncedas the temperature is increased.
Fig. 11 compares the transverse magnon half-widths
obtainedfromMD-SDandSDsimulationsforT =300K.
Although the difference between the half-widths is neg-
tween the MD-SD and MD half-widths outside the error ligible for small q values, for moderate to large q val-
bars, for all the temperatures considered. ues, half-widths for the MD-SD results are significantly
8
(a) Γ H P Γ N
500
0.0008 simulation 300 K (MD-SD)
300 K (SD)
Lorentzian fit
800 K (MD-SD)
400
800 K (SD)
ω = 270.57 ± 0.07 meV 1000 K (MD-SD)
0.0006 0 1000 K (SD)
Tω)S(q, ss0.0004 IΓ0 == (47.4.473 ±29 0 ±.0 02. 0m3e) Vx 10-4 Γ I0 ω (meV)230000
I/2
0
0.0002
100
2020 240 260 ω 280 300 320 0
0 0 0.5 1 0.5 0 0.5
ω (meV) [q 0 0] [q q q] [q q 0]
q (2π/a)
(b)
0.0008
simulation FIG. 9. Transverse magnon dispersion curvesfor T =300K,
curve fit T = 800K, and T = 1000K obtained from MD-SD and SD
0.0006 simulations for L=16.
ω = 104.5 ± 0.4 meV
0
ω) Γ = 25 ± 1 meV
TS(q, ss0.0004 I0 = (6.4 ± 0.2) x 10-4 limacaragngetnrosthnhoasrcntaettnhtienargtinfoogfr.tthheemSDagrneosnulltisfe.tTimheissidnudeictaotepshsoignnoinf--
ω = 93 ± 1 meV
c
0.0002 I = (7.9 ± 0.1) x 10-5
c
2. Longitudinal magnon modes
0
0 100 200 300 Our results for the longitudinal spin-spin dynamic
ω (meV) structure factor SL(q,ω) obtained from both MD-SD
ss
andSDsimulationsshowmanyverylow-intensityexcita-
FIG. 7. Transverse spin-spin dynamic structure factor ob- tionspeaks,forallwavevectorsconsidered. Fig.12shows
tained from MD-SD simulations for L = 16. (a) T = 300K SL(q,ω) for a small system size L = 8 at T = 300K,
and q = (1.1˚A-1, 0,0), (b) T = 1000K and q = (0.82˚A-1, ss
where we compare the SD results [panel (a)] with the
0,0). The symbols represent simulation data while the solid MD-SD results [panel (b)] for q= 2π(1,0,0).
lines are fits to functional forms presented by Eq. (12) [for La
In the context of classical Heisenberg models, Bunker
panel (a)] and Eq. (13) [for panel (b)].
et al. [18] showed that the excitation peaks observed in
160 SsLs(q,ω)aretwo-spin-wavecreationand/orannihilation
peakswhichresultfromthepairwiseinteractionsbetween
140 transverse magnon modes. For ferromagnetic systems,
MD-SD [1 0 0]
only spin wave annihilation peaks are present, and their
120 MD-SD [1 1 0]
frequencies are given by
MD-SD [1 1 1]
100
V) Exp [1 1 0] (Lynn) ω−(q ±q )=ω(q )−ω(q ), (14)
me 80 Exp [1 1 0] (Collins) ij i j i j
ω ( where q and q are the wave vectors of the two trans-
60 i j
verse magnon modes which comprise the two-spin-wave
40 excitation. Since the set of allowable wave vectors {qi}
depends on the system size L, the resultant two-spin-
20
wave spectrum also varies with L. For a real magnetic
crystal where L is practically infinite, the two-spin-wave
0
0 0.2 0.4 0.6 0.8 spectrum would become continuous.
-1
|q| (Å ) ToverifywhetherthepeaksweobserveinSL(q,ω)are
ss
two-spin-wave peaks, we chose a relatively small system
FIG.8. TransversemagnondispersionrelationsatT =300K size (L = 8) so that the set of allowable wave vectors is
obtained from MD-SD simulations for L = 16. The experi- reduced to a manageable size. Then, using MD-SD and
mental results reported by Lynn [61] and Collins [62] for the
SDsimulations,we separatelydeterminedthe transverse
[1 1 0] direction are also plotted for comparison.
magnon frequencies that correspond to the first few n
q
9
(a) (a)
10
%) 4 TT == 380000 KK [100] 8 MSDD-SD
(D T = 1000 K V)
S e
ω 2 m 6
) / D M (
ωS H 4 [100]
- D 0 HW
S
D- 2
M
ω -2
(
0
0.0 0.2 0.4 0.6 0.8 1.0
0.0 0.2 0.4 0.6 0.8 1.0
(b)
(b) 7
4
[110] 6 MD-SD
ωω) / (%)SDSD 231 TTT === 381000000 0KK K WHM (meV) 435 [110]SD
- D 0 H 2
S
D- 1
M
ω -1
( 0
0.0 0.1 0.2 0.3 0.4 0.5
-2
0.0 0.1 0.2 0.3 0.4 0.5
(c)
8
(c)
5 MD-SD
[111] 6 SD
%) 4 T = 300 K V)
ω (SD 3 TT == 810000 0K K M (me 4
) / D 2 WH [111]
ωS 1 H 2
- D
S 0
D-
M
ω -1 0
( 0.0 0.2 0.4 0.6 0.8 1.0
-2 q (2π/a)
0.0 0.2 0.4 0.6 0.8 1.0
q (2π/a)
FIG.11. Half-widthathalfmaximum(HWHM)ofthetrans-
verse magnons at T = 300K obtained from MD-SD and SD
FIG.10. Thefractionalshiftintransversemagnonfrequencies
simulations for L=16 in the (a) [1 0 0], (b) [1 1 0], and (c)
due to phonons for L = 16 at T = 300K, T = 800K, and
[1 1 1] lattice directions.
T =1000Kinthe(a)[100],(b)[110],and(c)[111]lattice
directions.
observedpeaksandthepredictedtwo-spin-wavepeakpo-
sitions, with the exception of the particular sharp peak
values along all possible lattice directions. With this in-
at ω ≈ 10meV which only appears in panel (b). Sur-
formationathand,wecanpredicttheexpectedpositions
prisingly,thepositionofthispeakcoincideswiththe fre-
ofallthetwo-spin-waveannihilationpeaksusingEq.(14)
quency of the longitudinal phonon mode for the same
forbothSDandMD-SDcase. Asanexample,letuscon- q as determined from the peak position of S (q,ω) or
nn
siderthewavevectorpairqi =(1,1,1)andqj =(1,1,0). SL(q,ω). Similar excitation peaks were observed for
vv
Sinceq −q =(0,0,1),theyproduceaspinwaveannihi-
i j all wave vectors, for all system sizes and temperatures
lation peak in SL(q,ω) for q=(0,0,1)at the frequency
ss considered. The origin of these coupled phonon-magnon
ω− = ω(q )− ω(q ). (Note that we have ignored the
i j modes can be explained as follows. Unlike transverse
common pre-factor 2π/La from the wave vectors.)
phonons,whenalongitudinalphononpropagatesalonga
In Fig. 12 (a) and (b), we have superimposed the pre- certain lattice direction, it generates fluctuations in the
dictedspinwaveannihilationpeakpositionscorrespond- local atom density along that direction with the corre-
ing to eachcase. We see an excellent match between the sponding phonon frequency. This, in turn, leads to fluc-
10
(a)
-2
10
-3
ω) 10
q,
L (
Sss
-4
10
0 10 20 30 40 50 60
ω (meV)
(b)
-2
10
0.002
ω)
q,
S(nn
ω) 10-3 08.7 9.2 9.8
q, ω (meV)
L (
Sss
-4
10
8.7 9.2 9.8
0 10 20 30 40 50 60
FIG. 12. The longitudinal component of the spin-spin dy- ω (meV)
namic structure factor SsLs(q,ω) for q = L2πa(1,0,0) obtained
from (a) SD and (b) MD-SD simulations for L = 8 at
FIG. 13. The longitudinal component of the spin-spin dy-
Thil=ati3o0n0pKe.aTkshearpereinddicitceadtepdobsiytiothnesovfertthiecatlwdoa-ssphiend-wlianvees.aTnnhie- namic structure factor SsLs(q,ω) for q = L2πa(1,0,0) obtained
from (a) SD and (b) MD-SD simulations for L = 8 at
dottedline[LA(100)]marksthefrequencyofthelongitudinal
T = 800K. The inset of (b) shows the density-density dy-
phonon mode for the same q.
namic structurefactor for thesame q.
tuations in the local density of the longitudinal compo-
nents of the spins (i.e. components of the spin vectors
parallel to the net magnetization). These longitudinal withincreasingtemperature,asthediffusivecentralpeak
spinfluctuationspropagatealongwiththephonon,yield- becomes more pronounced. At T = 1000K, the inten-
ingasharp,coupledmodeinthethelongitudinalmagnon sity of the peak is very low and is barely recognizable.
spectrum. Above the Curie temperature, spins are randomly ori-
Fig. 13 and Fig. 14 show SL(q,ω) for q = 2π(1,0,0) ented and the vector sum of spins per unit volume will
ss La
at T = 800K and T = 1000K, respectively, where we be zero on average. Hence, the coupled phonon-magnon
compare the SD results [panel (a)] with the MD-SD re- mode should entirely disappear; however, it is already
sults [panel (b)]. In each figure, the inset of panel (b) so faint at T = 1000K that we clearly would not have
showsthelongitudinaldensity-densitydynamicstructure sufficientresolutiontotestthisbehaviorabovethe Curie
factor for the same wave vector. In comparison to the temperature.
results for T = 300K, we observe that the diffusive cen-
tral peak becomes more pronounced as the temperature We would like to point out that the existence of these
rises,andmanyofthe low-intensitytwo-spin-wavepeaks coupledphonon-magnonmodes isa phenomenonthatso
broaden and disappear into its tail. These observations farhasn’tbeendiscoveredexperimentally. Infact,thisis
are in qualitative agreement with previous SD studies not surprising since the experimental detection of these
of the ferromagnetic Heisenberg model [18]. The cou- peaks would be extremely challenging due to their very
pledphonon-magnonpeak alsobecomes lesspronounced low intensities.