Many-body Green's function theory of electrons and nuclei beyond the Born-Oppenheimer approximation
aa r X i v : . [ c ond - m a t . o t h e r] J un Many-body Green’s function theory of electrons and nuclei beyondBorn-Oppenheimer approximation
Ville J. H¨ark¨onen, ∗ Robert van Leeuwen, and Eberhard K. U. Gross
1, 3 Max Planck Institute of Microstructure Physics, Weinberg 2, D-06112 Halle, Germany Department of Physics, Nanoscience Center, P.O. Box 35 FI-40014, University of Jyv¨askyl¨a, Jyv¨askyl¨a, Finland Fritz Haber Center for Molecular Dynamics, Institute of Chemistry,The Hebrew University of Jerusalem, Jerusalem 91904, Israel (Dated: June 18, 2019)The method of many body Green’s functions is used to describe an arbitrary system of electronsand nuclei in a rigorous manner given the Hamiltonian of Coulombic interactions and kinetic ener-gies. The theory given resolves the problem arising from the translational and rotational invarianceof the Hamiltonian afflicting the existing theory based on the same technique. As a result, we derivea coupled set of exact equations for the electron and nuclei Green’s functions giving a systematicway to potentially compute various properties of a rather arbitrary many-body systems of electronsand nuclei beyond Born-Oppenheimer approximation, including molecules and solids. We discuss aspecial case of crystalline solids in more detail.
I. INTRODUCTION
The Born-Oppenheimer (BO) approximation is arather central part in solving the many-body problemof electrons and nuclei in crystalline solids. This approx-imation has made the study of several physical propertiescomputationally feasible and our current understandingabout crystalline solids, for instance, is rather extensivelybased on it and its validity. In the BO approximation,the electron and nuclei problems are treated separately ina way that the nuclei coordinates are parameters in theelectronic problem and the nuclei Hamiltonian containsno operators acting on electrons . Only the BO energyfrom the electronic problem enters the nuclei Hamilto-nian. We can use various approaches, such as the den-sity functional theory or the method of many-bodyGreen’s functions , to solve the electronic problemalone. The nuclei part of the problem can be treated, forinstance, by using the technique of many-body Green’sfunctions to describe various properties of crystals aris-ing from nuclei . Usually the harmonic part of thenuclei Hamiltonian is diagonalized by a set of coordinatetransformations and the phonon Green’s function the-ory is written in terms of the transformed coordinates .Therefore we have a closed form exact solution of theharmonic problem given the BO energy and the secondorder derivatives of it with respect to the nuclei coor-dinates. In addition to the electronic structure, the BOapproximation has been successfully used to compute thephonon spectrum , thermal properties and ther-mal conductivities of various materials in some casesreproducing the experimental results with reasonable ac-curacy. In some systems, however, the BO approximationmay not be sufficiently accurate and more accurate meth-ods are needed. The validity of the BO approximationis weaker for instance in some metals and materialssuch as graphene . Further, a theory beyond the BO ap-proximation is needed in order to describe photovoltaics(solar cells, for example) , chemical reactions and even vision , just to list a few examples. For more examples,see Ref. and references therein.There are several first principle methods of goingbeyond the BO approximation. Among these arethe wave function methods , multicomponent densityfunctional theory , exact factorization , path in-tegrals with Green’s functions and the Green’s functionmethods . In the Green’s function approach westart with a general many-body Hamiltonian (essentiallycomprising the Coulombic potentials and kinetic ener-gies) and write the equations of motion (EOM) for theGreen’s functions involved in order to formally obtain thegeneral solution of the problem and then impose approx-imations. By starting with a generic Hamiltonian men-tioned above, a general theory, in most cases, within theharmonic approximation to the nuclei, has been obtainedfor the electron and nuclei Green’s functions . Inorder to calculate desired observables within this theory,we may solve a coupled set of equations for the involvedGreen’s functions: Hedin’s equations for the electronGreen’s function and related quantities and in a similarway the equations for the nuclei Green’s functions . Inmost cases, one has to impose further approximations inorder to solve the equations for the Green’s functions in-volved, but in principle we have a general theory of many-body systems of electrons and nuclei, usually within theharmonic approximation .There is, however, a known issue which prevents theuse of the existing form of the exact theory to correctlydescribe solids and molecules, for example. As has beenpointed out , the Hamiltonian used in the deriva-tion is not suitable for the description as such. This issuearises from the translational and rotational invariance ofthe original Hamiltonian we start with and lead to a con-stant electron and nuclei densities in position space forthe eigenstates of the Hamiltonian and spherically sym-metric densities for the ground state with vanishing totalangular momentum. In order to solve this issue, one canestablish a set of coordinate transformations toobtain a Hamiltonian suitable for the description and thishave already been done within the many-body Green’sfunction approach . However, no determining equationfor the nuclei Green’s functions (or for the nuclei density-density correlation function) were given in order to finda coupled set of self-consistent equations for an arbitrarysystem of electrons and nuclei. Moreover, the earlierwork was considering only crystalline solids excludingsome terms important when the theory is imposed inthe study of molecules, for example. The purpose of thepresent work is to derive such a set of equations with thegeneral Hamiltonian without any approximations estab-lished. As a special case, we discuss what is the form ofthe theory when applied to crystalline solids and compareour results with those obtained earlier.This paper is organized as follows. In Sec. II, thetransformed Hamiltonian suitable for the description ofa system of electrons and nuclei is given. The EOM forthe electrons are considered in Sec. III and for the nu-clei Green’s functions in Sec. IV. The Hedin like equa-tions for electrons are derived in Sec. V. We discuss howto determine some parameters related to the coordinatetransformations in Sec. VI. The general normal modefrequencies are derived in Sec. VII A and an expressionfor the nuclei self-energy (SE) due to the Coulombic in-teractions to the lowest order is given in Sec. VII B andis compared with the earlier results. The special case ofcrystalline solids is considered in Sec. VIII. Phonons andtheir interactions are discussed in Sec. VIII A. II. HAMILTONIAN
We denote the position coordinates of the system of N e electrons and N n nuclei in the laboratory frame as r i , i = 1 , . . . , N e , R k , k = 1 , . . . , N n . (1)The Hamiltonian of this system in position representa-tion can be written as H = T n + T e + V ee + V en + V nn , (2)where T n = − N n X k =1 ¯ h M k ∇ R k , T e = − ¯ h m e N e X i =1 ∇ r i ,V ee = N e X i,i ′ =1 ′ v ( r i , r i ′ ) , V en = N e X i =1 N n X k =1 v ( r i , R k ) ,V nn = N n X k,k ′ =1 ′ v ( R k , R k ′ ) . (3) Here V ee V en and V nn are the Coulombic potentials suchthat v ( r i , R ′ k ) = − Z k ς | r i − R ′ k | ,v ( r i , r i ′ ) = 12 ς | r i − r i ′ | ,v ( R k , R k ′ ) = 12 Z k Z k ′ ς | R k − R k ′ | , (4)where ς = e / (4 πǫ ). The Hamiltonian H given by Eq. 2is translationally and rotationally invariant and the one-body densities for the eigenstates are constants as a func-tion of position and spherically symmetric for the groundstate with vanishing total angular momentum . Inorder to derive a Hamiltonian describing correctly theeigenstates of crystalline solids, we establish a coordinatetransformation and write r ′′ i = R ( θ ) ( r i − R cmn ) , i = 1 , . . . , N e , R ′′ k = R k − R cmn , k = 1 , . . . , N n − , R ′′′ N n = R cm , (5)where R = R ( θ ) is the rotation matrix, θ = ( θ , θ , θ )is the vector of Euler angles and Euler angles θ α are as-sumed to be functions of all the N n − R ′′ . Moreover, R cm is the center-of-mass co-ordinate of the whole system of electrons and nuclei and R cmn is the nuclei center-of-mass given by R cmn = 1 M nuc N n X k =1 M k R k , M nuc ≡ N n X k =1 M k . (6)We make the coordinate transformation in two stages,first we establish r ′ i = r i − R cmn , i = 1 , . . . , N e , R ′ k = R k − R cmn , k = 1 , . . . , N n − , R ′′′ N n = R cm , (7)and then the rotation of the electronic coordinates shownin Eq. 5. Such body-fixed frame transformations are usedin the BO approximation as well when it is formulated formolecules . After transforming the Hamiltonian underthe transformation given by Eq. 7, we have H = T tot + V ′ ee + V ′ en + V ′ nn , (8)where V ′ ee V ′ en and V ′ nn are the Coulombic potentials writ-ten in terms of the primed variables and the kinetic en-ergy T tot = T e + T n is of the form T tot = − ¯ h M ∇ R cm − N e X i =1 ¯ h m i ∇ r ′ i − N n − X k =1 ¯ h M k ∇ R ′ k − ¯ h M nuc N e X i,i ′ =1 ∇ r ′ i · ∇ r ′ i ′ + ¯ h M nuc N n − X k,k ′ =1 ∇ R ′ k · ∇ R ′ k ′ . (9)After establishing the second transformation r ′′ i = R ( θ ) r ′ i and R ′′ k = R ′ k , we obtain H = T cm + T ′′ e + T ′′ n + T ′′ mpe + T ′′ mpn + T cvr + T ′ cvr + V ′′ ee + V ′′ en + V ′′ nn + V ext , (10)where V ext is the external potential part added to theHamiltonian (not originating from the coordinate trans-formation) and it is added in order to use the functionalderivative techniques in deriving the EOM. After theEOM are derived we put V ext to zero. The transformedkinetic energies are T cm = − ¯ h M ∇ R cm ,T ′′ e = − ¯ h m e N e X i =1 ∇ r ′′ i ,T ′′ n = − N n − X k =1 ¯ h M k ∇ R ′′ k ,T ′′ mpe = − ¯ h M nuc N e X i,j =1 ∇ r ′′ i · ∇ r ′′ j ,T ′′ mpn = ¯ h M nuc N n − X k,k ′ =1 ∇ R ′′ k · ∇ R ′′ k ′ , (11)where the total mass is M = M nuc + m e N e and thetransformed Coulombic potentials are V ′′ ee = N e X i,j =1 ′ v (cid:0) r ′′ i , r ′′ j (cid:1) ,V ′′ nn = N n X k,k ′ =1 ′ v ( R ′′ k , R ′′ k ′ ) ,V ′′ en = N e X i =1 N n X k =1 v ( r ′′ i , R R ′′ k ) , R ′′ N n ≡ − M N n N n − X k =1 M k R ′′ k . (12)Here, T cvr and T ′ cvr include the Coriolis and rotational-vibrational coupling terms and an explicit form of thesequantities is given in Appendix C (given in the abstractspace). The center-of-mass kinetic energy T cm commuteswith all the other terms in the Hamiltonian and doesnot enter to the EOM and is neglected from now on.There are only N n − V ′′ nn and V ′′ en involving R ′′ N n break the translational symmetry. The mass-polarizationterms T ′′ mpe and T ′′ mpn are expected to be rather small inthe case of crystals and larger molecules since they areproportional to the inverse of the total nuclei mass. So far the Euler angles in Eq. 5 are assumed to be generic func-tions of the nuclei coordinates R ′′ k ′ and we have not in-troduced any defining relations to them. There are manyways to choose these angles and here we assume that theEuler angles are defined through an implicit equation ofthe form F ( θ , R ′′ ) = . (13)For example, one possible form of Eq. 13 is the Eckartcondition and it can be written as N n − X k =1 M k x k × R R ′′ k = . (14)A relation connecting the quantities R ′′ k and x k is givenby Eq. 15. The relations given by Eqs. 13 and 14 definethe Euler angles and thus the rotation matrix as a func-tion of the N n − R = R [ θ ( R ′′ )]. Someimplicit conditions, like the Eckart condition, changesthe permutation symmetry of the Hamiltonian withrespect to our new nuclei variables leading to a morecomplicated commutators between the nuclei field oper-ators of identical nuclei and therefore changes the sym-metry of the wave functions written in terms of the trans-formed nuclei coordinates as well. However, the implicitcondition does not change canonical commutation rela-tions (Eq. 22) between the nuclei position operators.This means that the EOM for the electron and nucleiGreen’s functions remain the same whether or not theimplicit condition changes the permutation symmetry ofthe Hamiltonian with respect to nuclei variables. In thecase in which the implicit condition changes the afore-mentioned permutation symmetry and it turns out to bedifficult to apply the theory because of the complicatedsymmetry, we may assume, as in the previous works ,that the nuclei are distinguishable. By doing so we avoidthe problems potentially caused by a more complicatedpermutation symmetry. Some implicit conditions are notsuitable for describing arbitrary systems, for example, inthe case of linear molecules, the Eckart conditions as suchare not well defined . Thus, despite our theory is inprinciple general, the choice of F ( θ , R ′′ ) in Eq. 13 mayrestrict the range of validity of the theory.Next we establish yet another coordinate transforma-tion of the nuclei coordinates and write the position co-ordinates as a sum of rest positions x k (parameters) anddisplacements u k (quantum variables), namely R ′′ k = x k + u k , k = 1 , . . . , N n − . (15)We consider the determination of the rest positions x k in more detail in Sec. VI and for now these are takenas arbitrary real parameters. Under this transformation,the position variables in the Hamiltonian H transform asEq. 15 shows and in the kinetic energy terms (Eq. 11),the derivatives included transform as ∇ R ′′ k → ∇ u k . Wehave all the necessary steps made in order to write ourHamiltonian. We write the electronic parts of the Hamil-tonian in the second quantization. This can be done in asimilar way as in the laboratory frame formulation sincethe transformed Hamiltonian has the same permutationsymmetry with respect to the electronic variables as theoriginal one. The nuclei part is written in the first quan-tization. After using Eq. 15, we write the Hamiltonianoperator as (the center-of-mass kinetic energy neglectedsince it does not enter to the EOM)ˆ H = ˆ T tot + ˆ V ′′ ee + ˆ V ′′ en + ˆ V ′′ nn + ˆ V ext , (16)whereˆ T tot = ˆ T ′′ e + ˆ T ′′ n + ˆ T ′′ mpe + ˆ T ′′ mpn + ˆ T cvr + ˆ T ′ cvr . (17)Here ˆ T cvr and ˆ T ′ cvr are given in Appendix D and the otherkinetic energy terms areˆ T ′′ n = N n − X k =1 ˆ p k · ˆ p k M k , ˆ T ′′ mpn = − N n − X k,k ′ =1 ˆ p k · ˆ p k ′ M nuc , ˆ T ′′ e = − ¯ h m e Z d r ˆ ψ † ( r ) ∇ r ˆ ψ ( r ) , ˆ T ′′ mpe = − ¯ h M nuc Z d r Z d r ′ ˆ ψ † ( r ) ˆ ψ † ( r ′ ) ∇ r · ∇ r ′ × ˆ ψ ( r ′ ) ˆ ψ ( r ) . (18)The potetial energy terms in turn areˆ V ′′ ee = 12 Z d r Z d r ′ v ( r , r ′ ) ˆ ψ † ( r ) ˆ ψ † ( r ′ ) ˆ ψ ( r ′ ) ˆ ψ ( r ) , ˆ V ′′ en = Z d r Z d r ′ v ( r , r ′ ) ˆ n n ( r ′ ) ˆ n e ( r )= Z d r ˆ V en ( r ) ˆ n e ( r ) , ˆ V ′′ nn = N n X k,k ′ =1 ′ v ( x k + ˆ u k , x k ′ + ˆ u k ′ ) , (19)andˆ V ext = Z d r U ( r t ) ˆ n n ( r ) + Z d r U ( r t ) ˆ n e ( r )+ X p Z d r U p ( r t ) ˆ n ( p ) n ( r ) + N n − X k =1 J ( kt ) · ˆ u k + Z d r ϕ ( r t ) ˆ n e ( r ) , (20)where we use ˆ n e ( r ) = ˆ ψ † ( r ) ˆ ψ ( r ). In these equations v ( r , r ′ ) = ς/ | r − r ′ | andˆ n n ( r ) ≡ − N n X k =1 Z k δ (cid:16) r − R h θ (cid:16) ˆ R ′′ (cid:17)i ˆ R ′′ k (cid:17) , ˆ V en ( r ) = N n X k =1 − Z k ς (cid:12)(cid:12)(cid:12) r − R h θ (cid:16) ˆ R ′′ (cid:17)i ˆ R ′′ k (cid:12)(cid:12)(cid:12) . (21) The operators ˆ u α ( k ), ˆ p α ′ ( k ′ ), ˆ ψ † ( r ) and ˆ ψ ( r ) satisfy thefollowing commutation and anti-commutation relations[ˆ u α ( k ) , ˆ p α ′ ( k ′ )] − = i ¯ hδ αα ′ δ kk ′ , h ˆ ψ ( r ) , ˆ ψ † ( r ′ ) i + = δ ( r − r ′ ) . (22)When necessary (see Sec. VIII), we use the notation k = lκ and so on , which is perhaps more suitablein the case of crystalline solids, see Appendix A for thedescription of this notation. In the case of solids, we alsoimpose the periodic boundary conditions . So far wehave not established any approximations. We have nowset up the Hamiltonian to work with and derive the EOMfor the electron and nuclei Green’s functions by using itin Secs. III and IV. III. EQUATIONS OF MOTION FORELECTRONS
Here we derive the EOM for the electron Green’s func-tion. We denote an ensemble average for an operatorˆ o ( t ) (in the Heisenberg picture) as h ˆ o ( t ) i = X n h φ n | ˆ ρ ˆ o ( t ) | φ n i = Tr [ˆ ρ ˆ o ( t )] , (23)where {| φ n i} is a set of orthonormal basis kets. Thedensity operator is of the formˆ ρ = e − β ˆ H M Z , Z = Tr h e − β ˆ H M i , (24)where ˆ H M = ˆ H − µ e ˆ N e , ˆ N e = Z d r ˆ n e ( r ) , (25)and µ e is the chemical potential of the electrons. Thus,the density operator in Eq. 24 is chosen to describe sys-tems with a fixed number of nuclei . Here the notationis such that the time variable t could be taken as a realvariable or variable on the Keldysh contour . This isjustifiable since in both cases the EOM are the same . Ifthe time variables are taken to be variables in the Keldyshcontour, the time integrals are taken to be integrals onthe contour, see Ref. for details.We start by writing the EOM for the field operatorˆ ψ ( r t ) (operator in the Heisenberg picture), namely i ¯ h ∂∂t ˆ ψ ( r t ) = ˆ D ( r t ) ˆ ψ ( r t )+ Z d r ′ v ( r , r ′ ) ˆ ψ † ( r ′ t ) ˆ ψ ( r ′ t ) ˆ ψ ( r t ) − ¯ h M nuc Z d r ′ ˆ ψ † ( r ′ t ) ∇ r · ∇ r ′ ˆ ψ ( r ′ t ) ˆ ψ ( r t )+2 X β,σ,γ,γ ′ ˆ L (3) γσγ ′ β Z d r ′ × ˆ ψ † ( r ′ t ) r γ ′ r ′ γ ∂ ∂r ′ σ ∂r β ˆ ψ ( r ′ t ) ˆ ψ ( r t ) , (26)whereˆ D ( r t ) = − ¯ h m ∇ r + U ( r t ) + ϕ ( r t )+ Z d r ′ v ( r , r ′ ) ˆ n n ( r ′ ) + X β,γ ˆ L (1) βγ r γ ∂∂r β + N n − X k =1 X α,β,γ ˆ L (2) αβγ ( k ) ˆ p α ( k ) r γ ∂∂r β , (27)and ˆ L (1) βγ ≡ ˆ T (1) βγ + ˆ M (1) βγ , ˆ L (2) αβγ ( k ) ≡ ˆ T (2) αβγ ( k ) + ˆ M (2) αβγ ( k ) , ˆ L (3) γσγ ′ β ≡ ˆ M (3) γσγ ′ β + ˆ T (3) γσγ ′ β . (28)The quantities in Eq. 28 originate from the Coriolis andvibrational-rotational coupling terms, see Appendix C.Define the electron Green’s function as G ( r t, r ′ t ′ ) ≡ − i ¯ h Tr h e − β ˆ H M T n ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oi Tr h e − β ˆ H M i , (29)where D T n ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE = θ ( t − t ′ ) D ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) E − θ ( t ′ − t ) D ˆ ψ † ( r ′ t ′ ) ˆ ψ ( r t ) E , (30)and G ( r t, r ′ t ′ ) satisfies the Kubo-Martin-Schwingerboundary conditions . The EOM for the Green’s func-tion read i ¯ h ∂∂t G ( r t, r ′ t ′ ) = (cid:28) T (cid:26)(cid:20) ∂∂t ˆ ψ ( r t ) (cid:21) ˆ ψ † ( r ′ t ′ ) (cid:27)(cid:29) + δ ( t − t ′ ) δ ( r − r ′ ) . (31)After using Eq. 26 we find that (cid:20) i ¯ h ∂∂t − D ( r t ) (cid:21) G ( r t, r ′ t ′ )= δ ( t − t ′ ) δ ( r − r ′ ) + S ( r t, r ′ t ′ ) + S ′ ( r t, r ′ t ′ ) , (32)where the total density is defined by ˆ n ( r ′′ t ) ≡ ˆ n e ( r ′′ t ) +ˆ n n ( r ′′ t ) and S ′ ( r t, r ′ t ′ ) ≡ i ¯ h Z d r ′′ v ( r , r ′′ ) × D T n ˆ n ( r ′′ t ) ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE ,D ( r t ) ≡ − ¯ h m e ∇ r + U ( r t ) + ϕ ( r t ) . (33)By substracting the quantity (we use the shorthand no-tation given by Eq. A3 of Appendix A) Z d v (1 , h ˆ n (3) i G (1 , , (34) from both sides of Eq. 32 we find that (cid:20) i ¯ h ∂∂t + ¯ h m e ∇ r − V tot (1) (cid:21) G (1 , δ (1 −
2) + S (1 ,
2) + ˜ S (1 , , (35)where˜ S (1 , ≡ S ′ (1 , − Z d v (1 , h ˆ n (3) i G (1 , ,V tot (1) ≡ ϕ (1) + U (1) + Z d v (1 , h ˆ n (3) i . (36)All the terms originating from the transformed kineticenergies ˆ T ′′ mpe , ˆ T cvr and ˆ T ′ cvr are included in the quantity S (1 ,
2) given by S (1 ,
2) = X c =1 S c (1 , , (37)where S c (1 ,
2) for each c is given in Sec. V. Next wewrite the quantities S (1 ,
2) and ˜ S (1 ,
2) in terms of thecorresponding SE’s, define (here c = 1 , , , ≡ Z d S (1 , G − (2 , , Σ c (1 , ≡ Z d S c (1 , G − (2 , . (38)By inverting (38), we find that˜ S (1 ,
2) = Z d
3Σ (1 , G (3 , ,S c (1 ,
2) = Z d c (1 , G (3 , . (39)By using Eq. 39 in Eq. 35 δ (1 −
2) = (cid:20) i ¯ h ∂∂t + ¯ h m e ∇ r − V tot (1) (cid:21) G (1 , − Z d t (1 , G (3 , , (40)where we definedΣ t (1 , ≡ Σ (1 ,
2) + X c =1 Σ c (1 , . (41)We can also write Eq. 32 in the form of Dyson equation G (1 ,
2) = G (1 , Z d Z d G (1 ,
3) Σ t (3 , G (4 , , (42)the function G (1 ,
2) being a solution of (cid:20) i ¯ h ∂∂t + ¯ h m e ∇ r − V tot (1) (cid:21) G (1 ,
2) = δ (1 − . (43)The results obtained in the literature are similar to thoseobtained here, but here the equations are written in thebody-fixed frame in contrast to the laboratory frameformulation and the SE also contains the contribu-tions from the kinetic energy terms neglected in the ear-lier body-fixed frame formulation . By neglecting thekinetic energy terms ˆ T ′′ mpe , ˆ T cvr and ˆ T ′ cvr in the first placeleads to Eqs. 40 and 42 with Σ c (1 ,
3) = 0 for c = 1 , , t (1 ,
3) = Σ (1 , . Thus all the rather complicated kinetic en-ergy terms appearing in the theory, originating from thebody-fixed frame transformation, are hidden in the SE’sΣ c (1 , c (1 ,
3) in Eqs. 94 and 95 of Sec. V. Here we obtainedthe necessary EOM for electrons and we further write theHedin like equations for electrons in Sec. V, which maybe in some cases useful in deriving expressions for the SEoriginating from the Coulombic interaction, Σ (1 , IV. EQUATIONS OF MOTION FOR NUCLEI
Next we set out to derive the EOM for nuclei. The con-nection of the nuclei and electron equations is discussedin more detail in Sec. V where we give some justificationfor the choice of quantities for which the nuclei EOM arewritten. We start by writing the Heisenberg EOM forthe displacements, namely i ¯ h ∂∂t ˆ u α ( kt ) = h ˆ u α ( kt ) , ˆ H ( t ) i − . (44)After computing the commutators by using Eq. 16 [seealso Eq. 28] M k ∂∂t ˆ u α ( kt ) = ˆ p α ( kt ) − M k M nuc N n − X k ′′ =1 ˆ p α ( k ′′ t )+ M k X β,γ Z d r r γ ∂∂r β ˆ L (2) αβγ ( k ) ˆ n e ( r t ) . (45)After differentiating Eq. 45 with respect to time andtaking an ensemble average M k ∂ ∂t h ˆ u α ( kt ) i = K α ( kt ) + X c =1 K ( c ) α ( kt ) − J α ( kt ) , (46) where (we use ˆ T ′′ cvr ≡ ˆ T cvr + ˆ T ′ cvr ) K α ( kt ) ≡ i ¯ h (cid:28)h ˆ p α ( kt ) , ˆ V ′′ en + ˆ V ′′ nn i − (cid:29) + · · · ,K (1) α ( kt ) ≡ iM k ¯ hM nuc N n − X k ′ =1 (cid:28)h ˆ p α ( k ′ t ) , ˆ V ′′ en + ˆ V ′′ nn i − (cid:29) + · · · ,K (2) α ( kt ) ≡ i ¯ h (cid:28)h ˆ p α ( kt ) , ˆ T ′′ cvr i − (cid:29) ,K (3) α ( kt ) ≡ iM k ¯ hM nuc N n − X k ′ =1 (cid:28)h ˆ p α ( k ′ t ) , ˆ T ′′ cvr i − (cid:29) ,K (4) α ( kt ) ≡ i ¯ h X β,γ Z d r r γ ∂∂r β × (cid:28)h ˆ L (2) αβγ ( k ) ˆ n e ( r t ) , ˆ H ( t ) i − (cid:29) . (47)In the first two quantities of Eq. 47, all the other termsexcept those originating from the external potentials areexplicitly shown.In order to obtain the EOM for the nuclei Green’s func-tion we take a functional derivative of Eq. 45 with respectto J β ( k ′ t ′ ) and write M k ∂ ∂t D αβ ( kt, k ′ t ′ ) = − δ αβ δ kk ′ δ ( t − t ′ ) + K αβ ( kt, k ′ t ′ )+ X c =1 K ( c ) αβ ( kt, k ′ t ′ ) , (48)where δ h ˆ u α ( kt ) i δJ β ( k ′ t ′ ) = 1 i ¯ h hT { ˆ u α ( kt ) ˆ u β ( k ′ t ′ ) }i− i ¯ h h ˆ u α ( kt ) i h ˆ u β ( k ′ t ′ ) i≡ D αβ ( kt, k ′ t ′ ) ,K αβ ( kt, k ′ t ′ ) ≡ δK α ( kt ) δJ β ( k ′ t ′ ) ,K ( c ) αβ ( kt, k ′ t ′ ) ≡ δK ( c ) α ( kt ) δJ β ( k ′ t ′ ) . (49)We write yet another form of Eq. 48 by assuming that D αβ ( kt, k ′ t ′ ) is invertible, namely, Eq. 48 can be writtenin terms of the nuclei SE’sΠ αα ′ (cid:0) kt, ¯ k ¯ t (cid:1) ≡ − X k ′ ,β Z dt ′ K αβ ( kt, k ′ t ′ ) D − βα ′ (cid:0) k ′ t ′ , ¯ k ¯ t (cid:1) , Π ( c ) αα ′ (cid:0) kt, ¯ k ¯ t (cid:1) ≡ − X k ′ ,β Z dt ′ K ( c ) αβ ( kt, k ′ t ′ ) D − βα ′ (cid:0) k ′ t ′ , ¯ k ¯ t (cid:1) , (50)as M k ∂ ∂t D αβ ( kt, k ′ t ′ )= − X k ′′ ,α ′ Z dt ′′ Π tαα ′ ( kt, k ′′ t ′′ ) D α ′ β ( k ′′ t ′′ , k ′ t ′ ) − δ αβ δ kk ′ δ ( t − t ′ ) , (51)whereΠ tαα ′ ( kt, k ′ t ′ ) ≡ Π αα ′ ( kt, k ′ t ′ ) + X c =1 Π ( c ) αα ′ ( kt, k ′ t ′ ) . (52)Neglecting the kinetic energy terms ˆ T ′′ mpn , ˆ T cvr and ˆ T ′ cvr in the Hamiltonian given by Eq. 16 leads to Eq. (51)with Π ( c ) αα ′ ( kt, k ′ t ′ ) = 0 for all c and thus Π tαα ′ ( kt, k ′ t ′ ) =Π αα ′ ( kt, k ′ t ′ ). In other words, all the effects of the masspolarization, Coriolis and vibrational-rotational couplingterms on the EOM of the function D αβ ( kt, k ′ t ′ ) are hid-den in the SE’s Π ( c ) αα ′ ( kt, k ′ t ′ ).We also write the EOM for the momentum-displacement and momentum-momentum Green’s func-tions since these functions are in general needed in thecalculation of the total energy, see Appendix D. In or-der to use the functional derivative technique we add thefollowing potential to the external potential part of theHamiltonian N n − X k =1 3 X β =1 P β ( kt ) ˆ p β ( kt ) . (53)We note that this term would appear in the EOM of h ˆ u α ( kt ) i , but would not appear in the correspondingequations of D αβ ( kt, k ′ t ′ ). We start by writing the EOMfor momentum ensemble average, namely ∂∂t h ˆ p α ( kt ) i = K α ( kt ) + K (2) α ( kt ) − J α ( kt ) . (54)By taking a functional derivative of Eq. 54 with respectto P β ( k ′ t ′ ) we find that ∂∂t D ppαβ ( kt, k ′ t ′ ) = K ′ αβ ( kt, k ′ t ′ ) + K ′′ αβ ( kt, k ′ t ′ ) , (55)where D ppαβ ( kt, k ′ t ′ ) ≡ δ h ˆ p α ( kt ) i δP β ( k ′ t ′ ) ,K ′ αβ ( kt, k ′ t ′ ) ≡ δK α ( kt ) δP β ( k ′ t ′ ) ,K ′′ αβ ( kt, k ′ t ′ ) ≡ δK (2) α ( kt ) δP β ( k ′ t ′ ) . (56)Next we take a functional derivative of Eq. 54 with re-spect to J β ( k ′ t ′ ) and write ∂∂t D puαβ ( kt, k ′ t ′ ) = − δ αβ δ kk ′ δ ( t − t ′ ) + K αβ ( kt, k ′ t ′ )+ K (2) αβ ( kt, k ′ t ′ ) , (57) where Eq. 49 was used and D puαβ ( kt, k ′ t ′ ) ≡ δ h ˆ p α ( kt ) i /δJ β ( k ′ t ′ ). We can use Eq. 50 and write Eq.57 as ∂∂t D puαβ ( kt, k ′ t ′ )= − X k ′′ ,α ′ Z dt ′′ ˜Π tαα ′ ( kt, k ′′ t ′′ ) D α ′ β ( k ′′ t ′′ , k ′ t ′ ) − δ αβ δ kk ′ δ ( t − t ′ ) , (58)where˜Π tαα ′ ( kt, k ′ t ′ ) ≡ Π αα ′ ( kt, k ′ t ′ ) + Π (2) αα ′ ( kt, k ′ t ′ ) . (59)Moreover, we have D puαβ ( kt, k ′ t ′ ) = δ h ˆ p α ( kt ) i δJ β ( k ′ t ′ ) = δ h ˆ u β ( k ′ t ′ ) i δP α ( kt ) ≡ D upβα ( k ′ t ′ , kt ) . (60)From Eq. 58 we see that provided we have the neces-sary SE’s and the solution for D αβ ( kt, k ′ t ′ ), then we canobtain D puαβ ( kt, k ′ t ′ ) without solving the EOM for it.So far our derivation is general and we have not madeany simplifying assumptions. We still use the full Hamil-tonian without expanding in displacements ˆ u , which isthe usual procedure in the existing formulations .We can actually write the exact ensemble average of theHamiltonian in terms of the quantities appearing in theEOM, see Appendix D. With the Hamiltonian being ofgeneral form, we leave the opportunity to use other meth-ods, rather than expanding the Hamiltonian to a Taylorseries, in writing the quantities involved such as estab-lished in Ref. . We have obtained an exact Green’s func-tion theory of an arbitrary system of electrons and nucleigiven the Hamiltonian of kinetic energies and Coulombicinteractions. With some choices of F ( θ , R ′′ ) we may as-sume the nuclei to be distinguishable, see Sec. II. Thefull solution is probably rather difficult to obtain andapproximations are needed. We discuss a special caseof crystalline solids in Sec. VIII and give an explicitapproximate expression for the Fourier transformed SEΠ αα ′ ( k, k ′ , ω ) in Sec. VII B. V. HEDIN’S EQUATIONS
In order to derive the Hedin like equations, we rewriteEq. (32) by using the result D T n ˆ n ( r ′′ t ) ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE = i ¯ h δ D T n ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE δU ( r ′′ t )+ D T n ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE h ˆ n ( r ′′ t ) i , (61)and after using Eq. (61) in Eq. (32) we can write (cid:20) i ¯ h ∂∂t + ¯ h m e ∇ r − V tot (1) − i ¯ h Z d v (1 , δδU (3) (cid:21) × G (1 ,
2) = δ (1 −
2) + S (1 , . (62)Next we write Eq. (62) in terms of electron SE by usingthe following result for the electron Green’s function δG (1 , δU (3) = Z d Z d Z d G (1 ,
4) Γ (4 , , × ǫ − (6 , G (5 , , (63)where the vertex function is defined asΓ (4 , , ≡ − δG − (4 , δV tot (6) , (64)and the inverse of the dielectric function is ǫ − (1 ,
2) = δ (1 −
2) + Z d v (1 , δ h ˆ n (3) i δU (2) . (65)After using Eqs. 63 and 39 in Eq. 62, we obtain againthe EOM given by Eq. 40, but this time the SE, Σ (1 , ,
4) = i ¯ h Z d Z d W (1 , G (1 ,
3) Γ (3 , , , (66)where the screened Coulombic interaction is defined as W (1 , ≡ Z d v (1 , ǫ − (6 , . (67)Next we set out to derive the remaining Hedin’sequations in order to obtain determiningequations for the quantities like Γ (1 , ,
3) and W (1 , h ′ (1) ≡ i ¯ h ∂∂t + ¯ h m e ∇ r − V tot (1) , (68)and then invert Eq. 40 such that we have G − (1 ,
2) = h ′ (1) δ (1 − − Σ t (1 , . (69)By taking a functional derivative of Eq. 69 with respectto V tot (3) we find thatΓ (1 , ,
3) = δ (1 − δ (1 −
3) + δ Σ t (1 , δV tot (3) , (70)where the definition of Γ (1 , ,
3) was used. Next we use δ Σ t (1 , δV tot (3) = Z d Z d δ Σ t (1 , δG (4 , δG (4 , δV tot (3) , (71)and δG (4 , δV tot (3) = Z d Z d G (4 , G (7 ,
5) Γ (6 , , , (72) in Eq. 70 and writeΓ (1 , ,
3) = δ (1 − δ (1 −
3) + Z d Z d Z d Z d × δ Σ t (1 , δG (4 , G (4 , G (7 ,
5) Γ (6 , , . (73)Next we consider the screened Coulombic interactiongiven by Eq. 67 written as W (1 ,
2) = v (1 ,
2) + Z d Z d v (2 , P e (3 , W (1 , Z d Z d v (1 , δ h ˆ n n (4) i δU (3) v (2 , , (74)where the electronic polarization is defined as P e (1 , ≡ δ h ˆ n e (1) i δV tot (2) = − i ¯ h δG (1 , + ) δV tot (2) . (75)By making use of the result given by Eq. (72) we findthat the electronic polarization satisfies P e (1 ,
2) = − i ¯ h Z d Z d G (1 , G (cid:0) , + (cid:1) Γ (3 , , . (76)We rewrite now W (1 ,
2) in terms of electron and nucleiparts. In order to achieve this we write δ h ˆ n n (2) i δU (1) = − i ¯ h hT { ∆ˆ n n (2) ∆ˆ n e (1) }i + D (2 , , (77)where ∆ˆ n n (2) ≡ ˆ n n (2) − h ˆ n n (2) i , ∆ˆ n e (1) ≡ ˆ n e (1) −h ˆ n e (1) i and the so-called nuclear density-density corre-lation function is defined as D n (1 , ≡ − i ¯ h hT { ∆ˆ n n (1) ∆ˆ n n (2) }i . (78)Suppose that ˆ n ( p ) n (1) = ˆ n n (1) for some p in the externalpotential part of the Hamiltonian, then δ h ˆ n (1) i δU p (2) = D n (1 , Z d Z d P e (1 , v (3 , δ h ˆ n (4) i δU p (2) . (79)By using˜ ǫ e (1 , ≡ δ (1 − − Z d P e (1 , v (3 , , (80)we rewrite Eq. 79 as δ h ˆ n (4) i δU p (2) = δ h ˆ n n (2) i δU (4) = Z d ǫ − e (4 , D n (5 , . (81)By using Eq. 81 in Eq. 74 W (1 ,
2) = v (1 ,
2) + Z d Z d v (2 , P e (3 , W (1 , Z d Z d W e (1 , D n (5 , v (2 , , (82)where we used˜ W e (1 , ≡ Z d v (1 ,
3) ˜ ǫ − e (3 , . (83)Finally, after some re-arranging, Eq. 82 can be writtenas W (1 ,
2) = W e (1 ,
2) + W ph (1 , , (84)where W ph (1 , ≡ Z d Z d W e (1 , D n (3 , W e (4 , , (85)and W e (1 , ≡ Z d ǫ − e (2 , v (3 , ,ǫ e (1 , ≡ δ (1 − − Z d v (1 , P e (3 , . (86)By using Eq. 84 in Eq. 66, we see that the SE can bewritten as Σ (1 ,
4) = Σ e (1 ,
4) + Σ ph (1 , e (1 ,
4) = i ¯ h Z d Z d W e (1 , G (1 ,
3) Γ (3 , , , Σ ph (1 ,
4) = i ¯ h Z d Z d W ph (1 , G (1 ,
3) Γ (3 , , . (87)To summarize, the set of equations for the electronGreen’s function and the related quantities in the formof Hedin’s equations are (Eqs. 42, 66, 73, 76, 84 and 87) G (1 ,
2) = G (1 , Z d Z d G (1 ,
3) Σ t (3 , G (4 , , Γ (1 , ,
3) = δ (1 − δ (1 −
3) + Z d Z d Z d Z d × δ Σ t (1 , δG (4 , G (4 , G (7 ,
5) Γ (6 , , ,W (1 ,
2) = W e (1 ,
2) + W ph (1 , ,P e (1 ,
2) = − i ¯ h Z d Z d G (1 , G (cid:0) , + (cid:1) Γ (3 , , , Σ (1 ,
4) = Σ e (1 ,
4) + Σ ph (1 , . (88)Formally, the present theory is similar to the con-ventional ones and thus many of the approxima-tions, like the GW-approximation or other suit-able approximations , can be established in a similarway than in existing Green’s function theory for elec-trons within the BO approximation. Namely, in theGW-approximation Γ (1 , , ≈ δ (1 − δ (1 −
3) andthe electronic polarization and SE become P e (1 , ≈− i ¯ hG (1 , G (2 , + ) and Σ (1 , ≈ i ¯ hW (1 , G (1 , D n (1 , n n (1) and thus of D n (1 ,
2) by ex-panding the nuclei densities to a Taylor series in displace-ments and taking into account the terms of various order.It thus follows that in order to evaluate the expanded h ˆ n n (1) i and D n (1 , h ˆ u α ( kt ) ˆ u β ( k ′ t ) i . Similar terms appear if weexpand the nuclei terms included to the SE’s Σ c (1 , S c ( r t, r ′ t ′ ) appearing in Eq. 37 are defined as S ( r t, r ′ t ′ ) ≡ − i ¯ h ¯ h M nuc Z d r ′′ ∇ r · ∇ r ′′ × D T n ˆ n e ( r ′′ t ) ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE ,S ( r t, r ′ t ′ ) ≡ i ¯ h X β,γ r γ ∂∂r β × D T n ˆ D ′ βγ ( t ) ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE ,S ( r t, r ′ t ′ ) ≡ i ¯ h X β,σ,γ,γ ′ Z d r ′′ r γ ′ r ′′ γ ∂ ∂r ′′ σ ∂r β × D T n ˆ L (3) γσγ ′ β ˆ n e ( r ′′ t ) ˆ ψ ( r t ) ˆ ψ † ( r ′ t ′ ) oE , (89)whereˆ D ′ βγ ( t ) ≡ ˆ L (1) βγ + N n − X k =1 X α ˆ L (2) αβγ ( k ) ˆ p α ( k ) . (90)Perhaps among the simplest approximations to the SE’sΣ c ( r t, r ′ t ′ ) can be obtained if we write S ( r t, r ′ t ′ ) = − i ¯ h M nuc Z d r ′′ ∇ r · ∇ r ′′ × G ( r ′′ t, r t ; r ′ t ′ , r ′′ t ) ,S ( r t, r ′ t ′ ) ≈ X β,γ r γ ∂∂r β D ˆ D ′ βγ ( t ) E G ( r t, r ′ t ′ ) ,S ( r t, r ′ t ′ ) ≈ − hi X β,σ,γ,γ ′ Z d r ′′ r γ ′ r ′′ γ ∂ ∂r ′′ σ ∂r β × D ˆ L (3) γσγ ′ β E G ( r ′′ t, r t ; r ′ t ′ , r ′′ t ) , (91)where G (1 ,
2; 1 ′ , ′ ) ≡ − h D T n ˆ ψ (1) ˆ ψ (2) ˆ ψ † (2 ′ ) ˆ ψ † (1 ′ ) oE . (92)Similar approximation is established in Sec. VII B whenwe approximate h ˆ u α ( k t ) ˆ n e ( r t ) i ≈ h ˆ u α ( k t ) i h ˆ n e ( r t ) i .If we further make the Hartree-Fock approximation G (1 ,
2; 1 ′ , ′ ) ≈ G (1 , ′ ) G (2 , ′ ) − G (1 , ′ ) G (2 , ′ ) , (93)0we find by using Eq. 38 thatΣ ( r t, r ′′ t ′′ ) = ¯ h iM nuc Z d ¯ r ∇ r · ∇ ¯ r [ δ (¯ r − r ′′ ) G ( r t, ¯ r t ) − G (¯ r t, ¯ r t ) δ ( r − r ′′ )] δ ( t − t ′′ ) , Σ ( r t, r ′′ t ′′ ) = X β,γ D ˆ D ′ βγ ( t ) E r γ ∂∂r β δ ( r − r ′′ ) δ ( t − t ′′ ) , (94)andΣ ( r t, r ′′ t ′′ )= 2 i ¯ h X β,σ,γ,γ ′ D ˆ L (3) γσγ ′ β E Z d ¯ r r γ ′ ¯ r γ ∂ ∂ ¯ r σ ∂r β × [ G ( r t, ¯ r t ) δ (¯ r − r ′′ ) − G (¯ r t, ¯ r t ) δ ( r − r ′′ )] δ ( t − t ′′ ) . (95)After Taylor expanding ˆ D ′ βγ ( t ) and ˆ L (3) γσγ ′ β in displace-ments ˆ u , Σ (1 ,
2) and Σ (1 ,
2) become functions of h ˆ u α ( k ) i , D T αβ ( kt, k ′ t ′ ) and so on.Before continuing we give an overview of a work flowhow the coupled set of electrons and nuclei could besolved. First we give an initial guess for the rest positions x as described in Sec. VI. Given the rest positions we ex-pand all the quantities in the Hedin’s equations, whichare functions of the nuclei variables, to a Taylor series inˆ u . In the first iteration, we retain only the zeroth orderterms meaning that W ph (1 ,
2) = 0 from which followsthat W (1 ,
2) = W e (1 ,
2) and Σ (1 ,
2) = Σ e (1 ,
2) (see Eq.88). With the preceding choices, the Hedin’s equationsare actually similar to those in the BO approximation,but here written in the body-fixed frame. The zeroth or-der approximations for the SE’s Σ c (1 ,
2) can be obtainedby using the relations given by Eqs. 38 and 91. After theHedin’s equations are solved, we can calculate the nucleiSE by using the quantities obtained from the solution ofthe Hedin’s equation, see for instance Eq. 115. Given thenuclei SE, the EOM for the nuclei Green’s function canbe solved. Provided we have obtained the nuclei Green’sfunction, we can write the Hedin’s equation with non-zero W ph (1 ,
2) and solve the Hedin’s equation with theSE Σ (1 ,
2) = Σ e (1 ,
2) + Σ ph (1 , x are equilibrium positions (see Sec. VI). We iter-ate as long as we have found the equilibrium positions.Provided we have found the equilibrium positions, weseek the solution of the coupled equations such that thegrand potential, or in the case of zero temperature for-malism, the total energy is minimized. In the case whenwe have several configurations for which the rest posi-tions are equilibrium positions, the aim is to find thoseequilibrium positions which minimize the grand potential(or the total energy). VI. DETERMINATION OF REST POSITIONS
We have not yet discussed in detail how to determinethe rest positions x and so far these quantities have beenconsidered as arbitrary parameters. Consider, for in-stance, the total energy of the system, E tot , defined asan ensemble average of the Hamiltonian E tot = D ˆ H E . (96)The value of E tot must be independent of x provided ourexpansion (in displacements ˆ u ) of the SE’s, or the Hamil-tonian ensemble averege itself, converges for a given x .Hence, in principle we should have the value of E tot thesame for any given and reasonable x . However, we areprobably not able to make such an expansion up to arbi-trary order in ˆ u and thus the value of E tot may becomedependent on x due to the approximations established.For instance, if the positions x are chosen far away fromthe positions around which the nuclei would vibrate, thenin order to calculate the reasonable values of E tot , the dis-placements u would be rather large. This implies that wewould need a rather high order Green’s functions for nu-clei or alternatively a rather high order expression for thenuclei SE to describe the system properly. If we choosethe positions x rather poorly and at the same time takeonly the lowest order quantities in ˆ u into account, we ex-pect to find rather unrealistic results. Even though ourtheory is generally valid and we do not even need to havethe nuclei to vibrate around some particular region ofspace, it is beneficial to have some consistent strategy toobtain the parameters x such that the systems in whichthe nuclei perform such rather small vibrations aroundsome regions of space we can describe the system with areasonable accuracy by using the approximations of lowerorders in ˆ u .For the aforementioned systems, like some stable crys-tal lattices, we use the following method to determinethe positions x . The initial value of the x is obtained inthe same way as within the BO approximation, that is,from the experimental data of known structures or in thecase of hypothetical structures by using the methods ofstructural chemistry , for instance. After the initialguess, we aim to find the positions x such that they areequal to the most probable positions of the nuclei, thatis, we seek the positions such that x k = D ˆ R ′′ k E = x k + h ˆ u k i . (97)We note that h ˆ u k i is a function of x since the Hamil-tonian is and if the external potentials vanish, D ˆ R ′′ k E isindependent of time. If Eq. 97 holds, then it followsthat h ˆ u k i = . That is, if we are able to choose thepositions x as in Eq. 97, then the expected values ofthe displacements vanish. Next we seek a way how tofind the displacements satisfying h ˆ u k i = and thus therest positions x satisfying Eq. 97. For some systems, like1Bravais crystals, this condition maybe automatically sat-isfied due to symmetry . The determining equation for h ˆ u k i is given by Eq. 46. If we put the external potentialsto zero, then h ˆ u α ( kt ) i = h ˆ u α ( k ) i and thus M k ∂ ∂t h ˆ u α ( kt ) i = 0 , (98)meaning that the expected value of the total force on eachnuclei k vanish for all α . Thus, if the external potentialsvanish, Eq. 46 becomes0 = K α ( kt ) + X c =1 K ( c ) α ( kt ) . (99)Without the external potentials, K α ( kt ) and K (1) α ( kt )are of the form given by Eq. 47 with the visible termstaken into account while the rest of the terms K ( c ) α ( kt )remain the same. We write F α ( kt ) ≡ K α ( kt ) + X c =1 K ( c ) α ( kt ) , (100)and further F α ( kt ) = F ′ α ( kt ) + F ′′ α ( kt ), where F ′ α ( kt )is the part of F α ( kt ) which is a function of h ˆ u α ( k ) i and F ′′ α ( kt ) is the remaining part (not a function of h ˆ u α ( k ) i ).We can therefore write Eq. 99 as0 = F ′ α ( kt ) + F ′′ α ( kt ) . (101)This is the determining equation for h ˆ u α ( k ) i . If we nowwant to find x such that h ˆ u α ( k ) i = 0, then Eq. (101)becomes 0 = F ′′ α ( kt ) . (102)We now give an example such that we find an explicitand approximate form of Eq. 102. Suppose that the onlynon-zero term in Eq. 99 is K α ( kt ) meaning that we takeinto account terms which originate from the Coulom-bic potentials and neglect all the other terms. This isexpected to be relatively good approximation for somecrystalline solids provided we define the Euler angles suchthat the Coriolios and vibrational-rotational terms in theHamiltonian are small, see Sec. VIII. We approximate K α ( kt ) by expanding ˆ V ′′ en and ˆ V ′′ nn in displacements (seeAppendix D) and retain only the lowest orders. After cal-culating the commutators, we obtain to the lowest ordersin displacements K α ( kt ) ≈ − ∂V nn ( x ) ∂x α ( k ) − Z d r ∂V en ( r , x ) ∂x α ( k ) h ˆ n e ( r t ) i− X k ,α Z d r ∂ V en ( r , x ) ∂x α ( k ) ∂x α ( k ) × h ˆ u α ( k t ) ˆ n e ( r t ) i− X k ,α ∂ V nn ( x ) ∂x α ( k ) ∂x α ( k ) h ˆ u α ( k t ) i . (103) By using Eq. 113 we find that hT { ˆ u β ( k ′ t ′ ) ˆ n e ( r t ) }i≈ h ˆ n e ( r t ) i h ˆ u β ( k ′ t ′ ) i + i ¯ h Z d r Z dt × Z d r Z dt P e ( r t, r t ) ˜ W e ( r t , r t ) × X k ,α ∂n (0) n ( x , r ) ∂x α ( k ) D α β ( k t , k ′ t ′ ) . (104)Here, the quantity n (0) n ( x , r ) is given by Eq. D2 and re-lated to this we consider how to obtain the rotation ma-trix R [ θ ( x )] and its derivatives in Appendix B. When weuse Eq. 104 in Eq. 103 in approximating h ˆ u α ( kt ) ˆ n e ( r t ) i ,it can be seen that if we seek the solution such that h ˆ u k i = , meaning that F ′ α ( kt ) vanishes, Eq. 102 isof the following form0 ≈ ∂V nn ( x ) ∂x α ( k ) + Z d r ∂V en ( r , x ) ∂x α ( k ) h ˆ n e ( r t ) i + i ¯ h X k ,α Z d r ∂ V en ( r , x ) ∂x α ( k ) ∂x α ( k ) Z d r Z dt × Z d r Z dt P e ( r t, r t ) ˜ W e ( r t , r t ) × X k ,α ∂n (0) n ( x , r ) ∂x α ( k ) D T α β ( k t , k t ) , (105)where (Eq. 49) D T αβ ( kt, k ′ t ′ ) ≡ i ¯ h hT { ˆ u α ( kt ) ˆ u β ( k ′ t ′ ) }i . (106)Since the external potentials vanish, then D T αβ ( kt, k ′ t ′ )can be obtained from Eq. 51 due to the fact that h ˆ u α ( kt ) i = h ˆ u α ( k ) i . If we further approximate and ne-glect the last term of Eq. 105 we obtain0 = ∂V nn ( x ) ∂x α ( k ) + Z d r ∂V en ( r , x ) ∂x α ( k ) h ˆ n e ( r ) i . (107)Equations 105 and 107 are examples of approximate con-ditions for choosing the parameters x such that Eq. 97holds. We call the positions x the equilibrium positionsif Eq. 97 holds. With these approximations established,we have formally obtained the same condition (Eq. 107)to choose the parameters than in the BO approximationin the case of crystals , when the Hellmann-Feynmantheorem is used to calculate the total force. Thedifference is that our parameters x are in the body-fixed frame and the electron density h ˆ n e ( r ) i is an en-semble average taken with respect to the full many-bodyHamiltonian instead of an expected value within the BOapproximation taken with respect to an eigenstate ofthe electronic BO Hamiltonian. The electron densitycan be obtained from the electron Green’s function as h ˆ n e ( r t ) i = − i ¯ hG ( r t, r t + ). Thus, for a given x one solvesthe coupled set of equations for the electrons and nuclei2and from the solutions of these equations we obtain thequantities like G ( r t, r ′ t ′ ), P e ( r t, r ′ t ′ ), ˜ W e ( r t, r ′ t ′ ) and D T αβ ( kt, k ′ t ′ ) and we can check whether or not Eq. 105,or the corresponding general expression given by Eq. 102,holds and indicates that we have found the positions suchthat Eq. 97 is satisfied. VII. NUCLEI RELATED CONSIDERATIONSA. Nuclei vibrations
We start by writing the nuclei EOM in frequency space,but before Fourier transforming, all the time-dependentpotentials are set equal to zero. It follows that the nucleiGreen’s function and the nuclei SE are then functions oftime differences only and we can write Eq. 51 in terms ofFourier transformed quantities as (the convention used isgiven in Appendix A) δ αβ δ kk ′ = X k ′′ ,α ′ (cid:2) M k ω δ αα ′ δ kk ′′ − Π tαα ′ ( k, k ′′ , ω ) (cid:3) × D α ′ β ( k ′′ , k ′ , ω ) . (108)We write Eq. 108 in matrix form and re-arrange suchthat ˜ D ( ω ) = (cid:2) ω I − C t ( ω ) (cid:3) − , (109)where ˜ D ( ω ) ≡ M / D ( ω ) M / and C t ( ω ) ≡ M − / Π t ( ω ) M − / . The components of these matri-ces are denoted by αk and so on. Next we write the SEas C t ( ω ) = C A + C NA ( ω ) and assume that C A is theHermitian part of C t ( ω ), which we assume to be inde-pendent of frequency and C NA ( ω ) the non-Hermitianpart of C t ( ω ). One possible choise is C A = C ′ t (0)and C NA ( ω ) = C t ( ω ) − C ′ t (0) such that we collect allthe Hermitian zero frequency contributions of C t ( ω ) to C ′ t (0). For a Hermitian matrix, we can find a completeand orthonormal set of eigenvectors satisfying˜ ω j v α ( k | j ) = X k ′ ,β C Aαβ ( kk ′ ) v β ( k ′ | j ) ,δ jj ′ = X k,α v α ( k | j ′ ) v ∗ α ( k | j ) ,δ αβ δ kk ′ = X j v α ( k | j ) v ∗ β ( k ′ | j ) . (110)We call the quantities ˜ ω j the adiabatic normal mode fre-quencies. The eigenvalues ˜ ω j are real since C A is Her-mitian and ˜ ω j are real if C A is positive definite. Wetransform Eq. 109 by using the eigenvectors of C A suchthat D ( ω ) = h ω I − ˜ ω − C NA ( ω ) i − , (111)where D ( ω ) ≡ v † ˜ D ( ω ) v and C NA ( ω ) ≡ v † C NA ( ω ) v .If C NA ( ω ) is sufficiently small, the quasi-particle pic-ture holds and C NA ( ω ) can be pictured in generating interactions between adiabatic normal modes by shiftsof their values and finite lifetimes [imaginary part of C NA ( ω )]. Correct treatment requires the continua-tion of the frequency variable to the complex frequencyplane. The usual procedure is to consider the Mat-subara Green’s functions continued to the complex fre-quency plane . However, we can also obtain theGreen’s function of complex frequency by just replacingthe real frequency ω by the complex one in Eq. 111.We point out that the present discussion is general with-out any approximations made. We establish a similarconsideration of vibrations by using phonon basis in Sec.VIII A. B. Nuclei self-energy
In order to have an approximate form for the SEΠ αα ′ ( k, k ′ , ω ), we start by giving explicit form for some ofthe terms included in K αβ ′ ( kt, k ′ t ′ ) given by Eq. 49. Wehave already calculated the commutators included to thelowest orders in ˆ u α ( kt ) and we can obtain K αβ ′ ( kt, k ′ t ′ )by taking a functional derivative of K α ( kt ) given by Eq.103 with respect to J β ( k ′ t ′ ), namely K αβ ( kt, k ′ t ′ )= − X k ′′ ,α ′′ ∂ V nn ( x ) ∂x α ( k ) ∂x α ′′ ( k ′′ ) D α ′′ β ( k ′′ t, k ′ t ′ ) − X k α Z d r ∂ V en ( r , x ) ∂x α ( k ) ∂x α ( k ) δ h ˆ u α ( k t ) ˆ n e ( r t ) i δJ β ( k ′ t ′ ) − Z d r ∂V en ( r , x ) ∂x α ( k ) δ h ˆ n e ( r t ) i δJ β ( k ′ t ′ ) + O (cid:0) ˆ u (cid:1) . (112)Actually, the terms visible in Eq. 112 are all the termsincluded in K αβ ′ ( kt, k ′ t ′ ) if the Hamiltonian is expandedin displacements up to second order before writing theEOM. This is the usual procedure in the earlier labora-tory frame formulations . Next we use the followingresult δ h ˆ n e ( r t ) i δJ β ( k ′ t ′ ) = Z d r Z dt Z d r Z dt P e ( r t, r t ) × ˜ W e ( r t , r t ) δ h ˆ n n ( r t ) i δJ β ( k ′ t ′ ) , (113)where ˜ ǫ e ( r t , r t ) and ˜ W e ( r t , r t ) are given by Eqs.80 and 83, respectively. These quantities can be obtainedfrom the solution of the Hedin’s equations. We also writefor the quantity in the second term on the right hand sideof Eq. 112 δ h ˆ u α ( k t ) ˆ n e ( r t ) i δJ β ( k ′ t ′ ) = i ¯ h δ h ˆ n e ( r t ) i δJ β ( k ′ t ′ ) δJ α ( k t )+ h ˆ u α ( k t ) i δ h ˆ n e ( r t ) i δJ β ( k ′ t ′ )+ h ˆ n e ( r t ) i D α β ( k t, k ′ t ′ ) . (114)3Finally, by using Eqs. 112-114 in Eq. 50, we can writefor the nuclei SE under consideration in frequency spaceΠ αβ ( k, k ′ , ω )= ∂ V nn ( x ) ∂x α ( k ) ∂x β ( k ′ ) + Z d r ∂ V en ( r , x ) ∂x α ( k ) ∂x β ( k ′ ) h ˆ n e ( r ) i + Z d r Z d r ′ ∂n (0) n ( x , r ) ∂x α ( k ) h ˜ W e ( r , r ′ , ω ) − v ( r , r ′ ) i × ∂n (0) n ( x , r ′ ) ∂x β ( k ′ ) + X k ,α ∂ n (0) n ( x , r ′ ) ∂x α ( k ) ∂x β ( k ′ ) h u α ( k ) i + X k ,α Z d r Z d r ′ ∂ n (0) n ( x , r ) ∂x α ( k ) ∂x α ( k ) × h ˜ W e ( r , r ′ , ω ) − v ( r , r ′ ) i h ˆ u α ( k ) i× ∂n (0) n ( x , r ′ ) ∂x β ( k ′ ) + X k ,α ∂ n (0) n ( x , r ′ ) ∂x α ( k ) ∂x β ( k ′ ) h u α ( k ) i + · · · . (115)In Eq. 115, ˜ W e ( r , r ′ , ω ) is the Fourier transform of Eq.83 in the relative time variable. Terms explicitly shownin Eq. 115 are already included in the harmonic approx-imation, which means that the Hamiltonian is expandedup to second order in displacements (before writing theEOM). Actually not all the terms in the harmonic ap-proximation are even visible in Eq. 115. Namely, thefirst term on the right hand side of Eq. 114 is not shownand in practise this means that we have approximated h ˆ u α ( k t ) ˆ n e ( r t ) i ≈ h ˆ u α ( k t ) i h ˆ n e ( r t ) i . Further, onlysome of the second order terms of the right hand side ofEq. 113 are included. Sometimes it is convenient to writeΠ αβ ( k, k ′ , ω ) = Π Aαβ ( k, k ′ ) + Π NAαβ ( k, k ′ , ω ). The adia-batic SE Π Aαβ ( k, k ′ ) = Π αβ ( k, k ′ ,
0) is real (at least whenterms visible in Eq. 115 are taken into account) whilethe non-adiabatic part Π
NAαβ ( k, k ′ , ω ) = Π αβ ( k, k ′ , ω ) − Π αβ ( k, k ′ ,
0) is complex, in general. The adiabatic SE isreal since ˜ W e ( r t, r ′ t ′ ) is. In the present case (Eq. 115),the quantity Π nnαβ ( kt, k ′ t ′ ) appearing in Eq. D19 isΠ nnαβ ( kt, k ′ t ′ ) = ∂ V nn ( x ) ∂x α ( k ) ∂x β ( k ′ ) δ ( t − t ′ ) , (116)while the rest of the Fourier transforms of the terms visi-ble in Eq. 115 are included to Π enαβ ( kt, k ′ t ′ ). For a bettercomparison with the existing theory, we use Eq. 115 andwrite the non-adiabatic contribution asΠ NAαβ ( k, k ′ , ω ) = Z d r Z d r ′ ∂n (0) n ( x , r ) ∂x α ( k ) × h ˜ W e ( r , r ′ , ω ) − ˜ W e ( r , r ′ , i × ∂n (0) n ( x , r ′ ) ∂x β ( k ′ ) + · · · , (117) where terms including the quantity h ˆ u α ( k ) i are not ex-plicitly shown. Provided we can choose the parameters x such that Eq. 97 holds, then terms involving the quan-tities h ˆ u α ( k ) i vanish, see Sec. VI. The relation givenby Eq. 117 is formally analogous with the result ob-tained earlier . However, the difference is in the densi-ties n (0) n ( x , r ′ ) (see Eq. D2) and here all the variables arein the body-fixed frame. We note that if necessary, thehigher order expressions for the nuclei and phonon SE’scan be obtained by systematically including higher orderterms in the expansion of K αβ ( kt, k ′ t ′ ). Related to this,some results of a recent study on generic electron-bosonHamiltonians may be useful within the present theory aswell.Even though we are not giving explicit expressionsfor the SE’s Π ( c ) αα ′ ( kt, k ′ t ′ ), these quantities can be ob-tained with the same procedure as established here forΠ αα ′ ( kt, k ′ t ′ ). Namely, calculate the commutators con-tained in K ( c ) α ( kt ) (Eq. 47) by expanding the quantitiesinvolved in ˆ u , take the functional derivative with respectto J β ( k ′ t ′ ) in order to obtain K ( c ) αβ ( kt, k ′ t ′ ) and then useEq. 50 to obtain Π ( c ) αα ′ ( kt, k ′ t ′ ). In obtaining these ex-pressions, similar approximations may be established aswe did in the case of Π αα ′ ( kt, k ′ t ′ ). VIII. CRYSTALLINE SOLIDS
Here we apply the general theory to crystalline solids.First of all, the electron and nuclei mass polarizationterms ˆ T ′′ mpe , ˆ T ′′ mpn and also ˆ T ′ cvr appearing in the trans-formed kinetic energy are proportional to the inverse ofthe total nuclei mass and are thus small for crystals. Fur-ther, we assume that one can find an implicit conditionof the form (Eq. 13) such that the contributions from thekinetic energy ˆ T cvr become small. Some justifications fortheir neglect have been given , when the implicit con-dition is chosen to be the Eckart condition given by Eq.14. By neglecting the aforementioned terms the EOMfor electrons remain otherwise the same, except that wehave Σ t (1 ,
3) = Σ (1 ,
3) in Eqs. 40 and 42. In turn, thenuclei EOM remain otherwise the same, but in Eq. 51the SE is Π tαα ′ ( kt, k ′ t ′ ) = Π αα ′ ( kt, k ′ t ′ ). A. Phonons and their interactions
Here we use the notation k = lκ and so on, see Ap-pendix A. The aim is to define phonon frequencies beyondthe BO approximation. In the present case, Eq. 108 canbe written as X l ′′ ,κ ′′ ,α ′ (cid:2) M κ ω δ αα ′ δ ll ′′ δ κκ ′′ − Π αα ′ ( lκ, l ′′ κ ′′ , ω ) (cid:3) × D α ′ β ( l ′′ κ ′′ , l ′ κ ′ , ω ) = δ αβ δ ll ′ δ κκ ′ . (118)We also note that the quantities like D αβ ( lκ, l ′ κ ′ , ω ) andΠ αβ ( lκ, l ′ κ ′ , ω ) are dependent only on the difference the4cell indices l − l ′ . Therefore, we write D αβ ( lκ, l ′ κ ′ , ω )as a discrete Fourier transform in the relative coordinateallowed by the periodic boundary conditions , namely D αα ′ ( κκ ′ , l, ω ) = 1 N N X q D αα ′ ( κκ ′ , q , ω ) e i q · x l , (119)and in a similar way for Π αα ′ ( κκ ′ , l, ω ). In Eq. 119, N is the number of q points and thus number of unitcells in the Born-von Karman cell . By using Eq. 119and similar expression for Π αα ′ ( κκ ′ , l ′′ , ω ) in Eq. 118 weobtain X κ ′′ ,α ′ (cid:2) M κ ω δ αα ′ δ κκ ′′ − Π αα ′ ( κκ ′′ , q , ω ) (cid:3) × D α ′ β ( κ ′′ κ ′ , q , ω ) = δ αβ δ κκ ′ , (120)and after some re-arranging and by using matrix notation˜ D ( q , ω ) = (cid:2) ω I − C ( q , ω ) (cid:3) − , (121)where C ( q , ω ) ≡ M − / Π ( q , ω ) M − / , ˜ D ( q , ω ) ≡ M / D ( q , ω ) M / . (122)The components of the matrix like C ( q , ω ) are labeledby ακ and βκ ′ , for instance. We establish a similar pro-cedure as in Sec. VII A, namely, we write the SE as C ( q , ω ) = C A ( q )+ C NA ( q , ω ), where C A ( q ) = C ( q , C NA ( q , ω ) = C ( q , ω ) − C ( q , C A ( q ) can be written as ω q j e α ( κ | q j ) = X κ ′ ,β C Aαβ ( κκ ′ | q ) e β ( κ ′ | q j ) . (123)We call the quantities ω q j the adiabatic phonon fre-quencies. The eigenvalue equation given by Eq. 123is analogous to the eigenvalue equation written for thedynamical matrix in the conventional theory of latticedynamics . As the adiabatic phonon frequencies ω q j are defined by Eq. 123, in contrast to the conventionaltheory, already the non-interacting phonons potentiallycontain the terms of the Hamiltonian of arbitrary orderin ˆ u . In order for ω q j to be positive and thus ω q j tobe real, the matrix C A ( q ) have to be positive definite,which is not in general the case for a given x . This is alsothe case in the BO theory of lattice dynamics , wherethe positive definiteness of the dynamical matrix impliesthe minimum of the BO energy surface and the stabilityof the crystal lattice . Provided the matrix C A ( q ) isHermitian (as it is, for example, if Eq. 115 is used), thecomponents of the eigenvector e α ( κ | q j ) can be chosen tosatisfy the orthonormality and completeness conditions X κ,α e α ( κ | q j ′ ) e ∗ α ( κ | q j ) = δ jj ′ , X j e α ( κ | q j ) e ∗ β ( κ ′ | q j ) = δ αβ δ κκ ′ . (124) We use the eigenvectors of the adiabatic SE to transformEq. 121 to equation in phonon basis and after we estab-lish the transformation D ( q , ξ ) = h(cid:0) ξ − ω q (cid:1) I − C NA ( q , ξ ) i − , (125)where D ( q , ξ ) ≡ e † ( q ) ˜ D ( q , ξ ) e ( q ) and C NA ( q , ξ ) ≡ e † ( q ) C NA ( q , ξ ) e ( q ). We call D ( q , ξ ) the phononGreen’s function and C NA ( q , ξ ) the non-adiabaticphonon SE. Here, ξ is the complex frequency variable .If the non-adiabatic SE vanishes, we obtain the adiabaticphonon Green’s function D Ajj ( q , ω ) = 1 / (cid:0) ω − ω q j (cid:1) which is of the usual form appearing in the BO theoryof lattice dynamics . The adiabatic phonons have infi-nite lifetimes if the non-adiabatic SE vanishes. The non-adiabatic part can be pictured in generating interactionsbetween the adiabatic phonons which appear as shiftsto the adiabatic eigenvalues and finite lifetimes of thesequasiparticles. This picture is valid if the non-adiabaticSE is sufficiently small in comparison to the adiabaticpart. If the non-diagonal terms of C NA ( q , ξ ) are rela-tively small such that C NAjj ′ ( q , ξ ) ≈ j = j ′ , Eq.125 becomes D jj ( q , ξ ) = 1 ξ − ω q j − C NAjj ( q , ξ ) , (126)from which the shifts to the adiabatic phonon frequen-cies and finite lifetimes of adiabatic phonons are usuallydeduced .The non-adiabatic phonon SE in its general form canbe written as C NAjj ′ ( q , ξ ) = X l,κ,α X κ ′ ,α ′ e ∗ α ( κ | q j ) Π NAαα ′ ( κκ ′ , l, ξ ) √ M κ M κ ′ × e α ′ ( κ ′ | q j ′ ) e − i q · x l . (127)By using Eqs. 80 and 83 in Eq. 117 and then Fouriertransforming and changing the representation, we ob-tain the following approximate form for the non-adiabaticphonon SE C NAjj ′ ( q , ξ ) = 1 N Z d r Z d r ′ (cid:2) g q j ( r , ξ ) P e ( r , r ′ , ξ ) ˜ g ∗ q j ′ ( r ′ ) − g q j ( r , P e ( r , r ′ ,
0) ˜ g ∗ q j ′ ( r ′ ) (cid:3) . (128)In Eq. 128˜ g ∗ q j ( r ) = X l,κ,α M − / κ e α ( κ | q j ) e i q · x l ∂V en ( x , r ) ∂x α ( κl ) ,g q j ( r , ξ ) = Z d r ′ ˜ g q j ( r ′ ) ˜ ǫ − e ( r ′ , r , ξ ) , (129)where ˜ g q j ( r ) is the complex conjugate of ˜ g ∗ q j ( r ). Byre-arranging Eq. 129 we find that g q j ( r , ξ ) satisfies thefollowing equation g q j ( r , ξ ) = ˜ g q j ( r ) + Z d r ′ Z d r ′′ g q j ( r ′ , ξ ) P e ( r ′ , r ′′ , ξ ) × v ( r ′′ , r ) . (130)5The diagrams corresponding to D ( q , ξ ), C NA ( q , ξ ) and g q j ( r , ξ ) are of a similar form as in the laboratory frameformulation and are given in Ref. . B. Momentum functions
By establishing the approximations discussed at thebeginning of this section and then comparing Eqs. 51and 58 we see that ∂∂t D puαβ ( kt, k ′ t ′ ) = M k ∂ ∂t D αβ ( kt, k ′ t ′ ) , (131)and thus for the Fourier transforms D puαβ ( k, k ′ , ω ) = − iM k ωD αβ ( k, k ′ , ω ) . (132)Therefore, after we have obtained a solution for D αβ ( lκ, l ′ κ ′ , ω ) we can also obtain D puαβ ( lκ, l ′ κ ′ , ω ). Thelast function we need is the momentum Green’s functionand in the present case Eq. 55 becomes ∂∂t D ppαβ ( kt, k ′ t ′ ) = K ′ αβ ( kt, k ′ t ′ ) , (133)where we approximate K ′ αβ ( kt, k ′ t ′ ) ≈ − X k ,α ∂ V nn ( x ) ∂x α ( k ) ∂x α ( k ) δ h ˆ u α ( k t ) i δP β ( k ′ t ′ ) − Z d r ∂V en ( r , x ) ∂x α ( k ) δ h ˆ n e ( r t ) i δP β ( k ′ t ′ ) − X k ,α Z d r ∂ V en ( r , x ) ∂x α ( k ) ∂x α ( k ) × δ h ˆ u α ( k t ) ˆ n e ( r t ) i δP β ( k ′ t ′ ) + O (cid:0) ˆ u (cid:1) , (134)which can be obtained directly from Eq. 103. By us-ing the same approximations for K ′ αβ ( kt, k ′ t ′ ), as in Sec.VII B were established in writing K αβ ( kt, k ′ t ′ ), we findthat Eq. 133 can be written as ∂∂t D ppαβ ( kt, k ′ t ′ ) = − X k ′′ ,α ′ Z dt ′′ Π αα ′ ( kt, k ′′ t ′′ ) × D upα ′ β ( k ′′ t ′′ , k ′ t ′ ) . (135)The Fourier transform of this equation is iωD ppαβ ( k, k ′ , ω ) = X k ′′ ,α ′ Π αα ′ ( k, k ′′ , ω ) D upα ′ β ( k ′′ , k ′ , ω ) . (136)Since by Eq. 60, D upβα ( k ′ t ′ , kt ) = D puαβ ( kt, k ′ t ′ ), we canwrite D ppαβ ( k, k ′ , ω ) = − M k X k ′′ ,α ′ Π αα ′ ( k, k ′′ , ω ) × D βα ′ ( k ′ , k ′′ , − ω ) , (137)where Eq. 132 was used. Now we have all the necessaryequations in order to calculate the total energy of thesystem, for example. IX. CONCLUSIONS
In this work, we derived a coupled and self-consistentset of exact equations for the electron and nuclei Green’sfunctions following from the Hamiltonian of Coulombicinteractions and kinetic energies as a starting point. Thepresent theory, when applied to crystalline solids, resem-bles in many parts with the previous ones . However,we take into account an issue arising from the transla-tional and rotational symmetry of the Hamiltonian pre-venting the use of the existing Green’s function theoriesto describe systems other than those with constant den-sity eigenstates. The present theory is formally exact andit is not limited to the harmonic approximation. More-over, the theory developed can be applied to arbitrarymolecules, which in general, were out of the range of va-lidity of the previous theories based on the same method.The complexity of the given system of EOM is of the sameorder with the existing ones when applied to crystallinesolids.In addition to the general EOM, we considered the nor-mal modes of arbitrary molecular systems, a special caseof crystalline solids and reproduced the analogs of someof the results appearing in the existing theories. For in-stance, phonons and their interactions were defined in arather general way. While it is probably not a realisticgoal to obtain the solution of the EOM in general form inarbitrary systems, there is some work left to derive com-putationally accessible approximations to be used in theactual calculations. Our main emphasis in this work wason crystalline solids and we leave a more detailed treat-ment of molecules for future work. The present theoryallows one to go beyond BO approximation in a system-atic and exact way and therefore it may turn out to beuseful in understanding the behaviour of non-relativisticquantum mechanical many-body systems of interactingelectrons and nuclei.
Appendix A: Notation
The following notations are usedˆ u α ¯ m ( k ¯ m ) ≡ ˆ u α ( k ) · · · ˆ u α m ( k m ) , X k ¯ m ,α ¯ m = X k ,α · · · X k m ,α m , (A1)and so on. The following notation, suitable for the de-scription of crystalline solids, is essentially the same asRefs. and . We denote the nuclei with the labels k = lκ , and write x k = x lκ = x l + x κ , where x ( κ ) isthe position vector of the nuclei κ within the unit celland x l ≡ x l l l = l a + l a + l a , (A2)is the lattice translational vector of the l th unit cell andthe vectors a j are called the primitive translational vec-tors of the lattice. In Eq. A2, l , l , l are integers. With6this notation, we write R ′′ k = R ′′ lκ = x lκ + u lκ , where lκ go over N n − x k and u k , we sometimes use the following collectivenotation x ≡ x , . . . , x N n − and u ≡ u , . . . , u N n − .We use the following shorthand notation for the posi-tion and time variables i ≡ r i t i , δ ( i − j ) ≡ δ ( t i − t j ) δ ( r i − r j ) ,v (1 , ≡ δ ( t − t ) v ( r , r ) , Z di ≡ Z d r i Z dt i . (A3)In this work, we use the following convention for theFourier transforms f ( ω ) = Z ∞−∞ dtf ( t ) e iωt ,f ( t ) = 12 π Z ∞−∞ dωf ( ω ) e − iωt . (A4) Appendix B: Euler angles
The condition given by Eq. 13 defines the Euler anglesby an implicit equation of the form F ( θ , x + u ) = . Ingeneral, no explicit solution of such equations exists ,however. That is, it might not be possible to write θ asan explicit function of the nuclei variables x , u , but fora given x , u and the nuclei masses, a numerical solutioncould be sometimes obtained . We are in particular in-terested how to obtain the rotation matrix R [ θ ( x )] andits derivatives since these quantities appear in the Hamil-tonian if we establish the Taylor expansions of some partsof the Hamiltonian (see Appendix D) and thus they willeventually appear in the EOM.The rotation matrix R [ θ ( x + u )] can be obtainedfrom the generic conditions given by Eq. 13. In turn, R [ θ ( x )] can be obtained by solving Eq. 13 with u = 0.Thus, R [ θ ( x )] is obtained from the implicit equation F ( θ , x ) = . (B1)For instance, in the case of Eckart conditions given byEq. 14 N n − X k =1 M k x k × R ( θ ) x k = , (B2)with the Euler angles θ ( x ), such that 0 < θ β ( x ) < π for β = 1 , ,
3. What we need next are the derivatives of R [ θ ( x )] appearing in the Hamiltonian (see Eqs. D1 and D2). By the chain rule, we write ∂ R ∂x α ( k ) = X β =1 ∂ R ∂θ β ∂θ β ∂x α ( k ) ,∂ R ∂x β ( k ′ ) ∂x α ( k ) = X β ′ ,β ′′ =1 ∂ R ∂θ β ′ ∂θ β ′′ ∂θ β ′′ ∂x β ( k ′ ) ∂θ β ′ ∂x α ( k )+ X β ′ =1 ∂ R ∂θ β ′ ∂ θ β ′ ∂x α ( k ) ∂x β ( k ′ ) . (B3)Given the explicit form for the rotation matrix as a func-tion of the Euler angles we can calculate the derivativeslike ∂ R /∂θ β ′ so next we find a way how to calculate thederivatives like ∂θ β ′ /∂x α ( k ) and ∂ θ β ′ /∂x α ( k ) ∂x β ( k ′ ).After taking the total derivative of F σ ( θ , x , u ) with re-spect to x α ( k ) we write dF σ dx α ( k ) = X β ′ =1 A σβ ′ ∂θ β ′ ∂x α ( k ) + ∂F σ ∂x α ( k ) = 0 , (B4)where A σβ ′ ≡ ∂F σ /∂θ β ′ . Suppose that the matrix A isinvertible such that X σ =1 A − γσ A σβ ′ = δ γβ ′ . (B5)By multiplying Eq. B4 with A − γσ , taking a sum over σ and re-arranging, we obtain ∂θ γ ∂x α ( k ) = − X σ =1 A − γσ ∂F σ ∂x α ( k ) . (B6)The result given by Eq. B6 is essentially included tothe implicit function theorem . For the second orderderivative, we take the total derivative of Eq. B6 withrespect to x β ( k ′ ) and after some algebra find that ∂ θ γ ∂x α ( k ) ∂x β ( k ′ )= − X σ =1 A − γσ ∂ F σ ∂x α ( k ) ∂x β ( k ′ ) − X β ′ ,σ =1 A − γσ ∂θ β ′ ∂x α ( k ) × X β ′′ =1 ∂ F σ ∂θ β ′ ∂θ β ′′ ∂θ β ′′ ∂x β ( k ′ ) + ∂ F σ ∂θ β ′ ∂x β ( k ′ ) − X β ′′ ,σ =1 A − γσ ∂ F σ ∂x α ( k ) ∂θ β ′′ ∂θ β ′′ ∂x β ( k ′ ) . (B7)All the quantities in Eqs. B6 and B7 are to be evalu-ated at θ = θ ( x ) , u = 0 when applied in solving therelations of Eq. B3 aiming to evaluate R [ θ ( x )]. For in-stance, when the condition F σ ( θ , x , u ) is assumed to be7the Eckart condition given by Eq. 14, we have ∂F σ ∂θ β ′ = N n − X k =1 M k X η,ν =1 ǫ σην x η ( k ) X β =1 ∂ R νβ ∂θ β ′ R ′′ β ( k ) ,∂F σ ∂x α ( k ) = M k X ν =1 ǫ σαν X β =1 R νβ R ′′ β ( k )+ M k X η,ν =1 ǫ σην x η ( k ) R να , (B8)and so on. Here, ǫ σην is the Levi-Civita symbol. Appendix C: Coriolis and vibrational-rotationalcoupling terms
The kinetic energy contributions ˆ T cvr and ˆ T ′ cvr areˆ T cvr = X β,γ ˆ T (1) βγ Z d r ˆ ψ † ( r ) r γ ∂∂r β ˆ ψ ( r )+ N n − X k =1 X α,β,γ ˆ T (2) αβγ ( k ) ˆ p α ( k ) × Z d r ˆ ψ † ( r ) r γ ∂∂r β ˆ ψ ( r )+ X β,σ,γ,γ ′ ˆ T (3) γσγ ′ β Z d r Z d r ′ × ˆ ψ † ( r ) ˆ ψ † ( r ′ ) r γ ′ r ′ γ ∂ ∂r ′ σ ∂r β ˆ ψ ( r ′ ) ˆ ψ ( r ) , ˆ T ′ cvr = X σ,γ ˆ M (1) σγ Z d r ˆ ψ † ( r ) r γ ∂∂r σ ˆ ψ ( r )+ N n − X k =1 X α,σ,γ ˆ M (2) ασγ ( k ) ˆ p α ( k ) × Z d r ˆ ψ † ( r ) r γ ∂∂r σ ˆ ψ ( r )+ X β,σ,γ,γ ′ ˆ M (3) γσγ ′ β Z d r Z d r ′ × ˆ ψ † ( r ) ˆ ψ † ( r ′ ) r γ ′ r ′ γ ∂ ∂r ′ σ ∂r β ˆ ψ ( r ′ ) ˆ ψ ( r ) . (C1)In Eq. C1ˆ T (1) βγ ≡ − N n − X k =1 ¯ h M k X α,β ′ =1 ∂ ˆ R ββ ′ ∂u α ( k ) ˆ R Tβ ′ γ , ˆ T (2) αβγ ( k ) ≡ − i ¯ hM k X β ′ =1 ∂ ˆ R ββ ′ ∂u α ( k ) ˆ R Tβ ′ γ , ˆ T (3) γσγ ′ β ≡ − N n − X k =1 ¯ h M k X α,β ′ ,σ ′ =1 ∂ ˆ R σσ ′ ∂u α ( k ) ˆ R Tσ ′ γ × ∂ ˆ R ββ ′ ∂u α ( k ) ˆ R Tβ ′ γ ′ , (C2) and ˆ M (1) σγ ≡ ¯ h M nuc N n − X k,k ′ =1 3 X α,σ ′ =1 ∂ ˆ R σσ ′ ∂u α ( k ) ∂u α ( k ′ ) ˆ R Tσ ′ γ , ˆ M (2) ασγ ( k ) ≡ i ¯ hM nuc N n − X k ′ =1 3 X σ ′ =1 ∂ ˆ R σσ ′ ∂u α ( k ′ ) ˆ R Tσ ′ γ , ˆ M (3) γσγ ′ β ≡ ¯ h M nuc N n − X k,k ′ =1 3 X α,β ′ ,σ ′ =1 ∂ ˆ R σσ ′ ∂u α ( k ′ ) ˆ R Tσ ′ γ × ∂ ˆ R ββ ′ ∂u α ( k ) ˆ R Tβ ′ γ ′ . (C3)Here we divided the Coriolis and vibrational-rotationalcoupling terms into two parts, the other part ˆ T ′ cvr , is pro-portional to the inverse of the total nuclei mass while theother, ˆ T cvr , is not. All the quantities defined by Eqs. C2and C3 are functions of x and ˆ u . Strictly speaking, thederivatives like ∂ ˆ R σσ ′ /∂u α ( k ′ ), are not so well definedobjects from a notational point of view. By this nota-tion we mean that the quantities like ∂ R σσ ′ /∂u α ( k ′ ) and ∂ R σσ ′ /∂u α ( k ) ∂u α ( k ′ ), are some functions of x , u afterdifferentiation and we can then find the correspondingfunctions of operators in the abstract space. For exam-ple if ∂ R /∂u α ( k ) = G α ( k, x , u ), then in our notation ∂ ˆ R /∂u α ( k ) = G α ( k, x , ˆ u ).We also use the following form of the Coriolis andvibrational-rotational coupling termsˆ T cvr + ˆ T ′ cvr = X c =1 ˆ T ( c ) cvr , (C4)whereˆ T (1) cvr ≡ X β,γ Z d r r γ ∂∂r β ˆ L (1) βγ ˆ n e ( r t ) , ˆ T (2) cvr ≡ N n − X k =1 X α,β,γ Z d r r γ ∂∂r β ˆ L (2) αβγ ( k ) ˆ p α ( k ) ˆ n e ( r t ) , ˆ T (3) cvr ≡ X β,σ,γ,γ ′ Z d r Z d r ′ r γ ′ r ′ γ ∂ ∂r ′ σ ∂r β × ˆ L (3) γσγ ′ β ˆ ψ † ( r t ) ˆ ψ † ( r ′ t ) ˆ ψ ( r ′ t ) ˆ ψ ( r t ) , (C5)and the quantities ˆ L (1) βγ , ˆ L (2) αβγ ( k ) and ˆ L (3) γσγ ′ β are definedby Eq. 28. Appendix D: Useful relations
We expand ˆ n n ( r ), ˆ V en ( r ) and ˆ V ′′ nn to Taylor series indisplacements and write (see Appendix A for the nota-8tion) ˆ V ′ nn = ∞ X m =0 m ! X k ¯ m ,α ¯ m ∂ m V nn ( x ) ∂x α ¯ m ( k ¯ m ) ˆ u α ¯ m ( k ¯ m ) ,V nn ( x ) = N n X k,k ′ =1 ′ v ( x k , x k ′ ) , x N n ≡ − M N n N n − X k =1 M k x k , ˆ V en ( r ) = ∞ X m =0 m ! X k ¯ m ,α ¯ m ∂ m V en ( r , x ) ∂x α ¯ m ( k ¯ m ) ˆ u α ¯ m ( k ¯ m ) ,V en ( r , x ) = N n X k =1 − Z k ς | r − R [ θ ( x )] x k | , (D1)and ˆ n n ( r ) = ∞ X m =0 m ! X k ¯ m ,α ¯ m ∂ m n (0) n ( x , r ) ∂x α ¯ m ( k ¯ m ) ˆ u α ¯ m ( k ¯ m ) ,n (0) n ( x , r ) ≡ − N n X k =1 Z k δ ( r − R [ θ ( x )] x k ) . (D2)In Appendix B we discuss how to actually obtain therotation matrix R [ θ ( x )] and the necessary derivativesof it.Next we set out to derive an exact form of the totalenergy of the system. We start by writingˆ V ′′ nn = ∞ X m =0 ˆ V ( m ) nn , ˆ V en ( r ) = ∞ X m =0 ˆ V ( m ) en ( r ) , (D3)where ˆ V ( m ) nn ≡ m ! X k ¯ m ,α ¯ m ∂ m V nn ( x ) ∂x α ¯ m ( k ¯ m ) ˆ u α ¯ m ( k ¯ m ) , ˆ V ( m ) en ≡ m ! X k ¯ m ,α ¯ m ∂ m V en ( r , x ) ∂x α ¯ m ( k ¯ m ) ˆ u α ¯ m ( k ¯ m ) , (D4)such that ˆ V ′′ en = Z d r ∞ X m =0 ˆ V ( m ) en ( r ) ˆ n e ( r ) . (D5)We note that from these expansions it follows that D ˆ V ′′ en E = − i ¯ h Z d r V en ( r , x ) G (cid:0) r t, r t + (cid:1) − i ¯ h ∞ X m =1 m X k,α Z d r × (cid:28)h ˆ p α ( k ) , ˆ V ( m ) en ( r ) i − ˆ u α ( k ) ˆ n e ( r ) (cid:29) , (D6) and D ˆ V ′′ nn E = V nn ( x ) − i ¯ h ∞ X m =1 m X k,α × (cid:28)h ˆ p α ( k ) , ˆ V ( m ) nn i − ˆ u α ( k ) (cid:29) , (D7)where we used h ˆ n e ( r t ) i = − i ¯ hG ( r t, r t + ). We write thequantity defined by Eq. 47 as K α ( kt ) = K enα ( kt ) + K nnα ( kt ) such that K yα ( kt ) ≡ i ¯ h (cid:28)h ˆ p α ( kt ) , ˆ V ′′ y i − (cid:29) , (D8)where y = en, nn . In a similar way we write the quantitydefined by Eq. 49 as K αβ ( kt, k ′ t ′ ) = K enαβ ( kt, k ′ t ′ ) + K nnαβ ( kt, k ′ t ′ ), where K yαβ ( kt, k ′ t ′ ) ≡ δK yα ( kt ) δJ β ( k ′ t ′ ) . (D9)We defineΠ yαα ′ (cid:0) kt, ¯ k ¯ t (cid:1) ≡ − X k ′ ,β Z dt ′ K yαβ ( kt, k ′ t ′ ) D − βα ′ (cid:0) k ′ t ′ , ¯ k ¯ t (cid:1) . (D10)Next we use Eq. D3 and write K enα ( kt ) and K nnα ( kt )such that K yα ( kt ) = ∞ X m =1 K y, ( m ) α ( kt ) , (D11)and in a similar way for all the quantities derived byusing K enα ( kt ) and K nnα ( kt ). By using Eq. D11 in Eqs.D9 and D10 and these results with Eqs. D6 and D7 weeventually find that D ˆ V ′′ nn E = V nn ( x ) − ∞ X m =1 m X k,α K nn, ( m ) α ( kt ) h ˆ u α ( kt ) i + i ¯ h ∞ X m =1 m X k,α X k ′′ ,α ′ Z dt ′′ Π nn, ( m ) αα ′ ( kt, k ′′ t ′′ ) × D α ′ α ( k ′′ t ′′ , kt ) , (D12)and D ˆ V ′′ en E = − i ¯ h Z d r V en ( r , x ) G (cid:0) r t, r t + (cid:1) − ∞ X m =1 m X k,α K en, ( m ) α ( kt ) h ˆ u α ( kt ) i + i ¯ h ∞ X m =1 m X k,α X k ′′ ,α ′ Z dt ′′ Π en, ( m ) αα ′ ( kt, k ′′ t ′′ ) × D α ′ α ( k ′′ t ′′ , kt ) . (D13)We use these results to write the total energy E tot interms of the Green’s functions and related quantities. Wewrite E tot = D ˆ T tot E + D ˆ V ′′ ee E + D ˆ V ′′ en E + D ˆ V ′′ nn E , (D14)9where the external potentials are put to zero and ˆ T tot isgiven by Eq. 17. We have D ˆ T ′′ e E = i ¯ h m e Z d r ∇ r G (cid:0) r t, r t + (cid:1) , D ˆ T ′′ n E = N n − X k =1 3 X α =1 i ¯ h M k D pp, T αα (cid:0) kt, kt + (cid:1) , D ˆ T ′′ mpn E = − i ¯ h M nuc N n − X k =1 3 X α =1 D pp, T αα (cid:0) kt, k ′ t + (cid:1) . (D15)Here D pp, T αβ ( kt, k ′ t ′ ) ≡ i ¯ h hT { ˆ p α ( kt ) ˆ p β ( k ′ t ′ ) }i , (D16)and this function can be obtained from the solution ofEq. 55 since with the external potentials put to zero h ˆ p α ( kt ) i = h ˆ p α ( k ) i . By using the EOM for the fieldoperator given by Eq. 26 without the external potentialswe find that i ¯ h Z d r (cid:20) − i ¯ h ∂∂t − ¯ h m e ∇ r (cid:21) G (cid:0) r t, r t + (cid:1) = 12 hD ˆ V ′′ en E + D ˆ T (1) cvr E + D ˆ T (2) cvr Ei + D ˆ T (3) cvr E + D ˆ V ′′ ee E + D ˆ T ′′ mpe E , (D17)where ˆ T ( c ) cvr are given by Eqs. C4 and C5. This resultcan be found by multiplying Eq. 26 from the left withˆ ψ † ( r t ) /
2, integrating over r , taking an ensemble averageand then establishing some rearranging. From Eqs. 38and from Eqs. C5 and 89 of Appendix C, it follows that D ˆ T (1) cvr E + D ˆ T (2) cvr E = − i ¯ h Z d r Z d r ′ Z dt ′ Σ ( r t, r ′ t ′ ) × G ( r ′ t ′ , r t ) . (D18) By using Eqs. D12, D13, D17 and D18, the total energycan be written as E tot = i ¯ h N n − X k,k ′ =1 3 X α =1 (cid:0) M − k δ kk ′ − M − nuc (cid:1) D pp, T αα (cid:0) kt, k ′ t + (cid:1) − ∞ X m =1 m X k,α h K en, ( m ) α ( kt ) + 2 K nn, ( m ) α ( kt ) i h ˆ u α ( kt ) i + V nn ( x ) + i ¯ h ∞ X m =1 m X k,α X k ′′ ,α ′ Z dt ′′ × h Π en, ( m ) αα ′ ( kt, k ′′ t ′′ ) + 2Π nn, ( m ) αα ′ ( kt, k ′′ t ′′ ) i × D α ′ α ( k ′′ t ′′ , kt ) − i ¯ h Z d r (cid:20) i ¯ h ∂∂t − ¯ h m ∇ r + V en ( r , x ) (cid:21) G (cid:0) r t, r t + (cid:1) − i ¯ h Z d r Z d r ′ Z dt ′ Σ ( r t, r ′ t ′ ) G ( r ′ t ′ , r t ) . (D19)This is an exact form of E tot . Provided x are the equi-librium positions, then the second term on the righthand side of Eq. D19 involving h ˆ u α ( kt ) i vanishes and D βα ( k ′ t ′ , kt ) = D T βα ( k ′ t ′ , kt ). ACKNOWLEDGMENTS
V. J. H. thanks Dr. Ivan Gonoskov for useful discus-sions on various aspects of the present work. ∗ [email protected] M. Born and R. Oppenheimer,Ann. Phys. (Leipzig) , 457 (1927). K. Huang and M. Born,
Dynamical Theory of Crystal Lat-tices (Clarendon Press Oxford, 1954). P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). R. G. Parr and W. Yang,
Density Functional Theory ofAtoms and Molecules (Oxford University Press, 1989). R. M. Dreizler and E. K. U. Gross,
Density FunctionalTheory: An Approach to the Quantum Many-Body Prob-lem (Springer-Verlag, 1990). P. C. Martin and J. Schwinger,Phys. Rev. , 1342 (1959). G. Baym, Ann. Phys. , 1 (1961). E. G. Maksimov, Zh. Eksp. Teor. Fiz , 2236 (1975). L. Hedin, Phys. Rev. , A796 (1965). G. Mahan,
Many-Particle Physics (Plenum Press, 1990). A. Fetter and J. Walecka,
Quantum Theory of Many-Particle Systems (McGraw-Hill, 1971). G. Stefanucci and R. van Leeuwen,
Nonequilibrium Many-Body Theory of Quantum Systems (Cambridge UniversityPress, 2013). A. A. Maradudin and A. E. Fein,Phys. Rev. , 2589 (1962). R. A. Cowley, Adv. Phys. , 421 (1963). R. Shukla and E. Cowley, Phys. Rev. B , 4055 (1971). G. Rickayzen,
Green’s Functions and Condensed Matter,Acad (Academic Press, 1980). A. A. Maradudin,
Elements of The Theory of Lattice Dy-namics , Vol. 1 (North-Holland Publishing Company, 1974)pp. 1–82. P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni,Phys. Rev. B , 7231 (1991). S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi,Rev. Mod. Phys. , 515 (2001). A. Togo and I. Tanaka, Scr. Mater. , 1 (2015). G. A. S. Ribeiro, L. Paulatto, R. Bianco,I. Errea, F. Mauri, and M. Calandra,Phys. Rev. B , 014306 (2018). S. Baroni, P. Giannozzi, and E. Isaev, Rev. Mineral.Geochem. , 39 (2010). V. J. H¨ark¨onen and A. J. Karttunen, Phys. Rev. B ,024305 (2014). R. Mittal, M. Gupta, and S. Chaplot, Prog. Mater. Sci. , 360 (2018). V. J. H¨ark¨onen and A. J. Karttunen,Phys. Rev. B , 024307 (2016). J. Linnera and A. J. Karttunen,Phys. Rev. B , 014304 (2017). V. J. H¨ark¨onen and A. J. Karttunen,Phys. Rev. B , 054310 (2016). P. Norouzzadeh, J. S. Krasinski, and T. Tadano,Phys. Rev. B , 245201 (2017). H. Euchner, S. Pailh`es, V. M. Giordano, andM. de Boissieu, Phys. Rev. B , 014304 (2018). T. Tadano and S. Tsuneyuki,Phys. Rev. Lett. , 105901 (2018). G. V. Chester, Adv. Phys. , 357 (1961). S. Pisana, M. Lazzeri, C. Casiraghi, K. S. Novoselov, A. K.Geim, A. C. Ferrari, and F. Mauri, Nat. Mater. , 198(2007). J. C. Tully, J. Chem. Phys. , 22A301 (2012). M. E. Casida, B. Natarajan, and T. Deutsch,
Funda-mentals of Time-Dependent Density Functional Theory (Springer, 2012) pp. 279–299. D. Polli, P. Alto`e, O. Weingart, K. M. Spillane, C. Man-zoni, D. Brida, G. Tomasello, G. Orlandi, P. Kukura, R. A.Mathies, et al. , Nature , 440 (2010). L. Kantorovich, Phys. Rev. B , 014307 (2018). M. Baer,
Beyond Born-Oppenheimer: electronic nonadia-batic coupling terms and conical intersections (John Wiley& Sons, 2006). T. Kreibich and E. K. U. Gross,Phys. Rev. Lett. , 2984 (2001). T. Kreibich, R. van Leeuwen, and E. K. U. Gross,Phys. Rev. A , 022501 (2008). M. Marques,
Fundamentals of Time-dependent densityfunctional theory , Vol. 837 (Springer, 2012) p. 249. N. I. Gidopoulos and E. K. U. Gross, Phil. Trans. R. Soc.London , 20130059 (2014). A. Abedi, N. T. Maitra, and E. K. U. Gross,Phys. Rev. Lett. , 123002 (2010). A. Abedi, N. T. Maitra, and E. Gross, J. Chem. Phys. , 22A530 (2012). R. Requist and E. K. U. Gross,Phys. Rev. Lett. , 193001 (2016). N. S. Gillis, Phys. Rev. B , 1872 (1970). L. Hedin and S. Lundqvist,
Solid State Phys. , , 1 (1970). F. Giustino, Rev. Mod. Phys. , 015003 (2017). B. Sutcliffe, Adv. Chem. Phys. , 1 (2000). R. van Leeuwen, Phys. Rev. B , 115110 (2004). R. G. Littlejohn and M. Reinsch,Rev. Mod. Phys. , 213 (1997). B. T. Sutcliffe, Handbook of Molecular Physics and Quan-tum Chemistry , 485 (2003). I. Scivetti, J. Kohanoff, and N. I. Gidopoulos,Phys. Rev. A , 032516 (2009). J. Schwinger, Proc. Nat. Acad. Sci. U. S. A , 452 (1951). C. Eckart, Phys. Rev. , 552 (1935). E. B. Wilson, J. C. Decius, and P. C. Cross,
Molecularvibrations: the theory of infrared and Raman vibrationalspectra (McGraw-Hill Book Company, 1955). A. Sayvetz, J. Chem. Phys. , 383 (1939). A. Maradudin, E. W. Montroll, G. H. Weiss, and I. P.Ipatova,
Theory of The Lattice Dynamics in The HarmonicApproximation , Vol. Supplement 3 (Academic Press, 1971)pp. 6–81. L. V. Keldysh, Sov. Phys. JETP , 1018 (1965). J. Schwinger, Proc. Nat. Acad. Sci. U. S. A , 455 (1951). F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. ,237 (1998). A. Stan, N. E. Dahlen, and R. Van Leeuwen, J. Chem.Phys. , 114105 (2009). C. Rostgaard, K. W. Jacobsen, and K. S. Thygesen,Phys. Rev. B , 085103 (2010). P. Koval, D. Foerster, and D. S´anchez-Portal,Phys. Rev. B , 155417 (2014). G. Lani, P. Romaniello, and L. Reining, New J. Phys. ,013056 (2012). U. M¨uller,
Inorganic Structural Chemistry (John Wiley &Sons, 2006). A. J. Karttunen, T. F. Fassler, M. Linnolahti, and T. A.Pakkanen, Inorg. Chem. , 1733 (2010). H. Hellmann,
Einf¨uhrung in die Quantenchemie (Deuticke,Leipzig) (1937). R. P. Feynman, Phys. Rev. , 340 (1939). G. Baym and N. D. Mermin, J. Math. Phys. , 232 (1961). D. Karlsson and R. van Leeuwen,Phys. Rev. B , 125124 (2016). B. Farid,
Electron Correlation in the Solid State , editedby N. H. March (World Scientific Publishing, 1999) pp.103–261. A. Marini and Y. Pavlyukh,Phys. Rev. B , 075105 (2018). M. Born,
Math. Proc. Cambridge Philos. Soc. , , 160(1940). P. B. Allen, Phys. Rev. B , 2577 (1972). G. Grimvall,
The Electron-Phonon Interaction in Metals (North-Holland, 1981). A. Guzman,
Derivatives and Integrals of MultivariableFunctions (Springer Science & Business Media, 2003). S. G. Krantz and H. R. Parks,
The Implicit Function Theo-rem: History, Theory, and Applications (Springer Science& Business Media, 2012). S. V. Krasnoshchekov, E. V. Isayeva, and N. F. Stepanov,J. Chem. Phys. , 154104 (2014). O. de Oliveira et al. , Real Analysis Exch.39