Classical nuclear motion coupled to electronic non-adiabatic transitions

Classical nuclear motion coupled to electronic non-adiabatic transitions

Federica Agostini Max-Planck Institut für Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany    Ali Abedi Max-Planck Institut für Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany    E. K. U. Gross Max-Planck Institut für Mikrostrukturphysik, Weinberg 2, D-06120 Halle, Germany
August 21, 2019

We present a detailed derivation and numerical tests of a new mixed quantum-classical scheme to deal with non-adiabatic processes. The method is presented as the zero-th order approximation to the exact coupled dynamics of electrons and nuclei offered by the factorization of the electron-nuclear wave function [A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Rev. Lett. 105 (2010)]. Numerical results are presented for a model system for non-adiabatic charge transfer in order to test the performance of the method and to validate the underlying approximations.

31.15.-p, 31.50.-x, 31.15.xg, 31.50.Gh

I Introduction

Among the ultimate goals of condensed matter physics and theoretical chemistry is the atomistic description of phenomena such as vision Polli et al. (2010); Hayashi et al. (2009); Chung et al. (2012), photo-synthesis Tapavicza et al. (2011); Brixner et al. (2005), photo-voltaic processes Rozzi et al. (2013); Silva (2013); Jailaubekov et al. (2013), proton-transfer and hydrogen storage Sobolewski et al. (2002); do N. Varella et al. (2006); Fang and Hammes-Schiffer (1997); Marx (2006). These phenomena involve the coupled dynamics of electrons and nuclei beyond the Born-Oppenheimer (BO), or adiabatic, regime and therefore require the explicit treatment of excited states dynamics. Being the exact solution of the full dynamical problem unachievable for realistic molecular systems, as the numerical cost for solving the time-dependent Schrödinger equation (TDSE) scales exponentially with the number of degrees of freedom, approximations need to be introduced. Usually, a quantum-classical (QC) description of the full system is adopted, where only a small number of degrees of freedom, e.g. electrons or protons, are treated quantum mechanically, while the remaining degrees of freedom are considered as classical particles, e.g. nuclei or ions. In this context, the challenge resides in determining the force that generates the classical trajectory as effect of the quantum subsystem. As the BO approximation breaks down, this electronic quantum effect on the classical nuclei cannot be expressed by the single adiabatic potential energy surface (PES) corresponding to the occupied eigenstate of the BO Hamiltonian. An exact treatment would require to take into account several adiabatic PESs that are coupled via electronic non-adiabatic transitions in regions of strong coupling, as avoided crossings or conical intersections. In the approximate QC description, however, the concept of single PES and single force that drives the classical motion is lost. In order to provide an answer to the question “What is the classical force that generates a classical trajectory in a quantum environment subject to non-adiabatic transitions?”, different approaches to QC non-adiabatic dynamics have been proposed for the past 50 years Ehrenfest (1927); McLachlan (1964); Shalashilin (2009); Tully (1990); Kapral and Ciccotti (1999); Agostini et al. (2007); Curchod et al. (2011); Bonella and Coker (2005); Doltsinis and Marx (2002); Ben-Nun et al. (2000); Maddox (1995); Aleksandrov (1981); Pechukas (1969); Tully and Preston (1971); Meyer and Miller (1980); Subotnik (2010); Wu and Herman (2005); Brown and Heller (1981); Coker and Xiao (1995); Prezhdo and Rossky (1997); Kosloff and Hammerich (1991); Antoniou and Schwartz (1996); Stock (1995), but the problem still remains a challenge.

This paper investigates an alternative point of view on this longstanding problem. The recently proposed exact factorization of the electron-nuclear wave function Abedi et al. (2010, 2012) allows to decompose the coupled dynamics of electrons and nuclei such that a time-dependent vector potential and a time-dependent scalar potential generate the nuclear evolution in a Schrödinger-like equation. The time-dependent potentials represent the exact effect of the electrons on the nuclei, beyond the adiabatic regime. This framework offers the same advantages of the BO approximation, since the single force that generates the classical trajectory in the QC description can be determined from the potentials. We have extensively investigated the properties of the scalar potential Abedi et al. (2013a); Agostini et al. (2013, 2014) in situations where it can be calculated exactly, by setting the vector potential to zero with an appropriate choice of the gauge, also in the context of the QC approximation. Here, we describe a procedure Abedi et al. (2014) to derive an approximation to the vector and scalar potentials, that leads to a new mixed QC (MQC) approach to the coupled non-adiabatic dynamics of electrons and nuclei. This MQC scheme is presented as zero-th order approximation to the exact electronic and nuclear equations and intends to adequately describe those situations where nuclear quantum effects associated with zero-point energy, tunnelling or interference are not relevant. It is worth noting that the factorization, and the consequent decomposition of the dynamical problem, is irrespective of the specific properties, e.g. the masses, of the two sets of particles. Therefore, it can be generalized to any two-component system.

The paper is organized as follows. In Section II we recall the exact factorization Abedi et al. (2010, 2012). The procedure to derive the classical limit of the nuclear equation is described in Section III, with reference to the most common approaches to QC dynamics. The new MQC scheme is presented in Section IV and applied to a model system for non-adiabatic charge transfer in Section V. Here, apart from testing the performance of the new algorithm, the hypothesis underlying the classical approximation are validated by the numerical analysis. Our conclusions are presented in Section VI.

Ii Theoretical background

In the absence of an external field, the non-relativistic Hamiltonian


describes a system of interacting nuclei and electrons. Here, denotes the nuclear kinetic energy and is the BO Hamiltonian, containing the electronic kinetic energy and all interactions . As recently proven Abedi et al. (2010, 2012), the full wave function, , solution of the TDSE


can be written as the product


of the nuclear wave function, , and the electronic wave function, , which parametrically depends of the nuclear configuration Hunter (1975). Throughout the paper the symbols indicate the coordinates of the electrons and nuclei, respectively. Eq. (3) is unique under the partial normalization condition (PNC)


up to within a gauge-like phase transformation


The evolution equations for and ,


are derived by applying Frenkel’s action principle Kreibich et al. (2004); Frenkel (1934); McLachlan (1964) with respect to the two wave functions and are exactly equivalent Abedi et al. (2010, 2012) to the TDSE (2). Eqs. (6) and (7) are obtained by imposing the PNC Alonso et al. (2013); Abedi et al. (2013b) by means of Lagrange multipliers.

The electronic equation (6) contains the electronic Hamiltonian


which is the sum of the BO Hamiltonian and the electron-nuclear coupling operator ,


In Eq. (6), is the time-dependent PES (TDPES), defined as


and , along with the vector potential ,


mediate the coupling between electrons and nuclei in a formally exact way. Here, the symbol stands for an integration over electronic coordinates.

The nuclear evolution is generated by the Hamiltonian


according to the time-dependent Schrödinger-like equation (7). The nuclear wave function, by virtue of the PNC, reproduces the exact nuclear N-body density


that is obtained from the full wave function. Moreover, the nuclear N-body current density can be directly obtained from


thus allowing the interpretation of as a proper nuclear wave function.

The scalar and vector potentials are uniquely determined up to within the gauge transformations


The uniqueness can be straightforwardly proved by following the steps of the current density version Ghosh and Dhara (1988) of the Runge-Gross theorem Runge and Gross (1984). In this paper, as a choice of gauge, we introduce the additional constraint , on the gauge-dependent component of the TDPES. Therefore, this scalar potential will be only expressed in terms of its gauge-invariant Agostini et al. (2013, 2014) part, . It follows that the explicit expression of the TDPES is


where the second line is obtained from the action of the operator in Eq. (10) on the electronic wave function.

In the following, the electronic equation will be represented in the adiabatic basis. Therefore, it is worth introducing here the set of eigenstates of the BO Hamiltonian with eigenvalues . We expand the electronic wave function in this basis


as well as the full wave function


The coefficients of Eqs. (17) and (18) are related by


which follows from Eq. (3). is interpreted as the amount of nuclear density that evolves “on” the -th BO surface, as the nuclear density can be written as


by using the PNC and the orthonormality of the adiabatic states. In this basis, the PNC reads


Iii The classical limit

The nuclear wave function, without loss of generality, can be written as van Vleck (1928)


with a complex function. We suppose now that this function can be expanded as an asymptotic series in powers of , namely . Inserting this expression in Eq. (7), the lowest order term, , satisfies the equation


with defined as


Eq. (23) is obtained by considering only terms up to and is formally identical to the Hamilton-Jacobi equation, if is identified with the classical action and, consequently, with the th nuclear momentum,


Therefore, is a real function and Eq. (24) is the classical Hamiltonian corresponding to the quantum operator introduced in Eq. (12). It is worth noting that the canonical momentum derived from the classical Hamiltonian, as in the case of a classical charge moving in an electromagnetic field, is


A more intuitive, less rigorous, step in the process of approximating the nuclei with classical particles is the identification of the nuclear density with a -function Tully (1997), namely


where the symbol indicates the classical positions at time . is the classical trajectory.

Eqs. (25) and (27) are the conditions for the nuclear degrees of freedom to behave classically. They can be obtained by performing the following limit operations


Eq. (28a) follows from the fact that we consider only terms in Eq. (22). In Eq. (28b), indicates the variance of a Gaussian-shaped nuclear density, centered at the classical positions , that becomes infinitely localized as in Eq. (27), when the limit operation is performed. The effect of Eqs. (28a) and  (28b) is


where the term appears explicitly in the definition of the electron-nuclear coupling operator given in Eq. (9). In order to prove Eq. (29), we replace the nuclear wave function with its -expansion and we then take the limit . Only the zero-th order term survives, leading to Eq. (29), equivalent to


It will appear clear later that such term in the electronic equation is responsible for the non-adiabatic transitions induced by the coupling to the nuclear motion, as other MQC techniques, like the Ehrenfest method or the trajectory surface hopping Tully (1997); Barbatti (2011); Drukker (1999), also suggested. Here, we show that this term can be derived as the limit in the exact equations, but it represents only the lowest order contribution in a -expansion. Moreover, this coupling is expressed via , that is not the canonical momentum appearing in the classical Hamiltonian (whose expression is given in Eq. (26)).

Eq. (30) implies that the moduli of the coefficients in the expansion (17) become constant functions of when the nuclear density is infinitely localized at the classical positions. Indeed, when the classical approximation strictly applies, the delocalization or the splitting of a nuclear wave packet is negligible (). Therefore, any -dependence can be ignored and only the instantaneous classical position becomes relevant. It is worth underlining that this same hypothesis is at the basis of several MQC approaches Tully (1990, 1997); Barbatti (2011) and applies to the moduli and to the phases of the coefficients . Here, we try to give a rigorous explanation of this assumption for the moduli, , and we comment on the spatial dependence of the phases, , in the following section.

Let us first suppose that is a normalized Gaussian wave packet, namely


centered at with variance . In the classical limit, reduces to a -function at , consequently at each point where is zero, all terms on the right-hand-side of Eq. (20) have to be zero, since they are all non-negative. Therefore, should become -functions at . Since we are interested in this limit, we represent each by a not-normalized Gaussian ( is not normalized), centered at different positions, , than . Using this hypothesis, Eq. (20) becomes


where accounts for the normalization of . The pre-factors have been introduced because each is not normalized, is instead a normalized Gaussian centered at with variance . The variances and are allowed to be time-dependent even if we are not explicitly indicating this dependence. We will prove the following statements

(i) (34)
(ii) (35)

To this end, we compare the behavior of both sides of Eq. (33) for : we need to show that


where we wrote explicitly the expressions of the Gaussian functions in Eq. (33). The leading term on the right-hand-side of Eq. (36) has to satisfy the condition


or, equivalently,


A similar argument is applied to the Fourier Transform () of both sides of Eq. (33)




(and similarly for ) with . The Gaussian transforms to another Gaussian with inverse variance, and if we calculate the limit , we obtain a relation similar to Eq. (38), namely


Eqs. (38) and (41) are simultaneously satisfied if


and this proves statement (34). In order to prove statement (35), we study the behavior of Eq. (33) at , namely


Since the pre-factors sum up to one, the relation


must hold, since if . Using Eq. (19), we can now show that


or in other words is only a function of time and is constant in space. It is worth noting that Eqs. (42) and (44) have to be valid at all times.

As a consequence of the discussion presented so far, we obtain that the exact (quantum mechanical) population of the BO states as function of time


must equal the population calculated as


when the classical limit is performed.

iii.1 Spatial dependence of the phases

We will discuss in this section why the phases of the coefficients in Eq. (17) will be considered constant functions of , similarly to the moduli as shown in the previous section. We mention again that this hypothesis is usually introduced in the derivation of other MQC approaches Tully (1990, 1997); Barbatti (2011) and we try here to justify this approximation.

First of all, we introduce the expression of the vector potential in the adiabatic basis


and we define the quantities


We used here the definition of the first order non-adiabatic coupling (NAC) . Moreover, the symbol will be used to indicate the phase of the coefficient , of the expansion (18).

We will show now that the dependence on the index of the phases can be dropped. This is done by relating the quantity to the rate of displacement of the mean value(s) of and . We are going to use again the hypothesis of a Gaussian-shaped nuclear density, which infinitely localizes at in the classical limit.

Two equivalent exact expressions for the expectation value of the nuclear momentum are obtained from Eqs. (3) and (18), namely


where is the phase of . The identity of the terms in square brackets under the integral signs follows,


where the approximations used here are (i) replacing with from Eq. (23) (then using Eq. (25) for the nuclear momentum) and (ii) neglecting the spatial dependence of (as proven, in the previous section, to be consistent with the classical treatment of the nuclei). As will be derived in Appendix A, the term on the left-hand-side of Eq. (52) can be also written as


where is the displacement rate of , the mean value of the nuclear density. If we define the quantity , Eq. (52) can be re-written as


where the PNC has been used. By analogy with Eq. (53), the term in parenthesis can be defined as


the momentum associated to the motion of the mean value () of the BO-projected wave packet . This is also consistent with our previous statement Agostini et al. (2014) on the connection between the phase and the propagation velocity of based on semi-classical arguments. Therefore, we obtain