Many-body Green's function theory of electrons and nuclei beyond Born-Oppenheimer approximation

The method of many body Green's functions is used to describe an arbitrary system of electrons and nuclei in a rigorous manner given the Hamiltonian of Coulombic interactions and kinetic energies. The theory given resolves the problem arising from the translational and rotational invariance of the Hamiltonian afflicting the existing theory based on the same technique. As a result, we derive a coupled set of exact equations for the electron and nuclei Green's functions giving a systematic way to potentially compute various properties of a rather arbitrary many-body systems of electrons and nuclei beyond Born-Oppenheimer approximation, including molecules and solids. We discuss a special case of crystalline solids in more detail.


I. INTRODUCTION
The Born-Oppenheimer (BO) approximation 1,2 is a rather central part in solving the many-body problem of electrons and nuclei in crystalline solids. This approximation has made the study of several physical properties computationally feasible and our current understanding about crystalline solids, for instance, is rather extensively based on it and its validity. In the BO approximation, the electron and nuclei problems are treated separately in a way that the nuclei coordinates are parameters in the electronic problem and the nuclei Hamiltonian contains no operators acting on electrons 2 . Only the BO energy from the electronic problem enters the nuclei Hamiltonian. We can use various approaches, such as the density functional theory [3][4][5][6] or the method of many-body Green's functions [7][8][9][10][11][12][13] , to solve the electronic problem alone. The nuclei part of the problem can be treated, for instance, by using the technique of many-body Green's functions to describe various properties of crystals arising from nuclei 8, [14][15][16][17] . Usually the harmonic part of the nuclei Hamiltonian is diagonalized by a set of coordinate transformations 2,18 and the phonon Green's function theory is written in terms of the transformed coordinates 14 . Therefore we have a closed form exact solution of the harmonic problem given the BO energy and the second order derivatives of it with respect to the nuclei coordinates. In addition to the electronic structure, the BO approximation has been successfully used to compute the phonon spectrum [19][20][21][22] , thermal properties [23][24][25] and thermal conductivities [26][27][28][29][30][31] of various materials in some cases reproducing the experimental results with reasonable accuracy. In some systems, however, the BO approximation may not be sufficiently accurate and more accurate methods are needed. The validity of the BO approximation is weaker for instance in some metals 18,32 and materials such as graphene 33 . Further, a theory beyond the BO approximation is needed in order to describe photovoltaics (solar cells, for example) 34 , chemical reactions 35  There are several first principle methods of going beyond the BO approximation.
Among these are the wave function methods 2,38 , multicomponent density functional theory [39][40][41] , exact factorization [42][43][44][45] , path integrals with Green's functions 37 and the Green's function methods 8,9,[46][47][48] . In the Green's function approach we start with a general many-body Hamiltonian (essentially comprising the Coulombic potentials and kinetic energies) and write the equations of motion (EOM) for the Green's functions involved in order to formally obtain the general solution of the problem and then impose approximations. By starting with a generic Hamiltonian mentioned above, a general theory, in most cases, within the harmonic approximation to the nuclei, has been obtained for the electron and nuclei Green's functions 8,9,[46][47][48] . In order to calculate desired observables within this theory, we may solve a coupled set of equations for the involved Green's functions: Hedin's equations 10 for the electron Green's function and related quantities and in a similar way the equations for the nuclei Green's functions 8,48 . In most cases, one has to impose further approximations in order to solve the equations for the Green's functions involved, but in principle we have a general theory of manybody systems of electrons and nuclei, usually within the harmonic approximation 8,9,48 .
There is, however, a known issue which prevents the use of the existing form of the exact theory to correctly describe solids and molecules, for example. As has been pointed out 39,40,49,50 , the Hamiltonian used in the derivation is not suitable for the description as such. This issue arises from the translational and rotational invariance of the original Hamiltonian we start with and lead to a constant electron and nuclei densities in position space for the eigenstates of the Hamiltonian and spherically symmetric densities for the ground state with vanishing total angular momentum. In order to solve this issue, one can establish a set of coordinate transformations 39,40,49,51 to obtain a Hamiltonian suitable for the description and this have already been done within the many-body Green's function approach 50 . However, no determining equation for the nuclei Green's functions (or for the nuclei densitydensity correlation function) were given 50 in order to find a coupled set of self-consistent equations for an arbitrary system of electrons and nuclei. Moreover, the earlier work was considering only crystalline solids excluding some terms important when the theory is imposed in the study of molecules, for example. The purpose of the present work is to derive such a set of equations with the general Hamiltonian without any approximations established. As a special case, we discuss what is the form of the theory when applied to crystalline solids and compare our results with those obtained earlier.
This paper is organized as follows. In Sec. II, the transformed Hamiltonian suitable for the description of a system of electrons and nuclei is given. The EOM for the electrons are considered in Sec. III and for the nuclei Green's functions in Sec. IV. The Hedin like equations for electrons are derived in Sec. V. We discuss how to determine some parameters related to the coordinate transformations in Sec. VI. The general normal mode frequencies are derived in Sec. VII A and an expression for the nuclei self-energy (SE) due to the Coulombic interactions to the lowest order is given in Sec. VII B and is compared with the earlier results. The special case of crystalline solids is considered in Sec. VIII. Phonons and their 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 The Hamiltonian of this system in position representation can be written as where Here V ee V en and V nn are the Coulombic potentials such that where ς = e 2 / (4πǫ 0 ). The Hamiltonian H given by Eq. 2 is translationally and rotationally invariant and the onebody densities for the eigenstates are constants as a function of position and spherically symmetric for the ground state with vanishing total angular momentum 39,40,50 . In order to derive a Hamiltonian describing correctly the eigenstates of crystalline solids, we establish a coordinate transformation 39,40,49,50,52 and write where R = R (θ) is the rotation matrix, θ = (θ 1 , θ 2 , θ 3 ) is the vector of Euler angles and Euler angles θ α are assumed to be functions of all the N n − 1 nuclei variables denoted by R ′′ . Moreover, R cm is the center-of-mass coordinate of the whole system of electrons and nuclei and R cmn is the nuclei center-of-mass given by We make the coordinate transformation in two stages, first we establish and then the rotation of the electronic coordinates shown in Eq. 5. Such body-fixed frame transformations are used in the BO approximation as well when it is formulated for molecules 1,53 . After transforming the Hamiltonian under the transformation given by Eq. 7, we have where V ′ ee V ′ en and V ′ nn are the Coulombic potentials written in terms of the primed variables and the kinetic energy T tot = T e + T n is of the form After establishing the second transformation r ′′ where V ext is the external potential part added to the Hamiltonian (not originating from the coordinate transformation) and it is added in order to use the functional derivative techniques 54 in deriving the EOM. After the EOM are derived we put V ext to zero. The transformed kinetic energies are where the total mass is M = M nuc + m e N e and the transformed Coulombic potentials are Here, T cvr and T ′ cvr include the Coriolis and rotationalvibrational coupling terms and an explicit form of these quantities is given in Appendix C (given in the abstract space). The center-of-mass kinetic energy T cm commutes with all the other terms in the Hamiltonian and does not enter to the EOM and is neglected from now on. There are only N n − 1 primed nuclei coordinates appearing in the Hamiltonian. However, the number of degrees of freedom is still the same as before since one of the coordinates is the total center-of-mass coordinate of the system. The potential terms V ′′ nn and V ′′ en involving R ′′ Nn break the translational symmetry. The mass-polarization terms T ′′ mpe and T ′′ mpn are expected to be rather small in the case of crystals and larger molecules since they are proportional to the inverse of the total nuclei mass. So far the Euler angles in Eq. 5 are assumed to be generic functions of the nuclei coordinates R ′′ k ′ and we have not introduced any defining relations to them. There are many ways to choose these angles and here we assume that the Euler angles are defined through an implicit equation of the form For example, one possible form of Eq. 13 is the Eckart condition and it can be written as 49-51,55,56 A relation connecting the quantities R ′′ k and x k is given by Eq. 15. The relations given by Eqs. 13 and 14 define the Euler angles and thus the rotation matrix as a function of the N n −1 nuclei variables, R = R [θ (R ′′ )]. Some implicit conditions, like the Eckart condition, changes the permutation symmetry of the Hamiltonian 49,50 with respect to our new nuclei variables leading to a more complicated commutators between the nuclei field operators of identical nuclei and therefore changes the symmetry of the wave functions written in terms of the transformed nuclei coordinates as well. However, the implicit condition does not change canonical commutation relations (Eq. 22) between the nuclei position operators. This means that the EOM for the electron and nuclei Green's functions remain the same whether or not the implicit condition changes the permutation symmetry of the Hamiltonian with respect to nuclei variables. In the case in which the implicit condition changes the aforementioned permutation symmetry and it turns out to be difficult to apply the theory because of the complicated symmetry, we may assume, as in the previous works 8, 48 , that the nuclei are distinguishable. By doing so we avoid the problems potentially caused by a more complicated permutation symmetry. Some implicit conditions are not suitable for describing arbitrary systems, for example, in the case of linear molecules, the Eckart conditions as such are not well defined 55,57 . Thus, despite our theory is in principle general, the choice of F (θ, R ′′ ) in Eq. 13 may restrict the range of validity of the theory.
Next we establish yet another coordinate transformation of the nuclei coordinates and write the position coordinates as a sum of rest positions x k (parameters) and displacements u k (quantum variables), namely We consider the determination of the rest positions x k in more detail in Sec. VI and for now these are taken as arbitrary real parameters. Under this transformation, the position variables in the Hamiltonian H transform as Eq. 15 shows and in the kinetic energy terms (Eq. 11), the derivatives included transform as ∇ R ′′ k → ∇ u k . We have all the necessary steps made in order to write our Hamiltonian. We write the electronic parts of the Hamiltonian in the second quantization. This can be done in a similar way as in the laboratory frame formulation since the transformed Hamiltonian has the same permutation symmetry with respect to the electronic variables as the original one. The nuclei part is written in the first quantization. After using Eq. 15, we write the Hamiltonian operator as (the center-of-mass kinetic energy neglected since it does not enter to the EOM) wherê HereT cvr andT ′ cvr are given in Appendix D and the other kinetic energy terms arê The potetial energy terms in turn arê = drV en (r)n e (r) , andV ext = drU (rt)n n (r) + drU (rt)n e (r) where we usen e (r) =ψ † (r)ψ (r). In these equations v (r, r ′ ) = ς/ |r − r ′ | and The operatorsû α (k),p α ′ (k ′ ),ψ † (r) andψ (r) satisfy the following commutation and anti-commutation relations When necessary (see Sec. VIII), we use the notation k = lκ and so on 2,18,58 , which is perhaps more suitable in the case of crystalline solids, see Appendix A for the description of this notation. In the case of solids, we also impose the periodic boundary conditions 2,58 . So far we have not established any approximations. We have now set up the Hamiltonian to work with and derive the EOM for the electron and nuclei Green's functions by using it in Secs. III and IV.

III. EQUATIONS OF MOTION FOR ELECTRONS
Here we derive the EOM for the electron Green's function. We denote an ensemble average for an operator o (t) (in the Heisenberg picture) as where {|φ n } is a set of orthonormal basis kets. The density operator is of the form whereĤ M =Ĥ − µ eNe ,N e = drn e (r) , and µ e is the chemical potential of the electrons. Thus, the density operator in Eq. 24 is chosen to describe systems with a fixed number of nuclei 8 . Here the notation is such that the time variable t could be taken as a real variable or variable on the Keldysh contour 13,59 . This is justifiable since in both cases the EOM are the same 13 . If the time variables are taken to be variables in the Keldysh contour, the time integrals are taken to be integrals on the contour, see Ref. 13 for details. We start by writing the EOM for the field operator ψ (rt) (operator in the Heisenberg picture), namely wherê The quantities in Eq. 28 originate from the Coriolis and vibrational-rotational coupling terms, see Appendix C.
Define the electron Green's function as 7,10-13 where and G (rt, r ′ t ′ ) satisfies the Kubo-Martin-Schwinger boundary conditions 7,13 . The EOM for the Green's function read After using Eq. 26 we find that where the total density is defined byn (r ′′ t) ≡n e (r ′′ t) + n n (r ′′ t) and By substracting the quantity (we use the shorthand notation given by Eq. A3 of Appendix A) from both sides of Eq. 32 we find that wherẽ All the terms originating from the transformed kinetic energiesT ′′ mpe ,T cvr andT ′ cvr are included in the quantity S (1, 2) given by where S c (1, 2) for each c is given in Sec. V. Next we write the quantities S (1, 2) andS (1, 2) in terms of the corresponding SE's, define (here c = 1, 2, 3) By inverting (38), we find that By using Eq. 39 in Eq. 35 where we defined We can also write Eq. 32 in the form of Dyson equation the function G 0 (1, 2) being a solution of The results obtained in the literature are similar to those obtained here, but here the equations are written in the body-fixed frame in contrast to the laboratory frame formulation 10,48 and the SE also contains the contributions from the kinetic energy terms neglected in the earlier body-fixed frame formulation 50 . By neglecting the kinetic energy termsT ′′ mpe ,T cvr andT ′ cvr in the first place leads to Eqs. 40 and 42 with Σ c (1, 3) = 0 for c = 1, 2, 3 and thus Σ t (1, 3) = Σ (1,3). In this case, we obtain similar equations as in the earlier work on the body-fixed frame 50 . Thus all the rather complicated kinetic energy terms appearing in the theory, originating from the body-fixed frame transformation, are hidden in the SE's Σ c (1, 3). We give an explicit and approximative form of Σ c (1, 3) in Eqs. 94 and 95 of Sec. V. Here we obtained the necessary EOM for electrons and we further write the Hedin like equations for electrons in Sec. V, which may be in some cases useful in deriving expressions for the SE originating from the Coulombic interaction, Σ (1, 2).

IV. EQUATIONS OF MOTION FOR NUCLEI
Next we set out to derive the EOM for nuclei. The connection of the nuclei and electron equations is discussed in more detail in Sec. V where we give some justification for the choice of quantities for which the nuclei EOM are written. We start by writing the Heisenberg EOM for the displacements, namely After computing the commutators by using Eq. 16 [see also Eq. 28] After differentiating Eq. 45 with respect to time and taking an ensemble average where (we useT ′′ cvr ≡T cvr +T ′ cvr ) In the first two quantities of Eq. 47, all the other terms except those originating from the external potentials are explicitly shown.
In order to obtain the EOM for the nuclei Green's function we take a functional derivative of Eq. 45 with respect to J β (k ′ t ′ ) and write where .
We write yet another form of Eq. 48 by assuming that D αβ (kt, k ′ t ′ ) is invertible, namely, Eq. 48 can be written in terms of the nuclei SE's as where . In other words, all the effects of the mass polarization, Coriolis and vibrational-rotational coupling terms on the EOM of the function D αβ (kt, We also write the EOM for the momentumdisplacement and momentum-momentum Green's functions since these functions are in general needed in the calculation of the total energy, see Appendix D. In order to use the functional derivative technique we add the following potential to the external potential part of the Hamiltonian We note that this term would appear in the EOM of û α (kt) , but would not appear in the corresponding equations of D αβ (kt, k ′ t ′ ). We start by writing the EOM for momentum ensemble average, namely By taking a functional derivative of Eq. 54 with respect to P β (k ′ t ′ ) we find that where .
Next we take a functional derivative of Eq. 54 with respect to J β (k ′ t ′ ) and write where Eq. 49 was used and D pu αβ (kt, k ′ t ′ ) ≡ δ p α (kt) /δJ β (k ′ t ′ ). We can use Eq. 50 and write Eq. 57 as From Eq. 58 we see that provided we have the necessary SE's and the solution for D αβ (kt, k ′ t ′ ), then we can obtain D pu αβ (kt, k ′ t ′ ) without solving the EOM for it. So far our derivation is general and we have not made any simplifying assumptions. We still use the full Hamiltonian without expanding in displacementsû, which is the usual procedure in the existing formulations 8,9,48 . We can actually write the exact ensemble average of the Hamiltonian in terms of the quantities appearing in the EOM, see Appendix D. With the Hamiltonian being of general form, we leave the opportunity to use other methods, rather than expanding the Hamiltonian to a Taylor series, in writing the quantities involved such as established in Ref. 46 . We have obtained an exact Green's function theory of an arbitrary system of electrons and nuclei given the Hamiltonian of kinetic energies and Coulombic interactions. With some choices of F (θ, R ′′ ) we may assume the nuclei to be distinguishable, see Sec. II. The full solution is probably rather difficult to obtain and approximations are needed. We discuss a special case of crystalline solids in Sec. VIII and give an explicit approximate 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 rewrite Eq. (32) by using the result δU (r ′′ t) and after using Eq. (61) in Eq. (32) we can write Next we write Eq. (62) in terms of electron SE by using the following result for the electron Green's function where the vertex function is defined as and the inverse of the dielectric function is After using Eqs. 63 and 39 in Eq. 62, we obtain again the EOM given by Eq. 40, but this time the SE, Σ (1, 2), is given by (which is equivalent to the SE appearing in Eq. 38) where the screened Coulombic interaction is defined as Next we set out to derive the remaining Hedin's equations 10,48,50,60,61 in order to obtain determining equations for the quantities like Γ (1, 2, 3) and W (1, 2). We start with the vertex function defined by Eq. 64. We denote and then invert Eq. 40 such that we have By taking a functional derivative of Eq. 69 with respect to V tot (3) we find that where the definition of Γ (1, 2, 3) was used. Next we use and δG (4, 5) δV tot (3) = d6 d7G (4, 6) G (7, 5) Γ (6, 7, 3) , (72) in Eq. 70 and write δG (4, 5) G (4, 6) G (7, 5) Γ (6,7,3) . (73) Next we consider the screened Coulombic interaction given by Eq. 67 written as where the electronic polarization is defined as By making use of the result given by Eq. (72) we find that the electronic polarization satisfies We rewrite now W (1, 2) in terms of electron and nuclei parts. In order to achieve this we write where ∆n n (2) ≡n n (2) − n n (2) , ∆n e (1) ≡n e (1) − n e (1) and the so-called nuclear density-density correlation function is defined as Suppose thatn (79) By using we rewrite Eq. 79 as where we used Finally, after some re-arranging, Eq. 82 can be written as where and Formally, the present theory is similar to the conventional ones 10,48 and thus many of the approximations, like the GW-approximation 10,61-64 or other suitable approximations 65 , can be established in a similar way than in existing Green's function theory for electrons within the BO approximation. Namely, in the GW-approximation Γ (1, 2, 3) ≈ δ (1 − 2) δ (1 − 3) and the electronic polarization and SE become P e (1, 2) ≈ −ihG (1, 2) G (2, 1 + ) and Σ (1, 2) ≈ ihW (1, 2) G (1, 2). The nuclei problem enters to the electronic equations, for instance, through the SE's and the nuclei densitydensity correlation function D n (1, 2). We can obtain an approximate form ofn n (1) and thus of D n (1, 2) by expanding the nuclei densities to a Taylor series in displacements and taking into account the terms of various order. It thus follows that in order to evaluate the expanded n n (1) and D n (1, 2), we need to evaluate ensemble averages like û α (kt)û β (k ′ t) . Similar terms appear if we expand the nuclei terms included to the SE's Σ c (1,2). In order to determine these quantities, we can use the EOM for the nuclei Green's functions derived in Sec. IV. Next we give explicit and approximate expressions for the corresponding SE's arising from the Coriolis and vibrational-rotational coupling terms. The quantities S c (rt, r ′ t ′ ) appearing in Eq. 37 are defined as Perhaps among the simplest approximations to the SE's Σ c (rt, r ′ t ′ ) can be obtained if we write where (92) Similar approximation is established in Sec. VII B when we approximate û α1 (k 1 t)n e (rt) ≈ û α1 (k 1 t) n e (rt) . If we further make the Hartree-Fock approximation 13 we find by using Eq. 38 that and After Taylor expandingD ′ βγ (t) andL (3) γσγ ′ β in displacementsû, Σ 2 (1, 2) and Σ 3 (1, 2) become functions of û α (k) , D T αβ (kt, k ′ t ′ ) and so on. Before continuing we give an overview of a work flow how the coupled set of electrons and nuclei could be solved. First we give an initial guess for the rest positions x as described in Sec. VI. Given the rest positions we expand all the quantities in the Hedin's equations, which are functions of the nuclei variables, to a Taylor series in u. In the first iteration, we retain only the zeroth order terms meaning that W ph (1, 2) = 0 from which follows that W (1, 2) = W e (1, 2) and Σ (1, 2) = Σ e (1, 2) (see Eq. 88). With the preceding choices, the Hedin's equations are actually similar to those in the BO approximation, but here written in the body-fixed frame. The zeroth order approximations for the SE's Σ c (1, 2) can be obtained by using the relations given by Eqs. 38 and 91. After the Hedin's equations are solved, we can calculate the nuclei SE by using the quantities obtained from the solution of the Hedin's equation, see for instance Eq. 115. Given the nuclei SE, the EOM for the nuclei Green's function can be solved. Provided we have obtained the nuclei Green's function, we can write the Hedin's equation with nonzero W ph (1, 2) and solve the Hedin's equation with the SE Σ (1, 2) = Σ e (1, 2) + Σ ph (1, 2). After the new solutions of the electron and nuclei Green's functions have been obtained, we can test whether or not the rest positions x are equilibrium positions (see Sec. VI). We iterate as long as we have found the equilibrium positions. Provided we have found the equilibrium positions, we seek the solution of the coupled equations such that the grand potential, or in the case of zero temperature formalism, the total energy is minimized. In the case when we have several configurations for which the rest positions are equilibrium positions, the aim is to find those equilibrium positions which minimize the grand potential (or the total energy).

VI. DETERMINATION OF REST POSITIONS
We have not yet discussed in detail how to determine the rest positions x and so far these quantities have been considered as arbitrary parameters. Consider, for instance, the total energy of the system, E tot , defined as an ensemble average of the Hamiltonian The value of E tot must be independent of x provided our expansion (in displacementsû) of the SE's, or the Hamiltonian ensemble averege itself, converges for a given x.
Hence, in principle we should have the value of E tot the same for any given and reasonable x. However, we are probably not able to make such an expansion up to arbitrary order inû and thus the value of E tot may become dependent on x due to the approximations established. For instance, if the positions x are chosen far away from the positions around which the nuclei would vibrate, then in order to calculate the reasonable values of E tot , the displacements u would be rather large. This implies that we would need a rather high order Green's functions for nuclei or alternatively a rather high order expression for the nuclei SE to describe the system properly. If we choose the positions x rather poorly and at the same time take only the lowest order quantities inû into account, we expect to find rather unrealistic results. Even though our theory is generally valid and we do not even need to have the nuclei to vibrate around some particular region of space, it is beneficial to have some consistent strategy to obtain the parameters x such that the systems in which the nuclei perform such rather small vibrations around some regions of space we can describe the system with a reasonable accuracy by using the approximations of lower orders inû. For the aforementioned systems, like some stable crystal lattices, we use the following method to determine the positions x. The initial value of the x is obtained in the same way as within the BO approximation, that is, from the experimental data of known structures or in the case of hypothetical structures by using the methods of structural chemistry 66,67 , for instance. After the initial guess, we aim to find the positions x such that they are equal to the most probable positions of the nuclei, that is, we seek the positions such that We note that û k is a function of x since the Hamiltonian is and if the external potentials vanish, R ′′ k is independent of time. If Eq. 97 holds, then it follows that û k = 0. That is, if we are able to choose the positions x as in Eq. 97, then the expected values of the displacements vanish. Next we seek a way how to find the displacements satisfying û k = 0 and thus the rest positions x satisfying Eq. 97. For some systems, like Bravais crystals, this condition maybe automatically satisfied due to symmetry 8 . The determining equation for û k is given by Eq. 46. If we put the external potentials to zero, then û α (kt) = û α (k) and thus meaning that the expected value of the total force on each nuclei k vanish for all α. Thus, if the external potentials vanish, Eq. 46 becomes Without the external potentials, K α (kt) and K (1) α (kt) are of the form given by Eq. 47 with the visible terms taken into account while the rest of the terms K (c) α (kt) remain the same. We write which is a function of û α (k) and F ′′ α (kt) is the remaining part (not a function of û α (k) ). We can therefore write Eq. 99 as This is the determining equation for û α (k) . If we now want to find x such that û α (k) = 0, then Eq. (101) becomes We now give an example such that we find an explicit and approximate form of Eq. 102. Suppose that the only non-zero term in Eq. 99 is K α (kt) meaning that we take into account terms which originate from the Coulombic potentials and neglect all the other terms. This is expected to be relatively good approximation for some crystalline solids provided we define the Euler angles such that the Coriolios and vibrational-rotational terms in the Hamiltonian are small, see Sec. VIII. We approximate K α (kt) by expandingV ′′ en andV ′′ nn in displacements (see Appendix D) and retain only the lowest orders. After calculating the commutators, we obtain to the lowest orders in displacements By using Eq. 113 we find that Here, the quantity n (0) n (x, r) is given by Eq. D2 and related to this we consider how to obtain the rotation matrix R [θ (x)] and its derivatives in Appendix B. When we use Eq. 104 in Eq. 103 in approximating û α (kt)n e (rt) , it can be seen that if we seek the solution such that û k = 0, meaning that F ′ α (kt) vanishes, Eq. 102 is of the following form where (Eq. 49) Since the external potentials vanish, then D T αβ (kt, k ′ t ′ ) can be obtained from Eq. 51 due to the fact that û α (kt) = û α (k) . If we further approximate and neglect the last term of Eq. 105 we obtain Equations 105 and 107 are examples of approximate conditions for choosing the parameters x such that Eq. 97 holds. We call the positions x the equilibrium positions if Eq. 97 holds. With these approximations established, we have formally obtained the same condition (Eq. 107) to choose the parameters than in the BO approximation in the case of crystals 20 , when the Hellmann-Feynman theorem 68,69 is used to calculate the total force. The difference is that our parameters x are in the bodyfixed frame and the electron density n e (r) is an ensemble average taken with respect to the full many-body Hamiltonian instead of an expected value within the BO approximation taken with respect to an eigenstate of the electronic BO Hamiltonian. The electron density can be obtained from the electron Green's function as n e (rt) = −ihG (rt, rt + ). Thus, for a given x one solves the coupled set of equations for the electrons and nuclei and from the solutions of these equations we obtain the quantities like G (rt, r ′ t ′ ), P e (rt, r ′ t ′ ),W e (rt, 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 such that Eq. 97 is satisfied.

A. Nuclei vibrations
We start by writing the nuclei EOM in frequency space, but before Fourier transforming, all the time-dependent potentials are set equal to zero. It follows that the nuclei Green's function and the nuclei SE are then functions of time differences only and we can write Eq. 51 in terms of Fourier transformed quantities as (the convention used is given in Appendix A) We write Eq. 108 in matrix form and re-arrange such thatD The components of these matrices are denoted by αk and so on. Next we write the SE as C t (ω) = C A + C N A (ω) and assume that C A is the Hermitian part of C t (ω), which we assume to be independent of frequency and C N A (ω) the non-Hermitian part of C t (ω). One possible choise is C A = C ′t (0) and C N A (ω) = C t (ω) − C ′t (0) such that we collect all the Hermitian zero frequency contributions of C t (ω) to C ′t (0). For a Hermitian matrix, we can find a complete and orthonormal set of eigenvectors satisfying We call the quantitiesω j the adiabatic normal mode frequencies. The eigenvaluesω 2 j are real since C A is Hermitian andω j are real if C A is positive definite. We transform Eq. 109 by using the eigenvectors of C A such that where If C N A (ω) is sufficiently small, the quasi-particle picture holds and C N A (ω) can be pictured in generating interactions between adiabatic normal modes by shifts of their values and finite lifetimes [imaginary part of C N A (ω)]. Correct treatment requires the continuation of the frequency variable to the complex frequency plane. The usual procedure is to consider the Matsubara Green's functions continued to the complex frequency plane 8,14,17,70,71 . However, we can also obtain the Green's function of complex frequency by just replacing the real frequency ω by the complex one 72 in Eq. 111. We point out that the present discussion is general without any approximations made. We establish a similar consideration 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 of the terms included in K αβ ′ (kt, k ′ t ′ ) given by Eq. 49. We have already calculated the commutators included to the lowest orders inû α (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 Actually, the terms visible in Eq. 112 are all the terms included in K αβ ′ (kt, k ′ t ′ ) if the Hamiltonian is expanded in displacements up to second order before writing the EOM. This is the usual procedure in the earlier laboratory frame formulations 8,48 . Next we use the following result δ n e (rt) δJ β (k ′ t ′ ) = dr 2 dt 2 dr 3 dt 3 P e (rt, whereǫ e (r 1 t 1 , r 2 t 2 ) andW e (r 3 t 3 , r 2 t 2 ) are given by Eqs. 80 and 83, respectively. These quantities can be obtained from the solution of the Hedin's equations. We also write for the quantity in the second term on the right hand side of Eq. 112 Finally, by using Eqs. 112-114 in Eq. 50, we can write for the nuclei SE under consideration in frequency space n e (r) In Eq. 115,W e (r, r ′ , ω) is the Fourier transform of Eq. 83 in the relative time variable. Terms explicitly shown in Eq. 115 are already included in the harmonic approximation, which means that the Hamiltonian is expanded up to second order in displacements (before writing the EOM). Actually not all the terms in the harmonic approximation are even visible in Eq. 115. Namely, the first term on the right hand side of Eq. 114 is not shown and in practise this means that we have approximated û α1 (k 1 t)n e (rt) ≈ û α1 (k 1 t) n e (rt) . Further, only some of the second order terms of the right hand side of Eq. 113 are included. Sometimes it is convenient to write Π αβ (k, k ′ , ω) = Π A αβ (k, k ′ ) + Π N A αβ (k, k ′ , ω). The adiabatic SE Π A αβ (k, k ′ ) = Π αβ (k, k ′ , 0) is real (at least when terms visible in Eq. 115 are taken into account) while the non-adiabatic part Π N A αβ (k, k ′ , ω) = Π αβ (k, k ′ , ω) − Π αβ (k, k ′ , 0) is complex, in general. The adiabatic SE is real sinceW e (rt, r ′ t ′ ) is. In the present case (Eq. 115), the quantity Π nn αβ (kt, k ′ t ′ ) appearing in Eq. D19 is while the rest of the Fourier transforms of the terms visible in Eq. 115 are included to Π en αβ (kt, k ′ t ′ ). For a better comparison with the existing theory, we use Eq. 115 and write the non-adiabatic contribution as where terms including the quantity û α (k) are not explicitly shown. Provided we can choose the parameters x such that Eq. 97 holds, then terms involving the quantities û α (k) vanish, see Sec. VI. The relation given by Eq. 117 is formally analogous with the result obtained earlier 48 . However, the difference is in the densities n (0) n (x, r ′ ) (see Eq. D2) and here all the variables are in the body-fixed frame. We note that if necessary, the higher order expressions for the nuclei and phonon SE's can be obtained by systematically including higher order terms in the expansion of K αβ (kt, k ′ t ′ ). Related to this, some results of a recent study 73 on generic electron-boson Hamiltonians may be useful within the present theory as well.
Even though we are not giving explicit expressions for the SE's Π (c) αα ′ (kt, k ′ t ′ ), these quantities can be obtained with the same procedure as established here for Π αα ′ (kt, k ′ t ′ ). Namely, calculate the commutators contained in K (c) α (kt) (Eq. 47) by expanding the quantities involved inû, take the functional derivative with respect to J β (k ′ t ′ ) in order to obtain K (c) αβ (kt, k ′ t ′ ) and then use Eq. 50 to obtain Π (c) αα ′ (kt, k ′ t ′ ). In obtaining these expressions, similar approximations may be established as we 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 polarization termsT ′′ mpe ,T ′′ mpn and alsoT ′ cvr appearing in the transformed kinetic energy are proportional to the inverse of the total nuclei mass and are thus small for crystals. Further, we assume that one can find an implicit condition of the form (Eq. 13) such that the contributions from the kinetic energyT cvr become small. Some justifications for their neglect have been given 50 , when the implicit condition is chosen to be the Eckart condition given by Eq. 14. By neglecting the aforementioned terms the EOM for electrons remain otherwise the same, except that we have Σ t (1, 3) = Σ (1, 3) in Eqs. 40 and 42. In turn, the nuclei EOM remain otherwise the same, but in Eq. 51

A. Phonons and their interactions
Here we use the notation k = lκ and so on, see Appendix A. The aim is to define phonon frequencies beyond the BO approximation. In the present case, Eq. 108 can be written as We also note that the quantities like D αβ (lκ, l ′ κ ′ , ω) and Π αβ (lκ, l ′ κ ′ , ω) are dependent only on the difference the cell indices l − l ′ . Therefore, we write D αβ (lκ, l ′ κ ′ , ω) as a discrete Fourier transform in the relative coordinate allowed by the periodic boundary conditions 2 , namely and in a similar way for Π αα ′ (κκ ′ , l, ω). In Eq. 119, N is the number of q points and thus number of unit cells in the Born-von Karman cell 2,58 . By using Eq. 119 and similar expression for Π αα ′ (κκ ′ , l ′′ , ω) in Eq. 118 we obtain and after some re-arranging and by using matrix notatioñ where The components of the matrix like C (q, ω) are labeled by ακ and βκ ′ , for instance. We establish a similar procedure as in Sec. VII A, namely, we write the SE as 0). The eigenvalue equation for C A (q) can be written as We call the quantities ω qj the adiabatic phonon frequencies. The eigenvalue equation given by Eq. 123 is analogous to the eigenvalue equation written for the dynamical matrix in the conventional theory of lattice dynamics 2,58 . As the adiabatic phonon frequencies ω qj are defined by Eq. 123, in contrast to the conventional theory, already the non-interacting phonons potentially contain the terms of the Hamiltonian of arbitrary order inû. In order for ω 2 qj to be positive and thus ω qj to be real, the matrix C A (q) have to be positive definite, which is not in general the case for a given x. This is also the case in the BO theory of lattice dynamics 2 , where the positive definiteness of the dynamical matrix implies the minimum of the BO energy surface and the stability of the crystal lattice 74 . Provided the matrix C A (q) is Hermitian (as it is, for example, if Eq. 115 is used), the components of the eigenvector e α (κ|qj) can be chosen to satisfy the orthonormality and completeness conditions We use the eigenvectors of the adiabatic SE to transform Eq. 121 to equation in phonon basis and after we establish the transformation where D (q, ξ) ≡ e † (q)D (q, ξ) e (q) and C N A (q, ξ) ≡ e † (q) C N A (q, ξ) e (q). We call D (q, ξ) the phonon Green's function and C N A (q, ξ) the non-adiabatic phonon SE. Here, ξ is the complex frequency variable 72 .
If the non-adiabatic SE vanishes, we obtain the adiabatic phonon Green's function D A jj (q, ω) = 1/ ω 2 − ω 2 qj which is of the usual form appearing in the BO theory of lattice dynamics 14 . The adiabatic phonons have infinite lifetimes if the non-adiabatic SE vanishes. The nonadiabatic part can be pictured in generating interactions between the adiabatic phonons which appear as shifts to the adiabatic eigenvalues and finite lifetimes of these quasiparticles. This picture is valid if the non-adiabatic SE is sufficiently small in comparison to the adiabatic part. If the non-diagonal terms of C N A (q, ξ) are relatively small such that C N A jj ′ (q, ξ) ≈ 0 for all j = j ′ , Eq. 125 becomes from which the shifts to the adiabatic phonon frequencies and finite lifetimes of adiabatic phonons are usually deduced 48,75,76 . The non-adiabatic phonon SE in its general form can be written as By using Eqs. 80 and 83 in Eq. 117 and then Fourier transforming and changing the representation, we obtain the following approximate form for the non-adiabatic phonon SE In Eq. 128 whereg qj (r) is the complex conjugate ofg * qj (r). By re-arranging Eq. 129 we find that g qj (r, ξ) satisfies the following equation g qj (r, ξ) =g qj (r) + dr ′ dr ′′ g qj (r ′ , ξ) P e (r ′ , r ′′ , ξ) ×v (r ′′ , r) .
The diagrams corresponding to D (q, ξ), C N A (q, ξ) and g qj (r, ξ) are of a similar form as in the laboratory frame formulation and are given in Ref. 48 .

B. Momentum functions
By establishing the approximations discussed at the beginning of this section and then comparing Eqs. 51 and 58 we see that and thus for the Fourier transforms Therefore, after we have obtained a solution for D αβ (lκ, l ′ κ ′ , ω) we can also obtain D pu αβ (lκ, l ′ κ ′ , ω). The last function we need is the momentum Green's function and in the present case Eq. 55 becomes where we approximate which can be obtained directly from Eq. 103. By using the same approximations for K ′ αβ (kt, k ′ t ′ ), as in Sec. VII B were established in writing K αβ (kt, k ′ t ′ ), we find that Eq. 133 can be written as The Fourier transform of this equation is (136) Since by Eq. 60, D up βα (k ′ t ′ , kt) = D pu αβ (kt, k ′ t ′ ), we can write where Eq. 132 was used. Now we have all the necessary equations in order to calculate the total energy of the system, for example.

IX. CONCLUSIONS
In this work, we derived a coupled and self-consistent set of exact equations for the electron and nuclei Green's functions following from the Hamiltonian of Coulombic interactions and kinetic energies as a starting point. The present theory, when applied to crystalline solids, resembles in many parts with the previous ones 8, 9,48 . However, we take into account an issue arising from the translational and rotational symmetry of the Hamiltonian preventing the use of the existing Green's function theories to describe systems other than those with constant density eigenstates. The present theory is formally exact and it is not limited to the harmonic approximation. Moreover, the theory developed can be applied to arbitrary molecules, which in general, were out of the range of validity of the previous theories based on the same method. The complexity of the given system of EOM is of the same order with the existing ones when applied to crystalline solids.
In addition to the general EOM, we considered the normal modes of arbitrary molecular systems, a special case of crystalline solids and reproduced the analogs of some of the results appearing in the existing theories. For instance, phonons and their interactions were defined in a rather general way. While it is probably not a realistic goal to obtain the solution of the EOM in general form in arbitrary systems, there is some work left to derive computationally accessible approximations to be used in the actual calculations. Our main emphasis in this work was on crystalline solids and we leave a more detailed treatment of molecules for future work. The present theory allows one to go beyond BO approximation in a systematic and exact way and therefore it may turn out to be useful in understanding the behaviour of non-relativistic quantum mechanical many-body systems of interacting electrons and nuclei. this notation, we write R ′′ k = R ′′ lκ = x lκ + u lκ , where lκ go over N n − 1 values in total. In the case of the quantities x k and u k , we sometimes use the following collective notation x ≡ x 1 , . . . , x Nn−1 and u ≡ u 1 , . . . , u Nn−1 .
We use the following shorthand notation for the position and time variables In this work, we use the following convention for the Fourier transforms Appendix B: Euler angles The condition given by Eq. 13 defines the Euler angles by an implicit equation of the form F (θ, x + u) = 0. In general, no explicit solution of such equations exists 77-79 , however. That is, it might not be possible to write θ as an explicit function of the nuclei variables x, u, but for a given x, u and the nuclei masses, a numerical solution could be sometimes obtained 79 . We are in particular interested how to obtain the rotation matrix R [θ (x)] and its derivatives since these quantities appear in the Hamiltonian if we establish the Taylor expansions of some parts of the Hamiltonian (see Appendix D) and thus they will eventually appear in the EOM.
The rotation matrix R [θ (x + u)] can be obtained from 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 For instance, in the case of Eckart conditions given by Eq. 14 with the Euler angles θ (x), such that 0 < θ β (x) < 2π for β = 1, 2, 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) = 3 β=1 ∂R ∂θ β ∂θ β ∂x α (k) , . (B3) Given the explicit form for the rotation matrix as a function of the Euler angles we can calculate the derivatives like ∂R/∂θ β ′ so next we find a way how to calculate the derivatives like ∂θ β ′ /∂x α (k) and ∂ 2 θ β ′ /∂x α (k)∂x β (k ′ ). After taking the total derivative of F σ (θ, x, u) with respect to x α (k) we write where A σβ ′ ≡ ∂F σ /∂θ β ′ . Suppose that the matrix A is invertible such that By multiplying Eq. B4 with A −1 γσ , taking a sum over σ and re-arranging, we obtain .
The result given by Eq. B6 is essentially included to the implicit function theorem 80 . For the second order derivative, we take the total derivative of Eq. B6 with respect to x β (k ′ ) and after some algebra find that All the quantities in Eqs. B6 and B7 are to be evaluated at θ = θ (x) , u = 0 when applied in solving the relations of Eq. B3 aiming to evaluate R [θ (x)]. For instance, when the condition F σ (θ, x, u) is assumed to be the Eckart condition given by Eq. 14, we have and so on. Here, ǫ σην is the Levi-Civita symbol. n (x, r) ∂x αm (km)û αm (km) , In Appendix B we discuss how to actually obtain the rotation matrix R [θ (x)] and the necessary derivatives of it. Next we set out to derive an exact form of the total energy of the system. We start by writinĝ and where we used n e (rt) = −ihG (rt, rt + ). We write the quantity defined by Eq. 47 as K α (kt) = K en α (kt) + K nn α (kt) such that where y = en, nn. In a similar way we write the quantity defined by Eq. 49 as K αβ (kt, k ′ t ′ ) = K en αβ (kt, k ′ t ′ ) + K nn αβ (kt, k ′ t ′ ), where .
We define (D10) Next we use Eq. D3 and write K en α (kt) and K nn α (kt) such that and in a similar way for all the quantities derived by using K en α (kt) and K nn α (kt). By using Eq. D11 in Eqs. D9 and D10 and these results with Eqs. D6 and D7 we eventually find that ×D α ′ α (k ′′ t ′′ , kt) .
We use these results to write the total energy E tot in terms of the Green's functions and related quantities. We write