Collective Spin and Charge Excitations in Planar Aromatic Molecules
K. Haghighi Mood, S. A. Jafari, E. Adibi, G. Baskaran, M. R. Abolhassani
aa r X i v : . [ c ond - m a t . s t r- e l ] O c t Collective Spin and Charge Excitations in Planar Aromatic Molecules
K. Haghighi Mood, S. A. Jafari ∗ ,
2, 3
E. Adibi, G. Baskaran, and M. R. Abolhassani Department of Physics, Science and Research Branch (IAU), Tehran, Iran Department of Physics, Sharif University of Technology, Tehran 11155-9161, Iran School of Physics, Institute for Research in Fundamental Sciences, Tehran 19395-5531, Iran Department of Physics, Isfahan University of Technology, Isfahan 84156, Iran Institute of Mathematical Sciences, Chennai 600113, India
Employing high accuracy fixed node diffusion Monte Carlo (DMC) method we calculated the lowest tripletcollective excitation (spin gap), as well as an upper bound for the singlet excitations (charge gap) in a seriesof charge neutral planar non-ladder aromatic compounds. Both excitation energies lie below the continuum ofparticle-hole excitation energies obtained from Hartree-Fock orbitals. Hence they can be interpreted as genuinebound states in the particle-hole channel. Assuming a resonating valence bond (RVB) ground state which hasbeen recently suggested for sp bonded systems [ M. Marchi, et. al. , Phys. Rev. Lett. , 086807 (2011)],offers a unified description of both excited states as two-spinon and doublon-holon bound states. We corroborateour interpretation, by Exact diagonalization study of a minimal model on finite honeycomb clusters. PACS numbers: 73.21.-b, 71.10.-w, 75.10.Kt
Introduction : Correlation effects are characteristic of pi- π conjugated systems composed essentially of hexagonal ar-rangements of sp bonds [1]. Pauling took the initiative todescribe the bonding in benzene (C H ), a prototype of thesesystems, in terms of valence bonds (VB), focusing his atten-tion on the spin part of the bonding wave function, namely asinglet. Such singlet valence bonds can be formally describedas the ground state of an effective Heisenberg exchange in-teraction J S . S , where J is the exchange integral betweenthe overlapping atomic orbitals [2]. When the coordinationnumber is low, in the above effective (model) Hamiltonian,the condition is more conducive for superposition of valencebond singlets to constitute the ground state – a unique oppor-tunity provided by three-fold coordination in these aromaticsystems. Hence Pauling’s formulation of the energy levels ofmolecules in terms of quantum mechanical superposition ofvalence bond configurations, the so called resonating valencebond (RVB) [3] becomes an important alternative rout to un-derstand energy levels of molecules.Some recent works [4–6] have combined the power ofMonte Carlo methods with the basic notion of RVB [7], toconstruct variational wave functions in terms of geminals (de-terminants composed of ”pairing” wave functions). Optimiza-tion of such RVB based many-body wave-functions, deter-mine the properties of sp bonded systems with remarkableaccuracy [4–6]. More specifically this technique captures theKekule and Dewar contributions to the ground state of ben-zene [6]. Therefore, the notion of RVB in these systems is ca-pable of capturing interesting many-body effects in the groundstate, at much lower computational cost, compared to moreinvolved quantum chemical methods. The application of theabove method to undoped graphene indicates that the groundstate is a short range, gapped spin liquid [6] which agrees with ∗ Electronic address: [email protected] other proposals based on the Hubbard model [8–10].In view of the above mentioned evidences for a possiblespin liquid ground state in planar sp bonded systems, thenext natural question would be about the nature of excitedstates. As a simple prototype which demonstrates the inad-equacy of single-particle description of energy states in thisfamily of molecules, consider benzene, C H for which theMO would predict a singlet ground state for six p z electronson the hexagonal ring. Within MO picture, the first singletexcited state ( S ) and the first triplet excited state ( T ) areexpected to be degenerate. However, observation of a remark-able splitting between the low-lying triplet and singlet excitedstates [12] indicates the importance of correlation effects evenin the excited states of these molecules. Since such correla-tion effects are based on local interactions, one expects thesame picture to hold even in extremely extended members ofthis family, such as graphene [13] and carbon nano-tubes. Theweak coupling (itinerant) limit of the graphene is well knownto represent a Dirac liquid [16], and can be described by stan-dard single-particle approach [13]. However, ab-initio calcu-lations [14] show that the strength of short-range part of theCoulomb interaction in these materials is ∼ eV, which isremarkably high and comparable to the estimated values ofthese parameters in conjugated polymers [15]. For such largevalues of Hubbard parameter U in these systems, emergenceof a non Fermi liquid state, such as spin liquid [8, 9, 11] be-comes conceivable.In this paper, we investigate the nature of low-lying ex-cited states in small molecules belonging to the family of sp -bonded carbon systems. Here we employ the state ofthe art QMC method to investigate the nature of many-bodyexcitations in such hydrocarbons. This numerically accuratemethod suggests that the lowest excitation in such moleculesis a triplet state, separated by a substantial gap from the nextsinglet excited state, for which we obtain an upper bound. Weargue that these two lowest excited states, namely, T and S can be naturally understood in terms of a picture based onspin-charge separation. This suggests that the ground statecould be viewed as a resonating valence bond state, in agree-ment with a recent proposal by Marchi and coworkers [6]. Method : Considering computational cost and accuracy,Variational Monte Carlo (VMC) and Diffusion Monte Carlo(DMC) [20] algorithms are methods of choice for the cal-culation of many body properties of medium electronic sys-tems. These QMC methods can achieve chemical accuracywith a typical computational cost ranging from the secondto fourth power of the number of particles [24]. In this pa-per we use these methods as implemented in CASINO pack-age to calculate spin and charge gap of some aromatic com-pounds. The CASINO code employs important samplingDMC method [21, 22] to project out the many-body lowestenergy state. In this method, the important sampled imagi-nary time Schrodinger equation is of the following form: f ( R , t + ∆ τ ) = Z K ( R , t + ∆ τ ; R ′ , t ) f ( R , t ) d R ′ , (1)where f ( R , t + ∆ τ ) = Ψ t ( R ) ψ ( R , t + ∆ τ ) , Ψ t ( R ) is thetrial wave function and ψ ( R , t + ∆ τ ) is system wave func-tion. The kernel K ( R , t + ∆ τ ; R ′ , t ) is the propagator. As ∆ τ approaches to infinity, ψ ( R , t + ∆ τ ) tends to ground statein any sector corresponding to a definite set of quantum num-bers. For an efficient DMC calculation we need an optimizedtrial wave function. We used the multiplication of spin up anddown Slater determinants and a Jastrow factor as a trial wavefunction: Ψ t = e j ( R ) D ↿ ( r , ..., r N ) D ⇃ ( r , ..., r N ) . (2)Here R = ( r , r , ... r N ) denotes the spatial coordinates ofall the electrons. The single-particle orbitals employed inthe above Slater determinants have been constructed fromHartree-Fock (HF) mean field solutions which serve as a ref-erence basis for ”free” particle-hole excitations. Note that thisis not the exact Jastrow-Slater trial wave function form, as it isantisymmetric only with respect to the exchange of electronswith the same spin. Such wave functions can be used to obtainexpectation values with lower computational cost for any spinindependent operators [20]. CASINO uses Jastrow factors ofthe form proposed in Refs. [22, 23]. We have taken into ac-count the electron-electron terms u , electron-nucleus terms χ centered on the nuclei and 3 body electron-electron-nucleusterms f in our calculations: j ( r i , r j ) = N − X i =1 N X j = i +1 u ( r i,j ) + N ions X I =1 N X j = i +1 χ I ( r i,I )+ N ions X I =1 N − X i =1 N X j = i +1 f I ( r i,I , r j,I , r i,j ) . (3)Optimization with respect to the parameters contained in theJastrow factor was achieved by a VMC variance minimiza-tion procedure [20]. After VMC optimization we used the sooptimized wave function as a DMC trial wave function. Op-timization of Jastrow factors without optimizing orbitals did not affect the accuracy in our calculations. However, opti-mization of Jastrow factors provides a better trial wave func-tion for DMC calculation by making it more efficient. At thefinal stage of calculation, DMC projects out the ground statefrom this trial wave function.Using the above method, we calculate the many-bodyground state in a given sector corresponding to the conservedtotal S z and total number of particles N . To extract informa-tion about spin-charge splitting from total energies, we pro-ceed as follows: Let E ( N ↑ , N ↓ ) denote the ground state en-ergy for a system where N σ is the number of electrons, eachcarrying spin ~ σ/ with σ = ± corresponding to ↑ and ↓ spinorientations, respectively. N = P σ N σ = N ↑ + N ↓ , as wellas the total spin component, S z = ( N ↑ − N ↓ ) / are constantsof motion and hence do not change the numerical projectionby DMC procedure. Therefore quantum numbers ( N ↑ , N ↓ ) appropriately label various sectors of the spectrum.Let us define the the spin gap ( ∆ s ) and charge gap ( ∆ c ) as, ∆ s ≡ E ( N ↑ + 1 , N ↓ − − E ( N ↑ , N ↓ ) , , (4) c ≡ [ E ( N ↑ +1 , N ↓ )+ E ( N ↑ − , N ↓ )] − E ( N ↑ , N ↓ ) (5)where ( N ↑ , N ↓ ) correspond to neutral system. In all com-pounds considered here, the total number of electrons, N , iseven, so that the unpolarized configuration (i.e. the state withequal number of spin up and spin down electrons, N ↑ = N ↓ )turns out to be the ground state. The energy E ( N ↑ , N ↓ ) ofthis state can be calculated as follows: We generate a trialwave function from HF method with fixed total charge (neu-tral) and one spin multiplicity. Now to calculate the spin gap,Eq. (4) we flip one of the spins from, e.g. ↓ sector, withoutaltering the total charge. The lowest energy obtained by QMCprocedure in this sector will correspond to E ( N ↑ +1 , N ↓ − .Note that in this sector, the total charge is zero and spin mul-tiplicity is three. Note that since the spin and spatial sym-metries of the many-body Hamiltonian are not broken byHF solutions, the corresponding symmetry attributes are notchanged by QMC projection. This means that the energy ofthe ( N ↑ + 1 , N ↓ − state will represent any of the three de-generate states belonging to the triplet representation of the SU (2) group. Spin gap defined above, represents the exactvalue of triplet excitation energy.Now let us discuss the physical meaning of the charge gapdefined above: Imagine an infinitely large system, with equalnumber of ↑ and ↓ spin electrons. When an electron is movedfrom one point in the system to a distant location, the resultingexcitation will be a doublon-holon pair. ∆ c is half of the aver-age energy of a pair, and hence can be interpreted as an upperbound for the energy of a single holon. To calculate the en-ergy of such doublon-holon pair, an approximate scheme is toisolate two small sub-systems surrounding the holon, and thedoublon. In the absence of interactions, the doublon-holonenergy will be given by the first term in the right hand side ofEq. (5). However, in reality there will be an attractive inter-action between them which lowers their true energy. There- FIG. 1: Benzene C H , Phenanthrene C H , Pyrene C H ,Coronene C H TABLE I: Spin and charge excitations (eV)Compound Spin Gap ∆ c C H . . C H . . C H . . C H . . C H . . fore c defined above, is an upper bound for the energy ofdoublon-holon pair with respect to the neutral background.Because of the time-reversal symmetry of the Hamiltonianemployed here, for such excitations based on charge fluctua-tions the spin orientation of the added/removed electron doesnot matter. Results:
For five planar aromatic compounds depicted inFigs. 1,2 we have calculated the above charge and spin gapswithin the all electron fixed node DMC scheme. The re-sults are reported in Table I. The spin gaps obtained here arein good agreement with experimentally reported values [12].Also the charge gap ∆ c obtain here as an upper bound for the S state agrees with existing results. For example in case ofbenzene, ∆ c = 5 . eV represents a fair upper bound for thecalculated result E S = 4 . eV [12]. Therefore the methodprescribed here to calculate the spin gap does indeed give the”lowest” excited state, and also the doublon-holon interpreta-tion employed here does represent a true upper bound for theenergy of S state. For geometry optimization as well as trialwave function generation, we used 6-311G** Gaussian basiswhich has been done by Gaussian 03 code [19]. Note thatwe performed separate optimizations for ground, and excitedstates. All required energies are obtained with an accuracybetter than ∼ meV per atom. For each of the compoundsreported in Table I, and corresponding to each set of quantumnumbers ( N ↑ , N ↓ ) , we have optimized the geometry and thetrial wave function constructed based on HF method. Thenthe Jastrow parameters have been optimized using variance FIG. 2: Coronene with additional benzene ring C H (benzo-coronene) minimization VMC. All reported DMC results are all-electroncalculations, and we did not use any pseudo-potential. DMCtime step is taken to be . Hartree − . Optimized geome-tries have been verified to ensure they do not contain imagi-nary frequencies.To interpret the data in Table I let us represent them in adifferent way. In Fig. 3 we plot charge and spin gaps ver-sus the number of carbons. The new physical interpretationcome about, when we also plot the tower of particle-hole ex-citations obtained from the (weakly correlated) Hartree-Focktheory for single particle states. This tower is the molecularanalogue of the continuum of free particle-hole pairs. As canbe seen by increasing the system size, the tower of particle-hole excitations approaches to a continuum. Moreover, thecharge and spin gaps we obtain always remain below the thecontinuum of ”free” particle-hole pairs. Therefore, they canbe interpreted as the ”bound state” of underlying free particle-hole pairs which are caused by many-body effects. First im-portant point which is suggested by this figure is that, a verylarge energy difference between lower edge of the tower, andthe many-body states found here implies they are long-livedexcitations which do not decay into the tower. Therefore,they can be associated with new quasi-particles. Note thatthe blue circles (∆ c ) is an upper bound for the true energyof a duoblon (holon), so that the true energy of the doublonstate is even lower than ∆ c . The question is, what are thesequasi-particles?Consider the lowest excited state T , which is a tripletmany-body state for all system sizes considered here. The T state can be understood in terms of a simple RPA-like boundstate formation in the triplet channel of particle-hole pairs. Ashort range repulsion of Hubbard type translates into the at-traction in the triplet particle-hole channel, and binds themtogether [10]. However the S state whose exact location inFig. 3 is somewhere between the blue circle and red square cannot be understood in terms of simple RPA-like treatments, asthe RPA in singlet channel predicts an anti-bound-state above E xc i t a t i on E n e r g y ( e V ) N(Carbon) : Free particle−hole excitaitons : Triplet : Singlet (upper bound)
FIG. 3: Charge and spin gaps versus the number of carbon atoms inaromatic compounds studied here. To generate the ”free” particle-hole continuum, we have used the Hartree-Fock orbitals. For allstudied compounds, energy of spin excitation, T is below the sin-glet (charge) excitation S , and they both are below the continuumof free particle-hole excitations. the tower of free particle-hole states [10]. Therefore the sec-ond excited state S is a genuine many-body effect, much be-yond the simple RPA like treatments. The method used hereto obtain the upper bound for the singlet charge excitationssuggests that the S can be associated to an average energy ofa doublon and a holon. To corroborate this claim further, letus use a simplified model Hamiltonian, which can capture theessence of the present QMC calculation in a more transparentway. First of all note that the minimal model which captures T state is a Hubbard model. Moreover, our earlier study ofthe particle-hole excitation spectrum in 1D chains suggeststhat the singlet collective states below the particle-hole con-tinuum are controlled by the nearest neighbor Coulomb in-teraction [18]. Therefore the minimal effective model whichcaptures both states is an extended Hubbard model, H = − t X h i,j i σ c † iσ c jσ + h.c. + U X j n j ↑ n j ↓ + V X h ı ,j i n i n j . (6)Here i, j denote sites of a 2D honeycomb lattice and h i, j i indicates that they are nearest neighbors. c † jσ creates an elec-tron in the Wannier state corresponding to the p z orbital atsite j . Here U and V denote the strength of on-site and near-est neighbor Coulomb interactions. Estimates of these param-eters based on the ab-initio methods indicates that even thescreened of these parameters in graphene are substantial [14]and on the scale of U ∼ eV, which is comparable to corre-sponding estimates for smaller aromatic molecules [15].The result of the exact diagonalization for a -site hon-eycomb lattice is shown in Fig. 4. The values of U = 9 . eV and t = 2 . eV are adopted from Ref. [14]. For the con-sidered range of V , the ground state always remains a totalsinglet state ( S ). For small values of V , the first excited stateis the T triplet, followed by a singlet excited state, S . As V increases, the singlet excited state, S comes down and ap-proaches the energy of T excited state for V ≈ eV. Beyondthis point, the first excited state will be a singlet state. Think-ing from the limit of very large molecules, the S state will E xc i t a t i on E ne r g y V (eV)
16 site cluster
TripletSinglet
FIG. 4: (Color online) The excitation energy for a -site clustercorresponding to C H in the extended Hubbard model. Values ofHubbard U and t taken from Ref. [14] are in eV. The ground statealways remains a singlet. The first excited states for small values of V are triplet. By increasing V , the order of T (red) and S (blue)is switched, and beyond V ≈ eV, S will be the first excited state.All energies in this figure are in eV. have no analogue in terms of plasmon oscillations. Since, firstof all, plasmon oscillations require non-neutral system [25]Secondly, long range Coulomb interaction makes the singletbranch either an acoustic plasmon (for 2D coulomb repulsion)or a gapped pi-plasmon (3D coulomb repulsion) branch. Sothe S state can not be interpreted as molecular analogue ofplasmon mode. On the other hand, the decreasing behavior ofthe S energy with V is consistent with a doublon-holon inter-pretation: The repulsion V among the electrons will becomeattraction − V between the doublon-holon pair, and increasing V will lower their energy. Summary and discussions : We have used ab-initio
QMCmethod to obtain an accurate excited state T and an upperbound for the T . We then used exact diagonalization to studya minimal model which captures the same set of excitations.Assuming RVB ground state [6], offers a unified understand-ing of both states. In this scenario, the T can be understood asthe energy required to break a singlet in the RVB backgroundand render it triplet [17]. Moreover, the S can be attributed toa holon-doublon pair created by removing one electron fromone carbon site, and placing it in the p z orbital of another car-bon site. Such charge fluctuations are allowed because theon-site Coulomb energy U is finite. In this picture, the de-crease in S energy by increase in V becomes quite natural.This interpretation can be a possible description of the col-lective charge excitations observed in thick multi-wall carbonnano-tubes [26]. Acknowledgements : K.H. thanks Mehdi D. Davari for use-ful discussions. S.A.J. was supported by the National EliteFoundation (NEF) of Iran. We wish to thank Dr. M. Khazaeifor assistance in the calculations, and Dr. A. Vaezi and Prof.H. Fukuyama for useful discussions. [1] D. J. Klein, S. A. Alexander, W. A. Seitz, T. G. Schmalz and G.E. Hite, Theor. Chim. Acta (1986) 393.[2] For a very readable introduction with historical notes see: P. W.Anderson, Physics Today, (2008) 8.[3] L. Pauling, The Nature of the Chemical Bond , Cornell Univer-sity Press, 3rd Ed. (1960)[4] Michele Casula, Seiji Yunoki, Claudio Attaccalite, SandroSorella, Comp. Phys. Commun. (2005) 386;[5] Todd D. Beaudet, Michele Casula, Jeongnim Kim, SandroSorella, Richard M. Martin, J. Chem. Phys. , (2008) 164711[6] M. Marchi, S. Azadi, S. Sorella, Phys. Rev. Lett. , 086807(2011)[7] P. W. Anderson, Science (1987) 1196.[8] Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, A. Muramatsu,Nature (2010) 847.[9] A. Vaezi, X. G. Wen, arXiv:1010.5744 (2010).[10] G. Baskaran, S.A. Jafari, Phys. Rev. Lett. (2002) 016402; G.Baskaran, S.A. Jafari, Phys. Rev. Lett. , (2004) 199702.[11] H. Mosadeq, F. Shahbazi, S. A. Jafari, J. Phys. Condens. Matter, (2011) 226006.[12] E. C. da Silva, J. Gerratt, D. L. Cooper, M. Raimondi, J. Chem.Phys. (1994) 3866, and references therein.[13] For a review see: A. H. Castro Neto, F. Guinea, N. M. R. Peres,K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. , 109 (2009). [14] T. O. Wehling, E. Sasioglu, C. Friedrich, A. I. Lichtenstein, M.I. Katsnelson, and S. Blu¨ugel, arXiv:1101.4007 (2011).[15] H. Kiess, et al, Eds., Conjugated conducting polymers ,Springer, 1992.[16] S. A. Jafari, Eur. Phys. Jour. B , 537 (2009).[17] Z. Noorbakhsh, F. Shahbazi, S. A. Jafari and G. Baskaran, J.Phys. Soc. Jpn. (2009) 054701.[18] M. Hafez, S. A. Jafari, Eur. Phys. Jour. B, (2010) 323.[19] We thank Dr. M. Khazaei of IMR, Tohoku university for assis-tance with this part of calculations.[20] W. M. C. Foulkes, L. Mitas, R. J. Needs and G. Rajagopal Rev.Mod. Phys. (2001) 33.[21] C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem.Phys. (1993) 2865.[22] R. J. Needs, M. D. Towler, N. D. Drummond and P. L. opezRios, CASINO Manual, Ver. 2.4.0, Univ. of Cambridge (2009).[23] N. D. Drummond, M. D. Towler and R. J. Needs, Phys. Rev. B (2004) 235119.[24] R. J. Needs, M. D. Towler, N. D. Drummond and P. L. opezRios, J. Phys.: Condens. Matter (2010) 22 023201[25] E. H. Hwang, S. Das Sarma, Phys. Rev. B (2007) 205418.[26] Kramberger, C., Hambach, R., Giorgetti, C., R¨ummeli, M. H.,Knupfer, M., Fink, J., B¨uchner, B., Reining, Lucia, Einars-son, E., Maruyama, S., Sottile, F., Hannewald, K., Olevano,V., Marinopoulos, A. G. and Pichler, T., Phys. Rev. Lett.100