A simple derivation of the exact quasiparticle theory and its extension to arbitrary initial excited eigenstates
aa r X i v : . [ c ond - m a t . o t h e r] F e b A simple derivation of the exact quasiparticle theory and its extension toarbitrary initial excited eigenstates
Kaoru Ohno, ∗ Shota Ono,
1, 2 and Tomoharu Isobe Department of Physics, Yokohama National University,79-5 Tokiwadai, Hodogaya-ku, Yokohama 240-8501, Japan Department of Electrical, Electronic and Computer Engineering,Gifu University, 1-1 Yanagido, Gifu City 501-1193, Japan (Dated: October 23, 2018)
Abstract
The quasiparticle (QP) energies, which are minus of the energies required by removing or producedby adding one electron from / to the system, corresponding to the photoemission or inverse photoemission(PE / IPE) spectra, are determined together with the QP wave functions, which are not orthonormal and evennot linearly independent but somewhat similar to the normal spin orbitals in the theory of the configurationinteraction, by self-consistently solving the QP equation coupled with the equation for the self-energy. Theelectron density, kinetic and all interaction energies can be calculated using the QP wave functions. Weprove in a simple way that the PE / IPE spectroscopy and therefore this QP theory can be applied to anarbitrary initial excited eigenstate. In this proof, we show that the energy-dependence of the self-energyis not an essential di ffi culty, and the QP picture holds exactly if there is no relaxation mechanism in thesystem. The validity of the present theory for some initial excited eigenstates is tested using the one-shotGW approximation for several atoms and molecules. ∗ [email protected] . INTRODUCTION There are variety of excited states of atoms, molecules and any kind of materials, for exam-ple, irradiated by laser-, electron-, and ion-beams, and they show variety of fascinating behaviors.However, lack of a first-principles theory, which can treat arbitrary electronic excited states ofmaterials, has strongly prevented the progress of theoretical studies of such excited states of mate-rials. In 1964, Hohenberg and Kohn [1] invented density functional theory (DFT) on the basis ofthe variational principle in quantum mechanics. Next year, adopting DFT in the local density ap-proximation (LDA), Kohn and Sham [2] proposed to use the so-called Kohn–Sham (KS) equation,which has a simple form of a single particle Schr¨odinger equation with the Hartree and exchange-correlation potentials. If one could solve the KS equation for continuous fractional occupationnumbers f µ ( f ν ) for certain occupied (empty) KS levels µ ( ν ), one might obtain the quasiparticle(QP) energies, E Nf µ = − E N − f µ = = R ε µ ( f µ ) d f µ and E N + f ν = − E Nf ν = = R ε ν ( f ν ) d f ν , which should rigor-ously correspond to the experimental photoemission and inverse photoemission (PE / IPE) spectra,thanks to Janak’s theorem, ε λ ( f λ ) = ∂ E /∂ f λ for λ = µ or ν , by integration [3]. Moreover, in 1985,Almbladh and von Barth [4] demonstrated that the correct treatment of the exchange-correlationenergy produces the KS eigenvalue for the HOMO (highest-occupied molecular orbital) level thatis identical to the minus of the experimental, i.e. , exact, ionization potential (IP). However, even ifone could perform such calculations, the resulting information is just the usual PE / IPE energy orspectra of the system at the ground state only.In 1984, in order to treat the non-steady states and the electronic excited states of materials,the time-dependent density functional theory (TDDFT) was invented by Runge and Gross [5].Although this is a very important theory that may treat arbitrary excited states of materials, thereis a serious problem that no systematic way is given to determine the exchange-correlation kernelˆ µ xc , which is a functional of a whole temporal history of the electron density { ρ ( r , t ) } and the initialwave function Ψ ( { r i } , t ) at t = t . So, one usually adopts the adiabatic LDA (A-LDA), in whichthe whole temporal history is completely ignored in the time domain and the simplest LDA isused in the space domain in the functional. The reliability of the A-LDA has been discussed bymany authors [6–12]. Recently, it has been also emphasized to use accurate exchange-correlationfunctional in the study of excited states [13].Of course some of the electronic excited states may be treated with highly accurate quan-tum chemistry (QC) methods [14, 15] such as the configuration interaction (CI), complete active2pace (CAS), Møller–Plesset (MP), and coupled cluster (CC) methods. However, these wavefunction-based methods beyond the Hartree–Fock approximation (HFA) are computationally veryheavy and again restricted to the singly, doubly, or triply (perhaps at most quadruply) excitedneutral / charged states only.On the other hand, the Green’s function method on the basis of the many-body perturbationtheory (MBPT) can directly describe the PE / IPE energy spectra, which are usually called thequasiparticle (QP) energies. The number of the papers using this method such as the GW approx-imation has increased recent years. One di ffi culty of using this method is that the MBPT and theGreen’s function are not easy to understand for the researchers who are not familiar to the fieldtheory. Therefore, the Green’s function is not at all popular in the fields of quantum chemistry andmaterials science. Moreover, there is a serious problem related to the energy-dependence of theself-energy, which has been a long standing issue and the QP picture has sometimes been consid-ered as an approximation. And again, the main issue here is that the Green’s function has beenalways defined on the neutral ground state only, while the initial electronic state in the PE / IPEspectroscopy can be any of the neutral or charged excited states.The purpose of this paper is to present a simple derivation of the QP theory, which is applicableto any of the initial electronic excited eigenstates | Ψ M γ i of the M -electron system ( M can bedi ff erent from the total number of protons N in the system) without invoking the complicated fieldtheoretical MBPT. In the course of this derivation, we will show that the energy-dependence of theself-energy is not an essential di ffi culty, and the QP picture holds rigorously if there is no relaxationmechanism such as electron-photon and electron-phonon interactions or cascade damping in thecontinuum energy bands. (Even when there is such a relaxation mechanism in the system, the QPpicture should still hold within the lifetime [16, 17], or equivalently within the mean free path.)Moreover, it becomes clear that the QP wave functions φ λ ( r , s ), which are defined as the over-lap between the initial state | Ψ M γ i operated by a creation or annihilation operator ( ψ † s ( r ) or ψ s ( r ))and the final state h Ψ M ± λ | in a PE / IPE process, satisfy the completeness condition but are gen-erally non-orthonormal and linearly dependent. (A relation to the normal spin orbitals in the CItheory [14, 18] will be mentioned in Section II.) Although the QP theory looks like a one-electronapproximation, it is not an approximation at all. It should be emphasized that the QP energies andthe QP wave functions rigorously include all the necessary many-body information and formallyexact in a sense that these quantities yield exact electron density as well as the exact kinetic, localpotential, and electron-electron interaction energies. The present theory is exact for any initial3xcited eigenstate.The rest of this paper is organized as follows. In Section II, the QP energies and the QP wavefunctions are defined in connection to the PE / IPE spectroscopy, and the equations satisfied bythem are derived. The exact expressions for the electron density and the expectation values forthe kinetic, local potential, and interaction energies are derived. In Section III, the same equationsare rederived using the Green’s functions. In Section IV, a concrete method to calculate the self-energy is derived by invoking Brillouin–Winger’s perturbation theory. In Section V, several testcalculations of isolated atoms, ions and molecules using the one-shot GW approximation demon-strate the validity of the present theory. Section VI presents several important remarks: The virialtheorem, the Ward–Takahashi identity, which is equivalent to the local charge conservation (thecontinuity equation), and the macroscopic conservation laws hold, and the Luttinger–Ward func-tional exists for the total energy, for the excited eigenstate. Finally the DFT universal functionalof the density is extended to the excited eigenstate, for which the occupation numbers of all KSlevels become its additional dependence, and a variational principle still holds in a restricted Fockspace that are orthogonal to all lower lying excited eigenstates. Finally, Section VII is devoted tosummarizing this paper. II. BASIC THEORY
The present formulation can be simply extended to infinite crystal systems, although the sum-mation over k is not written explicitly here. The applicability to crystals is obvious in the case ofthe ground state because there have been so many papers using GW approximation for crystals.How about the case of the excited states? If only some finite number of electrons are excitedfrom the ground state, nothing changes in the infinite system. Therefore, the calculations are onlymeaningful when the finite density of electrons are excited from the ground state. Such calcula-tions would correspond to strong excitations induced by e.g. laser irradiations.We start from the definition of the QPs. The initial state can be a ground or any excited stateof the N -electron neutral system or any excited state of the M ( , N )-electron charged system.We write this state as | Ψ M γ i and its total energy as E M γ . They satisfy the eigenvalue equation H | Ψ M γ i = E M γ | Ψ M γ i . The electrons are interacting each other by the Coulomb interaction, andtherefore cannot be treated independently. However, if we introduce the QP picture, we have asimple way to treat | Ψ M γ i in terms of the QPs. The idea is based on Einstein’s photoelectronic4 ff ect, which is called in modern terminology the potoelectron or PE / IPE spectroscopy. Afterthe photoemission process, the system becomes a ( M − | Ψ M − µ i , whichsatisfies the eigenvalue equation H | Ψ M − µ i = E M − µ | Ψ M − µ i . The corresponding total energy E M − µ satisfies E M − µ − E M γ = h ν − K (1a)according to the energy conservation law, where h ν is the energy of the incident photon and K isthe kinetic energy of the emitted electron. Alternatively, after the inverse photoemission process,the system becomes a ( M + H | Ψ M + ν i = E M + ν | Ψ M + ν i . Inthis case the total energies satisfy E M γ − E M + ν = h ν − K . (1b)Therefore, the total energy di ff erences are the measurable quantities. We call these energies theQP energies and define them as ε µ = E M γ − E M − µ , (2a) ε ν = E M + ν − E M γ . (2b)The QP energy ε µ corresponds to the minus of the energy required to remove one electron from the | Ψ M γ i state and add it at at the vacuum level; see Fig. 1. Similarly, the QP energy ε µ correspondsto the minus of the energy gain when one electron is added to the | Ψ M γ i state from the vacuumlevel. Note that the final states includes many ( M ± M ± ff erent from the initial eigenstate are excluded or theiramplitudes are quite small even if they are included.The QP wave functions can be defined as an overlap between the initial state | Ψ M γ i with one-electron removed or added at a point ( r , s ) ( r and s are spacial and spin coordinates of the elec-tron) and the final state | Ψ M − µ i or | Ψ M + ν i (therefore, they are sometimes called the overlapamplitudes): φ µ ( r , s ) = h Ψ M − µ | ψ s ( r ) | Ψ M γ i , (3a) φ ∗ ν ( r , s ) = h Ψ M + ν | ψ † s ( r ) | Ψ M γ i . (3b)5 acuum level QP energyspectra
FIG. 1: Meaning of the QP energies φ ( r , s ) = < M -1, μ| ( r )| M, γ>ψ μ s |M -1, μ> | M, γ> rr | M, γ>>> r xxx r (a) φ ( r , s ) = < M +1, ν| ( r )| M, γ> ψ ν * s (cid:8433) |M +1, ν> | M, γ> | M, γ>>>>> xxx r rr (b) FIG. 2: Meaning of the QP wave functions (overlap amplitudes). Yellow and green parts represent the holeand electron distributions, respectively.
Here ψ † s ( r ) and ψ s ( r ) are the (creation and annihilation) field operators. These QP wave functionsrepresent the amplitude of the one-electron di ff erence between the initial and the final states (seeFig. 2). The total electron spin density of the | Ψ M γ i state is given by ρ s ( r ) = h Ψ M γ | ψ † s ( r ) ψ s ( r ) | Ψ M γ i = occ X µ h Ψ M γ | ψ † s ( r ) | Ψ M − µ ih Ψ M − µ | ψ s ( r ) | Ψ M γ i = occ X µ (cid:12)(cid:12)(cid:12) φ µ ( r , s ) (cid:12)(cid:12)(cid:12) , (4)6here the completeness condition for the eigenstates | Ψ M − µ i of the ( M − occ X µ | Ψ M − µ ih Ψ M − µ | = H (1) = ˆ T + ˆ v = X s Z ψ † s ( r ) h (1) s ( r ) ψ s ( r ) d r , h (1) s ( r ) = − ~ m ∇ + v ( r ) (6)( v ( r ) is the local potential such as the Coulomb potential caused by nuclei) are obtained by the QPwave functions as follows: h ˆ T i = − ~ m h Ψ M γ | X s Z ψ † s ( r ) ∇ ψ s ( r ) d r | Ψ M γ i = − ~ m occ X µ X s Z φ ∗ µ ( r , s ) ∇ φ µ ( r , s ) d r , (7) h ˆ v i = h Ψ M γ | X s Z ψ † s ( r ) v ( r ) ψ s ( r ) d r | Ψ M γ i = occ X µ X s Z v ( r ) (cid:12)(cid:12)(cid:12) φ µ ( r , s ) (cid:12)(cid:12)(cid:12) d r . (8)In fact, the QP wave functions diagonalize the density matrix ρ s ( r ′ , r ) = h Ψ M γ | ψ † s ( r ′ ) ψ s ( r ) | Ψ M γ i = occ X µ φ ∗ µ ( r ′ , s ) φ µ ( r , s ) . (9)The spin density ρ s ( r ) is given by ρ s ( r , r ), and the expectation values of the kinetic and localpotential energies in the | Ψ M γ i state are given by h ˆ T i = − ~ m X s Z lim r ′ → r ∇ ρ s ( r ′ , r ) d r , (10) h ˆ v i = X s Z v ( r ) ρ s ( r , r ) d r . (11)In this sense, the QP wave functions φ µ ( r , s ) have a property similar to the natural spin orbitalsintroduced by L¨owdin more than 60 years ago [18, 19] in connection with the CI theory [14].Here, it is beneficial to note that the completeness condition occ X µ φ µ ( r , s ) φ ∗ µ ( r ′ , s ′ ) + emp X ν φ ν ( r , s ) φ ∗ ν ( r ′ , s ′ ) = all X λ φ λ ( r , s ) φ ∗ λ ( r ′ , s ′ ) = δ ( r − r ′ ) δ ss ′ . (12)7olds for the QP states. This can be easily seen by sandwiching the anticommutation relation ofthe field operators n ψ s ( r ) , ψ † s ′ ( r ′ ) o = δ ( r − r ′ ) δ ss ′ , (13)with h Ψ M γ | and | Ψ M γ i and using the other completeness condition for the eigenstates | Ψ M + ν i ofthe ( M + emp X ν | Ψ M + ν ih Ψ M + ν | = i.e. , not nor-malized to unity, and even not linearly independent [20, 21], and therefore they are di ff erent fromthe normal spin orbitals mentioned above in an exact sense. However, one should notice that theoverlap matrix S between the QP wave functions satisfies the idempotency relation S = S , i.e. , S λλ ′ = X s Z φ ∗ λ ( r , s ) φ λ ′ ( r , s ) d r = X ss ′ Z φ ∗ λ ( r , s ) δ ( r − r ′ ) δ ss ′ φ λ ′ ( r ′ , s ′ ) d r d r ′ = X ss ′ Z φ ∗ λ ( r , s ) all X λ ′′ φ λ ′′ ( r , s ) φ ∗ λ ′′ ( r ′ , s ′ ) φ λ ′ ( r ′ , s ′ ) d r d r ′ = all X λ ′′ X s Z φ ∗ λ ( r , s ) φ λ ′′ ( r , s ) d r X s ′ Z φ ∗ λ ′′ ( r ′ , s ′ ) φ λ ′ ( r ′ , s ′ ) d r ′ = all X λ ′′ S λλ ′′ S λ ′′ λ , (15)where the completeness condition (12) was used in the third equality. So, if one diagonalizes theoverlap matrix S by a unitary transformation as e S = U † S U according to the canonical orthonor-malization procedure introduced by L¨owdin [18, 19, 21], it is realized that the resulting diagonalmatrix e S satisfies e S ( e S − = S satisfies S ( S − =
0. Therefore thediagonal element of e S , i.e. , the eigenvalues of S must be either zero or unity. Thus the restrictedset of the transformed functions e φ α ( r , s ) = all X λ φ λ ( r , s ) U λα , (16)which correspond to the eigenvalue unity only, becomes linearly independent and satisfies theorthonormality condition e S αα ′ = Z e φ ∗ α ( r , s ) e φ α ′ ( r , s ) d r = δ αα ′ . (17)8his restricted set of the transformed functions satisfies the idempotency relation e S = e S . Since e S is given by restrict X α ′′ e S αα ′′ e S α ′′ α = restrict X α ′′ X s Z e φ ∗ α ( r , s ) e φ α ′′ ( r , s ) d r X s ′ Z e φ ∗ α ′′ ( r ′ , s ′ ) e φ α ′ ( r ′ , s ′ ) d r ′ = X ss ′ Z e φ ∗ α ( r , s ) restrict X α ′′ e φ α ′′ ( r , s ) e φ ∗ α ′′ ( r ′ , s ′ ) e φ α ′ ( r ′ , s ′ ) d r d r ′ , (18)which must be equal to (17), the completeness condition X λ e φ λ ( r , s ) e φ ∗ λ ( r ′ , s ′ ) = δ ( r − r ′ ) δ ss ′ (19)is fulfilled also for the transformed functions [21]. In practical calculations, however, one usuallyassumes that the diagonal elements S λλ are unity. In such cases, one has to redetermine these normsafterwards, because they are generally not unity at all. In this procedure, the unitary transformed S ′ is diagonal but its diagonal elements are no more zero or unity. It has a form [21] S ′ = U ′† S U ′ = µ
00 0 , (20)where the submatrix µ is a diagonal matrix, whose elements are not necessarily unity. Then, theorthonormal set must be chosen as e φ α ( r , s ) = √ µ αα all X λ φ λ ( r , s ) U ′ λα (21)instead of (16). From this equation, the QP wave functions with the correct norm S λλ are found tobe φ λ ( r , s ) = p S λλ restrict X α √ µ αα e φ α ( r , s ) U ′† αλ . (22)From Eq. (16), they should be equal to φ λ ( r , s ) = restrict X α e φ α ( r , s ) U † αλ . (23)Comparing these two equations, we have √ µ αα U ′† αλ p S λλ = U † αλ . (24)9ecause U αλ is unitary, restrict X α U ′ λα µ αα U ′† αλ ′ = (1 / S λλ ) δ λλ ′ (25)is obtained. Therefore, one can determine the norm of each QP wave function S λλ using thisprocedure. Apart from this complexity, the non-orthogonality and linear dependence of the QPwave functions do not bring any essential di ffi culty in the present formulation.What we want to emphasize here is that there are vast amount of irrelevant ( M ± M − M ± H is given by H = H (1) + H (2) with the one-body part (6) and the two-body part, i.e. , the electron-electron Coulomb interaction, H (2) = X ss ′ Z ψ † s ( r ) ψ † s ′ ( r ′ ) h (2) ss ′ ( r , r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) d r d r ′ , h (2) ss ′ ( r , r ′ ) = e πε | r − r ′ | , (26)the commutation relation between ψ s ( r ) and H is given by[ ψ s ( r ) , H ] = h (1) s ψ s ( r ) + X s ′ Z ψ † s ( r ′ ) h (2) s ′ s ( r ′ , r ) ψ s ′ ( r ′ ) ψ s ( r ) d r ′ . (27)Then, by sandwiching this with h Ψ M − µ | and | Ψ M γ i or with h Ψ M γ | and | Ψ M + ν i , and using Eqs. (2)and (3), we derive ε µ φ µ ( r , s ) = h (1) s ( r ) φ µ ( r , s ) + X s ′ Z h (2) s ′ s ( r ′ , r ) h Ψ M − µ | ψ † s ′ ( r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) | Ψ M γ i d r ′ , (28a) ε ν φ ν ( r , s ) = h (1) s ( r ) φ ν ( r , s ) + X s ′ Z h (2) s ′ s ( r ′ , r ) h Ψ M γ | ψ † s ′ ( r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) | Ψ M + ν i d r ′ , (28b)10ere, the second term in the right-hand side (r.h.s.) of these equations are rewritten as X s ′ Z h (2) s ′ s ( r ′ , r ) h Ψ M − µ | ψ † s ′ ( r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) | Ψ M γ i d r ′ = Z Σ s ( r , r ′ , E M γ − E M − µ ) h Ψ M − µ | ψ s ( r ′ ) | Ψ M γ i d r ′ , (29a) X s ′ Z h (2) s ′ s ( r ′ , r ) h Ψ M γ | ψ † s ′ ( r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) | Ψ M + ν i d r ′ = Z Σ s ( r , r ′ , E M + ν − E M γ ) h Ψ M γ | ψ s ( r ′ ) | Ψ M + ν i d r ′ , (29b)by introducing the energy-dependent self-energy Σ s ( r , r ′ ; ε λ ) with ε λ = ε µ or ε ν given by Eq.(2). (Here and hereafter, the self-energy includes the Hartree term, which describes the classicalCoulomb interaction.) It is clear that both identities (29a) and (29b), which define the self-energy,are consistent with each other, in particular, in the excited-state formulation. In fact, since h Ψ M − µ | and | Ψ M γ i are arbitrary excited eigenstates in Eq. (29a), they can be changed to h Ψ M γ | and | Ψ M + ν i ,and then Eq. (29b) is readily obtained. In other words, in order to have the two identities (29a)and (29b) simultaneously, one must recognize that the ordinary formulation for the self-energy forthe ground state should hold for arbitrary excited eigenstates also.From (29a) and (29b), the QP equations h (1) s ( r ) φ µ ( r , s ) + Z Σ s ( r , r ′ , ε µ ) φ µ ( r ′ , s ) d r ′ = ε µ φ µ ( r , s ) , (30a) h (1) s ( r ) φ ν ( r , s ) + Z Σ s ( r , r ′ , ε ν ) φ ν ( r ′ , s ) d r ′ = ε ν φ ν ( r , s ) (30b)are derived. Since there is the energy-dependence in the self-energy, we have to solve theseequations level to level. Note that the self-energy has a state dependence only through the depen-dence on a particular QP energy. In the case of the ground state, it has been demonstrated that theenergy dependence can be approximately treated by linearization [22–24]. However, the energy-dependence is very important in the QP theory. Indeed, multiplying both sides of Eq. (29a) by h Ψ M γ | ψ † s ( r ) | Ψ M − µ i = φ ∗ µ ( r , s ), summing up with respect to µ and s , integrating with respect to r , and using the completeness condition (5) again, we obtain E int = X ss ′ Z h (2) ss ′ ( r , r ′ ) h Ψ M γ | ψ † s ( r ) ψ † s ′ ( r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) | Ψ M γ i d r d r ′ = occ X µ X s Z φ ∗ µ ( r , s ) Σ s ( r , r ′ ; ε µ ) φ µ ( r ′ , s ) d r d r ′ (31)11or the expectation value of the electron-electron Coulomb interaction for the | Ψ M γ i state, whichis known as the Galitskii and Migdal’s formula [25] E int = −
12 Tr Z Σ s ( r , r ′ ; ω ) G s ( r ′ , r ; ω ) d r ′ (32)(Notations appearing here will be explained in Section VI.) in the QP representation [26] for thecase of the ground state. The crucial points here are that the energy dependence of the self-energyappears in the same way in both Eqs. (30a) and (31) in a general scheme of the non-orthonormalQP wave functions and that Eq. (31) holds not only for the ground state but also for arbitraryexcited eigenstates. Equation (31) implies that the information of the self-energy is enough todetermine the total energy E M γ of the | Ψ M γ i state of the system, because the QP wave functions areuniquely obtained by solving the QP equation (30). In other words, if the self-energy is known, theQP energies and QP wave functions can be calculated by solving the QP equation (30) one by one,and then the electron spin density ρ s ( r ), the kinetic energy h ˆ T i , the local potential energy h ˆ v i , andthe interaction energy E int are obtained via Eqs. (4), (7), (8), and (31), respectively. Consequentlythe total energy E M γ is obtained from the knowledge of the self-energy only. How to determinethe self-energy will be described in Section IV. In next section (Sec. III), we will show that theequations derived in Section II can be rederived by using the Green’s function. III. RELATION TO THE GREEN’S FUNCTION FORMALISM
The equations (29) and (30) can be also derived from the Dyson equation for the Green’sfunction as follows. The Green’s function G on the eigenstate | Ψ M γ i is defined as G s ( x , x ′ ) = − i h Ψ M γ | T [ ψ s ( x ) ψ † s ( x ′ )] | Ψ M γ i , (33)where the convention to write x = ( r , t ) is used; T denotes the time-ordered product, and ψ † s ( x ) = e iHt ψ † s ( r ) e − iHt and ψ s ( x ) = e iHt ψ s ( r ) e − iHt are, respectively, the creation and annihilation operatorsin the Heisenberg representation. By means of the completeness conditions (5) and (14), G s ( x , x ′ ) = i occ X µ φ µ ( x , s ) φ ∗ µ ( x ′ , s ) θ ( t ′ − t ) − i emp X ν φ ν ( x , s ) φ ∗ ν ( x ′ , s ) θ ( t − t ′ ) (34)12s readily derived. Here we put φ µ ( x , s ) = h Ψ M − µ | ψ s ( x ) | Ψ M γ i = φ µ ( r , s ) e − i ε µ t , (35a) φ ν ( x , s ) = h Ψ M γ | ψ s ( x ) | Ψ M + ν i = φ ν ( r , s ) e − i ε ν t , (35b)and φ λ ( r , s ) and ε λ ( λ = µ or ν ) are the QP wave functions and the QP energies defined in Eqs.(2) and (3).The Heisenberg equation of motion for the annihilation operator is given by i ∂∂ t ψ s ( x ) = [ ψ s ( x ) , H ] = h (1) ( r ) ψ s ( x ) + X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h ψ † s ′′ ( x ′′ ) ψ s ′ ( x ′′ ) i t ′′ = t d r ′′ ψ s ( x ) . (36)Operating both sides ψ † s ( x ′ ) from the right and the time-ordered operator T from the left, andsandwiching them with h Ψ M γ | and | Ψ M γ i , we have " i ∂∂ t − h (1) s ( r ) G s ( x , x ′ ) = δ ( x − x ′ ) + i X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | T (cid:2) ψ s ( x ) ψ s ′′ ( x ′′ ) ψ † s ′′ ( x ′′ + ) ψ † s ( x ′ ) (cid:3) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ (37)where x ′′ + with t ′′ = t means ( r ′′ , t + ) with t + = t + + . The self-energy is introduced as Z Σ s ( x , x ′′ ) G s ( x ′′ , x ′ ) dx ′′ = i X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | T (cid:2) ψ s ( x ) ψ s ′′ ( x ′′ ) ψ † s ′′ ( x ′′ + ) ψ † s ( x ′ ) (cid:3) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ , (38)and then the Dyson equation " i ∂∂ t − h (1) s ( r ) G s ( x , x ′ ) = δ ( x − x ′ ) + Z Σ s ( x , x ′′ ) G s ( x ′′ , x ′ ) dx ′′ (39)holds for any eigenstate. Equation (39) is rewritten in the Fourier space as h ω − h (1) s ( r ) i G s ( r , r ′ ; ω ) = δ ( r − r ′′ ) + Z Σ s ( r , r ′′ ; ω ) G s ( r ′′ , r ′ ; ω ) d r ′′ . (40)The Fourier transform of the Green’s function (34) with (35) is given by G s ( r , r ′ ; ω ) = all X λ φ λ ( r , s ) φ ∗ λ ( r ′ , s ) ω − ε λ − i δ λ , (41)13here δ λ = + for occupied states and δ λ = − + for empty states, i.e. , δ µ = + and δ ν = − + .Inserting this into Eq. (40) and approaching ω to one of the ε λ ’s ( λ can be either ν (emp) or µ (occ)), we obtain the QP equation h ε λ − h (1) s ( r ) i φ λ ( r , s ) − Z Σ s ( r , r ′ ; ε λ ) φ λ ( r ′ , s ) d r ′ = , (42)which is equivalent to Eqs. (30a) and (30b). If the eigenvalues are m -fold degenerate, one canchoose the corresponding m QP wave functions to be orthogonal to each other in this degeneratesubspace. Then, substituting Eq. (41) with ω = ε λ into Eq. (40), multiplying the resulting equationby the complex conjugate of one of the m QP wave functions φ ∗ λ ( r ′ , s ), and integrating with respect r ′ , we derive Eq. (42).With the aid of Eqs. (34) and (35), the left-hand side (l.h.s.) of Eq. (38) can be written as i occ X µ Z Z t ′ −∞ e − i ε µ t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ µ ( r ′′ , s ) d r ′′ φ ∗ µ ( x ′ , s ) − i emp X ν Z Z ∞ t ′ e − i ε ν t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ ν ( r ′′ , s ) d r ′′ φ ∗ ν ( x ′ , s ) , (43)whereas the r.h.s. of Eq. (38) can be written as i θ ( t ′ − t ) X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | ψ † s ( x ′ ) ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ − i θ ( t − t ′ ) X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) ψ † s ( x ′ ) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ = i θ ( t ′ − t ) occ X µ X s ′′ h Ψ M γ | ψ † s ( x ′ ) | Ψ M − µ i Z h (2) s ′′ s ( r ′′ , r ) h Ψ M − µ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ − i θ ( t − t ′ ) emp X ν X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M + ν i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ h Ψ M + ν | ψ † s ( x ′ ) | Ψ M γ i = i θ ( t ′ − t ) occ X µ X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M − µ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ φ ∗ µ ( x ′ , s ) − i θ ( t − t ′ ) emp X ν X s ′′ Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M + ν i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ φ ∗ ν ( x ′ , s ) . (44)Comparing these two equations, we find Z Z t ′ −∞ e − i ε µ t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ µ ( r ′′ , s ) d r ′′ − f µ ( x , s ; t ′ ) = Z h (2) s ′′ s ( r ′′ , r ) h Ψ M − µ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ (45a)14or t ′ > t and f ν ( x , s ; t ′ ) − Z Z ∞ t ′ e − i ε ν t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ ν ( r ′′ , s ) d r ′′ = − Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M + ν i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ (45b)for t > t ′ . Here we put occ X µ f µ ( x , s ; t ′ ) φ ∗ µ ( x ′ , s ) = emp X ν Z Z ∞ t ′ e − i ε ν t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ ν ( r ′′ , s ) d r ′′ φ ∗ ν ( x ′ , s ) (46a)for t ′ > t and emp X ν f ν ( x , s ; t ′ ) φ ∗ ν ( x ′ , s ) = occ X µ Z Z t ′ −∞ e − i ε µ t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ µ ( r ′′ , s ) d r ′′ φ ∗ µ ( x ′ , s ) (46b)for t > t ′ . Note that f µ ( x , s ; t ′ ) and f ν ( x , s ; t ′ ) in Eq. (45) are implicitly defined by Eq. (46). Sincethe r.h.s. of Eqs. (45a) and (45b) does not depend on t ′ , the l.h.s. should not depend on t ′ also.Therefore, we can put t ′ = ∞ in Eqs. (45a) and (46a), and t ′ = −∞ in Eqs. (45b) and (46b). Since f µ ( x , s ; ∞ ) = f ν ( x , s ; −∞ ) =
0, we finally obtain
Z Z ∞−∞ e − i ε µ t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ µ ( r ′′ , s ) d r ′′ = Z h (2) s ′′ s ( r ′′ , r ) h Ψ M − µ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M γ i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ (47a)and Z Z ∞−∞ e − i ε ν t ′′ Σ s ( x , x ′′ ) dt ′′ ! φ ν ( r ′′ , s ) d r ′′ = Z h (2) s ′′ s ( r ′′ , r ) h Ψ M γ | ψ † s ′′ ( x ′′ ) ψ s ′′ ( x ′′ ) ψ s ( x ) | Ψ M + ν i (cid:12)(cid:12)(cid:12) t ′′ = t d r ′′ . (47b)These equations are definitely identical to Eqs. (29a) and (29b). IV. PERTURBATION SERIES FOR THE SELF-ENERGY
Next, we derive a concrete method to calculate the self-energy Σ s . Because here we discuss theenergies at an instant time only, even when the system is slowly varying in time, we simply dropall the time dependence in the following expressions.15irst, by the eqautions ψ s ( r ) | Ψ M γ i = occ X µ | Ψ M − µ ih Ψ M − µ | ψ s ( r ) | Ψ M γ i = occ X µ | Ψ M − µ i φ µ ( r , s ) = occ X µ a µ | Ψ M γ i φ µ ( r , s ) , (48a) ψ † s ( r ) | Ψ M γ i = emp X ν | Ψ M + ν ih Ψ M + ν | ψ † s ( r ) | Ψ M γ i = emp X ν | Ψ M + ν i φ ∗ ν ( r , s ) = emp X ν a † ν | Ψ M γ i φ ∗ ν ( r , s ) , (48b)and a ν | Ψ M γ i = a † µ | Ψ M γ i = ν and an occupied QP state µ , we introduce theannihilation and creation operators, a λ and a † λ , as ψ s ( r ) = X λ φ λ ( r , s ) a λ , (49a) ψ † s ′ ( r ′ ) = X λ ′ φ ∗ λ ′ ( r ′ , s ′ ) a † λ ′ , (49b)where λ , λ ′ includes both µ and ν . Then, the anticommutation relations of the field operators arewritten as n ψ s ( r ) , ψ † s ′ ( r ′ ) o = X λλ ′ φ λ ( r , s ) φ ∗ λ ′ ( r ′ , s ′ ) n a λ , a † λ ′ o = δ ( r − r ′ ) δ ss ′ , (50a) { ψ s ( r ) , ψ s ′ ( r ′ ) } = X λλ ′ φ λ ( r , s ) φ λ ′ ( r ′ , s ′ ) { a λ , a λ ′ } = , (50b) n ψ † s ( r ) , ψ † s ′ ( r ′ ) o = X λλ ′ φ ∗ λ ( r , s ) φ ∗ λ ′ ( r ′ , s ′ ) n a † λ , a † λ ′ o = . (50c)From these equations, we can assume that the QP annihilation and creation operators satisfy theanticommutation relations n a λ , a † λ ′ o = δ λλ ′ , { a λ , a λ ′ } = , n a † λ , a † λ ′ o = . (51)Note that these relations hold even when the QP states λ, λ ′ are not orthogonal to each other,because Eq. (50) with Eq. (51) satisfies the completeness condition (12) only.Second, we introduce a hypothetical non-interacting system, where the QPs are not interactingeach other, although in real systems the QPs generally interact with each other at least at shortdistances. This hypothetical system is described by the Hamiltonian H ′ = H (1) + H (2)self (52)16here H (2)self is given by (26) whose operation is restricted to construct the self-energy Σ s only andine ff ective for any other interaction. Its eigenvalue equation is given by H ′ | Φ M + n α i = E ′ M + n α | Φ M + n α i (53)where M + n indicates the total number of electrons in this system. When n = α = γ , weassume that | Φ M γ i has the same energy eigenvalue E ′ M γ = E M γ as the true system. All the otherstates of the hypothetical system are constructed by operating several a λ ’s and a † λ ’s on | Φ M γ i . Weassume H ′ as the unperturbed Hamiltonian and H ′′ = H (2) − H (2)self (54)as a perturbation, and adopt Brillouin–Wigner’s perturbation theory, which is well known to beapplicable not only to the ground state but also to any excited eigenstate of the full Hamiltonian.Multiplying both sides of (cid:16) E M γ − H ′ (cid:17) | Ψ M γ i = H ′′ | Ψ M γ i , (55)by the projection operator P = − | Φ M γ ih Φ M γ | , (56)and noting that H ′ commutes with P , we have (cid:16) E M γ − H ′ (cid:17) P | Ψ M γ i = PH ′′ | Ψ M γ i . (57)Equation (57) can be rewritten as | Ψ M γ i = | Φ M γ ih Φ M γ | Ψ M γ i + E M γ − H ′ PH ′′ | Ψ M γ i . (58)Then, by iteration, we have | Ψ M γ ih Φ M γ | Ψ M γ i = | Φ M γ i + E M γ − H ′ PH ′′ | Φ M γ i + E M γ − H ′ PH ′′ E M γ − H ′ PH ′′ | Φ M γ i + · · · . (59)So, the expectation value of A ( r , r ′ ) = ψ † s ( r ) ψ † s ′ ( r ′ ) ψ s ′ ( r ′ ) ψ s ( r ) can be evaluated as h A ( r , r ′ ) i = h Φ M γ | A ( r , r ′ ) | Φ M γ i + h Φ M γ | A ( r , r ′ ) 1 E M γ − H ′ PH (2) | Φ M γ i + h Φ M γ | A ( r , r ′ ) 1 E M γ − H ′ PH (2) E M γ − H ′ PH (2) | Φ M γ i + · · · , (60)17here we replaced H ′′ with H (2) without restriction because it is obvious that these interactionsare not related to the construction of the self-energy but simply related to the construction ofthe expectation value of A ( r , r ′ ). Then, multiplying the both sides of Eq. (60) by h (2) ss ′ ( r , r ′ ),summing over s , s ′ , and integrating with respect to r and r ′ , we have E int = h Φ M γ | H (2) | Φ M γ i + h Φ M γ | H (2) E M γ − H ′ PH (2) | Φ M γ i + h Φ M γ | H (2) E M γ − H ′ PH (2) E M γ − H ′ PH (2) | Φ M γ i + · · · . (61)In the k -th order term in this perturbation series, there are 4( k +
1) field operators sandwiched by h Φ M γ | and | Φ M γ i , which can be expanded as (49a) and (49b). For example, the first term in ther.h.s. of Eq. (61) is expressed as12 X ss ′ occ X µ µ µ µ Z h (2) ss ′ ( r , r ′ ) φ ∗ µ ( r , s ) φ ∗ µ ( r ′ , s ′ ) φ µ ( r ′ , s ′ ) φ µ ( r , s ) d r d r ′ h Φ M γ | a † µ a † µ a µ a µ | Φ M γ i . (62)In this matrix element, the QP states µ and µ are removed from and µ and µ are added to the M -particle state | Φ M γ i , and after these operations, the M -particle state must be restored to theoriginal state | Φ M γ i ; otherwise this matrix element becomes zero. That is, either ( µ = µ and µ = µ ) or ( µ = µ and µ = µ ) should be satisfied. Of course, this first order approximationrigorously corresponds to the Hartree–Fock energy in the extended HFA [27]. In fact, Eq. (62)becomes E HFint = X ss ′ occ X µµ ′ Z h (2) ss ′ ( r , r ′ ) (cid:12)(cid:12)(cid:12) φ µ ( r , s ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) φ µ ′ ( r ′ , s ′ ) (cid:12)(cid:12)(cid:12) d r d r ′ (63) − X s occ X µµ ′ Z h (2) ss ( r , r ′ ) φ ∗ µ ( r , s ) φ ∗ µ ′ ( r ′ , s ) φ µ ( r ′ , s ) φ µ ′ ( r , s ) d r d r ′ . (64)Its analogy holds for all orders. That is, all the intermediate QP states λ , λ , λ , λ , · · · must bepairwise ( λ i can be either occupied states µ i or empty states ν i ); each pair consists of one annihila-tion operator a λ i and one creation operator a † λ i of the same QP state λ i . All possible combinationsof the product of the pairs have to be summed up to evaluate the total contribution at this order.Because of the existence of the projection operator P in (61), the intermediate states cannot be theinitial state | Φ M γ i . Moreover, the energy denominator in (61) is just a simple addition or subtrac-tion of the QP energies ε λ ’s. If we draw a pair of the QP state by a thick sold line and H (2) bya dotted line, (61) is expressed by the diagrams shown in Fig. 3. This expression is exactly the18ame as the interaction energy expressed by the skeleton diagrams, i.e., the diagrams expressed bythe clothed Green’s function (a thick solid line in Fig. 3), in the usual MBPT [28–31]. They alsocorrespond to the same order contribution in the Møller–Plesset theory if the initial excited-stateconfiguration and the full self-consistency are imposed in the latter. Thus we derived exactly thesame formula only by using the true eigenstates of the full Hamiltonian and their derivatives suchas the QP energies and QP wave functions. An important point here is that we have never assumedthe orthogonality of the QP wave functions in these equations. E = + + + + (cid:567)(cid:567)(cid:567) int r r
1 21 2 r ’ r ’ λ s λ ’ λ s λ ’
1 1 1
2 2 FIG. 3: Diagrammatic representation of the interaction energy E int . Third, in order to obtain the expression for the self-energy Σ s , we compare the expression (61)with the Galitskii and Migdal’s formula of the interaction energy (31). Looking at Eq. (31), wefind that if we drop φ µ ( r , s ) φ ∗ µ ( r ′ , s ) from (61), i.e. , a thick solid line from Fig. 3, we should havethe self-energy Σ s ( r , r ′ ; ε µ ) equivalent to Fig. 4. In Fig. 4, we drew two external thin solid lines toindicate the dependence of the self-energy Σ s on the two positions ( r , s ) and ( r ′ , s ). This is againexactly the same skeleton diagrams for the self-energy in the MBPT [30, 31]. When we removea pair of φ µ ( r , s ) φ ∗ µ ( r ′ , s ) from the product of H (2) ’s in the expression of E int (61), we also have toremove the summation with respect to this special QP state, µ . Therefore, the corresponding QPenergy ε µ appearing in some parts of the energy denominators ( E M γ − H ′ ) − remains unsummedwith the special index µ after taking the summations with respect to all the intermediate QP states,resulting in the ε µ dependence in the final expression. Consequently, the resulting self-energy musthave the explicit ε µ dependence, where µ corresponds to the QP state of the external thin line inFig. 4.Last, we make a brief comment on the relationship between the present formulation using theresolvent operator ( E M γ − H ′ ) − and the formulation using the clothed Green’s function (41). Tosee this relationship, we insert the completeness condition X α | Φ M α ih Φ M α | = H (2) ’s in Eq. (61). Then H (2) is replaced by the sum over all possible inter-19 + + + + (cid:567)(cid:567)(cid:567) Σ FIG. 4: Diagrammatic representation of the self-energy Σ s . Two external thin solid lines are additionallydrawn to indicate the dependence of the self-energy on the two positions ( r , s ) and ( r ′ , s ). mediate QP states λ , λ , λ , λ of the product of four QP wave functions and the matrix element h Φ M λ | · · · | Φ M λ i replaced with h Φ M α | · · · | Φ M β i in a form (62). (Note that the number of electronsremains M in the intermediates states α and β , because H (2) has two annihilation and two creationoperators.) Then, all the QP states λ i are made pairwise, each of which has a form φ λ ( r , s ) φ ∗ λ ( r ′ , s )identical to the numerator of the Green’s function (41). Simultaneously, each resolvent operator isreplaced with its expectation value ( E M γ − E ′ M α ) − sandwiched by an intermediate state | Φ M α i , whichcorresponds to the denominator of the result of the ω -integration of the even-number product of theGreen’s functions. As an example, let us consider the third term of Fig. 3. It consists of two bub-ble diagrams, each of which corresponds to the polarization function P ( r i , r ′ i ; ω i ) ( i = ,
2) in therandom-phase approximation (RPA), connected by two dotted lines, h (2) s s ( r , r ) and h (2) s s ( r ′ , r ′ ).The polarization function is evaluated as P ( r , r ′ ; ω ) = − i X s Z ∞−∞ G s ( r , r ′ ; ω + ω ′ ) G s ( r ′ , r ; ω ′ ) d ω ′ π = X λλ ′ φ ∗ λ ( r , s ) φ λ ′ ( r , s ) φ ∗ λ ′ ( r ′ , s ) φ λ ( r ′ , s ) ω − ε λ + ε λ ′ − i δ λ [ f ( ε λ ′ ) − f ( ε λ )] (66)as usual, where f ( ε ) is the Fermi distribution function. Multiplying P ( r , r ′ ; ω ) by P ( r ′ , r ; ω ) = P ( r , r ′ ; − ω ) and i / π , and integrating it with respect to ω from −∞ to ∞ , we obtain a possibleterm, h Φ M γ | a † λ a † λ a λ ′ a λ ′ | Φ M α ih Φ M α | a † λ ′ a † λ ′ a λ a λ | Φ M γ i E M γ − E ′ M α = δ occ λ δ emp λ ′ δ occ λ δ emp λ ′ ε λ − ε λ ′ + ε λ − ε λ ′ (67)whose denominator is identical to the expectation values of the resolvent operator ( E M γ − H ′ ) − inthe second term of Eq. (61). Final result can be obtained by operating14 X λ X λ ′ X λ X λ ′ X s X s Z d r d r d r ′ d r ′ h (2) s s ( r , r ) h (2) s s ( r ′ , r ′ ) × φ ∗ λ ( r , s ) φ ∗ λ ( r , s ) φ λ ′ ( r , s ) φ λ ′ ( r , s ) × φ ∗ λ ′ ( r ′ , s ) φ ∗ λ ′ ( r ′ , s ) φ λ ( r ′ , s ) φ λ ( r ′ , s ) . (68)20he eight product of the QP wave functions comes from the expansions (49a) and (49b). All theseoperations correspond to taking the trace, Tr[ Ph (2) Ph (2) ], in the Green’s function method [23]. V. TEST CALCULATIONS
In order to verify the applicability of the QP theory for some initial excited eigenstates, letus first consider a hydrogen atom with N = M = M = M + = ffi nity (EA) of its cation A + or M + , providedthat the atomic geometries of M and M + are the same (This is not a bad approximation for typicalmolecules). Similarly, there are several other identities between the electron removal and absorp-tion energies of cationic, anionic and photoabsorbed states of A or M. These identities allow usto check the validity of the present theory. Here, we perform the perturbation expansion startingfrom the cationic, anionic, or photoabsorbed state as well as the neutral ground state of A or M tocalculate the EQP energies, which should coincide with each other.We use the standard one-shot GW approximation ( G W ) by beginning with the LDA [32] forboth the ground-state and excited-state configurations of isolated lithium, aluminum, and beryl-lium atoms and nitrogen, oxygen and lithium (diatomic) molecules with the bond lengths fixed at1.106 ˚A, 1.216 ˚A, and 2.723 ˚A, respectively, for N , O , and Li (calculated with B3LYP / , and Al, and that of 12 ˚A is used for Be, N , and O .Cuto ff energies for PWs, the Fock exchange and the correlation part of the self-energy (as well21 ABLE I: G W QP energies of neutral (A), cationic (A + ), anionic (A − ), and photoabsorbed (A ∗ ) statesof lithium, aluminum, and beryllium atoms in units of eV. The right atom / ion has one more electron thanthe left atom / ion; and the left and the right have N ( −
1) and N ( +
1) electron configurations, respectively( S + L states are also listed); In the same row, the QP energies of the forward, → , electron absorptionprocess and the reverse, ← , electron release process should coincide with each other. The values insidethe parentheses indicate the usual G W QP energy calculated for the neutral atom at the ground state. ECin the superscript in the Reference values means the value calculated by the EOM coupled cluster singlesand doubles (EOM-CCSD) with aug-cc-pV5Z basis set using Gaussian 09 [35]. All the other referencevalues a − g are experimental data [36–41].Atom N ( −
1) Electron S + L Atom N ( +
1) Electron S + L QP energies (eV) / Ion Configuration State / Ion Configuration State → ←
ReferenceLi + (1 s ) S ↔ Li (1 s ) (2 s ) S 5.1 (5.6) 5.39 a Li + (1 s ) S ↔ Li ∗ (1 s ) (2 p ) P 3.0 3.2 3.54 b Li + ∗ (1 s )(2 s ) S ↔ Li (1 s ) (2 s ) S 63.9 (63.4) 64.46 b Li + ∗ (1 s )(2 p ) P ↔ Li ∗ (1 s ) (2 p ) P 65.0 66.0 65.76 b Li (1 s ) (2 s ) S ↔ Li − (1 s ) (2 s ) S (0.6) 0.3 0.61 c Al + (3 s ) S ↔ Al (3 s ) (3 p ) S 5.7 (5.6) 5.98 d Al + (3 s ) S ↔ Al ∗ (3 s ) (4 s ) P 2.6 3.1 2.84 e Al + ∗ (3 s )(3 p ) S ↔ Al (3 s ) (3 p ) P 12.6 (10.2) 13.40 e Al (3 s ) (3 p ) P ↔ Al − (3 s ) (3 p ) S (0.3) - 0.43 f Al (3 s ) (3 p ) P ↔ Al −∗ (3 s ) (3 p )(4 s ) P (-0.1) 0.5 0.2 EC Al ∗ (3 s ) (4 s ) S ↔ Al −∗ (3 s ) (3 p )(4 s ) P 3.5 2.7 3.3 EC Be + (2 s ) S ↔ Be (2 s ) S 9.3 (9.3) 9.32 f Be + (2 s ) S ↔ Be ∗ (2 s )(2 p ) P 4.2 5.4 4.40 g Be + (2 s ) S ↔ Be ∗ (2 s )(2 p ) P 6.4 6.2 6.59 g a Reference [36]. b References [36, 37]. c References [38]. d Reference [39]. e References [39, 40]. f Reference [41]. h References [37, 41]. ABLE II: G W QP energies of neutral (M), cationic (M + ), and photoabsorbed (M ∗ ) states of nitrogen,oxygen, and lithium molecules in units of eV. The right molecule has one more electron than the leftmolecule; and the left and the right have N − N electron configurations, respectively (electronicstates are also listed); In the same row, the QP energies of the forward, → , electron absorption process andthe reverse, ← , electron release process should coincide with each other and also with the reference values[42–45]. The values inside the parentheses indicate the usual G W QP energy calculated for the neutralatom at the ground state.Molecule State Molecule State QP energies (eV) → ←
ReferenceN + Σ + g ↔ N Σ + g a N + Σ + g ↔ N ∗ Π g b O + Π g ↔ O Σ − g c O + Π g ↔ O ∗ Π g d Li + Σ + g ↔ Li Σ − g e Li + Σ + g ↔ Li ∗ Σ − u f a Reference [42]. b References [42, 43]. c Reference [44]. d References [43, 44]. e Reference [45]. f References [43, 45]. as the polarization part) are set as 11.05 Ry, 44.22 Ry, 11.05 Ry for Li, Li , and Al, 19.65 Ry,78.61 Ry, 19.65 Ry for Be and N , and 27.71 Ry, 122.83 Ry, 30.71 Ry for O , respectively. In thesummation in the correlation part, 400 levels are included for Li, Al, and Be; 3000, 6000, and 7000levels are included, respectively, for Li , N , and O . We use the generalized plasmon pole model[32] to avoid the ω integration. The energy dependence of the self-energy is dealt by introducingthe Z factor (linearization around the Kohn-Sham eigenvalues) as usual.The computational time scales as O ( N emp N occ N G ) where N emp , N occ , and N G are the numbers ofempty states, occupied states, and G vectors needed for the non-local Fourier components of thepolarization function and the correlation part of the self-energy.There is no essential di ff erence between the QP theories using the neutral ground state and23 IG. 5: QP energy spectra of Li + ∗ calculated by G W in units of eV. Experimental data are taken fromRefs. [36, 37, 46, 47]. Arrows represent electrons (spin) occupying the level. the other (neutral or charged) excited states as the initial state, and the GW calculations for thelatter excited-state case can be simply performed by just changing the order of levels after thediagonalization. However, in order to compare the results, we will call the former neutral ground-state case the ordinary QP (OQP) theory and the latter excited-state case the extended QP (EQP)theory.The calculated results of the QP energies of several QP states are listed in Table I for Li, Al,and Be, and in Table II for N , O , and Li , together with the reference values. EC in Table Imeans the value calculated by the EOM coupled cluster singles and doubles (EOM-CCSD) withaug-cc-pV5Z basis set using Gaussian 09 [35]. References [36–45] are all experimental data. Wefind good agreements between the EQP energies of the forward ( → ) electron absorption processand the backward ( ← ) electron release process. Maximum deviation from the reference valueis 1.0 eV in the case Be + ← Be ∗ for atoms, and 0.6 eV in the case Li + ← Li ∗ for molecules,which are smaller than the deviations in the OQP theory (values inside parentheses in Tables I andII); 3.2 eV in the case Al + ∗ ← Al. The average deviations are 0.4 eV and 0.3 eV, respectively,for atoms and molecules; they are the same orders for both the OQP and EQP theories (valuesoutside parentheses in Tables I and II). This means that there is no specific deviation in the EQPtheory, demonstrating the validity of the present theory. To obtain better agreement, it would benecessary to use the self-consistent GW( Γ ) methods [23, 24]. The resulting QP energy spectra ofLi + ∗ calculated by G W are shown in Fig. 5. It is seen from this figure that the agreement between24he theory and experiment is very good already in the present approximation. VI. SEVERAL IMPORTANT REMARKS
Before ending this paper, we will give several important remarks.(1) First of all, it is easy to see that the virial theorem holds for any excited eigenstate of the M -electron Hamiltonian H M { r i , R k } = − ~ m M X i ∇ i + e πε − M X i X k Z k | r i − R k | + M X i > j | r i − r j | + X k > l Z k Z l | R k − R l | , (69)where { r i } and { R k } denote the positions of electrons and nuclei. In fact, since the M -electron wavefunction Ψ M γ { r i } of this excited eigenstate γ is normalized as R Ψ M ∗ γ { α r i } Ψ M γ { α r i } Q Mi d ( α r i ) = α , the total energy of this excited eigenstate at the equilibriumatomic positions { R k } , E M γ = Z Ψ ∗ γ { α r i } H M { α r i , α R k } Ψ γ { α r i } M Y i d ( α r i ) , (70)has the α -dependence only in the form of h ˆ T i /α + h ˆ V total i /α , where H M { r i , R k } is given by (69), h ˆ T i and h ˆ V total i denote the kinetic and potential contributions to E M γ . If the atomic configurationis fully relaxed at this excited state, the total energy must be stationary with respect to the changein α . Therefore, di ff erentiating E M γ with respect to α and setting α at unity, the virial theorem2 h ˆ T i + h ˆ V total i = H (2) =
0) system, exchange therole of creation and annihilation operators for occupied and empty states of the excited eigenstate,and introduce the normal product as usual. When the interaction H (2) is gradually switched on,the QP levels may apparently change the order. However, this does not matter at all, because,in the case of an excited state, it is not necessary to require that the noninteracting ground statemoves onto the interacting ground state without level crossing. Even if the order of occupied andempty levels changes during the switch on process, we can keep track of each QP state; see Fig. 6.In the noninteracting system, the QP states refer to the one-particle states. There is a one-to-onecorrespondence between the one-particle states of the noninteracting system and the QP states of25he interacting systems, although the orders of these states may be di ff erent. There is no constrainton the orders of the one-particle states and the QP states, and they can be di ff erent to each other.Therefore, we should not use the standard words “the adiabatic switch on process” for this processbecause the word “adiabatic” involves the meaning of no level crossing. Instead, it would be betterto use the words “the gradual switch on process” to avoid confusion. The gradual switch on of theinteraction from time t = −∞ to t = t = −∞ ) intothe interacting excited state (at t = FIG. 6: (a) Adiabatic switch on process for the ground state: Occupied QP levels must be below empty QPlevels, i.e. , level crossing is not allowed between the occupied and empty QP levels. (b) Gradual switch onprocess for an excited state: QP levels can change the order, i.e. , occupied and empty QP levels can intersecteach other. (3) Third, according to the Klein’s derivation (see, for example, Section 5-6 in Ref. [31]), thetotal energy can be alternatively written as E M γ [ G ] − E M γ [ G (0) ] = Φ [ G ] + Tr " G s G (0) s − − ln G s G (0) s , (71)where G (0) s is the Green’s function of the non-interacting system, i.e. , the system described by theone-body part of the Hamiltonian H (1) only, Tr = i P s R d r R ∞−∞ e i ηω d ω/ π with η = + , and Φ [ G ]is the the Luttinger–Ward functional [23, 31, 48] Φ [ G ] = − ∞ X k = k Tr Z d r ′ Σ ( k ) s ( r , r ′ ; ω ) G s ( r ′ , r ; ω ) . (72)Here Σ ( k ) s ( r , r ′ ; ω ) represents the k -th order skeleton diagram of the self-energy. In Eq. (71),Tr[ G s / G (0) s −
1] is no other than Tr[ Σ s G s ] = − E int (see Galitsukii–Migdal’s formula (32) and the26yson equation (40) with G (0) − s = ω − h (1) s ), and − Tr ln G s / G (0) s can be evaluated as P occ µ ε µ − E M γ [ G (0) ] and E M γ [ G (0) ] = P occ µ ε (0) µ with the eigenvalue of the non-interacting system ε (0) µ . We havealso the variational principles δ Φ [ G ] /δ G s = − Σ s and δ E M γ [ G ] /δ G s = δ ( x − z ) G − s ( z , y ) − G − s ( x , z ) δ ( z − y ) = i X α = ∂∂ z α Γ ( α ) s ( x , y , z ) , (73)which is equivalent to the continuity equation for the electron density (gauge invariance), holdsalso for the excited eigenstate. Here Γ ( α ) s ( x , y , z ) are the scalar ( α =
0) and vector ( α = , , h Φ M γ | T [ j α ( z ) ψ s ( x ) ψ † s ( y )] | Φ M γ i = − Z G s ( x , x ′ ) Γ ( α ) s ( x ′ , y ′ , z ) G s ( y ′ , y ) dx ′ dy ′ , (74)where the four-current density j α ( x ) is defined as j ( x ) = ρ ( x ) = X s ψ † s ( x ) ψ s ( x ) , (75) j ( x ) = − ρ ( x ) A ( x ) − i ~ m X s n ψ † s ( x ) ∇ ψ s ( x ) − [ ∇ ψ † s ( x )] ψ s ( x ) o (76)with the vector potential A ( x ).(5) For the discussion of transport properties, it is important to guarantee the macroscopicconservation laws of the total energy, the total number of electrons, the total momentum, and soon. These macroscopic conservation laws have been rigorously discussed by Baym and Kadano ff [52, 53] for the ground state. Their argument may be generalized to the case of an arbitraryexcited eigenstate, and the macroscopic conservation laws hold for arbitrary excited states, whenthe necessary symmetries suggested by them are satisfied for the excited-state Green’s functions.In general, however, we expect that the total number of electrons is conserved in our excited-statetheory because the local charge conservation is guaranteed by the Ward–Takahashi identity (73).Moreover, we expect that the total energy and the total momentum of the electronically excitedmany-atom systems should be conserved, if the system is isolated. This is simply due to thetime-reversal symmetry and the translational symmetry in space of the total Hamiltonian of themany-atom system.(6) However, the fluctuation-discipation theorem is not guaranteed for arbitrary excited statesbecause it is based on the canonical ensemble, i.e. the thermal equilibrium, at a given temperature.Therefore, the linear response theory is not developed here and left for the future study.277) Last, we discuss how the DFT is modified for the excited state. Similar to the ordinary DFT,we introduce a universal function F defined as F = h Φ M γ | H eg | Φ M γ i , (77)where H eg is the Hamiltonian of the interacting electron gas. By definition, F is a unique functionalof the eigenstate | Φ M γ i of the full Hamiltonian H . This statement can be also straightforwardlydeduced from the Luttinger–Ward theory, because F is expressed as a unique functional of theGreen’s function, which is by definition (see Eq. (33)) a unique functional of the eigenstate | Φ M γ i .(Note that the rest contribution to the total energy is just the external potential v ( r ) coupled tothe electron density ρ ( r ) = P s ρ s ( r ).) Below we show that this eigenstate | Φ M γ i and therefore F are unique functionals of the electron density ρ ( r ) and the occupation numbers of all KS levels.For this purpose, we use the reductio ad absurdum along the same line as the Hohenberg–Kohn’stheory [1]. Let us assume that two di ff erent local potentials v ( r ) and v ′ ( r ) give the same electrondensity ρ ( r ), and show that this assumption leads to a contradiction. The Hamiltonians withthe local potentials v ( r ) and v ′ ( r ) are written as H and H ′ , respectively. Their eigenvalues andeigenstates are written as E M γ , E ′ M γ , | Ψ M γ i , and | Ψ ′ M γ i , which satisfy the eignvalue equations H | Ψ M γ i = E M γ | Ψ M γ i and H ′ | Ψ ′ M γ i = E ′ M γ | Ψ ′ M γ i . The excited eigenstate is orthogonal to alllower excited eigenstates as well as the ground state. Therefore, we can still use the variationalprinciple that E M γ (or E ′ M γ ) become minimum when | Ψ M γ i (or | Ψ ′ M γ i ) is the true excited eigenstateinside the Fock space orthogonal to all lower excited eigenstates and the ground state. Then, insidethis Fock space, we expect [12] E ′ M γ = h Ψ ′ M γ | H ′ | Ψ ′ M γ i < h Ψ M γ | H ′ | Ψ M γ i = h Ψ M γ | H | Ψ M γ i + Z [ v ′ ( r ) − v ( r )] ρ ( r ) d r = E M γ + Z [ v ′ ( r ) − v ( r )] ρ ( r ) d r (78a)and E M γ = h Ψ M γ | H | Ψ <γ i < h Ψ ′ M γ | H | Ψ ′ M γ i = h Ψ ′ M γ | H ′ | Ψ ′ M γ i + Z [ v ( r ) − v ′ ( r )] ρ ( r ) d r = E ′ M γ + Z [ v ( r ) − v ′ ( r )] ρ ( r ) d r . (78b)28dding these two inequalities, we obtain E M γ + E ′ M γ < E M γ + E ′ M γ . (79)Thus we find that the initial assumption was wrong and that v ( r ) is a unique functional of ρ ( r )when the corresponding state is constructed as the γ -th excited eignstate. Since v ( r ) determines theform of the Hamiltonian H uniquely, its eigenstate | Ψ M γ i itself is a unique functional of the density ρ ( r ) and the information that this is the excited level γ . Since the latter information is equivalentto the knowledge of the occupation numbers of all KS states, we thus conclude that | Ψ M γ i andtherefore F are unique functionals of the electron density ρ ( r ) and the occupation numbers ofall KS states only. The variational principles still holds within the Fock space orthogonal to alllower excited states and the ground state, which can be controlled by the choice of the occupationnumbers of KS states. Therefore only additional requirement in the extension of the ordinary DFTto the electronic excited state is that we have to impose the constraint on the occupation numbers ofall KS states. In other words, the ordinary DFT still holds for the excited state if we additionallygive the information on the occupation numbers of all KS states. This part only influences theexchange-correlation potential in the DFT. We expect that there is a close similarity between theDFT for excited states and the QP theory for initial excited eigenstates. This argument encouragesa direct use of the DFT to an electronic excited state such as in [54], if a reliable exchange-correlation functional is used. VII. SUMMARY
In this paper, we have shown that the quasiparticle (QP) picture, which correspond to thephotoemission or inverse potoemission (PE / IPE) spectroscopy, holds exactly for arbitrary excitedeigenstate | Ψ M γ i of the system. The QP energy spectra correspond to the PE / IPE energy spectra,and the QP wave functions yields the electron density and the expectation values of the kinetic andlocal potential energies. They are the solution of the QP equations coupled with the equation forthe energy-dependent self-energy. The interaction energy E int of | Ψ M γ i eigenstate is expressed bythe energy-dependent self-energy and the QP wave functions in Galitski–Migdal’s form. All theformulation is given either with or without introducing the Green’s function. We have tried to drivethe closed set of equations as simple as possible to be readable for general readers who are notfamiliar to the field theoretical MBPT. In this derivation, we showed that the energy-dependence29f the self-energy is not an essential di ffi culty. Some simple calculations of isolated atoms andmolecules using the one-shot GW approximation demonstrate the validity of the present theory.We also showed that the virial theorm, the Ward–Takahashi identity, and the conservation lawshold and the Luttinger–Ward functional exists for arbitrary excited eigenstates. Last, we gavea comment on the extension of density functional theory (DFT) to arbitrary excited eigenstate.We hope that the present theory makes a breakthrough toward the first-principles calculations ofarbitrary electronic excited states of materials. Acknowledgments
We thank Steven G. Louie for helpful discussions. This work has been supported by the grant-in-aid for Scientific Research on Innovative Areas (Grant No. 25104713) from MEXT and thegrant-in-aid for Scientific Research B (Grant No. 25289218) from JSPS, and also by the highperformance computer infrastructure (HPCI) (Project IDs. hp150231, hp150259, hp160234) andthe post ”K” projects (priority issue 7 and germinating issue 1) both promoted by MEXT. [1] P. Hohenberg, and W. Kohn, Phys. Rev. , B864 (1964).[2] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[3] J. F. Janak, Phys. Rev. B , 7165 (1978).[4] C.-O. Almbladh and U. von Barth, Phys. Rev. B , 3231 (1985).[5] E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984).[6] J. F. Dobson, Phys. Rev. Lett. , 2244 (1994).[7] G. Vignale, Phys. Rev. Lett. , 3233 (1995).[8] G. Vignale and W. Kohn, Phys. Rev. Lett. , 2037 (1996).[9] P. Hessler, J. Park, and K. Burke, Phys. Rev. Lett. , 378 (1999).[10] I. D’Amico and G. Vignale, Phys. Rev. B , 7876 (1999).[11] M. Levy and A. Nagy, Phys. Rev. Lett. , 4361 (1999).[12] K. Ohno, Mater. Trans. , 649 (2007).[13] A. Crawford-Uranga, U. De Giovannini, E. R¨as¨anen, M. J. T. Oliveira, D. J. Mowbray, G. M.Nikolopoulos, E. T. Karamatskos, D. Markellos, P. Lambropoulos, S. Kurth, and A. Rubio, Phys. ev. A , 033412 (2014).[14] A. Szabo, and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced ElectronicStructure Theory (Dover, New York, 1996).[15] I. Shavitt and R. J. Bartlett, Many-Body Methods in Chemistry and Physics (Cambridge Univ., Cam-bridge, U.K., 2009).[16] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloschinski, MethodsofQuantumFieldTheoryinStatis-tical Physics (Dover, New York, 1963) Chapter 2.[17] E. M. Lifshitz and L. P. Pitaevskii, “Statistical Physics Part 2” L. D. Landau and E. M. Lifshitz, Courseof Theoretical Physics, Volume 9 Translated from the Russian by J. B. Sykes and M. J. Kearsley,(Pergamon Press, Oxford, 1980) Chapter 1, Section 8.[18] P.-O. L ¨owdin, Phys. Rev , 1474 (1955).[19] P.-O. L ¨owdin, Adv. Phys. , 1 (1956).[20] L. Hedin and S. Lundqvist, E ff ectsof Electron-Electron and Electron-Phnon Interactions on the One-Electron States of Solids in Solid State Physics Vol. 23 edited by F. Seitz, D. Turnbull, and H. Ehren-reich (Academic Press, New York, 1969) pp. 1-180.[21] O. Goscinski and P. Linder, J. Math. Phys. , 1313 (1970).[22] M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. , 246403 (2007).[23] R. Kuwahara and K. Ohno, Phys. Rev. A , 032506 (2014).[24] R. Kuwahara, Y. Noguchi, and K. Ohno, Phys. Rev. B , 121116(R) (2016).[25] V. M. Galitskii and A. B. Migdal, Sov. Phys. JETP , 96 (1958); Russian original: ZhETF , 139(1958).[26] P. S´anchez-Friera and R. W. Godby, Phys. Rev. Lett. , 5611 (2000).[27] K. Morokuma and S. Iwata, Chem. Phys. Lett. , 192 (1972).[28] R. P. Feynman, Statistical Mechanics Advanced Book Classic (Westview Press, Oxford, 1998).[29] J. M. Luttinger, Phys. Rev. , 942 (1961).[30] D. Pines TheMany-body Problem Advanced Book Classic (Westview Press, Oxford, 1997).[31] P. Nozi`eres, Theory of Interacting Fermi Systems Advanced Book Classic (Westview Press, Oxford,1997).[32] M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. , 1418 (1985).[33] P. J. Linstrom and W. G. Mallard, NIST Chemistry Webbook NIST Standard Reference DatabaseNumber 69 (National Institute of Standards and Technology, 2015) http: // webbook.nist.gov
34] S. Ono, Y. Noguchi, R. Sahara, Y. Kawazoe, and K. Ohno, Computer Phys. Comm. , 20 (2015).[35] M. J. Frisch et al. , Gaussian 09, Revision A.02, Gaussian Inc., 2009.[36] B. A. Bushaw, W. N ¨ortersh¨auser, G. W. F. Drake, and H.-J. Kluge, Phys. Rev. A , 052503 (2007).[37] A. Kramida and W. C. Martin, J. Phys. Chem. Ref. Data , 1185 (1997).[38] G. Hae ffl er, D. Hanstorp, I. Kiyan, A. E. Klinkm¨uller, U. Ljungblad, and D. J. Pegg, Phys. Rev. A ,4127 (1996).[39] E. S. Chang, J. Phys. Chem. Ref. Data , 119 (1990).[40] W. C. Martin and R. Zalubas, J. Phys. Chem. Ref. Data , 817 (1979).[41] R. Belgang, D. Sxhmidt, and P. J. West, J. Phys. (Paris) Colloques , C7-229 (1983).[42] T. Trickl, E. F. Cromwell, Y. T. Lee, and A. H. Kung, J. Chem. Phys. , 6006 (1989).[43] K. P. Huber and G. Herzberg, Molecular Spectra and Molecular Structure Vol. IV. Constants of Di-atomic Molecules (Van Nostrand Reinhold, New York, 1979).[44] K. Kimura, S. Katsumata, Y. Achiba, T. Yamazaki, and S. Iwata, Ionization Energies, Ab Initio As-signments, andValenceElectronicStructurefor200MoleculesintheHandbookofHeIPhotoelectronSpectra ofFundamental Organic Compounds (Japan Scientific Society Press, Tokyo, 1981).[45] M. W. McGeogh and R. E. Schlier, Chem. Phys. Let. , 347 (1983).[46] G. W. Erickson, J. Phys. Chem. Ref. Data , 831 (1977).[47] G. W. F. Drake, Can. J. Phys. , 586 (1988).[48] J. M. Luttinger and J. C. Ward, Phys. Rev. , 1417 (1960).[49] J. C. Ward, Phys. Rev. , 182 (1950).[50] Y. Takahashi, Il Nuovo Cimento , 371 (1957).[51] J. R. Schrie ff er, TheoryofSuperconductivity Advanced Book Classic (Westview Press, Oxford, 1999)eq. (8-92), Chap. 8-5.[52] G. Baym and L. P. Kadano ff , Phys. Rev. , 287 (1961).[53] G. Baym, Phys. Rev. , 1391 (1962).[54] F. Mauri and R. Car, Phys. Rev. Lett. , 3166 (1995)., 3166 (1995).