Table Of Content

Spin 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 initial 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