Nonadiabaticity from first principles: Exact-factorization approach for solids

PHYSICAL REVIEW B 112, 075102 (2025) Editors’ Suggestion Nonadiabaticity from first principles: Exact-factorization approach for solids Galit Cohen ,1,*Rachel Steinitz-Eliyahu ,1,*E. K. U. Gross ,2,†Sivan Refaely-Abramson ,1,†and Ryan Requist 3,2,† 1Department of Molecular Chemistry and Materials Science, Weizmann Institute of Science, Rehovot 7610001, Israel 2Fritz Haber Center for Molecular Dynamics, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 91904, Israel 3Division of Theoretical Physics, Ruđer Boškovi´c Institute, Zagreb 10000, Croatia (Received 15 February 2025; accepted 2 June 2025; published 1 August 2025) The thorough treatment of electron-lattice interactions from first principles is one of the main goals in condensed matter physics. While the commonly applied adiabatic Born-Oppenheimer approximation is sufficient for describing many physical phenomena, it is limited in its ability to capture meaningful features originating from nonadiabatic coupling effects. The exact factorization method, starting from the full Hamiltonian of electrons and nuclei, provides a way to systematically account for nonadiabatic effects. This formalism was recently developed into an ab initio density functional theory framework. Within this framework we develop here a perturbative approach to the electronic states in solid state materials. We derive exact-factorization-based perturbations of the Kohn-Sham states up to second order in the nuclear displacements. These nonadiabatic features in the calculated energy and wave function corrections are expressed in terms of readily available density functional perturbation theory components. DOI: 10.1103/dmpv-zqdh I. INTRODUCTION Theoretical understanding of electron-phonon interactions is central to the realization and optimization of numerous physical phenomena, from the computation of fundamental ground and excited state properties in materials to optoelectronic phenomena and dynamical scattering processes dominating energy conversion and transport in materials [1–9]. Calculations of these interactions are typically made within the adiabatic Born-Oppenheimer (BO) approximation [10,11], in which it is assumed that the electrons adapt immediately to the nuclear configuration. As a result of this approximation, the wave function describing the system can be expressed as a single product of an electronic part and a nuclear part, where the electronic part depends parametrically on the nuclear positions. For many physical phenomena, this separation of the wave functions is extremely successful. However, it does not capture nonadiabatic (NA) effects [12–16]. Nonadiabatic effects are observed in the electronic degrees of freedom [16–19], e.g., in energy spectra, kinks in photoemission, charge transfer and localization, polaron formation [5–7,20,21], superconductivity [22,23], and latticeinduced decoherence [24,25], as well as in the phononic ones, for example, in frequency renormalization, Raman *These authors contributed equally to this work. †These authors share equal correspondence. Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. measurements, and x-ray spectra [12,26,27]. Accounting for the electron-phonon interaction in materials in a rigorous way, including nonadiabatic effects, is thus substantial for establishing a general theoretical framework. Yet, a full ab initio treatment of these phenomena remains a challenging task. The traditional and most direct way of treating electronphonon interactions starts from the BO approximation. An ab initio treatment of electrons and phonons within the BO approximation is achieved through density functional perturbation theory (DFPT) [19,28,29], by calculating the adiabatic response of the Kohn-Sham potential and BO potential energy surface to small displacements in the nuclear coordinates. It has been proposed [30] to extend DFPT by accounting for the electronic response to time-dependent nuclear displacements. A common starting point in many approaches to the calculation of electron-phonon interaction is the Fröhlich Hamiltonian, which comprises mean-field electrons, BO phonons, and the lowest order electron-phonon interaction. The key challenge is deducing a Fröhlich type Hamiltonian from the full ab initio Hamiltonian of electrons and nuclei. Several attempts have been made to overcome this difficulty [31,32]. A closed set of equations, called the Hedin-Baym equations, for electronic and phononic Green’s functions has been formulated by Guistino [19] starting from the seminal works of Baym and Hedin [33,34]. These equations have been further formulated by Stefanucci, van Leeuwen, and Perfetto for out-of-equilibrium systems, in terms of conserving diagrammatic expansions [9]. The work presented in this paper also starts from the full ab initio Hamiltonian. It then employs the exact factorization (EF) [35,36] of the many-electron many-nucleus wave function, which renders the formalism technically similar to the BO +DFPT treatment, but it remains formally exact, i.e., retains all nonadiabatic effects. This formalism was employed 2469-9950/2025/112(7)/075102(12) 075102-1 Published by the American Physical Society GALIT COHEN et al. PHYSICAL REVIEW B 112, 075102 (2025) to deduce a density functional framework for the complete system of electrons and nuclei [37,38]. For the case of small amplitude nuclear motion, this can be cast into a density functional theory (DFT) for electrons and phonons as recently derived [39]. In this approach the Kohn-Sham electrons are “dressed” by nonadiabatic phonon-induced interactions that appear through a nonadiabatic correction to the Kohn-Sham potential. The self-consistent density obtained is the conditional probability of finding an electron at position rgiven that the collection of nuclei resides at a certain point in nuclear configuration space. As a consequence, one can avoid transformation to the body-fixed coordinate frame [31,40], which is a big advantage in practical terms. In Ref. [39], the density functional formalism was applied to a simple Fröhlich model system as a test. We present here a full ab initio derivation of the EF-DFT equations for real materials by applying perturbation theory in the spirit of DFPT, including multiband transitions and higher-order electron-phonon interactions. We derive a theoretical framework within EF-DFT to account for nonadiabatic effects arising from these interactions. Our manuscript is organized as follows. We present the general EF-DFT theory for extended systems in Sec. II.We derive nonadiabatic DFPT-based corrections to the electronic energy bands and wave functions in Sec. III and discuss their implications in Sec. IV. II. EXACT FACTORIZATION APPROACH TO NONADIABATIC EFFECTS IN ELECTRON-LATTICE INTERACTIONS We begin with a general description of the full system via the many-body (MB) Hamiltonian ˆ H=ˆ Tn(R)+ˆ Wnn(R)+ˆ Te(r)+ˆ Wee(r)+ˆ Wen(R,r),(1) consisting of the kinetic and Coulomb interactions in the nuclear and electronic systems and the mutual interaction between electrons and nuclei. The sets of nuclear and electronic coordinates are denoted by Rand r, respectively. The complete many-electron many-nucleus wave function satisfies the time-independent Schrödinger equation (SE) ˆ H(r,R)=E(r,R),(2) where the eigenstates of the full Hamiltonian, the MB wave functions, are denoted by (r,R).1 Owing to the small ratio of electronic over nuclear mass, one intuitively expects that the electrons quickly adjust to the slow motion of the nuclear degrees of freedom. This intuitive picture suggests the adiabatic approximation for the electronnucleus wave function ad (r,R)=φBO α(r|R)χαλ(R).(3) Thus the MB wave function is approximated as a single product of separate electronic and nuclear wave functions; φBO αis the electronic wave function in the Born-Oppenheimer [10] approximation and χαλ(R) is the nuclear wave function. 1While the full Hamiltonian ˆ His translationally invariant, we assume eigenstates with spontaneously broken translational symmetry [41]. For each fixed nuclear configuration, the electronic BO wave function satisfies ˆ HBO(r,R)φBO α(r|R)=εα(R)φBO α(r|R),(4) with ˆ HBO(r,R)=ˆ Wnn(R)+ˆ Te(r)+ˆ Wee(r)+ˆ Wen(R,r)(5) being the Hamiltonian governing the correlated motion of electrons in the field of clamped nuclei. The nuclear wave functions χαλ(R) are the solutions of the nuclear SE (ˆ Tn(R)+εα(R))χαλ(R)=αλχαλ(R),(6) where the λindex denotes the vibrational states associated with each electronic eigenstate, α. The crucial simplification achieved by the adiabatic approximation is that nuclear and electronic degrees of freedom are treated in separate equations, Eqs. (4) and (6), where the solution of Eq. (4), namely the electronic eigenvalue εα(R), is used as input to the nuclear SE, Eq. (6). The nuclear kinetic energy reads ˆ Tn(R)=−μ· κ 1 2Mκ ∇2 Rκ,(7) where μ=m Mwith mbeing the electron mass and Mrepresenting a nuclear reference mass (typically chosen to be the proton mass) such that Mκ·Mare the true nuclear masses. Here and in the remainder of the article, all equations are written in atomic units (e2=¯h=m=1). While in Eq. (3) the nuclear wave function is the solution of Eq. (6), one may improve upon Eq. (3) by variationally optimizing the nuclear part. This leads to μ κ−i∇Rκ+AκφBO α(R)2 2Mκ +εα(R)+εgeoφBO α(R) טχαλ(R)=˜ αλ ˜χαλ(R),(8) where ˜χis the nuclear wave function that minimizes the total energy. In Eq. (8) the vector potential and geometric term are defined in terms of the general expressions Aκ[φα](R)=φα|−i∇Rκ|φα,(9) εgeo[φα](R)=μ κ 1 2Mκ [∇Rκφα|∇Rκφα−Aκ[φα](R)2] (10) evaluated for φα=φBO α. Note that εgeo carries a prefactor of the mass ratio. Recently, the exact factorization (EF) approach was introduced as a way to include nonadiabatic effects beyond the BO approximation [35,36,42]. This is achieved by expressing the MB wave function as a single product of electronic and nuclear exact wave functions (r,R)=φ(r|R)χ(R).(11) Note that this formally exact representation of the full wave function involves an electronic state which is different from the BO electronic wave function but shares with it the normalization condition |φ(r|R)|2dr =1. This differs from the Born-Huang expansion [43], in which multiple BO surfaces 075102-2 NONADIABATICITY FROM FIRST PRINCIPLES: … PHYSICAL REVIEW B 112, 075102 (2025) are directly accounted for, and from the adiabatic approximation in which a single BO surface is considered. Within the EF scheme, nonadiabatic effects are accounted for in a SE equivalent to Eq. (8), where the vector potential, Eq. (9), and the geometric term, Eq. (10), are functionals of the exact electronic wave function (rather than the BO wave function). The equations of motion of the two factors, φ(r|R) and χ(R), can be deduced from the Rayleigh-Ritz variational principle by minimizing the total energy functional E=|ˆ H| =χ∗(R)μ κ 1 2Mκ (−i∇Rκ+Aκ[φ](R))2χ(R) +|χ(R)|2(εBO[φ](R)+εgeo[φ](R)),(12) with εBO[φ](R)=dr φ∗(r|R)ˆ HBO(r,R)φ(r|R).(13) Using this variational principle, a density functional theory has been formulated in terms of the electronic conditional probability density, n(r|R)[37]. This leads to a Kohn-Shamlike equation for the electronic degrees of freedom, while the nuclear degrees of freedom are described by the N-body SE μ κ 1 2Mκ (−i∇Rκ+Aκ[φ](R))2+εBO[φ](R) +εgeo[φ](R)χ(R)=Eχ(R) (14) that follows directly from the EF. Hereafter, the vector potential is neglected, as justified for systems with vanishing nuclear current density in the ground state. Expanding εBO(R) and εgeo(R) to second order in the nuclear displacements from their equilibrium positions allows us to solve Eq. (14) in terms of harmonic phonon coordinates U={Uqλ}, with qand λbeing the phonon momentum and mode, respectively: μ qλ ˆ P−qλˆ Pqλ 2+ε(U)χ(U)=χ (U).(15) The electronic Kohn-Sham equation reads ˆp2 2m+ˆven(r,U)+ˆvhxc(r,U)+ˆvgeo(r,U)ψnkU(r) =nkUψnkU(r).(16) Equation (15) accounts for the nuclear interactions, where the first term is the nuclear kinetic energy, ε(U) is the exact potential energy surface, and =E−ε(R0) is the total energy of the phonons. R0is the equilibrium position corresponding to the bottom of the BO potential energy surface. The EF energy surface ε(U) is given by ε(U)=εBO(U)+εgeo(U),(17) where we emphasize that both εBO[φ] and εgeo[φ]inEqs.(13) and (10) are evaluated with the EF electronic wave function and therefore contain all powers of μ. In the present work, we approximate each term by its lowest order with respect to powers of μ. This corresponds to treating the εBO[φ]termby a standard BO-based DFT functional. For the εgeo term, which introduces the nonadiabatic “dressing” of the electron-phonon interaction, we adopt a determinantal approximation for φ. Returning to Eq. (16), we observe that it has the same form as the standard Kohn-Sham equation, where ˆp2/2mis the electronic kinetic energy, ˆven is the electron-lattice interaction potential, and ˆvhxc is the Hartree-exchange-correlation (HXC) potential. Nonadiabaticity is introduced through an additive potential, ˆvgeo, derived from the functional derivative of the εgeo term in the total energy, Eq. (12). In the original EF-DFT [37] framework ˆvgeo, like ˆvhxc, is a local multiplicative potential. Since in the present work we make the determinantal approximation for φ,ˆvgeo acts as a nonlocal potential, in the spirit of the generalized Kohn-Sham (GKS) approach [44], and is a functional of conditional single particle orbitals ψnkU(r,R); n,krepresent an electronic state with band index nand quasimomentum k. Hereafter, the subscript Uon the energies and orbitals is suppressed. Making the orbital-dependent energy functional stationary with respect to variations of the GKS orbitals, subject to orthonormality constraints imposed through a matrix of Lagrange multipliers, gives the matrix elements of ˆvgeo [39]as ψm,k+q|ˆvgeo|ψnk =ψm,k+q δεgeo δψ∗ nk−1 |χ|2 qλψm,k+q ∂ ∂U∗ qλ ×|χ|2δεgeo δ(∂ψ∗ nk/∂U∗ qλ).(18) Alternatively, the orbitals could be restricted to come from a local potential so that minimizing the total energy with respect to this local potential will lead to an optimized effective potential (OEP)-like approach [36]. It is useful to quantify nonadiabatic effects in terms of the powers of μ. The potential energy surface, being treated as a density functional, generates the GKS potential through functional differentiation. We see in Eqs. (17), (10), and (13) that the potential energy surface contains a zeroth-order term εBO that has no explicit μdependence and a first-order term εgeo that has a prefactor of μ. As a consequence, the GKS potential in Eq. (16) inherits zeroth-order terms ˆven and ˆvhxc and a first-order term given by ˆvgeo. To account for the mutual influence of electrons and phonons, in principle this set of equations must be solved self-consistently. As a first step, in the next section we adapt DFPT to the EF equations. III. NONADIABATIC PERTURBATIVE TREATMENT OF ELECTRON-PHONON INTERACTIONS In order to apply EF-DFT in realistic materials, we now turn to the derivation of the additive potential ˆvgeo. Here, we present a perturbative approach, in which corrections to the KS wave functions and the energy surface are introduced within firstand second-order perturbation theory. 075102-3 GALIT COHEN et al. PHYSICAL REVIEW B 112, 075102 (2025) The single-particle Hamiltonian can be written as the sum of a noninteracting part and an interaction term, ˆ H=ˆ H(0) s+ ˆ H. The noninteracting Hamiltonian corresponds to the KS Hamiltonian within BO ˆ H(0) s=ˆp2 2m+ˆv(0) en (r)+ˆv(0) hxc(r)≡ˆp2 2m+ˆv(0) s(r).(19) Here ˆv(0) en (r)=ˆven(r,U)U=0is the electron-lattice Coulomb interaction evaluated at R0,ˆv(0) hxc(r)=ˆvhxc(r,U)U=0corresponds to the standard HXC potential within BO, and ψ(0) nk(r) and (0) nk(r) are, respectively, the KS eigenstates and eigenenergies that solve the unperturbed KS equation. We define the perturbation ˆ Has ˆ H=ˆv s+ˆvgeo (20) to second order in the phonon displacements, Uqλ. Here, ˆv s=ˆv(1) s+ˆv(2) sis the sum of the firstand second-order perturbation terms implemented in standard DFPT [19,28,29] and ˆvgeo =ˆv(1) geo +ˆv(2) geo is the nonadiabatic correction. In the harmonic approximation, the nuclear wave functions are localized and have a zero-point amplitude of the order of μ1/4. Hence we can rationalize the orders of perturbation by recognizing that the characteristic nuclear displacement Uqλis of order μ1/4[10]. The superscripts correspond to powers in the displacements. Crucially, we have two separate sources of μdependence: on the one hand, the μdependence of the functional and, on the other hand, the μdependence arising from the powers in the displacements U. The latter is already present in standard DFPT. We follow the first-order corrections previously derived in Ref. [39] for a single band case and extend them for the manybands, many-mode case. Here, using the notations 1 ≡n1,k1 and 2 ≡n2,k2,ˆv(1) sis defined as ψ(0) 2ˆv(1) sψ(0) 1≡ λ gn2,n1,λ(k1,q)Uqλ Lqλ ,(21) where k2=k1+qas required by momentum conservation. Lqλis the amplitude of the zero-point motion, defined as Lqλ=¯h 2Mωqλ, and gn2,n1,λ(k1,q) is the first-order electronphonon coupling term representing interactions between two electrons through a phonon. We define the first-order correction for ˆv(1) geo using a Taylor expansion of the electronic wave functions. Following Eq. (18), this leads to ψ(0) 2ˆv(1) geoψ(0) 1= λ ¯hωqλUqλψ2 ∂ψ(1) 1 ∂Uqλ = λ gn2,n1,λ(k1,q)Uqλ Lqλ ¯hωqλ (0) 1−(0) 2±¯hωqλ , (22) where ¯hωqλis the phonon energy, accounting for both phonon emission and absorption channels. We note that these corrections are off-diagonal in the wave function basis, since they mix KS electron states upon phonon scattering. Since the first-order energy correction is defined as a diagonal term, and ˆv(1) sand ˆv(1) geo are off-diagonal by definition, this energy term vanishes, i.e., (1) nk=0. The first-order correction to the wave functions follows the form ψ(1) 1= 2=1ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 1 (0) 1−(0) 2ψ(0) 2 = λ 2=1 gn2,n1,λ(k1,q) (0) 1−(0) 2±¯hωqλ Uqλ Lqλψ(0) 2.(23) Since we consider systems at zero temperature, only the phonon emission case appears hereafter. These corrected wave functions can then be used to evaluate the changes in the density due to interaction with phonons. The electronic density, nU, is defined as nU(r)=nkfnkU|ψnkU|2, where fnkUis the Fermi-Dirac occupation function for zero temperature. Hence the first-order density correction is n(1)(r)=2Re λ 1,2=1 gn2,n1,λ(k1,q) (0) 1−(0) 2−¯hωqλ Uqλ Lqλ ψ(0)∗ 1(r) ×ψ(0) 2(r).(24) The general formula for the second-order energy correction is given by (2) 1= 2=1ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 1 2 (0) 1−(0) 2 +ψ(0) 1ˆv(2) s+ˆv(2) geoψ(0) 1.(25) Here, the two-phonon-induced potentials ˆv(2) sand ˆv(2) geo appear. While the second-order correction to the energy eigenvalue (2) 1only requires knowledge of the diagonal matrix elements of these potentials, since we ultimately seek second-order corrections to other observables such as the density, we proceed to evaluate all of the matrix elements of ˆv(2) sand ˆv(2) geo.In doing so, we need to take into account the gauge freedom in the choice of the electronic states; that is, the determinant of the occupied KS orbitals is invariant (up to a phase) to a unitary transformation among the occupied orbitals. Gonze (1995) [28] has found a convenient gauge choice, called the “parallel-transport” gauge, for carrying out higherorder perturbative calculations. Such a choice is specified through orthonormality constraints [as depicted in Eqs. (A4) and (A5)], which allow the calculations to be made with only the knowledge of the block off-diagonal components of the perturbations, i.e., the matrix elements between one occupied and one unoccupied state. As a second step, the results can be transformed back into the “diagonal gauge”—the familiar gauge choice in which the Lagrange multiplier matrix is diagonal with energy eigenvalues 1on the diagonal—producing corrections to the KS orbitals and eigenenergies which are consistent with the first-order corrections above. The matrix elements of the first-order perturbations ˆv(1) s,ˆv(1) geo were addressed above. The matrix elements of second-order perturbations consist of the phonon-induced deviations of the KS potential and the nonadiabatic effect accounted for by ˆvgeo. Starting from the diagonal elements of 075102-4 NONADIABATICITY FROM FIRST PRINCIPLES: … PHYSICAL REVIEW B 112, 075102 (2025) ˆv(2) sand ˆv(2) geo, required for the evaluation of Eq. (25), we have V(2) 1,1≡ψ(0) 1ˆv(2) sψ(0) 1= qλ g(2) n1,n1;λ,λ(k1,q,¯ q)|Uqλ|2 L2 qλ , (26) where V(2) 1,1is the potential perturbation that accounts for the two-phonon process originating from the Debye-Waller (DW) interaction g(2) n1,n1;λ,λ(k1,q,¯ q), which can be obtained in adiabatic DFPT calculations [1,19,45–49]. The expression for the ˆv(2) geo matrix element, which follows from Eq. (18) in the parallel-transport gauge (derivation elaborated in Appendix A), amounts to W(2) 1,1≡ψ(0) 1ˆv(2) geoψ(0) 1 = unocc  2 λλ g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ ×δλλ−U¯ qλ L¯ qλ Uqλ Lqλ − unocc  2 λλ g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ ×U¯ qλ L¯ qλ Uqλ Lqλ ,(27) where q=k2−k1and q=−q. When the quantity in Eq. (27) is multiplied by |χ(Uqλ)|2 and averaged over the phonon amplitudes, the first line cancels out. Using the definition for the matrix elements in Eqs. (21), (22), (26), and (27) and transforming the secondorder eigenenergy correction to the diagonal gauge, as shown in Appendix B, the energy correction results in (2) 1= 2=1 λλ|gn2,n1,λ(k1,q)|2 (0) 1−(0) 2−¯hωqλ |Uqλ|2 L2 qλ +g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ ×δλλ−U¯ qλ L¯ qλ Uqλ Lqλ + qλ g(2) n1,n1;λ,λ(k1,q,¯ q)|Uqλ|2 L2 qλ .(28) Note that conventionally the electronic band structure is evaluated for fixed nuclear positions (within the BO approximation); however, in the current framework it explicitly depends on the displacement. In the context of EF any purely electronic observable naturally comes out as a |χ(Uqλ)|2average. To make contact with the standard Fan-Migdal (FM) and DW result, we average over Uqλusing |χ(Uqλ)|2as a weighting function. As before, the middle term of Eq. (28) cancels out. The average energy renormalization is ¯(2) 1= 2=1 λ |gn2,n1,λ(k1,q)|2 (0) 1−(0) 2−¯hωqλ + qλ g(2) n1,n1;λ,λ(k1,q,¯ q). (29) The first term consists of the electron–one-phonon interaction and the associated nonadiabatic corrections. This reproduces the energy renormalization due to the well-known FM selfenergy [19,39,50], dressed with nonadiabatic interactions. The second term consists of the electron–two-phonon interaction. In Eqs. (28) and (29) there are second-order terms with nonadiabatic corrections entering through the dependence on phonon frequencies, as well as additional dependence on |Uqλ|2/L2 qλ. On the other hand, at this order there are no nonadiabatic corrections to the last term of Eq. (29) originating from the DW interaction. However, there will be corrections to the off-diagonal DW terms, as we discuss next. We continue to the evaluation of the off-diagonal elements of ˆv(2) sand ˆv(2) geo, which are required for the second-order wave function corrections. The expression for ˆv(2) selements follows: V(2) 2,1≡ψ(0) 2ˆv(2) sψ(0) 1 = qλλ g(2) n2,n1,λ,λ(k1,q,q)Uqλ Lqλ Uqλ Lqλ ,(30) where q=k2−k1−q,1∈occupied states subspace, and 2∈unoccupied states subspace. The expression for the ˆv(2) geo off-diagonal matrix element in the parallel-transport gauge allows for the construction of a differential equation for W(2) 2,1 (derivation elaborated in Appendix A), W(2) 2,1≡ψ(0) 2ˆv(2) geoψ(0) 1= qλ ¯hωqλUqλψ(0) 2 ∂ψ(2) 1 ∂Uqλ −δk2,k1 qλ ¯hωqλLqλL¯ qλψ(0) 2 ∂ψ(2) 1 ∂U¯ qλ Uqλ.(31) With the evaluation of ψ(0) 2|ψ(2) 1from the second-order Sternheimer equation with the additional nonadiabatic potential ˆv(2) geo and using the notation of V(2) 2,1,W(2) 2,1for the second-order matrix elements, the differential equation becomes W(2) 2,1= qλ ¯hωqλUqλ 1 (0) 1−(0) 2 ∂V(2) 2,1+W(2) 2,1 ∂Uqλ −δk2k1 qλ ¯hωqλLqλL¯ qλ 1 (0) 1−(0) 2 ∂V(2) 2,1+W(2) 2,1 ∂U¯ qλ Uqλ. (32) The derivation from here to the explicit form of the differential equation is given in Appendix A, leading to W(2) n2k2,n1k1= qλλ ¯hωqλ (0) n1k1−(0) n2k2Uqλ Lqλ Uqλ Lqλ −δk2k1δλλ ×g(2) n2,n1,λ,λ(k1,q,q)+g(2) n2,n1,λ,λ(k1,q,q) + qλ ¯hωqλ (0) n1k1−(0) n2k2 Uqλ ∂W(2) n2k2,n1k1 ∂Uqλ −δk2k1 qλ ¯hωqλ (0) n1k1−(0) n2k2 LqλL¯ qλ ∂2W(2) n2k2,n1k1 ∂U¯ qλ∂Uqλ , (33) 075102-5 GALIT COHEN et al. PHYSICAL REVIEW B 112, 075102 (2025) where q=k2−k1−q. The general solution of ˆv(2) geo can include nonanalytical terms or polynomials that reflect nonanalytic behavior and meaningful physical constraints. Here we introduce an analytical expression and reduce our solution to a particular case. A particular solution to the differential equation in terms of standard ab initio calculated quantities is W(2) n2k2,n1k1= qλλ 1 2g(2) n2,n1,λ,λ(k1,q,q)+g(2) n2,n1,λ,λ(k1,q,q) ׯhωqλ+¯hωq,λ (0) n1k1−(0) n2k2−¯hωqλ−¯hωqλ ×Uqλ Lqλ Uqλ Lqλ −δk1k2δλλ.(34) Finally, we demonstrate the orbital changes in this level of approximation. The second-order wave function correction in the parallel-transport gauge reads |ψ(2) 1,||=− 1 2 occ  3=1 unocc  2=1ψ(0) 3ˆv(1) s+ˆv(1) geoψ(0) 2 ((0) 3−(0) 2) ×ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 1 (0) 1−(0) 2ψ(0) 3 − unocc  3=1ψ(0) 3ˆv(2) s+ˆv(2) geoψ(0) 1 (0) 3−(0) 1ψ(0) 3 −1 2 unocc  2=1ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 1 2 (0) 2−(0) 12ψ(0) 1.(35) The correction to the single-particle KS states in the diagonal gauge is evaluated as ψ(2) 1,d=ψ(2) 1,||+ occ  3 U(1)∗ 13 ψ(1) 3,||+ occ  3 U(2)∗ 13 ψ(0) 3,(36) following Gonze [28] and the derivation in Appendix C. This shows that second-order wave function corrections account for modifications in the electronic KS states due to the two-phonon interaction. It also marks a step toward a self-consistent description of the solution, including mutual NA phonon-induced interactions between multiple electronic states. IV. DISCUSSION In the EF-DFPT method presented here we achieve explicit expressions correcting the KS wave functions and band structure to second order in the displacement due to nonadiabatic effects. The current method is rigorously derived from the full Hamiltonian of electrons and nuclei, in contrast to the commonly used model or Fröhlich Hamiltonian starting point. Thus in principle it accounts for all nonadiabatic effects. There is further sensitivity in the choice of the perturbative treatment applied to the Hamiltonian. Here we expand all quantities of the GKS Hamiltonian and the nuclear Hamiltonian systematically in powers of the electron over nuclei mass ratio. In principle, all orders of the mass ratio in the EF potential can be accounted for. The EF-KS Hamiltonian explicitly depends on the nuclear wave function through the introduction of ˆvgeo, which is a higher order term in the mass ratio compared to the standard KS potential. The corrected wave functions are conditional electronic states. For localized nuclear wave functions, the characteristic zero-point amplitude is of the order μ1/4, which implies that the nuclear displacements should be assigned to the same order in standard DFPT. This is a way for the mass dependence to enter in DFPT because the BO potential energy surface and the KS potential have no explicit dependence on the mass ratio. Therefore, DFPT introduces some nonadiabatic effects, whereas EF has systematically all orders of mass ratio. We emphasize that the corrections to the electronic wave functions derived above are in fact corrections to the conditional quantities, i.e., the KS orbitals for a given nuclear configuration. Consequently, the density evaluated from these corrected wave functions is also the conditional one. An important feature of the presented approach is the ability to quantify the nonadiabatic contribution in the corrections explicitly. This is reflected in the conditional wave functions and density as well as in any observable evaluated as a function of these quantities, such as the transition dipole moment in optical absorption, the total energy, electron and spin polarization, the dielectric screening function, and the electronic self-energy, extending beyond their evaluation at the DFT level, as well as in the update of the nuclear degrees of freedom. It is possible to relate the band structure corrections to self-energy corrections in Green’s function based methods and associate the energy corrections derived here with a nonadiabatic treatment of the Allen-Heine-Cardona (AHC) theory [47]. The relation of the wave function corrections to Green’s function based approaches is less straightforward. Although both the KS wave functions and Green’s functions can be used for the evaluation of the electronic density, it is crucial to stress that in the current context the density evaluated is the conditional one. Hence a potential mapping between the different approaches requires introduction of nuclear configuration dependence into the electronic Green’s function or potentially casting the EF formalism into a many-body Green’s function theory [51]. Such a mapping may benefit from the computation of the corrections presented above. We further note that, in the current derivation, temperature effects are not taken into account. These can be added through the electron Fermi-Dirac occupation function, fnk, and the phonon Bose-Einstein occupation function, Nqλ, e.g., in the density correction, Eq. (24)[52]. The theoretical framework presented here is designed in principle for ground-state properties and the introduction of thermal excitation as reflected in finite temperature should be addressed carefully. The extension of this approach to include thermal excitations and fluctuations in the occupation in a time-dependent formalism may provide an insight into nonadiabatic dynamical properties in the system. V. CONCLUSIONS We have developed an ab initio EF-based framework for the evaluation of the nonadiabatic signatures of electronlattice interaction in real materials. It allows for the explicit 075102-6 NONADIABATICITY FROM FIRST PRINCIPLES: … PHYSICAL REVIEW B 112, 075102 (2025) identification of these signatures in the electronic conditional band structure, wave function, density, and properties derived from them, framing them in terms of DFT/DFPT calculated quantities. Thus these EF-based perturbative corrections can be evaluated as a postprocessing step to standard DFT/DFPT codes and their code implementation will be presented in a follow-up work. Eventually, the perturbative corrections may be used iteratively to address the self-consistent solution of the nuclear and the electronic set of equations. Thus this work takes a first step towards the full EF treatment of the electron-lattice many-body problem in materials. ACKNOWLEDGMENTS G.C. acknowledges support from an Institute for Environmental Sustainability (IES) Fellowship. E.K.U.G. acknowledges support as a Mercator Fellow within SFB 1242 at the University Duisburg-Essen. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, Grant Agreements No. ERC-2017-AdG-788890 and No. ERC-2022-StG-101041159 and from the Israel Science Foundation, Grant Agreement No. 1208/19. DATA AVAILABILITY No data were created or analyzed in this study. APPENDIX A Below we derive the ˆv(2) geo matrix elements, using the notations 1 ≡n1,k1and 2 ≡n2,k2. The expression for ˆvgeo is given in Eq. (18) with the orbital-dependent functional [39], εgeo =¯h2 2M qλ occ  1∂ψ1 ∂Uqλ (1 −Pv) ∂ψ1 ∂Uqλ,(A1) where (1 −Pv)=Pc=unocc 3|ψ3ψ3|is the projector onto the unoccupied subspace. By taking the functional derivative of the εgeo term in the total energy, the ˆvgeo matrix element given in Eq. (18) becomes ψ2|ˆvgeo|ψ1 =− ¯h2 2M qλ occ  3ψ2 ∂ψ3 ∂Uqλ∂ψ3 ∂Uqλ ψ1(i) −¯h2 2M qλ occ  1 ∂ln |χ|2 ∂U∗ qλψ2 (1 −Pv) ∂ψ1 ∂Uqλ(ii) −¯h2 2M qλ occ  1ψ2 (1 −Pv) ∂2ψ1 ∂U∗ qλ∂Uqλ(iii) +¯h2 2M qλ∂ψ2 ∂Uqλ ∂ψ1 ∂Uqλ(iv) +¯h2 2M qλ occ  3ψ2 ∂ψ3 ∂U∗ qλψ3 ∂ψ1 ∂Uqλ(v).(A2) The second-order contribution to the matrix element of ˆvgeo is given by the sum of three terms. The first is the matrix element of the second-order component of ˆv(2) geo between the unperturbed wave functions. The second and third are the matrix elements of the first-order component ˆv(1) geo between an unperturbed state and the first-order corrected state: ψ2|ˆvgeo|ψ1(2) =ψ(0) 2ˆv(2) geoψ(0) 1+ψ(0) 2ˆv(1) geoψ(1) 1+ψ(1) 2ˆv(1) geoψ(0) 1. (A3) The perturbed KS states are subjected to orthonormalization conditions, ψ(0) 2ψ(i) 1+ψ(i) 2ψ(0) 1=−i−1 j=1ψ(j) 2ψ(i−j) 1,i>1, 0,i=1, (A4) where iis the order of the perturbation and 1,2 are either both occupied or both unoccupied states. The additional constraints defining the parallel-transport gauge [28]are ψ(0) 2ψ(i) 1−ψ(i) 2ψ(0) 1=0.(A5) 1. Diagonal elements of ˆ v(2) geo Using this gauge choice and expanding the terms in Eq. (A2)forψ1|ˆvgeo|ψ1, up to second order in the displacements Uqλ, it is possible to show that terms (i) and (v), and terms (ii), (iii), do not contribute at the second order. Thus this choice of gauge allows a convenient representation of ˆv(2) geo.To show this explicitly, we first focus on (1) which can be written as (i)=− ¯h2 2M qλ occ  3∂ψ1 ∂Uqλ ψ3ψ3 ∂ψ1 ∂Uqλ.(A6) Expanding the last factor we obtain ψ3 ∂ψ1 ∂Uqλ=ψ(0) 3 ∂ψ(0) 1 ∂Uqλ+ψ(0) 3 ∂ψ(1) 1 ∂Uqλ+O(U2). (A7) In this term 1,3∈occupied subspace. The derivative of the unperturbed wave function with respect to the displacement is zero; hence the first term vanishes. Due to the orthogonality of the basis and the parallel-transport gauge, Eqs. (A4) and (A5), it follows that the second term vanishes as well and this factor results in terms of O(U2). Therefore, since the first factor ∂ψ2 ∂Uqλ|ψ3is also O(U2), the (i) term is of order O(U4). Similarly, (v)isalsooforderO(U4) and therefore they do not contribute at second order. The diagonal elements in terms (ii) and (iii) consist of ψ1|(1 −Pv)for1∈occupied subspace. The projector, as well as ψ1|, contain all orders of the perturbation. Therefore, (ii) and (iii) vanish due to orthonormality. The only contribution to the second-order diagonal elements comes from (iv). After expanding the wave functions in (iv) and considering the gauge choice, the contributing elements are constructed from the first-order perturbed wave 075102-7 GALIT COHEN et al. PHYSICAL REVIEW B 112, 075102 (2025) function ψ(1) 1[expressed in Eq. (23)]. Thus ψ1|ˆvgeo|ψ1(2) =¯h2 2M qλ∂ψ(1) 1 ∂Uqλ ∂ψ(1) 1 ∂Uqλ = λ unocc  2 g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hωqλ2. (A8) Note that the degree of freedom of qis accounted for by the sum over 2 ≡(n2k2), with k2=k1+q. From here we can derive the diagonal matrix element W(2) 11 ≡ψ(0) 1|ˆv(2) geo|ψ(0) 1as W(2) 11 =ψ1|ˆvgeo|ψ1(2) −ψ(0) 1ˆv(1) geoψ(1) 1−ψ(1) 1ˆv(1) geoψ(0) 1 = λ unocc  2 g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hωqλ2 − λλ unocc  2 g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)(¯hω¯ qλ+¯hωqλ) (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ ×U¯ qλ L¯ qλ Uqλ Lqλ .(A9) Since the factor g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q) (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ U¯qλ L¯ qλ Uqλ Lqλ (A10) is symmetric with respect to interchanging qλwith ¯ qλwe finally get the diagonal elements of ˆv(2) geo as in Eq. (27). 2. Off-diagonal elements of ˆ v(2) geo To work out the off-diagonal terms of ˆv(2) geo we turn back to Eq. (A2) in its off-diagonal form and derive it to second order. From the same considerations elaborated above terms (i) and (v) do not contribute at this order. Term (iv), which contributes to the diagonal part of ˆv(2) geo, does not contribute to the off-diagonal elements due to the gauge choice. The remaining terms are (ii) and (iii). As before, these terms consist of ψ2|(1 −Pv), where here 2 ∈unoccupied subspace, in contrast to the diagonal case. Therefore, these terms do not vanish. In the current case ψ2|(1 −Pv) amounts to the unperturbed wave function ψ(0) 2|. The second-order contribution of (ii) results in −¯h2 2M qλ ∂ln |χ|2 ∂U∗ qλψ(0) 2 ∂ψ1 ∂Uqλ =¯h2 2M qλ Uqλ L2 qλ ∂ψ(0) 2ψ(2) 1 ∂Uqλ = qλ ¯hωqλUqλ ∂ψ(0) 2ψ(2) 1 ∂Uqλ .(A11) Similarly, the second-order contribution from (iii) follows −¯h2 2M qλ unocc  n2 δk2,k1 ∂ψ(0) 2ψ(2) 1 ∂U∗ qλ Uqλ.(A12) These contributions add up to form the off-diagonal ψ2|ˆvgeo|ψ1(2). We can deduce the matrix element of W(2) 21 following Eq. (A3) where ψ(0) 2|ˆv(1) geo|ψ(1) 1=ψ(1) 2|ˆv(1) geo|ψ(0) 1= 0 due to the block-off diagonality of the matrix elements in the parallel-transport gauge. Therefore, we get a differential equation, Eq. (31). Here we adopt a verbose indexing notation, denoting the orbitals explicitly by band and momenta, niki. To evaluate the first and second derivatives in Eq. (31) we first evaluate ψ(0) n2k2|ψ(2) n1k1from the nonadiabatic second-order Sternheimer equation Pch(0) s−(0) nkPcψ(2) nk =−Pcˆv(2) s+ˆv(2) geoψ(0) nk−Pcˆv(1) s+ˆv(1) geoψ(1) nk.(A13) For nk∈occupied subspace the consideration of the gauge choice eliminates the second term on the right hand side of the last equation. Projecting the equation on |ψ(0) n2k2, we get ψ(0) n2k2ψ(2) n1k1=V(2) n2k2,n1k1+W(2) n2k2,n1k1 (0) n1k1−(0) n2k2 ,(A14) with V(2) n2k2,n1k1,W(2) n2k2,n1k1defined above. Substitution of Eq. (A14)inEq.(31) leads to the differential equation determining W(2) n2k2,n1k1,Eq.(32). In order to get an explicit form of this partial differential equation, one needs to work out the derivatives of V(2) n2k2,n1k1. The V(2) n2k2,n1k1term is given in Eq. (30). The derivative with respect to the phonon displacements Uqλthen follows ∂V(2) n2k2,n1k1 ∂Uqλ = λg(2) n2,n1,λ,λ(k2,q,q)Uqλ Lqλ 1 Lqλ +g(2) n2,n1,λ,λ(k1,q,q)Uqλ Lqλ 1 Lqλ,(A15) where q=k2−k1−q. We follow with a second derivation with respect to U¯ qλ, where U¯ qλ=U∗ qλ=U−qλ. This requires k2=k1. Hence ∂V(2) n2k2,n1k1 ∂U¯ qλ Uqλ=δk1,k2 1 L−qλ 1 Lqλg(2) n2,n1,λ,λ(k1,q,¯ q) +g(2) n2,n1,λ,λ(k1,¯ q,q).(A16) In light of these derived expressions, Eqs. (A15) and (A16), the differential equation in Eq. (32) leads to Eq. (33). The solution to the differential equation (33) corresponds to the sum of a multiple of the homogeneous solution and the particular one. We distinguish two cases: one where k1=k2, leading to a second-order partial differential equation, and the other where k1= k2, amounting to a first-order differential equation. The particular solution for W(2) n2k2,n1k1in Eq. (34) contains both the k1=k2and k1= k2cases. APPENDIX B The transformation of the eigenenergies (Lagrange multipliers) from the parallel-transport gauge [Eqs. (A4) and (A5)] to the diagonal gauge is done via a unitary transformation matrix [28] and is shown below. The eigenvalues in the 075102-8 NONADIABATICITY FROM FIRST PRINCIPLES: … PHYSICAL REVIEW B 112, 075102 (2025) parallel-transport gauge for the firstand second-order corrections are (1) ||,2,1=ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 1(B1) and (2) ||,2,1=ψ(0) 2ˆv(2) s+ˆv(2) geoψ(0) 1 +ψ(1) 2ˆv(1) s+ˆv(1) geoψ(0) 1+ψ(0) 2ˆv(1) s+ˆv(1) geoψ(1) 1 +ψ(1) 2ˆ H(0) sψ(1) 1−1 2ψ(1) 2ψ(1) 1(0) 2+(0) 1, (B2) respectively. The transformation from the parallel-transport gauge to the diagonal gauge for the eigenvalues follows: (2) d,1=(2) ||,1,1− occ  3(1) ||,3,1 2 (0) 3−(0) 1 (B3) and by substituting Eqs. (B1) and (B2) it becomes (2) d,1= occ  2=1ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 1 2 (0) 1−(0) 2 (a) +ψ(0) 1ˆv(2) s+ˆv(2) geoψ(0) 1(b) +ψ(1) 1ˆv(1) s+ˆv(1) geoψ(0) 1(c) +ψ(0) 1ˆv(1) s+ˆv(1) geoψ(1) 1(d) +ψ(1) 1ˆ H(0) sψ(1) 1(e) −1 2ψ(1) 1ψ(1) 1(0) 1+(0) 1(f).(B4) In order to get an analytical expression for (2) d,1, we simplify the terms. Using the definitions above for Eq. (23)inthe parallel-transport gauge, the sum of terms (c)–( f) becomes (c)+(d)+(e)+(f)= unocc  2=1ψ(0) 1ˆv(1) s+ˆv(1) geoψ(0) 2 2 (0) 1−(0) 2 . (B5) The second-order correction to the energies in the diagonal gauge follows the addition of all the terms, resulting in Eq. (25). The first term in Eq. (25) consists of first-order matrix elements that are given in Eqs. (21) and (22). The second term in Eq. (25) consists of the diagonal elements of the second-order matrix elements as in Eq. (26) and Eq. (27). Hence the resulting term for the second-order correction to the eigenenergies in the diagonal gauge becomes (2) d,1= 2=1 λλ g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)(0) 1−(0) 2 (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ U¯ qλ L¯ qλ Uqλ Lqλ + qλ g(2) n1,n1,λ,λ(k1,q,¯ q)|Uqλ|2 L2 qλ + 2=1 λλ g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλδλλ−U¯ qλ L¯ qλ Uqλ Lqλ − 2=1 λλ g∗ n2,n1,λ(k1,q)gn2,n1,λ(k1,q)¯hωqλ (0) 1−(0) 2−¯hω¯ qλ(0) 1−(0) 2−¯hωqλ U¯ qλ L¯ qλ Uqλ Lqλ .(B6) The sum of these leads to the more compact expression in Eq. (28). APPENDIX C Following the gauge transformation [28] the second-order correction to the orbitals in the diagonal gauge is given in Eq. (36). The first term in Eq. (36) is the second-order correction to the wave function in the parallel-transport gauge, Eq. (35). The following terms consist of the unitary transformation matrix, U13, to first and second order in the perturbation, U(1)∗ 13 =0,1=3, −(1) ||,31 (0) 3−(0) 1 ,1= 3,(C1) U(2)∗ 13 =⎧ ⎨ ⎩ −1 2occ 2=3U(1) 23  2,1=3, −1 (0) 3−(0) 1occ 2U(1)∗ 12 (1) ||,32+(2) ||,31 −(1) d,1U(1)∗ 13 ,1= 3.(C2) These contain the firstand second-order Lagrange multipliers, as given in Eq. (B1) and Eq. (B2), respectively. We continue by evaluating the second term in Eq. (36) and substitute Eqs. (C1), (B1), and (23) in the parallel gauge: occ  3 U(1)∗ 13 ψ(1) 3,||= occ  3=1−(1) ||,31 (0) 3−(0) 1ψ(1) 3,||=− occ  3=1 unocc  2ψ(0) 2ˆv(1) s+ˆv(1) geoψ(0) 3 (0) 3−(0) 1ψ(0) 3ˆv(1) s+ˆv(1) geoψ(0) 1 (0) 3−(0) 2ψ(0) 2.(C3) 075102-9