Table Of Content

Ab Initio Calculation of Spin Gap Behavior in CaV O 4 9 C. Stephen Hellberg,1 W. E. Pickett,1,2 L. L. Boyer,1 Harold T. Stokes,3 and Michael J. Mehl1 1Center for Computational Materials Science, Naval Research Laboratory, Washington DC 20375 2Department of Physics, University of California, Davis CA 95616 3Department of Physics and Astronomy, Brigham Young University, Provo UT 84602 (February 1, 2008) 9 9 SecondneighbordominatedexchangecouplinginCaV4O9 hasbeenobtainedfromabinitiodensity 9 functional (DF) calculations. A DF-based self-consistent atomic deformation model reveals that 1 the nearest neighbor coupling is small dueto strong cancellation among thevarious superexchange n processes. ExactdiagonalizationofthepredictedHeisenbergmodelyieldsspin-gapbehavioringood a agreement with experiment. Themodelis refinedbyfittingtotheexperimentalsusceptibility. The J resulting model agrees very well with the experimental susceptibility and triplet dispersion. 4 2 PACS numbers: 75.10.Jm, 75.40.Cx, 75.50.Ee ] i c CaV O wasthefirsttwo-dimensionalsystemobserved 4 9 s to enter a low-temperature quantum-disordered phase 5 6 4 3 1 - l with a spin gap ∆ ≈ 110K. The gap was first apparent r t in its susceptibility, which vanishes at low temperatures m as χ(T → 0) ∼ exp(−∆/kT) [1], and was observed di- 8 7 5 6 4 . t rectly in the dispersion of triplet spin excitations (ΩQ) a measuredbyneutronscattering[2]. Thisunexpectedbe- m havior has stimulated considerable theoretical study of 1 2 8 7 - d the exchange couplings between S=1 spins on the V lat- 2 n tice using Heisenberg models [3–6]. o CaV4O9 is a layered compound—the interlayer dis- 6 4 3 1 2 c tanceissufficientlylargetomakeinterlayerV-Vcoupling [ negligible. Withinalayer,theVatomsforma 1-depleted 5 1 square lattice shown as the circles in Fig. 1 [7,8]. The 7 5 6 4 3 v lattice was originally viewed as an array of square “pla- 5 quettes” of V ions (e.g., 1-2-3-4 in Fig. 1) tending to- 5 2 ward singlet formation since isolated plaquettes have a J J’ 1 1 1 singlet ground state. Examination of the structure how- 0 ever suggestsintra- andinter-plaquette nearestneighbor J2 J2’ 9 9 V-V coupling should be similar, so the limit of isolated FIG. 1. Couplings in CaV4O9. The circles represent V plaquettes is not realistic. atoms, and the lines between V’s show the couplings. The / at Self-consistent electronic structure work [8] identified numberslabel the sites used in the LSDAcalculations to de- m the V4+ spin orbital as d , which implied that it was termine the couplings. Line thicknesses are proportional to xy a larger square of V ions, the “metaplaquette,” where the best set of couplings found by fitting to the experimen- - ′ d singlet formation arises. Fitting Heisenberg Hamiltoni- talsusceptibility. Thestrongest coupling,J2,isshownasthe n thick dot-dashed lines, forming metaplaquettes, e.g., 1-6-3-8. ans to the measured dispersion of the triplet excitations o confirmed that the dominant second neighbor exchange c : coupling is crucial to account for the shape of ΩQ [2]. sumsrunovernearest-neighborbondsandthennnsums v i ThecompleteHeisenbergHamiltonianforCaV4O9has run over next-nearest-neighbor bonds. Unprimed sums X four different coupling constants: nearest-neighbor (nn) connect V’s in the same plaquette, while primed sums r andnext-nearest-neighbor(nnn)couplings and, for each connect V’s in different plaquettes. The four couplings a ofthese,intra-andinter-plaquettecouplings. Innotation are drawn in different line styles in Fig. 1. of Gelfand et al. [3], the Hamiltonian is given by In this Letter we show that the spin gap behavior of CaV O , even considering its complex structure with H =J1 Si·Sj +J1′ Si·Sj eight ve4ry9low symmetry V4+ ions in the primitive cell, nn nn′ can be calculated in ab initio fashion. Our work has X X +J S ·S +J′ S ·S , (1) three separate aspects. 1) Local spin density approxi- 2 i j 2 i j nnn nnn′ mation (LSDA) calculations are used to obtain energies X X for various magnetic configurations. The resultant ex- where S denotes the spin 1 operator in site i. The nn i 2 1 TABLE I. Magnetic configurations of the eight V ions in TABLE II. Values for the four couplings (in meV). The the primitive cell for the states used to determine the ex- LSDA values are derived from the energies in Table I. The changeconstantsfromLSDA.Mostconfigurationsaredefined SCADresults arederivedfrom thelocal orbital method,and in the text. V ions are numbered as in Fig. 1. The final col- the Fit results come from fitting the experimental suscepti- umn shows the relative LSDAenergies. bility. Both are described later in thepaper. Also shown are thecouplings deducedfrom neutron scattering data[4,5]. 1 2 3 4 5 6 7 8 ∆E/8 FM + + + + + + + + 0.0 meV Method J1 J1′ J2 J2′ FMPL + + + + – – – – -95.2 meV LSDA 8.9 1.1 6.5 23.8 FiM + + + – + – + + -70.7 meV SCAD 9.7 12.5 3.9 19.3 AFMP + – + – + – + – -130.6 meV Fit 9.3 9.6 3.7 14.2 Néel + – + – – + – + -35.4 meV Neutron 6.8 6.8 1.7 14.0 STEP + – – + – + + – -74.8 meV STEP2 – + + – – + + – -81.6 meV conditions on the four Js, and a least-squares fit gives the values listed asLSDA in Table II, eachwith a fitting change interactions are obtained by fitting these ener- uncertaintyofabout1meV.Sincebothnearestandnext gies to the mean-fieldHeisenberg modelas describedbe- nearest couplings are AF in sign, there is a great deal of low. 2) An approximate but physically motivated local frustration in the magnetic system. The large value of orbital method called the self-consistent atomic defor- J′ indicates thatsingletformationonthe metaplaquette 2 mation (SCAD) method [9] is used to provide explicit is the driving force for the spin gap. local orbitals, eigenvalues, and hopping integrals for cal- To understand how these values of the exchange pa- culatingtheexchangeinteractionsfromperturbationthe- rameters arise, we evaluate the fourth-order expressions ory. This method reveals that the nn interactions are for the exchange constants, using an approximate but not intrinsically small, but the net value of the superex- parameter-freemethodbasedontheSCAD method. For changecouplingis smalldue to cancellationsamongvar- each coupling constant in CaV O , we focus on the rel- 4 9 ious fourth-order processes. It also indicates that direct evant clusters for each coupling. The nnn interactions V-V exchange coupling is important. 3) The Heisenberg requireaV OclusterwithtwoVions(eachwithonerel- 2 Hamiltonian is solved using exact diagonalization tech- evant orbital) and one O in between. The nn exchange niques on finite periodic clusters. Spin gap behavior is interactionsrequireaV O cluster. All three 2porbitals 2 2 obtained, and χ(T) is similar to the data. The Heisen- in each O are relevant, since the low symmetry makes berg couplings are refined by fitting to χ(T). The re- them non-degenerate and oriented in directions deter- sulting Hamiltonian agrees well with χ(T) and with the mined not by symmetry but by electronic interactions. triplet dispersion determined from neutron scattering. WeneglecttheHubbardU andHund’srulecouplingon The LSDA calculations of the energy for various mag- the Oions. Inwhatfollows,U is theVon-siterepulsion, netic configurations were more precise extensions of pre- ǫ and ǫ are site energies of the V and α-th O orbitals, V α vious work on CaV4O9 [8,10]. The magnetism of the V andtiα isthe hoppingamplitude betweenthei-thVand ion is found to be robust, allowing us to break the spin the α-th O orbital. Defining the energy denominators symmetry in any manner we choose and obtain the en- ∆ =U +ǫ −ǫ simplifies the expressions. α V α ergy from a self-consistent calculation. The symmetry The initial state has each O orbital doubly filled and of the non-magnetic state is initially broken as desired each V with one electron. The perturbation theory is by applying the necessary local magnetic fields to the V givenby three fourth-order terms and the direct second- ions. The seven configurations we have chosen include order V-V term: the ferromagnetic (FM) state, one ferrimagnetic (FiM) state, and five antiferromagnetic (AF) states with zero J =j1+j2+j3+jd net spin. These AF states include the Néelstate, a state 2 4 t t (t t )2 1α 2α 1α 2α in which FM plaquettes are antialigned (FMPL), and a = +4 U ∆ ∆3 state in which the metaplaquettes are aligned antiferro- α α ! α α X X magnetically (AFMP). The configurations, given explic- t t t t 1 1 2 4t2 itly in Table I, were chosen either because of their phys- +4 1α 2α 1β 2β + + 12 (2) ∆ +∆ ∆ ∆ U ical relevance (AFMP was anticipated to be lowest in αX<β α β (cid:18) α β(cid:19) energy, as found) or computational considerations such In the nnn case, α and β sum over the three orbitals in as retaining inversion symmetry. thesingleoxygenatom. Inthenncase,αandβ sumover The resulting energies were fit to the mean-field the six orbitals in both oxygen atoms. The first three Heisenberg model, which contains simply the Sz or Ising i terms in (2) can be categorized by their configurations terms of the full Hamiltonian (1), to determine the four after the second hop of the four-hop process: 1) One couplingconstants. Thesixenergydifferencesleadtosix vanadium empty; 2) One oxygen orbital empty; 3) Two 2 oxygen orbitals half filled. The last term has an extra 8 factoroftwobecauseitarisestwice: thetotalspinsinglet case is reduced in energy and the total spin triplet is increased by the same amount. The latter picks up a 6 minus sign due to electron exchange. ) g This expression is evaluated with the SCAD model, u/ m which expresses the total density n(r) as a sum over lo- e4 calized densities |φ(αi)(r−Ri)|2 centered at the atomic −60 Experiment sites Ri [9]. The orbitals φ(αi) are solutions to atom- χ (1 LSSCDAAD centeredone-electronHamiltoniansHi foreachsite. The 2 Neutron potentialsinH aredeterminedself-consistentlyfromthe i Fit expression for the functional derivative of the total energy. It includes a localapproximationfor exchange and 0 correlation energy [11] and the Thomas-Fermi function 0 100 200 300 400 500 600 700 for kinetic energy of overlapping densities. T (K) EachVionhasthe lowestofits five3dlevelsoccupied FIG. 2. Uniform magnetic susceptibilities calculated by by a single electron, giving the V4+, O2− ionic descrip- exact diagonalization of a 20-spin cluster. The theoretical tion. U ≈3.5eVwascomputedbyminimizingtheSCAD curves using the coupling constants from Table II are shown energy subject to the constraint that one V ion has its aslines,whilethecirclesshowtheexperimentalsusceptibility charge increased by unity. The electron comes mainly ofTaniguchi,etal.[1]. Thetheoreticalfittothesusceptibility from the other V ions with only a minor portion coming is thesolid curvethat lies overthe experimental points. from the nearby O ions. The matrix elements, t = hψ |H|ψ i require the full ij i j where n is the number of V atoms per gram and N is Hamiltonian H and orthogonalized orbitals ψ. The ψ’s the number of sites in the cluster. We take g = 1.67 for are obtained from the SCAD orbitals using Lo¨wdin’s all plots. This was determined from the fit to the exper- method [12], andH is determined fromthe site centered imental magnetic susceptibility data described below. SCADHamiltoniansbyremovingthekineticenergyover- To evaluate (3), we calculate all eigenvalues of the lap contributions from the latter. This gives expressions Hamiltonian—eigenvectors are not required. We block- for H that differ in the site selected for spherical har- diagonalize the Hamiltonian with all possible symme- monic expansionofthe potential. We findthe twopossi- tries: translations, rotations, S, and S [13]. The blocks bilities, t andt ,maydifferby ∼20%,whichleadsto a z ij ji are left with no degeneracies, so the eigenvalues are muchlargeruncertaintyinthefourth-orderJ’s. Sincethe calculated very efficiently using the Lanczos algorithm vanadiumsitesofagivenpairofVionsareequivalentby with no reorthogonalization developed by Cullum and symmetry, the direct interaction, j , has no such uncer- d Willoughby [14]. This allows χ to be calculated exactly tainty. To be consistent with the direct interaction cal- at all temperatures using one Lanczos run for each sym- culation, we use the vanadium-site-expanded potentials metrysector. TheHamiltonianforthe20-spinclusterhas forevaluatingmatrixelementsbetweenoxygen-vanadium blocks as large as 36950. Within each block at least the pairs. The net values obtained (labelled SCAD in Table 400 lowest and highest eigenvalues are calculated, and II) agree rather well with those derived from LSDA en- ′ an analytic density of states is assumed for the middle ergies for J , J , and J . The close agreement may be 1 2 2 eigenvalues. This technique will be described elsewhere. fortuitous in view of the uncertainties mentioned above The susceptibility of the full Hamiltonian (1) calcu- and the approximations inherent in the SCAD method. lated with each set of coupling constants in Table II Nevertheless, we believe certain qualitative features of ′ is shown in Fig. 2. The experimental susceptibility of the SCAD results are real: 1) The values for J and J 1 1 Taniguchi, et al. [1] is shown for comparison. All curves result primarily from j , with relatively small contribu- d exhibitaspingap,asevidencedbytheirlowtemperature tions from fourth-order terms due to cancellation within ′ behavior,χ(T →0)∼exp(−∆/kT),where∆is the gap. j andbetweenj andj . 2)ThevalueforJ ,thelargest 1 2 3 2 Both the LSDA and SCAD approaches overestimate the coupling, is dominated by a single term in j , resulting 1 gap,indicatingthatthecalculatedcouplingconstantsare from V overlap with the middle O 2p level. too large. The coupling constants deduced from neutron For each set of four coupling constants, we calculated scattering are also shown [4,5]. the uniform susceptibility of the Hamiltonian (1) on pe- Also shown in Fig. 2 is a curve generated using the riodic 20-spin clusters. The susceptibility is given by: coupling constants obtained from a least-squares fit of n(gµ )2 the susceptibility to the experimental results. In the fit- χ(T)= B hSzSzi, (3) Nk T i j ting procedure, we allow the g-value in eq. (3) and all B Xij four J’s to vary. At the best fit, we obtain the cou- 3 nalization technique, which shows the Hamiltonians de- 20 termined from both ab initio approaches have quantum- ) V disorderedphases. TheHamiltonianthatbestfitstheex- e m perimentalsusceptibilityiscalculated,andtheagreement n ( 15 isremarkable. Finallythetripletdispersionoftheabini- o si tio and best susceptibility-fit Hamiltonians are shown to r e p 10 agree well with the neutron scattering data. s Di We thank Z. Weihong for the code to calculate the et Experiment curves in Fig. 3 and N.E. Bonesteel, J.L. Feldman, R.E. Tripl 5 LFSit DtoA s sucsacleepdt ibbyil it0y.58 Rudd, M. Sato, R.R.P. Singh, and C.C. Wan for stim- ulating conversations. This work was supported by the SCAD scaled by 0.65 Office of Naval Research. C.S.H was supported by the 0 (0,0) (π,0) (π,π) (0,0) National Research Council, and W.E.P. by NSF Grant DMR-9802076. Computations were done at the Arctic FIG.3. The triplet dispersion ΩQ in CaV4O9, calculated Region Supercomputing Center and at the DoD Major from the fifth-order metaplaquette series expansion of Wei- SharedResourceCentersatNAVOCEANOandCEWES. hong, Oitmaa, and Hamer [5]. The circles are the neutron scatteringdataofKodama,et al.[2]. Theab initiocouplings havebeen rescaled so their minimum gaps match the experimental minimum. The solid line shows the dispersion of the (unrescaled)couplingsdeterminedbyfittingtheexperimental susceptibility. [1] S.Taniguchiet al.,J. Phys.Soc.Japan 64,2758 (1995). [2] K. Kodama et al.,J. Phys. Soc. Japan 66, 793 (1997). pling constants listed as “Fit” in Table II and shown as [3] M. P. Gelfand et al., Phys.Rev.Lett. 77, 2794 (1996). the line thicknesses in Fig. 1. We find g = 1.67, which [4] Y. Fukumoto and A. Oguchi, J. Phys. Soc. Japan 67, is smaller than the g-value indicated by ESR measure- 2205 (1998). ments [15]. Near the minimum, the fitting function is [5] Z. Weihong, J. Oitmaa, and C. J. Hamer, Phys. Rev. B quadratic. The eigenvalues of the Hessian (scaled by 58, 14147 (1998). an arbitrary constant) are 1, 0.046, 0.013, and 0.00039. [6] N. Katoh and M. Imada, J. Phys. Soc. Japan 64, 4105 The smallness of the last eigenvalue indicates that in (1995); K. Sano and K. Takano, ibid. 65, 46 (1996); K. ′ ′ the δ{J ,J ,J ,J } = {0.09,−0.57,0.81,−0.09} direc- Ueda, H.Koutani,M. Sigrist, and P.A.Lee, Phys.Rev. 1 1 2 2 tionfromthe minimum, the least-squaresfit is verysoft. Lett. 76, 1932 (1996); M. Troyer, H. Kontani, and K. The 20-spincluster is sufficiently large comparedwith Ueda, ibid. 76, 3822 (1996); O. A. Starykh et al., ibid. the correlation length to describe the infinite system ac- 77,2558(1996): M.Albrecht,F.Mila,andD.Poilblanc, Phys. Rev. B 54, 15856 (1996); Z. Weihong et al., ibid. curately. The minimum triplet gap hardly varies be- 55, 11377 (1997); K. Takano and K. Sano, cond-mat tween 20 and 32-spin clusters: ∆ = 9.92 meV while 20 (9805153); M. A. Korotin et al.,cond-mat (9901214). ∆ =10.02 meV for the Fit Hamiltonian. 32 [7] J.-C. Bouloux and J. Galy, Acta Cryst. B 29, 1335 Fig. 3 shows the triplet dispersion Ω of the LSDA, Q (1973). SCAD, and susceptibility-fit coupling constants calcu- [8] W. E. Pickett, Phys. Rev.Lett. 79, 1746 (1997). lated with the expansion in Ref. [5]. Since the LSDA [9] M.J.Mehl,H.T.Stokes,andL.L.Boyer,J.Phys.Chem. and SCAD coupling constants overestimate the gap, we Solids 57, 1405 (1996); L. L. Boyer, H. T. Stokes, and rescaledtheirJ’sby0.58and0.65,respectively. Boththe M. J. Mehl, Ferroelectrics 194, 173 (1997); L. L. Boyer, Fit and rescaled LSDA Ω agree with the neutron scat- H. T. Stokes, and M. J. Mehl, in First Principles Cal- Q tering data reasonably well; in particular, they correctly culations for Ferroelectrics, AIP Conf. Proc. edited by have minima at Q=(0,0). R. E. Cohen (AIP,Woodbury,NY,1998), No. 436. To conclude, we have shown that the quantum- [10] An increased basis cutoff of Emax=22.6 Ry and a fixed disordered phase in CaV O can be predicted in ab ini- 16 k-point mesh were used in the linearized augmented 4 9 plane wave calculations. tio fashion. We calculated the coupling constants of [11] L.HedinandB.I.Lundqvist,J.Phys.C4,2064(1971). the Heisenberg Hamiltonian for CaV4O9 in two very [12] P. O. Löwdin, J. Chem. Phys. 18, 365 (1950). different first-principles approaches. In both methods, [13] C. Gros, Z. Phys.B 86, 359 (1992). the strongest coupling is found between next-nearest- [14] J.K.CullumandR.A.Willoughby,Lanczos Algorithms neighborVatomsonmetaplaquettes—theweakcoupling (Birkhauser, Boston, 1985). between nearest-neighbor V’s results from the cancella- [15] S.Taniguchiet al.,J. Phys.Soc.Japan 66,3660 (1997). tion among superexchange processes. The uniform magnetic susceptibility for each set of coupling constants is calculated using a novel finite-temperature exact diago- 4