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

TL;DR

This paper develops a comprehensive many-body Green's function framework for electrons and nuclei that goes beyond the Born-Oppenheimer approximation, addressing invariance issues and enabling systematic property calculations.

Contribution

It introduces a coupled set of exact equations for electronic and nuclear Green's functions applicable to arbitrary systems beyond the Born-Oppenheimer approximation.

Findings

Addresses translational and rotational invariance problems.

Provides systematic approximation methods.

Applicable to crystalline solids.

Abstract

The method of many-body Green's functions is developed for arbitrary systems of electrons and nuclei starting from the full (beyond Born-Oppenheimer) Hamiltonian of Coulomb interactions and kinetic energies. The theory presented here resolves the problems arising from the translational and rotational invariance of this Hamiltonian that afflict the existing many-body Green's function theories. We derive a coupled set of exact equations for the electronic and nuclear Green's functions and provide a systematic way to approximately compute the properties of arbitrary many-body systems of electrons and nuclei beyond the Born-Oppenheimer approximation. The case of crystalline solids is discussed in detail.

Equations

$$H=T_n+T_e+V_{ee}+V_{en}+V_{nn},$$

$$T_n=-\sum _{k=1}^{N_n}\frac{{\hslash }^2}{2M_k}{\nabla }_{{\mathbf{R}}_k}^2,T_e=-\frac{{\hslash }^2}{2m_e}\sum _{i=1}^{N_e}{\nabla }_{{\mathbf{r}}_i}^2,$$

$$V_{ee}$$

$$V_{nn}$$

$$v_{en}\left({\mathbf{r}}_i,{\mathbf{R}}_k\right)$$

$$v_e\left({\mathbf{r}}_i,{\mathbf{r}}_{i^{\prime }}\right)$$

$$v_n\left({\mathbf{R}}_k,{\mathbf{R}}_{k^{\prime }}\right)$$

$${\mathbf{r}}_i,i=1,\dots ,N_e,{\mathbf{R}}_k,k=1,\dots ,N_n,$$

$$H{\Psi }_s=E_s{\Psi }_s,$$

$${\Psi }_s=e^{i{\mathbf{P}}_s\cdot {\mathbf{R}}_{cm}}{\Phi }_s\left({\tilde{\mathbf{r}}}_2,\dots ,{\tilde{\mathbf{r}}}_{N_e},{\tilde{\mathbf{R}}}_1,\dots ,{\tilde{\mathbf{R}}}_{N_n}\right),$$

$${\tilde{\mathbf{r}}}_i$$

$${\tilde{\mathbf{R}}}_k$$

$$n_s\left({\mathbf{r}}_1\right)=\int d{\mathbf{r}}_2\cdots \int d{\mathbf{r}}_{N_e}\int d\mathbf{R}{\left|{\Psi }_s\left({\mathbf{r}}_1,\dots ,{\mathbf{r}}_{N_e},\mathbf{R}\right)\right|}^2,$$

$$G_{gs}\left(\mathbf{r}t,{\mathbf{r}}^{\prime }{t}^{\prime }\right)=-\frac{i}{\hslash }⟨{\Psi }_0⟩\mathcal{T}\left\{\hat{\psi }\left(\mathbf{r}t\right){\hat{\psi }}^{†}\left({\mathbf{r}}^{\prime }{t}^{\prime }\right)\right\}{\Psi }_0.$$

$${\Psi }_s={\Phi }_{cm}\left({\mathbf{R}}_{cm}\right){\Phi }_s\left({\tilde{\mathbf{r}}}_2,\dots ,{\tilde{\mathbf{r}}}_{N_e},{\tilde{\mathbf{R}}}_1,\dots ,{\tilde{\mathbf{R}}}_{N_n}\right),$$

$$-\frac{{\hslash }^2}{2M}{\nabla }_{{\mathbf{R}}_{cm}}^2{\Phi }_{cm}\left({\mathbf{R}}_{cm}\right)=E_{cm}{\Phi }_{cm}\left({\mathbf{R}}_{cm}\right),$$

$${\Phi }_{cm}\left({\mathbf{R}}_{cm}\right)=j_l\left(R_{cm}\right)Y_{lm}\left(\vartheta ,\phi \right),$$

$${\mathbf{r}}_i^{\prime }$$

$${\mathbf{R}}_k^{\prime }$$

$${\mathbf{R}}_{N_n}^{\prime \prime }$$

$${\mathbf{R}}_{cmn}=\frac{1}{M_{nuc}}\sum _{k=1}^{N_n}{M}_k{\mathbf{R}}_k,M_{nuc}\equiv \sum _{k=1}^{N_n}{M}_k.$$

$$H$$

$$T_{cm}$$

$$T_e^{\prime }$$

$$T_n^{\prime }$$

$$T_{mpe}^{\prime }$$

$$T_{mpn}^{\prime }$$

$$V_{ee}^{\prime }$$

$$V_{nn}^{\prime }$$

$$V_{en}^{\prime }$$

$${\mathbf{R}}_{N_n}^{\prime }$$

$${\mathbf{R}}_k^{\prime }={\mathbf{x}}_k+{\mathbf{u}}_k,k=1,\dots ,N_n-1.$$

$$\hat{H}={\hat{T}}_{tot}+{\hat{V}}_{ee}^{\prime }+{\hat{V}}_{en}^{\prime }+{\hat{V}}_{nn}^{\prime }+{\hat{V}}_{ext},$$

$${\hat{T}}_{tot}={\hat{T}}_e^{\prime }+{\hat{T}}_n^{\prime }+{\hat{T}}_{mpe}^{\prime }+{\hat{T}}_{mpn}^{\prime }+{\hat{T}}_{cvr}+{\hat{T}}_{cvr}^{\prime }.$$

$${\hat{T}}_n^{\prime }$$

$${\hat{T}}_{mpn}^{\prime }$$

$${\hat{T}}_e^{\prime }$$

$${\hat{T}}_{mpe}^{\prime }$$

Full text

Many-body Green’s function, Born-Oppenheimer approximation, Lattice dynamics

I Introduction

The Born-Oppenheimer (BO) approximation Born and Oppenheimer (1927); Huang and Born (1954) is among the most fundamental ingredients of modern condensed-matter theory. Much of our current understanding of molecules and crystalline solids is heavily based on this approximation and its validity. The BO approximation not only makes calculations computationally feasible, the motion of nuclear wave packets in the lowest BO potential energy surface also provides us with an intuitive picture of many chemical reactions.

The numerical efficiency of the BO approximation comes from treating the electronic and nuclear problems separately in such a way that the nuclear coordinates only enter as parameters in the electronic many-body problem while the nuclear Hamiltonian contains, as scalar potential, the total electronic ground-state energy as a function of the nuclear coordinates. Then the total wave function of the system is a single product of a nuclear wave function (typically a vibrational state) and the many-electron wave function which parametrically depends on the nuclear coordinates.

To deal with the electronic-structure problem at fixed nuclear configuration, various approaches are available, such as density functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965); Parr and Yang (1989); Dreizler and Gross (1990) or the method of many-body Green’s functions Martin and Schwinger (1959); Baym (1961); Maksimov (1975); Hedin (1965); Mahan (1990); Fetter and Walecka (1971); Stefanucci and van Leeuwen (2013); Golze et al. (2019). The nuclear many-body Schrödinger equation is usually treated differently: One first rewrites the nuclear Hamiltonian in terms of normal coordinates which represent the nuclear vibrational and rotational degrees of freedom. In terms of these coordinates, the BO surface (i.e. the scalar potential in the nuclear Schrödinger equation) is expanded to second order around the equilibrium positions, so that the solution can be written analytically as a product of harmonic-oscillator wave functions (phonons). Higher orders of the expansion appear as phonon-phonon interactions in the nuclear Schrödinger equation and can be dealt with using bosonic many-body Green’s functions Baym (1961); Maradudin and Fein (1962); Cowley (1963); Shukla and Cowley (1971); Maradudin (1974); Rickayzen (1980). With this strategy, phonon spectra Giannozzi et al. (1991); Baroni et al. (2001); Togo and Tanaka (2015); Ribeiro et al. (2018), thermal conductivities and other thermal properties Baroni et al. (2010); Härkönen and Karttunen (2014); Mittal et al. (2018); Härkönen and Karttunen (2016a); Linnera and Karttunen (2017); Härkönen and Karttunen (2016b); Norouzzadeh et al. (2017); Euchner et al. (2018); Tadano and Tsuneyuki (2018); Shen et al. (2020) have been successfully calculated. In most cases, especially in solids, the results of such BO-based calculations are in very good agreement with experiment. Somewhat surprisingly, this is also true for metals. As there is no gap in metals, there is no clear separation of electronic and nuclear energy scales in the low-lying excitations, so the naive expectation would have been that the validity of the BO approximation is questionable. Still, BO-based calculations of phononic properties of metals are usually very successful. There are only a few notable exceptions such as $\text{MgB}_{2}$ Calandra et al. (2007), graphene Pisana et al. (2007) and, possibly, the recently discovered lanthanum hydride high-temperature superconductors Somayazulu et al. (2019); Flores-Livas and Arita (2019); Errea et al. (2020).

In spite of this tremendous success of the BO approximation, a number of fascinating phenomena in physics and chemistry appear in the so-called non-adiabatic regime where the coupled motion of electrons and nuclei beyond the Born-Oppenheimer approximation is essential. Examples are the process of vision Polli et al. (2010), exciton dynamics in photovoltaic systems Tully (2012); Nelson et al. (2014), the splitting of the nuclear wave packet in the Zewail Nobel prize experiments Dantus et al. (1987); Potter et al. (1992); Zewail (2000) and the occurrence of local electronic currents associated with nuclear motion Barth et al. (2009); Schild et al. (2016); Scherrer et al. (2017).

To go beyond the BO approximation is notoriously difficult because one has to start from scratch, that is from the complete many-body Hamiltonian of interacting electrons and nuclei. In the non-relativistic regime, this Hamiltonian is given by:

[TABLE]

where, in position representation, the kinetic energies are

[TABLE]

while the potential energy contributions are of the form

[TABLE]

Here the primed sums are over the values $i\neq i^{\prime}$, $k\neq k^{\prime}$ and $v_{en}$, $v_{e}$ and $v_{n}$ are the bare Coulomb interactions

[TABLE]

where $\varsigma=e^{2}/\left(4\pi\epsilon_{0}\right)$. The system consists of $N_{e}$ electrons and $N_{n}$ nuclei. The position coordinates of these particles

[TABLE]

as they appear in Eqs. 1-4, refer to the laboratory frame. The Hamiltonian $H$ given by Eq. 1 is translationally and rotationally invariant. This has a number of important consequences which will be discussed in the following. The eigenstates in position representation of the Hamiltonian $H$ are the solutions of the stationary Schrödinger equation

[TABLE]

where $s$ labels the eigenstates. The translational invariance of the Hamiltonian implies that it commutes with the total momentum operator. Hence the eigenstates $\Psi_{s}=\Psi_{s}\left(\mathbf{r}_{1},\ldots,\mathbf{r}_{N_{e}},\mathbf{R}_{1},\ldots,\mathbf{R}_{N_{n}}\right)$ can be chosen to be simultaneous eigenstates of the Hamiltonian and the total momentum operator:

[TABLE]

where $\mathbf{R}_{cm}$ is the center-of-mass position vector of the complete system, and $\mathbf{P}_{s}$ is the total momentum of state $\Psi_{s}$. The coordinates $\tilde{\mathbf{r}}_{i}$ and $\tilde{\mathbf{R}}_{k}$ are given by:

[TABLE]

We now calculate the standard probability density $n_{s}\left(\mathbf{r}_{1}\right)$ of finding an electron at point $\mathbf{r}_{1}$ (in the laboratory frame) if the system is in the eigenstate $\Psi_{s}$:

[TABLE]

where we denote all the nuclear variables by $\mathbf{R}$. Plugging Eq. 7 into Eq. 9 and substituting the integration variables $\mathbf{r}_{2},\ldots,\mathbf{r}_{N_{e}},\mathbf{R}_{1},\ldots,\mathbf{R}_{N_{n}}$ by the variables given by Eq. 8, we immediately realize that the density $n_{s}\left(\mathbf{r}_{1}\right)$ does not depend on $\mathbf{r}_{1}$. This is true for any atom, molecule or solid and is perfectly reasonable: The Hamiltonian $H$ does not localize the system anywhere in space and hence finding an electron is equally likely everywhere. A similar feature appears for the electronic one-body Green’s function when defined with respect to the ground state $\Psi_{0}$ of the Hamiltonian $H$:

[TABLE]

Here $\mathcal{T}$ is the usual time-ordering operator and $\hat{\psi}^{\dagger}\left(\mathbf{r}^{\prime}t^{\prime}\right),\hat{\psi}\left(\mathbf{r}t\right)$ are electron creation and annihilation operators in the Heisenberg picture, the latter referring to the full Hamiltonian given by Eq. 1. It is easy to verify that the so-defined Green’s function $G_{gs}\left(\mathbf{r}t,\mathbf{r}^{\prime}t^{\prime}\right)$ only depends on $\mathbf{r}-\mathbf{r}^{\prime}$. Again this feature holds true for all atoms, molecules and solids. This property of the Green’s function, and likewise the constancy of the probability density, are consequences of the translational invariance of the underlying Hamiltonian $H$ and are, as such, not at all surprising. If one wants to develop a density functional theory or a Green’s-function-based many-body theory for the complete system of electrons and nuclei then, clearly, the probability density of Eq. 9 and the Green’s function of Eq. 10 are not useful quantities because they do not reveal any features characteristic of the system, for example the molecule, to be described. We need to find other densities or Green’s functions on which the theory can be built, densities or Green’s functions that reflect the internal features of the molecule or crystal.

Eigenfunctions do not necessarily have the same symmetry as the Hamiltonian they come from, they may have lower symmetry. In our case, having the eigenstates $\Psi_{s}$ to be simultaneous eigenstates of the total momentum was a choice we made, and it is exactly this choice that produced the constant density. However, it is important to realize that we can always separate off the center-of-mass motion, and hence the eigenstates of the Hamiltonian given by Eq. 1 can always be written as

[TABLE]

where $\Phi_{cm}$ satisfies the equation

[TABLE]

where $M$ is the total mass of electrons and nuclei. Equation 7 was a particular solution of Eq. 12 that preserves the translational symmetry of the Hamiltonian. But other solutions can be chosen, for instance angular-momentum eigenstates (i.e. partial waves) of the form

[TABLE]

as well as linear combinations of these functions with the same absolute total momentum $P_{s}$, as these span the subspace of degenerate eigenfunctions with energy $E_{cm}=P^{2}_{cm}/2M$. Clearly, an eigenstate of the form 11 with $\Phi_{cm}$ being of the form 13 will not have a constant probability density $n_{s}\left(\mathbf{r}_{1}\right)$. There will be a genuine $\mathbf{r}_{1}$-dependence. However, this $\mathbf{r}_{1}$-dependence reflects the shape of the angular-momentum eigenstate 13 of the center of mass but, again, will not be characteristic of the internal structure of the molecule. So, for building a density functional theory or Green’s function theory, such choice will not be helpful. Alternatively, one might consider the possibility of localizing the whole system with an external confining potential. While this will avoid a constant density, the shape of $n_{s}\left(\mathbf{r}_{1}\right)$ will reflect the shape of the confining potential but will not reveal the internal structure of the molecule. So, once again, this choice is not suitable for constructing a density functional or Green’s function approach to the complete system of electrons and nuclei.

There are essentially two distinctly different choices one can make for the density or Green’s function that are useful in constructing a density-functional or Green’s function approach: One either refers the coordinates on which the density/Green’s function depends to a body-fixed coordinate frame, or one works with conditional probabilities. The second option, i.e. the choice of conditional probabilities, is taken in the BO framework: The standard BO-based Green’s function for electrons is defined in terms of an electronic many-body wave function $\Phi_{\mathbf{R}}\left(\mathbf{r}_{1},\ldots,\mathbf{r}_{N_{e}}\right)$ which is a conditional probability amplitude for finding the electrons at $\mathbf{r}_{1},\ldots,\mathbf{r}_{N_{e}}$, given the nuclei are at $\mathbf{R}$. Likewise, the standard Hohenberg-Kohn-Sham ground-state density functional theory is based on the conditional probability density $n_{\mathbf{R}}\left(\mathbf{r}_{1}\right)$ of finding an electron at $\mathbf{r}_{1}$, given the nuclei are at $\mathbf{R}$. Using the framework of the exact factorization Abedi et al. (2010, 2012), a density functional framework beyond the BO limit has been developed Gidopoulos and Gross (2014); Requist and Gross (2016); Li et al. (2018); Requist et al. (2019) on the basis of the exact (rather than BO) conditional density. A Green’s function approach based on the exact factorization has not been attempted so far.

In this article we follow the first option: Our goal is to develop a Green’s function-based many-body theory for the complete system of electrons and nuclei where the Green’s functions are defined in terms of coordinates that refer to a body-fixed coordinate frame. First steps in this direction were taken with the formulation of a multicomponent density-functional theory Kreibich and Gross (2001); Kreibich et al. (2008); Butriy et al. (2007) and with the derivation of a Green’s-function framework van Leeuwen (2004). In the latter, the electronic degrees of freedom are treated in a fully consistent way. However, no determining equation for the nuclear Green’s function was derived. A nuclear density-density correlation function must be imported from outside the theory. Furthermore, the earlier work considered only crystalline solids excluding some important terms when the theory is imposed in the study of molecules, for example. The purpose of the present article is to deduce, from the general Hamiltonian of Eq. 1, a coupled set of self-consistent equations for the electronic and nuclear Green’s functions that does not require any input from outside.

This paper is organized as follows. In Sec. II we discuss in detail the coordinate transformation to the body-fixed frame and how the Hamiltonian looks when expressed in terms of these coordinates. The equation of motion (EOM) for the electronic Green’s functions are presented in Sec III, and for the nuclear Green’s function is Sec. IV. Hedin-like equations for the electrons are deduced 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.1 and an expression for the nuclear self-energy (SE) due to the Coulombic interactions to the lowest order is given in Sec. VII.2. The special case of crystalline solids is considered in Sec. VIII. Phonons and their interactions are discussed in Sec. VIII.1.

II Coordinate transformation and the transformed Hamiltonian

II.1 Hamiltonian

As demonstrated above, the density $n_{s}\left(\mathbf{r}_{1}\right)$ is constant for the wave function given by Eq. 7 if $\mathbf{r}_{1}$ denotes a coordinate vector in the laboratory frame. By contrast, the density $n_{s}\left(\tilde{\mathbf{r}}_{1}\right)$ with $\tilde{\mathbf{r}}_{1}$ being the position measured from the center of mass is not a constant and reveals to some extent the internal structure of the system. However, the density $n_{s}\left(\tilde{\mathbf{r}}_{1}\right)$ is not good enough yet because it is spherical for all atoms, molecules and solids. Clearly, we still need to deal with the rotational degrees of freedom. This is achieved by the following transformation Sutcliffe (2000, 2003); Kreibich and Gross (2001); van Leeuwen (2004); Kreibich et al. (2008):

[TABLE]

where $\mathcal{R}=\mathcal{R}\left(\boldsymbol{\theta}\right)$ is a rotation matrix and $\boldsymbol{\theta}=\left(\theta_{1},\theta_{2},\theta_{3}\right)$ is the associated vector of Euler angles which, for now, are assumed to be functions of all the $N_{n}-1$ nuclear variables denoted by $\mathbf{R}^{\prime}$. We used $\mathbf{R}^{\prime\prime}_{N_{n}}$ for one of the transformed coordinates in order to keep our notation consistent in writing the transformed Coulomb potential terms. Furthermore, $\mathbf{R}_{cmn}$ is the nuclear center of mass given by

[TABLE]

We can think of the coordinate transformation given by Eq. 14 as two sequential steps. First, the coordinates are written relative to the nuclear center of mass and the center of mass of the total system is separated. Secondly, the rotation of the electronic coordinates shown in Eq. 14 is established to fix an orientation of the electronic subsystem relative to the nuclear center-of-mass frame. Such body-fixed-frame transformations are also used in the BO approximation when it is formulated for molecules Born and Oppenheimer (1927); Scivetti et al. (2009). In terms of the new coordinates given by Eq. 14, the Hamiltonian reads

[TABLE]

Here $V_{ext}$ is an external potential added to the Hamiltonian (not originating from the coordinate transformation) which is introduced in order to use the functional derivative techniques Schwinger (1951) when we derive the EOM. We will specify its explicit form later. After the EOM are obtained we put $V_{ext}\equiv 0$. The transformed kinetic energies in Eq. 16 are

[TABLE]

where the total mass is $M=M_{nuc}+m_{e}N_{e}$ and the transformed Coulomb potential terms are

[TABLE]

In Eq. 16, $T_{cvr}$ and $T^{\prime}_{cvr}$ include the Coriolis and rotational-vibrational coupling terms. The explicit form of these quantities is given in Appendix A. 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, hence, can be disregarded without loss of generality. There are only $N_{n}-1$ primed nuclear 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^{\prime}_{nn}$ and $V^{\prime}_{en}$ involving $\mathbf{R}^{\prime}_{N_{n}}$ are no longer translation invariant in the new coordinates. The mass-polarization terms $T^{\prime}_{mpe}$ and $T^{\prime}_{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 nuclear mass.

Next we perform yet another transformation of the nuclear coordinates and write the position coordinates as a sum of reference positions $\mathbf{x}_{k}$ (which act as parameters) and displacements $\mathbf{u}_{k}$ (quantum variables), namely

[TABLE]

We will consider the determination of the reference positions $\mathbf{x}_{k}$ in detail in Sec. VI. For now these are taken as arbitrary real parameters. Depending on the situation, we use the following notations interchangeably $\mathbf{u}_{k}=\mathbf{u}\left(k\right)$ for the displacements and, analogously, for $\mathbf{R}^{\prime}_{k}$ and $\mathbf{x}_{k}$. This allows us to denote the $\alpha$th Cartesian components of $\mathbf{x}_{k}$ and $\mathbf{u}_{k}$ by $x_{\alpha}\left(k\right)$ and $u_{\alpha}\left(k\right)$. Further, we sometimes denote all the $N_{n}-1$ nuclear reference positions and displacements by $\mathbf{x}$ and $\mathbf{u}$, respectively. We emphasize that the theory is still exact after the transformation given by Eq. 19. The position variables in the Hamiltonian $H$ transform as Eq. 19 shows and the derivatives in the kinetic energy terms (Eq. 17) transform as $\nabla_{\mathbf{R}^{\prime}_{k}}\rightarrow\nabla_{\mathbf{u}_{k}}$. We now have all the necessary prerequisites in place to write our Hamiltonian. We choose to write the electronic parts of the Hamiltonian in 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. With respect to the nuclear coordinates $\mathbf{R}^{\prime}$, the transformed Hamiltonian does not necessarily have the same permutation symmetry as the original one (for identical nuclei). Therefore, the nuclear part is kept in first quantization. Note that the transformed Hamiltonian only contains $N_{n}-1$ independent nuclear coordinate and the total (electron-nuclear) center of mass. After using Eq. 19 and ignoring the center-of-mass kinetic energy, we obtain the final form of the Hamiltonian

[TABLE]

where

[TABLE]

Here $\hat{T}_{cvr}$ and $\hat{T}^{\prime}_{cvr}$ are given in Appendix A and the other kinetic energy terms are

[TABLE]

The energy terms associated with the particle-particle interactions are

[TABLE]

and the external potential part is

[TABLE]

where we use $\hat{n}_{e}\left(\mathbf{r}\right)=\hat{\psi}^{\dagger}\left(\mathbf{r}\right)\hat{\psi}\left(\mathbf{r}\right)$. Furthermore, the quantities $U\left(\mathbf{r}t\right)$, $U^{\prime}\left(\mathbf{r}t\right)$, $\mathbf{J}\left(kt\right)$ and $\varphi\left(\mathbf{r}t\right)$ are time-dependent external potentials. In these equations $v\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\varsigma/\left|\mathbf{r}-\mathbf{r}^{\prime}\right|$ and

[TABLE]

The operators $\hat{u}_{\alpha}\left(k\right)$, $\hat{p}_{\alpha^{\prime}}\left(k^{\prime}\right)$, $\hat{\psi}^{\dagger}\left(\mathbf{r}\right)$ and $\hat{\psi}\left(\mathbf{r}\right)$ satisfy the following commutation and anti-commutation relations

[TABLE]

So far the Euler angles in Eq. 14 were assumed to be generic functions of the nuclear coordinates $\mathbf{R}^{\prime}$. We have not introduced any defining relations for them. There are many ways to choose these angles. Without loss of generality we assume that the Euler angles are defined through an implicit equation of the form

[TABLE]

For example, one possible form of Eq. 27 is the Eckart condition which can be written as Eckart (1935); Wilson et al. (1955); Littlejohn and Reinsch (1997); Sutcliffe (2000); van Leeuwen (2004)

[TABLE]

The relations given by Eqs. 27 and 28 define the Euler angles and thus the rotation matrix as a function of the $N_{n}-1$ nuclear variables, $\mathcal{R}=\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{R}^{\prime}\right)\right]$. Some implicit conditions, like the Eckart condition, can change the permutation symmetry present in the original Hamiltonian Sutcliffe (2000); van Leeuwen (2004) with respect to the old $\mathbf{R}$ and new nuclear variables $\mathbf{R}^{\prime}$ leading to possible complications when (anti)symmetrizing the nuclear subsystem(s) properly. However, the implicit condition does not change the canonical commutation relations (Eq. 26) between the nuclear position and momentum operators. This means that the EOM for the electron and nuclear Green’s functions remain the same whether or not the implicit condition changes the permutation symmetry of the Hamiltonian with respect to the nuclear variables. In the case in which the implicit condition changes the aforementioned permutation symmetries and it turns out to be difficult to apply the theory because of the complicated symmetry, we may assume, as in the previous works Baym (1961); Giustino (2017), that the nuclei are distinguishable. By doing so we avoid the problems potentially caused by a more complicated permutation symmetry. There are few physical effects in which the permutation symmetry of the nuclei is crucial. Some implicit conditions are not suitable for describing arbitrary systems, for example, in the case of linear molecules, the Eckart conditions are not well defined Eckart (1935); Sayvetz (1939). Thus, while our theory is general, the specific choice of $\mathbf{F}\left(\boldsymbol{\theta},\mathbf{R}^{\prime}\right)$ in Eq. 27 may restrict the range of applicability of the theory.

Up to this point, we have not introduced any approximations. We have set up the Hamiltonian to work with and derive the EOM for the electronic and nuclear Green’s functions by using it in Secs. III and IV. Before that, we discuss in more detail how to actually obtain the Euler angles appearing in the Hamiltonian.

II.2 Euler angles

The condition given by Eq. 27 defines the Euler angles by an implicit equation of the form $\mathbf{F}\left(\boldsymbol{\theta},\mathbf{x}+\mathbf{u}\right)=\mathbf{0}$. Explicit solutions of such equations do not always exist Guzman (2003); Krantz and Parks (2012); Krasnoshchekov et al. (2014), that is, it may be impossible to write $\boldsymbol{\theta}$ as an explicit function of the nuclear variables $\mathbf{x},\mathbf{u}$. However, for given $\mathbf{x},\mathbf{u}$ and given nuclear masses, a numerical solution can always be obtained Krasnoshchekov et al. (2014). We are particularly interested in how to obtain the rotation matrix $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$ and its derivatives. These quantities appear in the Hamiltonian when we do the Taylor expansions of some parts of the Hamiltonian (see Appendix B) and thus they will ultimately appear in the EOM.

The rotation matrix $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}+\mathbf{u}\right)\right]$ can be obtained from a generic condition given by Eq. 27. In turn, $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$ can be obtained by solving Eq. 27 with $\mathbf{u}_{k}=\mathbf{0}$ for all $k$. Thus, $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$ is obtained from the implicit equation $\mathbf{F}\left(\boldsymbol{\theta},\mathbf{x}\right)=\mathbf{0}$. For instance, in the case of Eckart conditions given by Eq. 28

[TABLE]

with the Euler angles $\boldsymbol{\theta}\left(\mathbf{x}\right)$, such that $0<\theta_{\beta}\left(\mathbf{x}\right)<2\pi$ for $\beta=1,2,3$. What we need next are the derivatives of $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$ appearing in the Hamiltonian (see Eqs. 131 and 132 of Appendix B). By the chain rule, we write

[TABLE]

Given the explicit form for the rotation matrix as function of the Euler angles we can calculate the derivatives like $\partial{\mathcal{R}}/\partial{\theta_{\beta^{\prime}}}$. Next we calculate the derivatives $\partial{\theta_{\beta^{\prime}}}/\partial{x_{\alpha}\left(k\right)}$ and $\partial^{2}{\theta_{\beta^{\prime}}}/\partial{x_{\alpha}\left(k\right)}\partial{x_{\beta}\left(k^{\prime}\right)}$, etc. After taking the total derivative of $F_{\sigma}\left(\boldsymbol{\theta},\mathbf{x},\mathbf{u}\right)$ with respect to $x_{\alpha}\left(k\right)$ we write

[TABLE]

where $A_{\sigma\beta^{\prime}}\equiv\partial{F_{\sigma}}/\partial{\theta_{\beta^{\prime}}}$. Suppose that the matrix $A$ is invertible such that

[TABLE]

By multiplying Eq. 31 with $A^{-1}_{\gamma\sigma}$, taking a sum over $\sigma$ and re-arranging, we obtain

[TABLE]

The result given by Eq. 33 is essentially included in the implicit function theorem de Oliveira et al. (2013). For the second order derivative, we take the total derivative of Eq. 33 with respect to $x_{\beta}\left(k^{\prime}\right)$ and after some algebra find that

[TABLE]

All the quantities in Eqs. 33 and 34 are to be evaluated at $\boldsymbol{\theta}=\boldsymbol{\theta}\left(\mathbf{x}\right),\mathbf{u}_{k}=\mathbf{0}$ for all $k$ when applied in solving the relations of Eq. 30 aiming to evaluate $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$. For instance, when the condition $F_{\sigma}\left(\boldsymbol{\theta},\mathbf{x},\mathbf{u}\right)$ is the Eckart condition given by Eq. 28, we have

[TABLE]

and analogously for the derivatives w.r.t. the other components. Here, $\epsilon_{\sigma\eta\nu}$ is the Levi-Civita symbol.

III Equation of motion for the electronic Green’s Function

Here we derive the exact EOM for the electronic Green’s function. We denote an ensemble average for an operator $\hat{o}\left(t\right)$ (in the Heisenberg picture) as

[TABLE]

where $\left\{\ket{\phi_{n}}\right\}$ is an orthonormal basis of the electron-nuclear Hilbert space. The density operator is the standard grand canonical statistical operator

[TABLE]

where

[TABLE]

and $\mu_{e}$ is the chemical potential of the electrons. One has to emphasize that the grand canonical Hamiltonian (Eq. 38) refers to a fixed number of nuclei, while for the electrons the chemical potential is fixed (rather than the particle number) Baym (1961). Here the time variable $t$ could be taken as a variable on a general time-contour such as the real time axis or on a more general time contour such as the Keldysh contour Keldysh (1965); Stefanucci and van Leeuwen (2013) allowing for a non-equilibrium formulation. This is justifiable since in both cases the EOM are the same Stefanucci and van Leeuwen (2013). In the following, we assume that the time variables are on the Keldysh contour $\gamma$ and denote the time variables on the contour with $z$.

We start by writing the EOM for the field operator $\hat{\psi}\left(\mathbf{r}z\right)$ (in the Heisenberg picture, that refers to the full electron-nuclear Hamiltonian), namely

[TABLE]

where

[TABLE]

and

[TABLE]

The quantities in Eq. 41 originate from the Coriolis and vibrational-rotational coupling terms, see Appendix A. Next we define the electron Green’s function as Martin and Schwinger (1959); Hedin (1965); Stefanucci and van Leeuwen (2013); Fetter and Walecka (1971); Mahan (1990)

[TABLE]

where $\mathcal{T}_{\gamma}$ is the contour time ordering operator such that

[TABLE]

and the Green’s function $G\left(1,2\right)$ satisfies the Kubo-Martin-Schwinger boundary conditions Martin and Schwinger (1959); Stefanucci and van Leeuwen (2013). Here and from now on, whenever convenient, we use the following shorthand notations

[TABLE]

The EOM for the Green’s function reads

[TABLE]

After using Eq. 39 in Eq. 45 we find that

[TABLE]

where

[TABLE]

and the total density is defined by $\hat{n}\left(1\right)\equiv\hat{n}_{e}\left(1\right)+\hat{n}_{n}\left(1\right)$. All the terms originating from the transformed kinetic energies $\hat{T}^{\prime}_{mpe},\hat{T}_{cvr}$ and $\hat{T}^{\prime}_{cvr}$ are included in the quantity $S\left(1,2\right)$ given by

[TABLE]

where $S_{c}\left(1,2\right)$ for each $c$ is given in Sec. V. By subtracting the quantity

[TABLE]

from both sides of Eq. 46 we find

[TABLE]

where

[TABLE]

Next we write the quantities $S\left(1,2\right)$ and $\tilde{S}\left(1,2\right)$ in terms of the corresponding SE’s, define (here $c=1,2,3$)

[TABLE]

By inverting (52), we find that

[TABLE]

By using Eq. 53 in Eq. 50

[TABLE]

where

[TABLE]

We can also write Eq. 46 in the form of a Dyson equation

[TABLE]

the function $G_{0}\left(1,2\right)$ being a solution of

[TABLE]

Most of the previous literature on the electron-nuclear many-body problem employed the laboratory frame formulation Baym (1961); Giustino (2017). By contrast the present article works with a body-fixed frame. This allows for an explicit inclusion of rotational and vibrational degrees of freedom and their coupling, in a consistent way. A first important step towards a consistent body fixed-frame formulation was taken in Ref. van Leeuwen (2004). In the present formulation, the self-energy $\Sigma_{t}$ also contains the kinetic energy contributions $T^{\prime}_{mpe},T_{cvr}$ and $T^{\prime}_{cvr}$. Neglecting these contributions leads to Eqs. 54 and 56 with $\Sigma_{c}\left(1,3\right)\equiv 0$ for $c=1,2,3$ and thus $\Sigma_{t}\left(1,3\right)=\Sigma\left(1,3\right)$. The resulting equations are analogous to those of Ref. van Leeuwen (2004). All the rather complicated kinetic energy terms, originating from the body-fixed frame transformation, are hidden in the SE’s $\Sigma_{c}\left(1,3\right)$. We will give an explicit, approximate form of each $\Sigma_{c}\left(1,3\right)$ in Eqs. LABEL:eq:CoriolisAndVibrationalRotCouplingTermsEq_16 and 87 of Sec. V. In the present chapter we have obtained the necessary EOM for the electrons. We will deduce the Hedin-like equations for the electrons in Sec. V.

IV Equation of motion for the nuclear Green’s function

Next we set out to derive the EOM for the nuclear Green’s function. The connection of the nuclear and electronic equations will be addressed in more detail in Sec. V. We start by writing the Heisenberg EOM for the displacements, namely

[TABLE]

After computing the commutators by using Eq. 20 [see also Eq. 41], we obtain

[TABLE]

Differentiating Eq. LABEL:eq:EquationsOfMotionForDisplacementsEq_2 with respect to time and taking an ensemble average leads to

[TABLE]

where (we use $\hat{T}^{\prime\prime}_{cvr}\equiv\hat{T}_{cvr}+\hat{T}^{\prime}_{cvr}$)

[TABLE]

In the first two quantities of Eq. 61, all the other contributions except those originating from the external potential terms are explicitly shown.

In order to obtain the EOM for the nuclear Green’s function $D_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$, we take a functional derivative of Eq. LABEL:eq:EquationsOfMotionForDisplacementsEq_2 with respect to $J_{\beta}\left(k^{\prime}z^{\prime}\right)$ and write

[TABLE]

where

[TABLE]

We write Eq. 62 in yet another form by using the inverse of $D_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ that respects the Kubo-Martin-Schwinger boundary conditions, namely, Eq. 62 can be written in terms of the nuclear SE’s

[TABLE]

as

[TABLE]

where

[TABLE]

By neglecting the kinetic energy terms $\hat{T}^{\prime}_{mpn},\hat{T}_{cvr}$ and $\hat{T}^{\prime}_{cvr}$ in the Hamiltonian given by Eq. 20 leads to Eq. 65 with $\Pi^{\left(c\right)}_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)=0$ for all $c$ and thus $\Pi^{t}_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)=\Pi_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$. In other words, all the effects of the mass-polarization, Coriolis and vibrational-rotational coupling terms on the EOM of the function $D_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ are hidden in the SE’s $\Pi^{\left(c\right)}_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$.

In the remaining part of this section, we deduce the EOM for the momentum-displacement and momentum-momentum Green’s functions since these functions are required in the calculation of the total energy, see Appendix B. In order to use the functional derivative technique we add the following potential to the external potential part of the Hamiltonian

[TABLE]

We note that this term would appear in the EOM of $\left\langle\hat{u}_{\alpha}\left(kz\right)\right\rangle$, but would not appear in the corresponding equations of $D_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$. We start by writing the EOM for momentum ensemble average, namely

[TABLE]

By taking a functional derivative of Eq. 68 with respect to $P_{\beta}\left(k^{\prime}z^{\prime}\right)$ we find that

[TABLE]

where

[TABLE]

Next we take a functional derivative of Eq. 68 with respect to $J_{\beta}\left(k^{\prime}z^{\prime}\right)$ and write

[TABLE]

where Eq. 63 was used and $D^{pu}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)\equiv\delta{\left\langle\hat{p}_{\alpha}\left(kz\right)\right\rangle}/\delta{J_{\beta}\left(k^{\prime}z^{\prime}\right)}$. We can use Eq. 64 and write Eq. 71 as

[TABLE]

where

[TABLE]

Furthermore, we have

[TABLE]

From Eq. 72 we see that provided we have the necessary SE’s and the solution for $D_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$, then we can obtain $D^{pu}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ without solving the EOM for it.

So far our derivation is exact. We have not made any approximations or simplifying assumptions. We still use the full Hamiltonian without restricting the expansion in displacements $\hat{\mathbf{u}}$ to some particular order, which is the usual procedure in the existing formulations Baym (1961); Maksimov (1975); Giustino (2017). We can actually write the exact total energy, defined as an ensemble average of the Hamiltonian, in terms of the quantities appearing in the EOM for the electrons and nuclei, see Appendix B. Hence we have deduced a formally exact Green’s function theory for arbitrary systems of electrons and nuclei given the Hamiltonian of kinetic energies and Coulombic interactions. Clearly the full solution will be hard to obtain for real-world systems and approximations are needed in practice. We discuss a special case of crystalline solids in Sec. VIII and give an explicit approximate expression for the Fourier transformed SE $\Pi_{\alpha\alpha^{\prime}}\left(k,k^{\prime},\omega\right)$ in Sec. VII.2.

V Hedin’s equations

Hedin-like equations for the complete system of electrons and nuclei can be derived in a similar way as has been done earlier van Leeuwen (2004). Namely, the equations can be written as

[TABLE]

Here, the vertex function $\Gamma\left(1,2,3\right)$ is almost the same as in Ref. van Leeuwen (2004), but contains contributions from the Coriolis, vibrational-rotational coupling and mass-polarization terms. The electronic polarization $P_{e}\left(1,2\right)$ is formally of the same form as in the earlier works, but also includes the aforementioned additional terms through the vertex function. The electronic part, $\Sigma_{e}\left(1,2\right)$, and the nuclear conribution, $\Sigma_{ph}\left(1,2\right)$ to the SE’s are

[TABLE]

The contributions to the screened Coulomb interaction can be written as

[TABLE]

and

[TABLE]

where

[TABLE]

Furthermore, the nuclear density-density correlation function appearing in Eq. (77) is

[TABLE]

where $\Delta\hat{n}_{n}\left(1\right)=\hat{n}_{n}\left(1\right)-\left\langle\hat{n}_{n}\left(1\right)\right\rangle$.

Formally, the present theory is similar to the conventional ones Hedin (1965); Giustino (2017) and thus many of the approximations, like the GW-approximation Hedin (1965); Aryasetiawan and Gunnarsson (1998); Stan et al. (2009); Rostgaard et al. (2010); Koval et al. (2014) or other suitable approximations Lani et al. (2012), can be established in a similar way as in the existing Green’s function theory for electrons within the BO approximation. Namely, in the GW-approximation $\Gamma\left(1,2,3\right)\approx\delta\left(1-2\right)\delta\left(1-3\right)$ and the electronic polarization and SE become $P_{e}\left(1,2\right)\approx-i\hbar G\left(1,2\right)G\left(2,1^{+}\right)$ and $\Sigma\left(1,2\right)\approx i\hbar W\left(1,2\right)G\left(1,2\right)$. The nuclear problem enters to the electronic equations, for instance, through the SE’s and the nuclear density-density correlation function $\rho_{n}\left(1,2\right)$. We can write $\hat{n}_{n}\left(1\right)$ and $\rho_{n}\left(1,2\right)$ by expanding them in a Taylor series with respect to the displacements and taking into account consistently all terms up to within a given order. It thus follows that in order to evaluate the expanded $\left\langle\hat{n}_{n}\left(1\right)\right\rangle$ and $\rho_{n}\left(1,2\right)$, we need to evaluate ensemble averages like $\left\langle\hat{u}_{\alpha}\left(kz\right)\right\rangle$, $\left\langle\hat{u}_{\alpha}\left(kz\right)\hat{u}_{\beta}\left(k^{\prime}z\right)\right\rangle$ and so on. Similar terms appear if we expand the nuclei terms included in the SE’s $\Sigma_{c}\left(1,2\right)$. In order to determine these quantities, we can use the EOM for the nuclear Green’s functions derived in Sec. IV.

Next we give exact and approximate expressions for the SE’s arising from the Coriolis, vibrational-rotational coupling and mass-polarization terms. The quantities $S_{c}\left(\mathbf{r}z,\mathbf{r}^{\prime}z^{\prime}\right)$ appearing in Eq. 48 are defined as

[TABLE]

where

[TABLE]

The perhaps simplest approximation to the SE’s $\Sigma_{c}\left(\mathbf{r}z,\mathbf{r}^{\prime}z^{\prime}\right)$ is obtained by introducing a mean-field-like factorization of the nuclear contributions:

[TABLE]

where

[TABLE]

A similar approximation is introduced in Sec. VII.2 when we approximate $\left\langle\hat{u}_{\alpha_{1}}\left(k_{1}z\right)\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle\approx\left\langle\hat{u}_{\alpha_{1}}\left(k_{1}z\right)\right\rangle\left\langle\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle$. If we further make the Hartree-Fock approximation Stefanucci and van Leeuwen (2013)

[TABLE]

we find by using Eq. 52

[TABLE]

and

[TABLE]

After Taylor expanding $\hat{\mathcal{D}}^{\prime}_{\beta\gamma}\left(t\right)$ and $\hat{L}^{\left(3\right)}_{\gamma\sigma\gamma^{\prime}\beta}$ in displacements $\hat{\mathbf{u}}$, $\Sigma_{2}\left(1,2\right)$ and $\Sigma_{3}\left(1,2\right)$ become functions of $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$, $D_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ and so on. One has to emphasize that, although derived in a BO-like manner, the terms in Eqs. (LABEL:eq:CoriolisAndVibrationalRotCouplingTermsEq_16) and (87) represent physical effects beyond the BO approximation.

VI Choice of reference positions

We have not yet discussed in detail how to choose the reference positions $\mathbf{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

[TABLE]

The value of $E_{tot}$ must be independent of $\mathbf{x}$ provided our expansion (in displacements $\hat{\mathbf{u}}$) of the SE’s, or the Hamiltonian ensemble average itself, converges for a given $\mathbf{x}$. However, in practice we are not able to make such an expansion up to arbitrary order in $\hat{\mathbf{u}}$ and solve the exact equations. Consequently, the value of $E_{tot}$ will become dependent on $\mathbf{x}$ due to the approximations made. For instance, if the positions $\mathbf{x}$ are chosen far away from the positions around which the nuclei would vibrate, then in order to calculate $E_{tot}$ accurately, the displacements $\mathbf{u}$ would be rather large. This implies that we would need a rather high order nuclear Green’s functions in order to describe the system properly and this is expected to be extremely difficult. If we choose the positions $\mathbf{x}$ rather poorly, i.e. far away from the equilibrium positions, and at the same time take only the lowest order quantities in $\hat{\mathbf{u}}$ into account, we expect to find rather unrealistic results. Our theory is formally exact and hence it also describes situations where the nuclei do not vibrate around some particular regions of space (like close to the melting point of a crystal, or beyond it in the liquid phase). However, it is beneficial to have a consistent strategy to obtain the parameters $\mathbf{x}$ such that those systems in which the nuclei perform such rather small vibrations, can be described with reasonable accuracy by using the approximations of lowest orders in $\hat{\mathbf{u}}$. We note that in these cases, the (anti)symmetrization with respect to those variables $\hat{\mathbf{R}}^{\prime}$ that refer to identical nuclei may not be necessary to obtain accurate results. If the lowest order approximations in $\hat{\mathbf{u}}$ describe the system in sufficient detail, the nuclei are well-localized near the positions $\mathbf{x}$, which also means that the effects of the (anti)symmetrization on the observables are expected to be rather small.

For the aforementioned systems, like some stable crystal lattices, we use the following method to determine the positions $\mathbf{x}$. The initial value of $\mathbf{x}$ is obtained in the same way as within the BO approximation, that is, either from experimental data of known structures or in the case of hypothetical structures by using the methods of structural chemistry Müller (2006); Karttunen et al. (2010). After the initial guess, which serves as starting point of an iteration, we aim to find the positions $\mathbf{x}$ such that they are equal to the expectation values of the nuclear positions (in the body-fixed frame), that is, we seek the positions such that

[TABLE]

We note that $\left\langle\hat{\mathbf{u}}_{k}\right\rangle$ is a function of $\mathbf{x}$ since the Hamiltonian is, as consequence of truncating the Taylor expansion in $\hat{\mathbf{u}}$ at a finite order. In the absence of the external potential terms, $\left\langle\hat{\mathbf{R}}^{\prime}_{k}\right\rangle$ is independent of time, and hence the equation determining the reference positions $\mathbf{x}$ becomes

[TABLE]

That is, if we are able to choose the positions $\mathbf{x}$ such that Eq. 89 holds, then the expectation values of the displacements vanish. Next we seek a way how to find the displacements satisfying Eq. 90 and thus the reference positions $\mathbf{x}$ satisfying Eq. 89. The determining equation for $\left\langle\hat{\mathbf{u}}_{k}\right\rangle$ is given by Eq. 60. If we put the external potential terms to zero, truncate the Taylor expansion (w.r.t. $\hat{\mathbf{u}}$) of the Hamiltonian at a finite order and choose as initial condition the equilibrium ensemble associated with this truncated Hamiltonian, then $\left\langle\hat{u}_{\alpha}\left(kz\right)\right\rangle=\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$ and thus

[TABLE]

meaning that the expectation value of the total force on each nucleus $k$ vanishes. Thus, if the external potential terms vanish, Eq. 60 becomes

[TABLE]

The quantities $K_{\alpha}\left(kz\right)$ and $K^{\left(c\right)}_{\alpha}\left(kz\right)$ are given by Eq. 61 (here without the external potential terms). We define the total force

[TABLE]

and decompose it further as $F_{\alpha}\left(kz\right)=F^{\prime}_{\alpha}\left(kz\right)+F^{\prime\prime}_{\alpha}\left(kz\right)$, where $F^{\prime}_{\alpha}\left(kz\right)$ is the part of $F_{\alpha}\left(kz\right)$ which is a function of $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$ and $F^{\prime\prime}_{\alpha}\left(kz\right)$ is the remaining part (not a function of $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$). We can therefore write Eq. 92 as

[TABLE]

This is the determining equation for $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$. If we now want to find $\mathbf{x}$ such that $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle=0$, then Eq. (94) becomes

[TABLE]

We now give an example where an explicit and approximate form of Eq. 95 is given. Suppose that the only non-zero term in Eq. 92 is $K_{\alpha}\left(kz\right)$ meaning that we take into account terms which originate from the Coulomb interactions and neglect all the other terms. This is expected to be a relatively good approximation for most 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_{\alpha}\left(kz\right)$ by expanding $\hat{V}^{\prime}_{en}$ and $\hat{V}^{\prime}_{nn}$ in displacements (see Appendix B) and retain only the lowest orders. After calculating the commutators, we obtain to the lowest orders in displacements

[TABLE]

Next we use the following result

[TABLE]

where $\tilde{W}_{e}\left(\mathbf{r}_{3}z_{3},\mathbf{r}_{2}z_{2}\right)$ and $\tilde{\epsilon}_{e}\left(\mathbf{r}_{1}z_{1},\mathbf{r}_{2}z_{2}\right)$ are given by Eqs. 78 and 79, respectively. These quantities are obtained from the solution of Hedin’s equations. By using Eq. 97 we find that

[TABLE]

Here, the quantity $n^{\left(0\right)}_{n}\left(\mathbf{x},\mathbf{r}\right)$ is given by Eq. 132 and related to this, we consider how to obtain the rotation matrix $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$ and its derivatives in Sec. II.2. When we use Eq. 98 in Eq. 96 to approximate $\left\langle\hat{u}_{\alpha}\left(kz\right)\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle$, it can be seen that if we seek the solution such that $\left\langle\hat{\mathbf{u}}_{k}\right\rangle=\mathbf{0}$, meaning that $F^{\prime}_{\alpha}\left(kz\right)$ vanishes, Eq. 95 is of the following form

[TABLE]

where (Eq. 63)

[TABLE]

If we further approximate and neglect the last term of Eq. 99 we obtain

[TABLE]

Equations 99 and 101 are examples of approximate conditions for choosing the parameters $\mathbf{x}$ such that Eq. 89 holds. We call the positions $\mathbf{x}$ the nuclear equilibrium positions. The approximate condition given by Eq. 101 certainly makes a lot of sense: The equilibrium position of the $k$-th nucleus is the location where the total electrostatic force on this nucleus (coming from the other nuclei and the electron cloud) vanishes. This is formally the same condition as in the BO approximation Baroni et al. (2001), when the Hellmann-Feynman theorem Hellmann (1937); Feynman (1939) is used to calculate the total force. The difference is that our parameters $\mathbf{x}$ refer to the body-fixed frame. Further, the electron density $\left\langle\hat{n}_{e}\left(\mathbf{r}\right)\right\rangle$ is an ensemble average taken with respect to the full electron-nuclear Hamiltonian rather than an expectation value with respect to an eigenstate of the electronic BO Hamiltonian. The electron density can be obtained from the electron Green’s function as $\left\langle\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle=-i\hbar G\left(\mathbf{r}z,\mathbf{r}z^{+}\right)$. Thus, for a given $\mathbf{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\left(\mathbf{r}z,\mathbf{r}^{\prime}z^{\prime}\right)$, $P_{e}\left(\mathbf{r}z,\mathbf{r}^{\prime}z^{\prime}\right)$, $\tilde{W}_{e}\left(\mathbf{r}z,\mathbf{r}^{\prime}z^{\prime}\right)$ and $D^{\mathcal{T}}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ and we can check whether or not Eq. 99, or the corresponding general expression given by Eq. 95, holds and indicates that we have found the positions such that Eq. 89 is satisfied.

Next we give an overview of the work flow how the coupled set of equations for the electrons and nuclei could be solved. The procedure is summarized in Fig. 1. First we give an initial guess for the reference positions $\mathbf{x}$ and solve $\boldsymbol{\theta}$ from Eq. 27 with $\mathbf{u}_{k}=\mathbf{0}$, see Sec. II.2. Given the reference positions we expand all the quantities in Hedin’s equations, which are functions of the nuclear variables, into a Taylor series in $\hat{\mathbf{u}}$. In the first iteration, we retain only the zeroth order terms meaning that $W_{ph}\left(1,2\right)=0$ from which it follows that $W\left(1,2\right)=W_{e}\left(1,2\right)$ and $\Sigma\left(1,2\right)=\Sigma_{e}\left(1,2\right)$ (see Eq. 75). With the preceding choices, Hedin’s equations are actually similar to those in the BO approximation apart from the SE’s $\Sigma_{c}\left(1,2\right)$, but here all quantities refer to the body-fixed frame. The zeroth order approximations for the SE’s $\Sigma_{c}\left(1,2\right)$ can be obtained by using the relations given by Eqs. 52 and 83. After the Hedin equations are solved, we calculate the nuclear SE, see for example Sec. VII.2. Given the nuclear SE, the EOM for the nuclear Green’s function can be solved. Once we have obtained the nuclear Green’s function, we can write Hedin’s equations with non-zero $W_{ph}\left(1,2\right)$, solve these equations with the SE $\Sigma\left(1,2\right)=\Sigma_{e}\left(1,2\right)+\Sigma_{ph}\left(1,2\right)$ and by including the nuclear contributions beyond the zeroth order in $\hat{\mathbf{u}}$ to the SE’s $\Sigma_{c}\left(1,2\right)$. We iterate the electronic and nuclear equations in order to minimize the grand potential, or in the case of zero temperature formalism, the total energy $E_{tot}$ (see Appendix B). After the convergence, we test whether or not the reference positions $\mathbf{x}$ are the equilibrium positions by checking if Eq. 95 is satisfied. We iterate this whole process until we have found the equilibrium positions.

VII Nuclear vibrations and self-energy

VII.1 Nuclear vibrations

We start by writing the nuclear EOM in frequency space, but before Fourier transforming, all the external potential terms are set equal to zero. It follows that the nuclear Green’s function and the nuclear SE’s are then functions of time differences only and we can write Eq. 65 in terms of the Fourier transformed quantities as

[TABLE]

Here the reader is free to choose, for example, the time-ordered or retarded components of these functions as the Eq. 65 is defined on a general time contour. We write Eq. 102 in matrix form and re-arrange such that

[TABLE]

where $\tilde{\mathbf{D}}\left(\omega\right)\equiv\mathbf{M}^{1/2}\mathbf{D}\left(\omega\right)\mathbf{M}^{1/2}$ and $\mathbf{C}^{t}\left(\omega\right)\equiv\mathbf{M}^{-1/2}\mathbf{\Pi}^{t}\left(\omega\right)\mathbf{M}^{-1/2}$. The components of these matrices are denoted by $\alpha k$ and so on. Next we write the SE as $\mathbf{C}^{t}\left(\omega\right)=\mathbf{C}^{A}+\mathbf{C}^{NA}\left(\omega\right)$ where we choose the ”adiabatic part” $\mathbf{C}^{A}$ as the zero-frequency-limit $\mathbf{C}^{A}=\mathbf{C}^{t}\left(0\right)$. The ”non-adiabatic” part $\mathbf{C}^{NA}$ is the remainder, i.e. $\mathbf{C}^{NA}\left(\omega\right)=\mathbf{C}^{t}\left(\omega\right)-\mathbf{C}^{t}\left(0\right)$. It is easy to see that $\mathbf{C}^{A}$ is a Hermitian matrix. Hence we can find a complete and orthonormal set of eigenvectors satisfying

[TABLE]

We call the quantities $\tilde{\omega}_{j}$ the adiabatic normal mode frequencies. The eigenvalues $\tilde{\omega}^{2}_{j}$ are real since $\mathbf{C}^{A}$ is Hermitian and $\tilde{\omega}_{j}$ are real if $\mathbf{C}^{A}$ is positive definite. We transform Eq. 103 by using the eigenvectors of $\mathbf{C}^{A}$ such that

[TABLE]

where $\boldsymbol{\mathcal{D}}\left(\omega\right)\equiv\mathbf{v}^{\dagger}\tilde{\mathbf{D}}\left(\omega\right)\mathbf{v}$ and $\boldsymbol{\mathcal{C}}^{NA}\left(\omega\right)\equiv\mathbf{v}^{\dagger}\mathbf{C}^{NA}\left(\omega\right)\mathbf{v}$. If $\boldsymbol{\mathcal{C}}^{NA}\left(\omega\right)$ is sufficiently small, the quasi-particle picture holds and $\boldsymbol{\mathcal{C}}^{NA}\left(\omega\right)$ can be pictured in generating interactions between adiabatic normal modes by shifts of their values and finite lifetimes [imaginary part of $\boldsymbol{\mathcal{C}}^{NA}\left(\omega\right)$]. We can obtain the Green’s function of complex frequency by just replacing the real frequency $\omega$ by the complex one Farid (1999) in Eq. 105. The usual procedure is to consider the Matsubara Green’s functions continued to the complex frequency plane Baym and Mermin (1961); Baym (1961); Maradudin and Fein (1962); Rickayzen (1980); Karlsson and van Leeuwen (2016). We point out that the present discussion is still completely general. In Sec. VIII.1 we will deduce a representation of $\mathbf{D}\left(\omega\right)$ in terms of the phonon basis of crystalline solids.

VII.2 Nuclear self-energy

In order to derive an approximate form for the SE $\Pi_{\alpha\alpha^{\prime}}\left(k,k^{\prime},\omega\right)$, we start from the explicit form of the terms included in $K_{\alpha\beta^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$ (see Eq. 63). We have already calculated the commutators included to the lowest orders in $\hat{u}_{\alpha}\left(kz\right)$ and we can obtain $K_{\alpha\beta^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$ by taking a functional derivative of $K_{\alpha}\left(kz\right)$ given by Eq. 96 with respect to $J_{\beta}\left(k^{\prime}z^{\prime}\right)$, namely

[TABLE]

The terms visible in Eq. 106 are all the terms included in $K_{\alpha\beta^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$ 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 Baym (1961); Giustino (2017). We further evaluate the quantity in the second term on the right hand side of Eq. 106

[TABLE]

Finally, by using Eqs. 97, 106 and LABEL:eq:PhononSelfEnergyEq_5 in Eq. 64, we can write the nuclear SE in frequency space as:

[TABLE]

In Eq. 108, $\tilde{W}_{e}\left(\mathbf{r},\mathbf{r}^{\prime},\omega\right)$ is obtained by Fourier transforming $\tilde{W}_{e}\left(\mathbf{r}z,\mathbf{r}^{\prime}z^{\prime}\right)$ (Eq. 78) in the relative time variable $z-z^{\prime}$. Terms explicitly shown in Eq. 108 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 found in the harmonic approximation appear in Eq. 108. Namely, the first term on the right hand side of Eq. LABEL:eq:PhononSelfEnergyEq_5 does not appear because we have made the approximation $\left\langle\hat{u}_{\alpha}\left(kz\right)\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle\approx\left\langle\hat{u}_{\alpha}\left(kz\right)\right\rangle\left\langle\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle$. Further, only some of the second order terms of the right hand side of Eq. 97 are included. The validity of this approximation, i.e. of making $\hat{u}_{\alpha}\left(kz\right)$ and $\hat{n}_{e}\left(\mathbf{r}z\right)$ uncorrelated in the sense that the expectation value of their product is the product of their expectation values has to be assessed carefully when the states with respect to which we define the Green’s functions are not the BO eigenstates.

Sometimes it is convenient to write Giustino (2017) $\Pi_{\alpha\beta}\left(k,k^{\prime},\omega\right)=\Pi^{A}_{\alpha\beta}\left(k,k^{\prime}\right)+\Pi^{NA}_{\alpha\beta}\left(k,k^{\prime},\omega\right)$ where the adiabatic SE $\Pi^{A}_{\alpha\beta}\left(k,k^{\prime}\right)=\Pi_{\alpha\beta}\left(k,k^{\prime},0\right)$ is Hermitian while the non-adiabatic part $\Pi^{NA}_{\alpha\beta}\left(k,k^{\prime},\omega\right)=\Pi_{\alpha\beta}\left(k,k^{\prime},\omega\right)-\Pi_{\alpha\beta}\left(k,k^{\prime},0\right)$ is non-Hermitian, in general. In the present case (Eq. 108), the contribution from the Coulombic nuclei-nuclei interaction to the SE is $\partial^{2}{V_{nn}\left(\mathbf{x}\right)}/\partial{x_{\alpha}\left(k\right)}\partial{x_{\beta}\left(k^{\prime}\right)}$, while the rest of the terms visible in Eq. 108 originates from the Coulomb electron-nuclei interaction. For a better comparison with existing theories, we use Eq. 108 to write the non-adiabatic contribution as

[TABLE]

where terms including the quantity $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$ are not explicitly shown. Provided we can choose the parameters $\mathbf{x}$ such that Eq. 89 holds, then terms involving the quantities $\left\langle\hat{u}_{\alpha}\left(k\right)\right\rangle$ vanish, see Sec. VI. The relation given by Eq. 109 is formally analogous with the result obtained earlier Giustino (2017). However, the difference is in the densities $n^{\left(0\right)}_{n}\left(\mathbf{x},\mathbf{r}^{\prime}\right)$ (see Eq. 132 of Appendix B) and here all the variables refer to the body-fixed frame. We note that if necessary, the higher order expressions for the nuclear and phonon SE’s can be obtained by systematically including higher order terms in the expansion of $K_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$. Related to this, some results of a recent study Marini and Pavlyukh (2018) on generic electron-boson Hamiltonians may be useful within the present theory as well.

While we do not give explicit expressions for the SE’s $\Pi^{\left(c\right)}_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$, it is obvious that these quantities can be obtained with the same procedure as established here for $\Pi_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$. Namely, calculate the commutators contained in $K^{\left(c\right)}_{\alpha}\left(kz\right)$ (Eq. 61) by expanding the quantities involved in $\hat{\mathbf{u}}$, take the functional derivative with respect to $J_{\beta}\left(k^{\prime}z^{\prime}\right)$ in order to obtain $K^{\left(c\right)}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ and then use Eq. 64 to calculate $\Pi^{\left(c\right)}_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$. In obtaining these expressions, similar approximations can be made as we did in the case of $\Pi_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$. We leave a more detailed discussions of these terms for future work.

VIII Crystalline solids

In this section, we apply the general theory to crystalline solids. First of all, the electronic and nuclear mass polarization terms $\hat{T}^{\prime}_{mpe}$, $\hat{T}^{\prime}_{mpn}$ and $\hat{T}^{\prime}_{cvr}$ appearing in the transformed kinetic energy are proportional to the inverse of the total nuclear mass and are thus small for crystals. Further, we assume that one can find an implicit condition of the form (Eq. 27) such that the contributions from the kinetic energy $\hat{T}_{cvr}$ become small. Some justifications for their neglect have been given van Leeuwen (2004), when the implicit condition is chosen to be the Eckart condition given by Eq. 28. By neglecting the aforementioned terms the EOM for electrons remain otherwise the same, except that we have $\Sigma_{t}\left(1,3\right)=\Sigma\left(1,3\right)$ in Eqs. 54 and 56 and in the Hedin’s equations given by Eq. 75. In turn, the nuclear EOM remains otherwise the same, but in Eq. 65 the SE is $\Pi^{t}_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)=\Pi_{\alpha\alpha^{\prime}}\left(kz,k^{\prime}z^{\prime}\right)$.

VIII.1 Phonons and their interactions

Here we use a notation suitable for the description of crystalline solids and impose periodic boundary conditions Huang and Born (1954); Maradudin et al. (1971); Maradudin (1974). We label the nucleus by $k=l\kappa$, and write $\mathbf{x}_{k}=\mathbf{x}_{l\kappa}=\mathbf{x}_{l}+\mathbf{x}_{\kappa}$, where $\mathbf{x}_{\kappa}$ is the position vector of the nuclei $\kappa$ within the unit cell and

[TABLE]

is the lattice translational vector of the $l$th unit cell with integers $l_{1}$, $l_{2}$, $l_{3}$, and the vectors $\mathbf{a}_{j}$ being the primitive translational vectors of the lattice. With this notation, we write $\mathbf{R}^{\prime}_{k}=\mathbf{R}^{\prime}_{l\kappa}=\mathbf{x}_{l\kappa}+\mathbf{u}_{l\kappa}$, where $l\kappa$ goes over $N_{n}-1$ values in total.

Our aim is to define phonon frequencies beyond the BO approximation as has been established earlier by using the Green’s function approach Baym (1961); Giustino (2017). However, we do not assume the harmonic approximation. In the present case, Eq. 102 can be written otherwise the same, but with $\Pi^{t}_{\alpha\alpha^{\prime}}\left(l\kappa,l^{\prime\prime}\kappa^{\prime\prime},\omega\right)=\Pi_{\alpha\alpha^{\prime}}\left(l\kappa,l^{\prime\prime}\kappa^{\prime\prime},\omega\right)$. We also note that the quantities like $D_{\alpha\beta}\left(l\kappa,l^{\prime}\kappa^{\prime},\omega\right)$ and $\Pi_{\alpha\beta}\left(l\kappa,l^{\prime}\kappa^{\prime},\omega\right)$ are dependent only on the difference of the cell indices $l-l^{\prime}$. Therefore, we write $D_{\alpha\beta}\left(l\kappa,l^{\prime}\kappa^{\prime},\omega\right)$ as a discrete Fourier transform in the relative coordinate allowed by the periodic boundary conditions Huang and Born (1954), namely

[TABLE]

and in a similar way for $\Pi_{\alpha\alpha^{\prime}}\left(\kappa\kappa^{\prime},l,\omega\right)$. In Eq. 111, $N$ is the number of $\mathbf{q}$ points and thus the number of unit cells in the Born-von Karman cell Huang and Born (1954); Maradudin et al. (1971). By using Eq. 111 and the analogous expression for $\Pi_{\alpha\alpha^{\prime}}\left(\kappa\kappa^{\prime},l^{\prime\prime},\omega\right)$ in Eq. 102, we obtain

[TABLE]

and after some re-arranging and by using matrix notation

[TABLE]

where $\mathbf{C}\left(\mathbf{q},\omega\right)\equiv\mathbf{M}^{-1/2}\mathbf{\Pi}\left(\mathbf{q},\omega\right)\mathbf{M}^{-1/2}$ and $\tilde{\mathbf{D}}\left(\mathbf{q},\omega\right)\equiv\mathbf{M}^{1/2}\mathbf{D}\left(\mathbf{q},\omega\right)\mathbf{M}^{1/2}$. The components of the matrix like $\mathbf{C}\left(\mathbf{q},\omega\right)$ are labeled by $\alpha\kappa$ and $\beta\kappa^{\prime}$, etc. Following the procedure of Sec. VII.1, namely, we write the SE as $\mathbf{C}\left(\mathbf{q},\omega\right)=\mathbf{C}^{A}\left(\mathbf{q}\right)+\mathbf{C}^{NA}\left(\mathbf{q},\omega\right)$, where $\mathbf{C}^{A}\left(\mathbf{q}\right)=\mathbf{C}\left(\mathbf{q},0\right)$ and $\mathbf{C}^{NA}\left(\mathbf{q},\omega\right)=\mathbf{C}\left(\mathbf{q},\omega\right)-\mathbf{C}\left(\mathbf{q},0\right)$. The eigenvalue equation for $\mathbf{C}^{A}\left(\mathbf{q}\right)$ can be written as

[TABLE]

We call the quantities $\omega_{\mathbf{q}j}$ the adiabatic phonon frequencies. The eigenvalue equation given by Eq. 114 is analogous to the eigenvalue equation written for the dynamical matrix in the conventional theory of lattice dynamics Huang and Born (1954); Maradudin et al. (1971). As the adiabatic phonon frequencies $\omega_{\mathbf{q}j}$ are defined by Eq. 114, in contrast to the conventional theory and existing theories beyond the BO approximation Baym (1961); Giustino (2017), already the non-interacting phonons potentially contain terms up to arbitrary order in $\hat{\mathbf{u}}$. In order for $\omega^{2}_{\mathbf{q}j}$ to be positive and thus $\omega_{\mathbf{q}j}$ to be real, the matrix $\mathbf{C}^{A}\left(\mathbf{q}\right)$ has to be positive definite, which is not in general the case for a given $\mathbf{x}$. This is also the case in the BO theory of lattice dynamics Huang and Born (1954), where the positive definiteness of the dynamical matrix implies the minimum of the BO energy surface and the stability of the crystal lattice Born (1940). Provided the matrix $\mathbf{C}^{A}\left(\mathbf{q}\right)$ is Hermitian (as it is, for example, if Eq. 108 is used), the components of the eigenvector $e_{\alpha}\left(\kappa|\mathbf{q}j\right)$ can be chosen to satisfy the orthonormality and completeness conditions

[TABLE]

We use the eigenvectors of the adiabatic SE to rewrite Eq. 113 in terms of the adiabatic phonon basis. After the transformation we obtain

[TABLE]

where $\boldsymbol{\mathcal{D}}\left(\mathbf{q},\xi\right)\equiv\mathbf{e}^{\dagger}\left(\mathbf{q}\right)\tilde{\mathbf{D}}\left(\mathbf{q},\xi\right)\mathbf{e}\left(\mathbf{q}\right)$ and $\boldsymbol{\mathcal{C}}^{NA}\left(\mathbf{q},\xi\right)\equiv\mathbf{e}^{\dagger}\left(\mathbf{q}\right)\mathbf{C}^{NA}\left(\mathbf{q},\xi\right)\mathbf{e}\left(\mathbf{q}\right)$. We call $\boldsymbol{\mathcal{D}}\left(\mathbf{q},\xi\right)$ the phonon Green’s function and $\boldsymbol{\mathcal{C}}^{NA}\left(\mathbf{q},\xi\right)$ the non-adiabatic phonon SE. Here, $\xi$ denotes a complex frequency variable Farid (1999). If the non-adiabatic SE vanishes, we recover the adiabatic phonon Green’s function $\mathcal{D}^{A}_{jj}\left(\mathbf{q},\omega\right)=1/\left(\omega^{2}-\omega^{2}_{\mathbf{q}j}\right)$ which is of the usual form and appears also in the BO theory of lattice dynamics Maradudin and Fein (1962). The adiabatic phonons have infinite lifetimes if the non-adiabatic SE vanishes. The non-adiabatic 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. The shifts to the adiabatic phonon frequencies and finite lifetimes of adiabatic phonons can be obtained as in the laboratory frame formulation Allen (1972); Grimvall (1981); Giustino (2017).

By using Eqs. 78 and 79 in Eq. 109 and then Fourier transforming and changing the representation, we obtain the following approximate form for the non-adiabatic phonon SE

[TABLE]

In Eq. 117

[TABLE]

where $\tilde{g}_{\mathbf{q}j}\left(\mathbf{r}\right)$ is the complex conjugate of $\tilde{g}^{\ast}_{\mathbf{q}j}\left(\mathbf{r}\right)$. The diagrams corresponding to $\boldsymbol{\mathcal{D}}\left(\mathbf{q},\xi\right)$, $\boldsymbol{\mathcal{C}}^{NA}\left(\mathbf{q},\xi\right)$ are of a similar form as in the laboratory frame formulation Giustino (2017, 2019).

VIII.2 Momentum functions

By making the approximations discussed at the beginning of this section and then comparing Eqs. 65 and 72 we see that

[TABLE]

and thus for the Fourier transforms

[TABLE]

Therefore, after we have obtained a solution for $D_{\alpha\beta}\left(l\kappa,l^{\prime}\kappa^{\prime},\omega\right)$ we can also obtain $D^{pu}_{\alpha\beta}\left(l\kappa,l^{\prime}\kappa^{\prime},\omega\right)$. The last function we need is the momentum-momentum Green’s function $D^{pp}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$ and in the present case Eq. 69 becomes

[TABLE]

where we approximate

[TABLE]

which can be obtained directly from Eq. 96. By employing the same approximations for $K^{\prime}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$, as used in Sec. VII.2 in writing $K_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$, we find that Eq. 121 can be written as

[TABLE]

The Fourier transform of this equation is

[TABLE]

Since by Eq. 74, $D^{up}_{\beta\alpha}\left(k^{\prime}z^{\prime},kz\right)=D^{pu}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$, we can write

[TABLE]

where Eq. 120 was used. Now we have all the necessary equations in place to calculate the total energy of the system, e.g.

IX Conclusions

In this work, we have derived a coupled and self-consistent set of exact equations for the electronic and nuclear Green’s functions following from the Hamiltonian of Coulomb interactions and kinetic energies as a starting point. The present theory, when applied to crystalline solids, resembles in some aspects the previous ones Baym (1961); Maksimov (1975); Giustino (2017). However, we resolve an issue arising from the translational and rotational symmetry of the Hamiltonian. This symmetry prevents the use of the existing many-body Green’s function theories beyond the BO approximation in describing systems other than those with constant density eigenstates. The present theory is formally exact and it is not limited to the harmonic approximation. The complexity of our system of EOM is of the same order as the existing ones when applied to crystalline solids.

In addition to the general EOM, we specifically consider the normal modes. For the special case of crystalline solids, phonons and their interactions beyond the BO approximation are rigorously defined and discussed in detail. 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. For instance, numerically tractable approximations to some parts of the nuclear SE originating from the Coriolis and vibrational-rotational couplings are yet to be derived. Our main emphasis in this work was on crystalline solids and we leave a more detailed treatment of molecules for future work. The practical implementations of the present method are under preparation.

To summarize, the present theory allows one to go beyond the BO approximation in a systematic and formally exact way by using the method of many-body Green’s functions. We expect it will become a useful tool in treating beyond-Born-Oppenheimer effects in solids and molecules.

Appendix A Coriolis and vibrational-rotational coupling terms

The kinetic energy contributions $\hat{T}_{cvr}$ and $\hat{T}^{\prime}_{cvr}$ are

[TABLE]

In Eq. 126

[TABLE]

and

[TABLE]

Here we divided the Coriolis and vibrational-rotational coupling terms into two parts: one part, $\hat{T}^{\prime}_{cvr}$, is proportional to the inverse of the total nuclear mass while the other, $\hat{T}_{cvr}$, is not. All the quantities defined by Eqs. 127 and 128 are functions of $\mathbf{x}$ and $\hat{\mathbf{u}}$. Strictly speaking, the derivatives like $\partial{\hat{\mathcal{R}}_{\sigma\sigma^{\prime}}}/\partial{u_{\alpha}\left(k^{\prime}\right)}$, are not so well defined objects from a notational point of view. By this notation we mean that the quantities like $\partial{\mathcal{R}_{\sigma\sigma^{\prime}}}/\partial{u_{\alpha}\left(k^{\prime}\right)}$ and $\partial^{2}{\mathcal{R}_{\sigma\sigma^{\prime}}}/\partial{u_{\alpha}\left(k\right)}\partial{u_{\alpha}\left(k^{\prime}\right)}$, are some functions of $\mathbf{x},\mathbf{u}$ after differentiation and we can then find the corresponding functions of operators $\hat{\mathbf{u}}$. For example if $\partial{\mathcal{R}}/\partial{u_{\alpha}\left(k\right)}=G_{\alpha}\left(k,\mathbf{x},\mathbf{u}\right)$, then in our notation $\partial{\hat{\mathcal{R}}}/\partial{u_{\alpha}\left(k\right)}=G_{\alpha}\left(k,\mathbf{x},\hat{\mathbf{u}}\right)$.

We also use the following form of the Coriolis and vibrational-rotational coupling terms

[TABLE]

where

[TABLE]

and the quantities $\hat{L}^{\left(1\right)}_{\beta\gamma}$, $\hat{L}^{\left(2\right)}_{\alpha\beta\gamma}\left(k\right)$ and $\hat{L}^{\left(3\right)}_{\gamma\sigma\gamma^{\prime}\beta}$ are defined by Eq. 41.

Appendix B Total Energy

We expand $\hat{n}_{n}\left(\mathbf{r}\right)$, $\hat{V}_{en}\left(\mathbf{r}\right)$ and $\hat{V}^{\prime}_{nn}$ into Taylor series with respect to the displacements:

[TABLE]

and

[TABLE]

Here, the following notations are used

[TABLE]

and so on. In Sec. II.2 we discuss how to actually obtain the rotation matrix $\mathcal{R}\left[\boldsymbol{\theta}\left(\mathbf{x}\right)\right]$ and the necessary derivatives of it.

Next we deduce an exact form of the total energy of the system. We start by writing

[TABLE]

where

[TABLE]

such that

[TABLE]

We note that from these expansions it follows that

[TABLE]

and

[TABLE]

where we used $\left\langle\hat{n}_{e}\left(\mathbf{r}z\right)\right\rangle=-i\hbar G\left(\mathbf{r}z,\mathbf{r}z^{+}\right)$. We write the quantity defined by Eq. 61 as $K_{\alpha}\left(kz\right)=K^{en}_{\alpha}\left(kz\right)+K^{nn}_{\alpha}\left(kz\right)$ such that

[TABLE]

where $y=en,nn$. In a similar way we write the quantity defined by Eq. 63 as $K_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)=K^{en}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)+K^{nn}_{\alpha\beta}\left(kz,k^{\prime}z^{\prime}\right)$, where

[TABLE]

We define

[TABLE]

Next we use Eq. 134 and write $K^{en}_{\alpha}\left(kz\right)$ and $K^{nn}_{\alpha}\left(kz\right)$ such that

[TABLE]

and in a similar way for all the quantities derived by using $K^{en}_{\alpha}\left(kz\right)$ and $K^{nn}_{\alpha}\left(kz\right)$. By using Eq. 142 in Eqs. 140 and 141 and these results with Eqs. LABEL:eq:UsefulRelationsEq_10_1 and 138 we eventually find that

[TABLE]

and

[TABLE]

We use these results to write the total energy $E_{tot}$ in terms of the Green’s functions and related quantities as

[TABLE]

where the external potential terms are put to zero and $\hat{T}_{tot}$ is given by Eq. 21. We have

[TABLE]

Here

[TABLE]

and this function can be obtained from the solution of Eq. 69 since with the external potentials put to zero $\left\langle\hat{p}_{\alpha}\left(kz\right)\right\rangle=\left\langle\hat{p}_{\alpha}\left(k\right)\right\rangle$. By using the EOM for the field operator given by Eq. 39 without the external potentials we find that

[TABLE]

where $\hat{T}^{\left(c\right)}_{cvr}$ are given by Eqs. 129 and 130. This result can be found Hedin and Lundqvist (1970) by multiplying Eq. 39 from the left with $\hat{\psi}^{\dagger}\left(\mathbf{r}z\right)/2$, integrating over $\mathbf{r}$, taking an ensemble average and then establishing some rearranging. From Eqs. 52, LABEL:eq:CoriolisAndVibrationalRotCouplingTermsEq_11 and 130, it follows that

[TABLE]

Now we have all the necessary ingredients in place to express the exact total energy in terms of the Green’s functions and SE’s. By using Eqs. 143, 144, 148 and 149, the total energy can be written as

[TABLE]

This is an exact form of $E_{tot}$. Provided $\mathbf{x}$ are the equilibrium positions, the second term on the right hand side of Eq. 150 involving $\left\langle\hat{u}_{\alpha}\left(kz\right)\right\rangle$ vanishes and $D_{\beta\alpha}\left(k^{\prime}z^{\prime},kz\right)=D^{\mathcal{T}}_{\beta\alpha}\left(k^{\prime}z^{\prime},kz\right)$.

Acknowledgements.

V.J.H. thanks Dr. Ivan Gonoskov for useful discussions on various aspects of the present work. E.K.U.G. acknowledges financial support by the European Research Council Advanced Grant FACT (ERC-2017-AdG-788890).