Table Of ContentSpin waves in paramagneticBCC iron: spin dynamics simulations
Xiuping Tao,1,∗ D. P. Landau,1 T. C. Schulthess,2 and G. M. Stocks3
1Center for Simulational Physics, University of Georgia, Athens, GA 30602
2ComputerScienceandMathematicsDivision, OakRidgeNationalLaboratory, OakRidge, TN37831-6164
3Metals and Ceramics Division, OakRidge National Laboratory, Oak Ridge, TN 37831-6114
5 (Dated: February2,2008)
0
0 Largescalecomputersimulationsareusedtoelucidatealongstandingcontroversyregardingtheexis-
2 tence,orotherwise,ofspinwavesinparamagneticBCCiron. Spindynamicssimulationsofthedynamic
structurefactorofaHeisenbergmodelofFewithfirstprinciplesinteractionsrevealthatwelldefinedpeaks
n
a persistfaraboveCurietemperatureTc. Atlargewavevectorsthesepeakscanbeascribedtopropagating
J spinwaves,atsmallwavevectorsthepeakscorrespondtoover-damped spinwaves. Paradoxically, spin
8 waveexcitationsexistdespiteonlylimitedmagneticshort-rangeorderatandaboveTc.
2
PACSnumbers:75.10.Hk,75.40.-s,75.40.Gb,75.50.Bb
Keywords:transitionmetal,spindynamics,Heisenbergmodel,short-rangeorder
]
i
c
s
- For over three decades, the nature of magnetic excita- propagatingmodes.
rl tions in ferromagneticmaterials abovethe Curie tempera- With new algorithmic and computational capabilities,
mt ture Tc has been a matter of controversy amongst experi- qualitativelymoreaccurateSDsimulationscannowbeper-
mentalistsandtheoristsalike. Earlyneutronscatteringex- formed. In particular, it can follow many more spins for
.
t periments on iron suggested that spin waves were renor- much longer integration time. We use these techniques
a
m malized to zero at Tc [1]; however, in 1975, using unpo- andamodeldesignedspecificallytoemulateBCCironand
larized neutron scattering techniques, Lynn at Oak Ridge havebeenable to unequivocallyidentify propagatingspin
-
d (ORNL) reported [2] that spin waves in iron persisted as wavemodesin theparamagneticstate, lendingsubstantial
n excitationsuptothehighesttemperaturemeasured(1.4T ), supportto Lynn’s[2]experimentalfindings. Interestingly,
c
o
and no further renormalization of the dispersion relation spinwavesarefounddespiteonlylimitedmagneticSRO.
c
wasobservedaboveT . To describe the high temperature dynamics we use a
[ c
1 Experimentally, it was challenged primarily by Shirane classicalHeisenbergmodelH = −(1/2) r6=r′Jr,r′Sr ·
v and collaborators at Brookhaven(BNL) [3]. Using polar- Sr′,forwhichtheexchangeinteractions,JrP,r′,areobtained
3 ized neutrons, they reported that spin wave modes were from first principles electronic structure calculations. For
1 not presentaboveT and suggested that the ORNL group Fe this is a reasonableapproximationsince thesize ofthe
c
7 neededpolarizedneutronstosubtractthebackgroundscat- magnetic moments associated with individual Fe-sites are
1 tering properly. Utilizing full polarization analysis tech- only weakly dependenton the magnetic state [14] and by
0
niquestheORNLgroupsubsequentlyconfirmedtheirear- including interactions up to fourth nearest neighbors it is
5
0 lier work and, in addition, they analyzed data from both possibletoobtainareasonablygoodTc.
/ groups and concludedthat their resolution was more than Large scale computer simulations using SD techniques
t
a an order of magnitude better than that employed by the to study the dynamic properties of Heisenberg ferromag-
m
BNL researchers [4]. Moreover, angle-resolved photoe- nets[15]andantiferromagnets[16]havebeenquiteeffec-
- missionstudies[5, 6]suggestedtheexistenceofmagnetic tive,andthedirectcomparisonofRbMnF SDsimulations
d 3
short-rangeorder(SRO)inparamagneticironandthatthis with experimentswas especially satisfying [16]. We have
n
could giverise to propagatingmodes. Theoretically, SRO adoptedthesetechniquesandusedL L LBCClattices
o
c ofratherlong lengthscales(25 A˚)was postulatedto exist withperiodicboundaryconditionsan×dL×= 32and40. At
: far above T [7, 8] and a more subtle kind was proposed eachlatticesite,thereisathree-dimensionalclassicalspin
v c
i later [9]. Contrarily, it was also suggested that above T , ofunitlength(weabsorbspinmomentsintothedefinition
X c
all thermal excitations are dissipative [10, 11]. To further of the interaction parameters)and each spin has a total of
r complicate matters, analytical calculations for a Heisen- 50 interacting neighbors. We use interaction parameters,
a
berg model of iron, with exchange interactions extending J ,fortheT = 0ferromagneticstateofBCCFecalculated
i
to fifth-nearest neighbors and a three pole approximation using the standard formulation [17] and the layer-KKR
[12], did not reproduce the line shape measured by either method [18]. The calculated values are J = 36.3386
1
experimentalgroupmentionedabove. Inaddition,Shastry meV, J = 20.6520 meV, J = 1.625962 meV, and
2 3
−
[13]performedspindynamics(SD)simulationsofanear- J = 2.39650meV.
4
−
est neighborHeisenberg model of paramagneticiron with In our simulations, a hybrid Monte Carlo method was
8192 spins and showed some plots of dynamic structure used to study the static properties and to generate equi-
factor S(q,ω) with a shoulder at nonzero ω for some q. librium configurations as initial states for integrating the
It was explained to be due to statistical errors instead of coupled equations of motion of SD [19]. At T and for
c
2
L = 32, the measured nonlinear relaxation time in the
8
equilibrating process and the linear relaxation time be-
0.95T
tween equilibrated states for the total energy and for the 1.00Tc
c
magnetization[20]arebothsmallerthan500 hybridsteps 1.10T
6 1.20Tc
per spin. We discarded 5000 hybrid steps (for equilibra- c
tion) and used every 5000th hybrid step’s state as an ini-
tial state for the SD simulations. For the Ji’s used here, ω)
T = 919(1)K, which is slightly smaller than the exper- q,4
c (
imental value Texp = 1043K. The equilibrium magne- S
c
tization m (1/N) S (1 T/T )1/3 in the
| | ≡ | r r| ∼ − c
vicinityofT andthisisPinagreementwithexperiments. 2
c
TheSDequationsofmotionare
dS
dtr = Heff ×Sr, (1) 00 100 200 300
ω (meV)
where Heff ≡ − r′Jr,r′Sr′ is an effective field at site
r due to its interacPting neighbors. The integration of the FIG. 1: Calculated energy dependence of S(q,ω) at q =
equationsdeterminesthetimedependenceofeachspinand π/a(1,0,0) and for T = 0.95Tc (700 runs), 1.0Tc (2000 runs),
was carriedoutusing analgorithm basedonsecond-order 1.1Tc (2240runs),and1.2Tc (2240runs)forL = 32. Thesolid
Suzuki-Trotterdecompositionsofexponentialoperatorsas linesarefitstothedataasexplainedintext.Errorbarsareshown
describedin[21]. Thealgorithmviewseachspinasunder- atafewtypicalpoints.
going Larmor precession around its effective field H ,
eff
which is itself changing with time. To deal with the fact
socalled,constant-qscansareforq= π/a(1,0,0)(q =
thatweareconsideringfourshellsofinteractingneighbors,
1.09 A˚−1), whichishalfwaytoBrillouinZoneboun|da|ry.
the BCC lattice is decomposed into sixteen sublattices.
At 0.95T , S(q,ω) already has a 3-peak structure: one
Thisalgorithmallowstimestepsaslargeasδt = 0.05(in c
units oft = J−1). Typically, the integrationwas carried weak central peak at zero energy and two symmetric spin
0 1 wavepeaks(weonlyshowdataforω 0sincethestruc-
outtot = 20000δt = 1000t .
max 0 ≥
ture factoris symmetric aboutω = 0). Note thatthe spin
The space- and time-displaced spin-spin correlation
functionCk(r r′,t)andtherelateddynamicalstructure wave peaks are already quite wide. As T goes to Tc and
factor, Sk(q,ω−), are fundamentalin the study ofspin dy- above,thecentralpeakbecomesmorepronounced. Inad-
dition,thespinwavepeaksshifttolowerenergies,broaden
namics[22]andaredefinedas
furtherandbecomelessobvious,howevertheystillpersist.
Ck(r r′,t) = Srk(t)Sr′k(0) Srk(t) Sr′k(0) , This3-peakstructureathightemperaturesisincontrastto
− h i−h ih i
(2)
the2-peakspinwavestructure foundatlow temperatures.
where k = x, y or z and the angle brackets denote In the neutron scattering from 54Fe(12%Si) experiments
h···i
theensembleaverage,and
[4], MookandLynnalsonoticedacentralpeak,butcould
Sk(q,ω) = eiq·(r−r′) +∞eiωtCk(r r′,t) dt , nooftaldleocyiidnegwofhseitlhiceornit.was intrinsic to pure iron or a result
Xr,r′ Z−∞ − √2π Ingeneral,constant-qscansareisotropicinthe(q,0,0),
(3) (q,q,0), and(q,q,q)directions. Forverysmall q, there
whereqandω aremomentumandenergy(E ω)trans- | |
is onlya centralpeakinthescans(asis expected)andthe
∝
ferrespectively. ItisSk(q,ω)thatwasprobedintheneu- 3-peakstructure onlydevelopsfor larger q. We fit the 3
tronscatteringexperimentsdiscussedearlier. peaksinS(q,ω)usingdifferentfittingfun|cti|onsandfound
By calculating partial spin sums ‘on the fly’ [15],
thebestresultswitheitheraGaussiancentralpeakplustwo
it is possible to calculate Sk(q,ω) without storing a
Lorentzianpeaksat ω :
0
huge amount of data associated with each spin config- ±
S(q,ω) = G+L +L , (4)
uration. Because L is finite, only a finite set of q + −
values are accessible: q = 2πnq/(La) with nq = or a Gaussian central peak plus two additional Gaussian
1, 2,..., L for the (q,0,0) and (q,q,q) directions peaksat ω :
a±nd n± = 1±, 2,..., L/2 for the (q,q,0) direction. ± 0
q ± ± ± S(q,ω) = G+G +G , (5)
(a is lattice constant.) ForT T , the ensembleaverage + −
c
≥
inEq.2wasperformedusingatleast2000startingconfig- whereG = I exp( ω2/ω2),L = I ω2/((ω ω )2 +
c − c ± 0 1 ∓ 0
urations. We averageSk(q,ω) overequivalentdirections ω2),andG = I exp( (ω ω )2/ω2). Formoderate q
1 ± 0 − ∓ 0 1 | |
andthisaveragedstructurefactorisdenotedasS(q,ω). theresultsarefitbestwithEq.4,whileEq.5worksbetter
InFig.1weshowthefrequencydependenceofS(q,ω) atlarger q. InFig.2weshow,forT = T ,theresultsof
c
obtainedforfourdifferenttemperaturesaroundT . These, fittingco|ns|tant-qscansat q = 0.48A˚−1 and q =1.16
c
| | | |
3
200 120
(a) simulation 0.3T (exp)
Eq.4 fitting c
T--1.4T (exp)
150 GL++L- 100 0.c3Tc(sicm)
1.0T (sim)
c
ω) 1.1T (sim NSW)
q,100 1.1Tc(sim)
S( 80 c
1.2T (sim NSW)
c
1.2T (sim)
50 ) c
V
me 60
(
0 E
0 20 40 60
ω (meV)
6 40
(b) simulation
Eq.5 fitting
5
G
G +G
+ - 20
4
)
ω
q,3
(
S 0
0 0.2 0.4 0.6 0.8 1
2
-1
|q| (Å )
1
FIG. 3: Comparison of dispersion curves obtained in our sim-
0
0 100 200 300 ulations (sim)withLynn’s experimental (exp)(Ref. [2])results
ω (meV) forthe(q,q,0)direction. Opensymbolsindicateexcitationswith
mixednatureandarenotduetospinwaves(NSW).
FIG. 2: Fits toS(q,ω)at T = Tc for two|q|-points along the
(q,q,0) direction for L = 32. (a) |q| = 0.48 A˚−1 fit to Eq. 4,
withIc =58.7,ωc =27.6meV,I0 =80.6,ω0 =21.4meV,and tail at slightly different ωmax to get an average ω0; these
ω1 = 11.3meV;(b)|q|= 1.16A˚−1 fittoEq.5withIc = 2.49, error bars are found to be no larger than symbols. In this
ωc = 193.3 meV, I0 = 3.02, ω0 = 175.3 meV, and ω1 = 61.2 figure,filledsymbolsindicatemodesthatareclearlypropa-
meV. The vertical scale in (b) is much smaller than that in (a). gating(ω /ω < 1)whileopensymbolsindicatethat,even
1 0
Errorbarsareshownatafewtypicalpoints.
though there are peaks at ω = 0, the peaks have widths
0
6
ω > ω . The calculated result for T = 0.3T is very
1 0 c
closetothatfromtheexperimentsandpropagatingmodes
A˚−1 in the (q,q,0) direction. The q = 0.48 A˚−1 result existforverysmall q. ForT T ,ourcurvesliebelow
| | | | ≥ c
fits well to Eq. 4 and has ω /ω < 1, i.e., the excitation theexperiments’sandsoftenwithincreasingtemperatures,
1 0
lifetimeislongerthanitsperiodandthusitcanberegarded apropertynotseenintheexperiments. Onepossibilityde-
as a spin wave excitation. It should be noted that this q servingoffurtherstudyis thatouruseoftemperatureand
valueisveryclosetothat(0.47A˚−1)forwhichLynnfou|nd| configurationindependentexchangeinteractions,inpartic-
propagating modes in contradiction to the findings of the ular those appropriate to the T = 0 ferromagnetic state,
BNL group. At q = 1.16 A˚−1, the structure factor has breaksdownathightemperatureswhenthespinmoments
| |
much weaker intensity and fits best to Eq. 5 with a ratio arehighlynon-collinear.
ω1/ω0thatisevensmallerthanat q = 0.48A˚−1. Thisis In our simulations we have equal access to constant-q
| |
illustrative of the general conclusion that the propagating scans and constant-E scans; however, this is not the case
natureoftheexcitationmodesismostpronouncedatlarge in neutronscattering experiments. Becausethe dispersion
q. curvesofFearegenerallyverysteep,experimentalistsusu-
| |
Figure3showsthedispersionrelationsobtainedbyplot- allyperformconstant-Escans.InFig.4weshowconstant-
ting the peak positions, ω , determined from the fits to E scansforseveralE valuesatT = 1.1T basedonsim-
0 c
S(q,ω) along the (q,q,0) direction. Calculated disper- ulations. Clearly, the constant-E scans have two peaks
sion curves are shown at several temperatures in the fer- (symmetricabout q = 0)thatbecomesmallerandwider
| |
romagneticandparamagneticphasestogetherwiththeex- andshifttohigher q asE increases. Peaksinconstant-E
| |
perimental results of Lynn [2]. To estimate errors, we fit- scansstronglysuggestthatSROpersistsaboveT [7].
c
ted each constant-q scan several times by cutting off the The degree of magnetic SRO can be obtained directly
4
25 performedatORNL-CCS(www.ccs.ornl.gov)andNERSC
(www.nersc.gov) using elements of the Ψ-Mag toolset
41.1 meV which may be obtained at http://mri-fre.ornl.gov/psimag.
20 54.8 meV
WorksupportedbyComputationalMaterials ScienceNet-
68.5 meV
96.0 meV work sponsored by DOE BES-DMSE (XT), BES-DMSE
)15 (GMS), NSF Grant No. DMR-0341874 (XT and DPL)
ω
q, andDARPA(XTandTCS)undercontractNo. DE-AC05-
(
S10 00OR22725withUT-BattelleLLC.
5
0 ∗ Electronicaddress: [email protected]
0 0.2 0.4 0.6 0.8 1
[1] M.F.Collins,V.J.Minkiewicz,R.Nathans,L.Passell,and
|q|/q
zb G.Shirane,Phys.Rev.179,417(1969).
[2] J.W.Lynn,Phys.Rev.B11,2624(1975).
FIG.4: T =1.1Tc constant-E scansalong(q,q,0)directionfor [3] J.P.Wicksted, G.Shirane,andO.Steinsvoll,Phys.Rev.B
E =41.1meV,54.8meV,68.5meV,and96.0meVwithL=40. 29, R488 (1984); G. Shirane, O. Steinsvoll, Y. J. Uemura,
BrillouinZoneboundaryqzb =1.55A˚−1inthedirection. andJ.Wicksted,J.Appl.Phys.55,1887(1984);J.P.Wick-
sted,P.Boni,andG.Shirane,Phys.Rev.B30,3655(1984).
[4] H.A.MookandJ.W.Lynn,J.Appl.Phys.57,3006(1985).
from the behavior of static correlation function Ck(r [5] E. M. Haines, V. Heine, and A. Ziegler, J. Phys. F: Metal
r′,0)(i.e. Eq.2witht =0),whichcanbecalculatedfro−m Phys. 15, 661 (1985); E. M. Haines, R. Clauberg, and
theMonteCarlo configurationsalone. ForT = 1.1T we R.Feder,Phys.Rev.Lett.54,932(1985).
c
find a correlationlengthofapproximately2a ( 6 neigh- [6] E. Kisker, R.Clauberg, and W. Gudat, Z.Phys. B 61, 453
bor shells), indicative of only limited SRO. Th∼us, in gen- (1985).
[7] V.Korenman,J.L.Murray,andR.E.Prange,Phys.Rev.B
eral, extensiveSRO is notrequiredto supportspin waves.
16,4032(1977);R.E.PrangeandV.Korenman,Phys.Rev.
Moreover, inspection of Fig. 3 for T = 1.1T shows that
c B19,4691(1978).
the point q & 0.77A˚−1, at which these peaks first cor-
[8] H.Capellmann,SolidStateCommun.30,7(1979);Z.Phys.
respond to propagating modes, is when their wavelength
B35,269(1979).
(λ 2a) first becomes the order of the static correlation [9] V.HeineandR.Joynt,Europhys.Lett.5,81(1988).
∼
length. [10] J.Hubbard,Phys.Rev.B20,4584(1979);Phys.Rev.B23,
Insummary,ourSDsimulationsclearlypointtotheex- 5974(1981).
istenceofspinwavesintheparamagneticstateofBCCFe [11] T.Moriya,J.Magn.Magn.Mat.100,261(1991).
and support the original conclusions of Lynn. Their sig- [12] B.S.Shastry,D.M.Edwards, andA.P.Young, J.Phys.C
nature is seen as spin wave peaks in dynamical structure 14,L665(1981).
factorinconstant-qandconstant-Escans. Detailedanaly- [13] B.S.Shastry,Phys.Rev.Lett.53,1104(1984).
[14] A. J. Pindor, J. Staunton, G. M. Stocks, and H. Winter, J.
sis ofthe constant-qscansshowsthatthe propagatingna-
Phys.F13,979(1983).
ture of these excitations is clearest at large q, in agree-
| | [15] K.ChenandD.P.Landau,Phys.Rev.B49,3266(1994).
mentwith experiment. This is alsoconsistentwith the re-
[16] S.-H.Tsai,A.Bunker, andD.P.Landau, Phys.Rev.B61,
quirementthattheirwavelengthbetheorderof,orshorter
333(2000).
than, the static correlation length. While the inclusion of [17] A. I. Liechenstein, M. I. Katsnelson, V. P. Antropov, and
four shells of first-principles-determined interactions into V.A.Gubanov,J.Mag.Mag.Mater.67,65(1987).
theHeisenbergmodelmakesourresultsspecificallyrelate [18] T. C.Schulthess and W. H. Butler, J. App. Phys. 83, 7225
to BCC Fe, we have also found spin waves in a Heisen- (1998).
berg model containing only nearest neighborinteractions. [19] P.PeczakandD.P.Landau,Phys.Rev.B47,14260(1993).
Inadditiontoelucidatingthelongstandingcontroversyre- [20] D.P.LandauandK.Binder,AGuidetoMonteCarloSim-
ulationsinStatisticalPhysics(CambridgeUniversityPress,
gardingtheexistenceofspinwavesaboveT ,thesesimula-
c Cambridge,2000).
tionsalsopointtotheimportantrolethatinelasticneutron
[21] M. Krech, A. Bunker, and D. P. Landau, Comput. Phys.
scatteringstudiesoftheparamagneticstatecanhaveinun-
Commun.111,1(1998).
derstandingthenatureofmagneticexcitations,particularly
[22] S.W.Lovesey,Theoryofneutronscatteringfromcondensed
whencoupledwithstate-of-the-artSDsimulations. matter(Clarendon,Oxford,1984).
WethankShan-HoTsai,H.K.Lee,V.P.Antropov,and
K.Binderforinformativediscussions. Computationswere
