On the mass of atoms in molecules: Beyond the Born-Oppenheimer approximation

Describing the dynamics of nuclei in molecules requires a potential energy surface, which is traditionally provided by the Born-Oppenheimer or adiabatic approximation. However, we also need to assign masses to the nuclei. There, the Born-Oppenheimer picture does not account for the inertia of the electrons and only bare nuclear masses are considered. Nowadays, experimental accuracy challenges the theoretical predictions of rotational and vibrational spectra and requires to include the participation of electrons in the internal motion of the molecule. More than 80 years after the original work of Born and Oppenheimer, this issue still is not solved in general. Here, we present a theoretical and numerical framework to address this problem in a general and rigorous way. Starting from the exact factorization of the electron-nuclear wave function, we include electronic effects beyond the Born-Oppenheimer regime in a perturbative way via position-dependent corrections to the bare nuclear masses. This maintains an adiabatic-like point of view: the nuclear degrees of freedom feel the presence of the electrons via a single potential energy surface, whereas the inertia of electrons is accounted for and the total mass of the system is recovered. This constitutes a general framework for describing the mass acquired by slow degrees of freedom due to the inertia of light, bounded particles. We illustrate it with a model of proton transfer, where the light particle is the proton, and with corrections to the vibrational spectra of molecules. Inclusion of the light particle inertia allows to gain orders of magnitude in accuracy.


I. INTRODUCTION
The Born-Oppenheimer (BO) [1], or adiabatic, treatment of the coupled motion of electrons and nuclei in molecular systems is among the most fundamental approximations in condensed matter physics and chemical physics. Based on the hypothesis that part of the system, usually electrons or protons, evolves on a much shorter time-scale than the rest, i.e. (heavy) nuclei or ions, the BO approximation allows one to visualize molecules as a set of nuclei moving on a single potential energy surface that represents the effect of the electrons in a given eigenstate. Yet, it is an approximation, yielding the correct dynamics only in the limit of infinite nuclear masses. Consequently when compared to highly accurate molecular spectroscopy measurements, theoretical predictions might deviate from experimentally observed behavior.
In those situations, the question of which masses [2][3][4][5] are to be considered when calculating rotational and vibrational spectra of light molecules, for instance hydrogenbased [6][7][8][9][10][11][12], often appears in the literature to rationalize this problem. In the BO approximation, the electrons appear only implicitly in the dynamics, as a potential energy contribution to the Hamiltonian driving the motion of the nuclei. The kinetic energy arising from the molecular motion then involves only the bare nuclear masses. However, electrons are carried along with the nuclei, thus how is their inertia accounted for? It has been proposed that more accurate results are obtained when employing atomic masses rather than bare nuclear masses [13].
The measured ro-vibrational spectrum of hot water in sunspots, for example, is very dense, with about 50 lines per wavenumber [14]. However, their assignment can not be performed at the BO level, either using nuclear or atomic masses, because of the lack of accuracy. Adding half an electron mass to the proton to effectively include non-adiabatic effects has been shown to lead to better results [15]. Such fractional masses account for the bond ionicity but there is no systematic way to include such corrections.
One solution to the problem is to perform a full non-adiabatic treatment of the coupled electron-nuclear problem, but the numerical cost is much larger than a BO calculation. Also, from a fundamental point of view, this does not answer the question of what is the mechanism by which the inertia of the electrons affects the mass of the heavy degrees of freedom. An alternative approach, pioneered by Bunker and Moss [6,7,16], is to treat perturbatively non-adiabatic effects, but applications are still limited to di-and tri-atomic molecules. In connection to the perturbation idea of Bunker and Moss, accurate numerical calculations have been performed on small molecules, like H 2 , D 2 , HD, H + 3 [17][18][19][20][21]. However, despite the effort to push forth the applications, it seems that the basic formalism still represents a major obstacle for the treatment of molecular systems comprising more than three atoms. The main reason is to be found in the use of internal coordinates, obtained after separation of the rotational and vibrational degrees of freedom of the center of mass of the molecule, as starting point for the application of the perturbation approach. In internal coordinates, the Hamiltonian of the molecular system is usually only handled numerically already for the tri-atomic case and difficulties are encountered when trying to rationalize the outcome of the computation.
In the present paper we examine this problem in the framework of the exact factorization of the electron-nuclear wave function [22]. This (non-adiabatic) reformulation of the quantum-mechanical problem is used as a starting point to develop a procedure that settles the issue described above in a rigorous way. The key point in the exact factorization is that the electronic effect on the nuclear system is taken into account by timedependent vector and scalar potentials. These concepts are the generalization of similar, but static, quantities appearing also within the BO approximation. We show that nonadiabatic effects can be accounted for, by formulating a theory that treats these effects as a perturbation to the BO problem. Such a framework has been discussed in previous work [23] to derive the nuclear velocity perturbation theory [24] for vibrational circular dichroism [25]. As we will show below, here we propose a new perspective on the nuclear velocity perturbation theory, which will allow us to access a broader class of both static, e.g. energetics, and dynamical, e.g. vibrational spectra, problems in quantum mechanics.
Within the nuclear velocity perturbation theory, non-adiabatic effects can be included by taking into account corrections to the BO approximation up to within linear order in the classical nuclear velocity. We show here that this is equivalent to a perturbation approach where the small parameter is the electron-nuclear mass ratio.
The major achievement of such formulation is presented in this paper: electronic nonadiabatic effects appear as a position-dependent mass correction to the bare nuclear mass, up to within linear order in the perturbation. From a fundamental perspective, we prove that it is possible to recover an adiabatic-like structure of the Hamiltonian governing the dynamics of the heavy degrees of freedom, with a kinetic energy contribution and a separate potential energy term. Since the mass correction can be fully identified with the electronic mass, totally missing in the BO approximation, we propose a theory able to restore a fundamental property, often overlooked, of the dynamical problem: the translational invariance of an isolated system with its physical mass, i.e. nuclear and electronic.
If in the BO approximation the nuclear masses are made position-dependent in the way proposed in this paper, the center of mass can be separated from rotations and internal vibrations and evolves as a free particle with mass equal to the total mass of the system (expected from the Galilean invariance of the problem [5]). This property enables us to apply the perturbation approach before moving to the molecular center of mass reference frame, with the formal advantage of a very simple and intuitive theory. From an algorithmic perspective, the corrections to the mass involve only ground-state properties and can be calculated as a response to the nuclear motion, within standard perturbation theory [26][27][28][29]. Therefore, we are able to perform numerical studies of molecular systems, easily pushing the applications beyond di-and tri-atomic molecules. The experimental implications are clear: the approach proposed here has the potential to predict and to describe ro-vibrational spectroscopic data for a large class of molecular systems when high accuracy is required.
The paper is organized as follows. First we show how, starting from the exact factorization, non-adiabatic effects are included by constructing a perturbative scheme based on the BO approach. Then, we prove that the vector potential of the theory can be expressed as a position-dependent correction to the bare nuclear mass. In the nuclear Hamiltonian, non-adiabatic effects are taken into account in an adiabatic-like picture, if the nuclear masses are corrected for the electronic contribution. We prove that (i) the position-dependent corrections sum up to the total electronic mass of the complete system and (ii) the Hamiltonian with position-dependent dressed masses is appropriate to compute rotational and vibrational spectra as it is possible to exactly separate the center of mass motion. Results are presented, discussing a model of a hydrogen bond and corrections to the vibrational frequencies of small molecular systems.

A. Exact factorization of the electron-nuclear wave function
The exact factorization of the electron-nuclear wave function has been presented [22] and discussed [30,31] in previous work. Therefore, we only introduce here the basic formalism and we refer to the above references for a detailed presentation.
A system of interacting particles, which will be taken as electrons of mass m e and nuclei of masses M ν , is described by the HamiltonianĤ =T n +Ĥ BO , withT n the nuclear kinetic energy andĤ BO the standard BO Hamiltonian. The evolution of the electronnuclear wave function Ψ(r, R, t), in the absence of an external time-dependent field, is described by the time-dependent Schrödinger equation (TDSE)ĤΨ = i ∂ t Ψ. The symbols r, R collectively indicate the Cartesian coordinates of N el electrons and N n nuclei, respectively, in a fix laboratory frame. When the exact factorization is employed, the solution of the TDSE is written as the product Ψ(r, is the nuclear wave function and Φ R (r, t) is an electronic factor parametrically depending on the nuclear configuration R. Φ R (r, t) satisfies the partial normalization condition dr|Φ R (r, t)| 2 = 1 ∀ R, t, which makes the factorization unique up to a gauge transformation. Starting from the TDSE for Ψ(r, R, t), Frenkel's action principle [32][33][34] and the partial normalization condition yield the evolution equations for Φ R (r, t) and χ(R, t), Here, the electronic and nuclear Hamiltonians areĤ respectively. The index ν is used to label the nuclei. The electron-nuclear coupling operator (ENCO), the time-dependent vector potential (TDVP), and the time-dependent potential energy surface (TDPES), mediate the exact coupling between the two subsystems, thus they include all effects beyond BO. The symbol . . . r indicates integration over the electronic coordinates. The TDVP and TDPES transform [22] as standard gauge potentials when the electronic and nuclear wave functions transform with a phase θ(R, t). The gauge, the only freedom in the definition of the electronic and nuclear wave functions, will be fixed below.

B. Large nuclear mass limit
Starting from the XF described above, we now consider the limit of large nuclear masses. The ENCO is inversely proportional to the nuclear masses M ν , then the BO limit [35] corresponds to the solution of Eqs. (1) setting the ENCO to zero [23]. Formally, however, approaching this limit of large but finite nuclear masses depends on the physical situation considered [36]. In the time-dependent case, keeping fixed the kinetic energy, it has been shown [36] that the BO limit is recovered asymptotically in terms of a small expansion parameter µ 4 used to scale the nuclear mass, M → M (µ) ≡ M/µ 4 . Making µ approach zero corresponds to the ratio of the nuclear mass over the electron mass M (µ) /m e going to infinity. This scaling factor will be used only to estimate perturbatively the order of the terms in the electronic equation, and will be set equal to unity to recover the values of the physical masses. The nuclear mass being made larger, the nuclear dynamics is slower such that time variable must then be scaled as well, by a factor µ 2 , i.e.
t → t/µ 2 [36], increasing the separation of time-scales between the light and heavy particles. Similarly, following a simple scaling argument, the nuclear momentum behaves as µ −2 in the semi-classical limit (see Appendix A). Then, the ENCO from Eq. (2) scales with where λ ν (R, t) = µ 2 −i ∇ν χ(R,t) . λ ν (R, t) tends towards a quantity independent of µ in the limit of small µ, since −i ∇ ν χ/χ is related to the nuclear momentum [23,31] and thus scales as µ −2 .
Using the definition in Eq. (4), we define the scaled TDPES, noting the second term in Eq. (5) does not contribute (by construction) to the TDPES.

C. Perturbative expansion
The electronic equation thus obtained, can be solved perturbatively in powers of µ 4 , with its solution of the form Φ R (r, t) = [1,37]. The time dependence appears only at order µ 2 , as it is clear from Eqs. (6) and (7).
Therefore the time dependence of Φ with (0) (R) the first term on the right-hand-side of Eq. (6). Here, ϕ The electronic equation at the next order yields where λ ν (R, t) = [λ ν (R, t) + µ 2 A ν (R, t)]/M ν from Eq. (2). We neglected the TDVP from the term in parenthesis since A ν (R, t) is O(µ 2 ). Furthermore, λ ν contains a term O(µ 2 ), which will be analyzed below along with the TDVP. Appendix B presents the connection between Eq. (9) and the nuclear velocity perturbation theory, thus providing a numerical scheme [23] to compute Φ (1) R (r, t) within perturbation theory [24]. The electronic wave function up to within O(µ 2 ) is where ϕ R,ν (r) is implicitly defined by Eq. (9). Eq. (10) is valid also as initial condition, i.e. the correction is included if at the initial time the nuclear velocity (the classical limit of λ ν (R, t)) is non-zero [38].
Φ R (r, t) is complex and can thus sustain an electronic current density [39,40] induced by the nuclear motion. The crucial point is that this current influences the nuclear motion through the TDVP.

D. Expression of the time-dependent vector potential
The TDVP becomes non-zero when inserting Eq. (10) in Eq. (3). As described in Ap- The singly underlined symbols A(R, t), λ (R, t) and ϕ (1) Since λ (R, t) depends on A(R, t), we find A(R, t) self-consistently, which amounts to include an infinite number of terms of order µ 2n .
Here, M ≡ M ν δ νi,ν j is the (3N n × 3N n ) diagonal mass matrix. If µ 4 = 1, expressions where the physical masses appear are recovered. From Eq. (11) it is evident that A(R) is a purely electronic quantity, which affects the nuclear momentum through the TDVP.
Such correction, however, also appears in the nuclear evolution equation (1).

E. Nuclear time-dependent Schrödinger equation
Using the expressions of the TDVP and of the TDPES, as described in Appendix D, we get an important result: the nuclear TDSE, in matrix form, becomes where the superscript T indicates the transpose vector and The second term is the diagonal BO correction (DBOC). The kinetic energy term now involves dressed nuclear masses. It is important to notice that such canonical form of the nuclear TDSE arises from the self-consistent solution for A.
The corresponding classical Hamiltonian [5] is simply The A-matrix is the new and fundamental quantity introduced in this study, for which we are able to provide a rigorous derivation, in the context of the exact factorization, an intuitive interpretation, in terms of electronic mass carried along by the motion of the nuclei, and an efficient computation scheme, based on perturbation theory [23].

III. PROPERTIES OF THE DRESSED POSITION-DEPENDENT MASS
When Cartesian coordinates are employed as done here, the A-matrix has the property of yielding the total electronic mass of the system when summed up over all nuclei, supporting its interpretation as a correction term to the nuclear mass (indices ν and ν run over the nuclei, i and j over the three spatial dimensions). It should also be noticed that the A-matrix is positive-definite in a ground-state dynamics [5]. The proof of Eq. (15) uses the property of the BO electronic wave function of being invariant under a translation of the reference system [5,41], and Eq. (9) (see Appendix E). This leads to where BO is the electronic contribution to the atomic polar tensor, defined as the variation with respect to nuclear positions of the electronic dipole moment (here averaged over the BO state) [42]. The second equality in Eq. (16) is obtained using the known property of the atomic polar tensor of yielding the total electronic charge of the system when summed over all nuclei [41,43].
It is common to separate the center of mass (CoM) motion before introducing the BO approximation. Within the molecular frame, the procedure presented here can be straightforwardly applied, by choosing coordinates in which the kinetic energy operator is the sum of two separated terms, i.e. nuclear and electronic. Using the above sum rule, Eq. (16), it is instead possible to separate of the CoM motion a posteriori and recover in that case the full mass of the system.
Starting from the Cartesian coordinates, we make the following change of coordinates with M tot = ν M ν + m e N el . From the sum rule (16), the nuclear Hamiltonian of Eq. (13) becomesĤ P CoM is the momentum (operator) associated to the center of mass (CoM) coordinate in Eq. (17), thus the first term accounts for the motion of the CoM as a free particle. The mass associated to the CoM is, correctly, the total mass of the system, i.e. nuclei and electrons, rather than the nuclear mass only, as in the BO approximation. The following terms in Eq. (18) are the kinetic and potential energies corresponding to the internal, rotational and vibrational, degrees of freedom (see Appendix F for a detailed derivation).

IV. APPLICATIONS
The formalism introduced above is employed to construct a numerical procedure that is (i) fundamentally adiabatic, namely only a single (static) potential energy surface is ex-

A. Proton transfer
As a first application, we consider a model of a proton involved in a one-dimensional [44], in which non-adiabatic effects are known to be important [45]. The light particle is the proton, assumed to be in its vibrational ground state. The mass ratio with the heavy particles, the two oxygens, is much larger than the electron-nuclear mass ratio, thus suggesting possible deviations from the BO approximation. We use an asymmetric potential mimicking a strong hydrogen bond (as shown in Fig. 1 Proton density corresponding to the BO ground state.
large distances we expect the effective mass of O − to be close to 17 a.m.u. as it carries along the proton. This is clear in Fig. 3, where it is shown that the element the A-matrix tends to a constant (equal to 1 a.m.u., the mass of the proton) at R > 3Å, whereas all other components are zero, as expected from the sum rule of Eq. (16). We show this schematically in Fig. 4 where we plot the proton density along the O−O bond.
We also report an estimate of the amount of electronic mass associated to each oxygen, as the sum over the columns of the A-matrix, e.g.
At short distances instead the proton is shared by the oxygens: the elements of the A- matrix are non-zero, but the O − diagonal contribution remains dominant. Notice that it is not surprising that the off-diagonal elements of the A-matrix are negative, as only two conditions are physically relevant: the diagonal elements must be non-negative, in a ground-state dynamics, and the sum of the elements must yield the electronic mass, in a translationally invariant system. As seen in Fig. 4, the two oxygens have then similar masses at very short distances.   Compared to Ehrenfest dynamics, the BO+M dynamics is much less computationally expensive, having a similar cost as the BO dynamics itself. Furthermore, the proposed Hamiltonian formulation allows for a full quantum treatment of the nuclear dynamics: it Overall, also in the static situation the mass correction leads to highly accurate results.
At a mass ratio µ −4 = 1600 an accuracy on the eigenvalues of about 10 −5 cm −1 is reached whereas it is only 0.5 cm −1 using the BO approximation. The results are shown as functions of the inverse mass ratio µ −4 .

B. Corrections to harmonic frequencies
Next, we consider non-adiabatic effects on vibrational frequencies and predict cor-  proposed [4][5][6]9] to cure some fundamental inconsistencies of the BO treatment. The novelty of the present study is thus to be found in the overall picture that our work conveys: the theory is developed based on a rigorous starting point, the exact factorization of the molecular wave function; the perturbation treatment is justified in terms of the electron-nuclear mass ratio, as in the seminal paper of Born and Oppenheimer; the algebraic procedure is very simple, easily allowing for applications not restricted to di-and tri-atomic molecules; the proposed numerical scheme requires standard electronic structure calculations to determine the mass corrections, as the expression of such corrections are explicitly given in terms of electronic properties. We expect that the theory will be able to provide solid information to predict and interpret highly accurate spectroscopy experiments on a large class of molecular systems.
Conceptually, we have resolved a well-known [54] fundamental inconsistency of the BO approximation. In a translationally invariant problem, the CoM moves as a free particle with mass that equals the total mass of the systems, i.e. nuclei and electrons, not only the nuclear mass. This feature is naturally built in the theory and corrects for a deficiency of the BO approximation, providing exactly the missing mass of the electrons. From a more practical point of view, our approach is very general and can be applied whenever a "factorization" of the underlying physical problem is possible, e.g. in the case of proton and oxygen atoms or in the case of electrons and nuclei.
Further applications are indeed envisaged, since the perturbative incorporation of non-adiabatic effects greatly reduces the complexity of the fully coupled problem. For instance, the approximations can be applied to nuclear wave packet methods for the calculation of highly accurate vibrational spectra beyond the BO approximation. The position dependent mass is also shown to be related to the ionicity of the bonds and may serve as a proxy to access electronic properties.

Appendix A: The adiabatic limit of the exact factorization
We will argue in this section that (i) the correct scaling of the time variable is µ 2 , when the parameter µ 4 is used to scale the nuclear mass M ν , and (ii) the term in the electronnuclear coupling operator of Eq. (2) containing the nuclear wave function scales as well with µ 2 .
Statement (i) is obtained by taking the large mass limit, or small µ 4 limit, as in [36]. In this situation, the dynamics of the heavy nuclei becomes semi-classical and our scaling argument will make the nuclear kinetic energy tend towards a constant. In the classical Following Ref. [36], the nuclear wave packet can be considered to be a Gaussian wave packet localized at the position R(t) , with momentum P (t): with σ(t) a (3N n × 3N n ) symmetric matrix yielding the spatial extension of the wave packet. From this expression, we see that statement (ii) holds: tends towards a quantity independent of µ.

Appendix B: Nuclear velocity perturbation theory
In this section we show the relation between the µ 4 −expansion proposed in the paper and the nuclear velocity perturbation theory (NVPT) of Ref. [23]. We recall here the definition of λ ν (R, t), In the framework of NVPT we have used λ ν (R, t)/M ν as the perturbation parameter that controls the degree of non-adiabaticity of the problem. The electronic equation (7) can be written using λ ν (R, t) as Also, as will be proved in Appendix C, the time-dependent vector potential (TDVP) is itself O(µ 2 ), thus it will be neglected from the term in square brackets on the right-handside. If we solve this equation order by order, Eqs. (8) and (9) are easily obtained. In particular, we recall here Eq. (9) whose solution yields Φ In Ref. [23] we started from the electronic Hamiltonian of the form and we have solved it perturbatively, usingĤ BO as the unperturbed Hamiltonian. It is clear, as stated above, that λ ν (R, t) is the small parameter that controls the strength of the perturbation and that −i ∇ ν is the (non-adiabatic) perturbation. We have looked for the eigenstates ofĤ el in the form as straightforwardly follows from the application of standard time-independent perturbation theory. The first order perturbation to the BO ground state can be written as with ω e0 (R) = ( R r , the non-adiabatic coupling vectors. This leads to a new expression of Φ R (r, t), which is exactly Eq. (10) when setting µ 2 = 1, to obtain the physical nuclear mass.
In the framework of NVPT, the perturbation parameter has been interpreted classically as the nuclear velocity [23,31,55]. It is worth mentioning here that, when performing a numerical simulation, such dependence on the nuclear velocity shall be correctly accounted for, also in the preparation of the initial electronic state. When using NVPT to perform the calculations, the electronic evolution is not explicit, in the sense that at each time the electronic wave function is simply reconstructed using ground state properties that are then inserted in Eq. (B8). However, when NVPT results are (or can be) compared with quantum-mechanical fully non-adiabatic results, the initial electronic state cannot be simply prepared in the ground state, unless the initial nuclear velocity is zero. If this is not the case, then the first order contribution in Eq. (B8), proportional to the finite value of the initial nuclear velocity, has to be included in the initial condition. Then NVPT and non-adiabatic results can be directly compared, as the same initial conditions are used in both.
Equating the first order corrections to the BO eigenstate, from the µ 4 − and the λ ν (R, t)− expansion, yields The comparison between the µ 4 −expansion and NVPT allows, first of all, to derive an ex- R (r) as given in Eq. (B7), and, second, to decompose the perturbed state as a sum of independent (linear) responses to the non-adiabatic perturbations, thus leading to As above, the index ν is used to label the nuclei and α labels the Cartesian components of the gradient. This equation can now be easily solved by employing density functional perturbation theory as described in Ref. [23].

Appendix C: Analysis of the perturbation parameter
The TDVP, defined in Eq. (3), is written using Eq. (10) as Up to within the linear order in µ 2 (or more precisely µ 2 λ ν (R, t)), this expression is where we can use Eq. (B10) to identify the A-matrix, We derive the following expression of the TDVP, namely Once again we keep the term O(µ 2 ) in λ , but we will show below how it will be included in the definition of the small parameter λ. A(R) is a matrix, thus the double-underlined notation, with (3N n × 3N n ) elements, whereas ϕ (1) R (r) is a vector with (3N n ) components. We have written also the TDVP and the parameter in matrix notation, with A(R, t) and λ (R, t) (3N n )−dimensional vectors. The elements of the A-matrix are with i, j labeling the Cartesian components and ν , ν the nuclei. When using Eq. (B7), the elements of the A-matrix can be written in terms of the non-adiabatic coupling vectors and of the BO eigenvalues as from which it follows that the A-matrix is symmetric. The A-matrix is also positive definite (i.e. for all non-zero real vectors v, the relation v T Av ≥ 0 holds) with non-negative diagonal elements, i.e.
This property is essential for the interpretation of the A-matrix as a position-dependent mass. The components of the TDVP can be expressed in terms of the components of the A-matrix, This expression is used in the definition of the parameter λ νi (R, t), given in Eq. (B1), where which, we recall, tends towards a quantity independent of µ if µ → 0.
Writing Eq. (C9) in matrix form and solving for λ(R, t) we obtain where M is a diagonal (3N n × 3N n ) matrix containing the masses of the nuclei and we have defined a position-dependent mass matrix M(R). This equation can be inverted to yielding the TDVP in the form given in Eq. (12) where only λ appears.
Eq. (C8) shows that the TDVP is at least first order in the perturbation parameter and this is the reason why it is not considered in the definition of the perturbed electronic Hamiltonian in Eq. (B5). Due to the explicit dependence of A ν (R, t) on λ ν (R, t), which is known via the A-matrix, we have been able to isolate the "actual" small parameter, i.e.
λ(R, t). In all expressions, however, we find λ (R, t), the matrix product of M −1 (R) and λ(R, t), which is a gauge-invariant quantity.
We use here a simplified notation, also using the property that the only time-dependence in the electronic wave function appears via λ ν (R, t). Using this form of the electronic wave function, we write and, by using the partial normalization condition up to within second order, since the normalization condition is already satisfied at zero-th order. We insert this result in Eq. (D5) to obtain In the second term on the right-hand-side we identify the A-matrix and we thus write where Eq. (D10) is a rewriting of Eq. (D9) in matrix form. Inserting the expression of λ (R, t) in terms of λ(R, t) given in Eq. (C12), we can express the second term of Eq. (D10) which exactly cancels the second term on the right-hand-side of Eq. (D3). The nuclear Hamiltonian of Eq. (13) is thus derived, The potential energy is time-independent and contains the BO energy, from the first term in Eq. (D10), and an additional contribution, according to It is worth noting that the first order contribution to the time-dependent potential (R, t) is zero, thus only (0) (R), the zeroth order term, appears as potential energy in the nuclear Hamiltonian of Eq. (13). This statement has been already proven in Ref. [23] using the definition in Eq. (4) and the expression of the electronic wave function up to within first order terms in the perturbation. The second term on the right-hand-side is referred to as Born-Huang diagonal correction in the applications proposed in the paper. Among the second order contributions to the potential energy (it appears at the order µ 4 in Eq. (6)), only this term beyond (0) BO (R) will be included in the calculations, due to the fact that at this stage the theory does not allow us to efficiently compute higher order terms.
The correspondence principle of quantum mechanics enables us to determine the classical nuclear Hamiltonian as where P = M(R)Ṙ is the nuclear momentum.

Appendix E: Electronic mass and the A-matrix
The derivation of Eq. (15) uses the property of the BO electronic wave function of being invariant under a translation of the coordinate reference system, namely ϕ which is valid for all values of ∆. Identifying ∇ k as the position representation of the momentum operatorp k corresponding to the k-th electron (divided by −i ), which can be written also asp and projecting the two terms in square brackets in Eq. (E1) onto ϕ (1) R,νi (r), from Eq. (B7), we identify the A-matrix on the left-hand-side and, for each Cartesian component j, we From the term on the right-hand-side we derive the expression of the atomic polar tensor (APT). First of all we write explicitly the commutator and we use Eq. (B10) to obtain then we identify the expectation value of the electronic dipole moment operator over the BO wave function in the following expression The derivative with-respect-to the i-th when we further sum over the index ν. This result states that when the A-matrix is summed up over all nuclei it yields the total electronic mass of the complete system. In Eqs. (13) and (D11) this means that the mass effect of the electrons is completely taken into account by the position-dependent mass corrections to the nuclear masses within the order of the perturbation considered here.
as we will now prove. First of all, we recall the expression of the position-dependent mass matrix,

Calculation of the A-matrix
We have computed the A-matrix using density functional perturbation theory [24,26,28,29,47,56] as described in Ref. [23] and checked that the sum rule of Eq. (16) is satisfied. The numerical scheme has been implemented in the electronic structure package CPMD [48]. Calculations have been performed using Troullier-Martins [49] pseudopotentials in the Becke-Lee-Yang-Parr [50,51] (BLYP) approximation of the exchangecorrelation kernel. The molecular geometry is the equilibrium geometry at the BLYP level, employing the aug-cc-pVTZ basis set [52] in the Gaussian electronic structure program [53].
and similarly for the other Cartesian components. This result is obtained by summing the entries of the matrix in Table II, and we find indeed the total electronic mass (m = 1, N el = 2) of the system as expected from Eq. (16). The evaluation of Eq. (E4) using only local pseudo potentials in the commutator in Eq. (E2) gives rise to the A-matrix contribution due to the local pseudo-potentials, in the following termed local part of the A-matrix. The local part of the A-matrix is indeed symmetric and has positive diagonal elements. However, it does not satisfy the sum rule of Eq. (16), as shown in Table III. In order to correct for this error, we can calculate the correction due to the full commutator where also the effect of the non-local pseudo-potential is included. For all 3N n nuclear coordinates, labeled by the indices i, ν, we obtain a commutator for each Cartesian component j. The correction hence gives rise to a (3N n × 3)−dimensional matrix ∆ ij ν , i.e. the non-local contribution to the electronic APT. However, the appropriate dimension of the matrix to be used to correct the A-matrix should be (3N n × 3N n ), as the A-matrix itself. Unfortunately, there is no protocol that allows us to match the dimensions of the two matrices, i.e. the A-matrix and the correction matrix, based on some physical properties. Therefore, we develop such protocol according to the following prescription. The correction matrix is denoted ∆A ij νν and is shown in Table IV. We add the symmetric part of the correction ∆ ij ν , in Table V,

Normal mode analysis
It is easy to prove that given a Lagrangian of the form we neglect the velocity-dependent term, we obtain with K the Hessian matrix computed from the ground state electronic potential. The term in square brackets is evaluated at the equilibrium geometry R 0 . The diagonalization of the matrix in square brackets yields corrected ν + ∆ν frequencies, as ∆ν includes the effect of electrons that follow the motion of the nuclei non-adiabatically, namely not instantaneously.